专栏首页育种数据分析之放飞自我表型相关系数与标准误

表型相关系数与标准误

一个同学发信给我:

“请问我想计算几个表型数据的相关系数和标准误,如何用R语言操作?”

我回答:“R中默认的函数有cor计算相关系数,标准误的话估计要用重抽样去操作?,但是很少有人会计算标准误这个数值。”

于是,该同学发过来一个截图:

上三角为表型相关系数,下三角为表型相关系数的标准误。

于是,我陷入了沉思,这里的表型相关系数,应该是和遗传相关系数对应的,属于数理遗传学的范畴,而我说的相关系数,应该是统计学的范畴。

然后我告诉该同学,你等着,我写篇公众号,解释一下这个误区。

统计学范畴的相关系数计算方法:

首先我们生成两个数据:

> set.seed(123)
> x=rnorm(1000)
> y= 0.3*x + rnorm(1000)
> dat=cbind(x,y)
> head(dat)
               x          y
[1,] -0.56047565 -1.1639414
[2,] -0.23017749 -1.1090083
[3,]  1.55870831  0.4496323
[4,]  0.07050839 -0.1110226
[5,]  0.12928774 -2.5105565
[6,]  1.71506499  1.5550930

然后使用cor函数计算x和y的相关系数:

> cor(dat)
          x         y
x 1.0000000 0.3573148
y 0.3573148 1.0000000

检验x和y的相关系数的显著性:

> cor.test(dat)
Error in cor.test.default(dat) : 缺少参数"y",也没有缺省值
> cor.test(dat$x,dat$y)

  Pearson's product-moment correlation

data:  dat$x and dat$y
t = 12.086, df = 998, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.3020116 0.4102211
sample estimates:
      cor 
0.3573148 

可以看出,两者显著性达到极显著。

数量遗传学范畴的相关系数计算方法:

一般是混合线性模型,有随机因子,比如父本,比如家系,比如个体动物模型,然后计算性状间的方差组分和协方差组分,最后计算遗传相关系数和表型相关系数。

> library(asreml)
> data("harvey")
> head(harvey)
  Calf   Sire Dam Line ageOfDam  y1  y2  y3
1  101 Sire_1   0    1        3 192 390 224
2  102 Sire_1   0    1        3 154 403 265
3  103 Sire_1   0    1        4 185 432 241
4  104 Sire_1   0    1        4 183 457 225
5  105 Sire_1   0    1        5 186 483 258
6  106 Sire_1   0    1        5 177 469 267

这里,将Sire作为随机因子,分析y1,y3两个性状的遗传相关和表型相关,以及他们的标准误。

> mod = asreml(cbind(y1,y3) ~ trait + trait:Line ,
+              random = ~ us(trait):Sire, residual = ~ units:us(trait), data=harvey)
Model fitted using the sigma parameterization.
ASReml 4.1.0 Sun Mar 29 17:46:47 2020
          LogLik        Sigma2     DF     wall    cpu
 1      -437.786           1.0    124 17:46:47    0.0
 2      -433.479           1.0    124 17:46:47    0.0
 3      -429.450           1.0    124 17:46:47    0.0
 4      -427.635           1.0    124 17:46:47    0.0
 5      -427.223           1.0    124 17:46:47    0.0
 6      -427.206           1.0    124 17:46:47    0.0
 7      -427.206           1.0    124 17:46:47    0.0
> summary(mod)$varcomp
                        component std.error    z.ratio bound %ch
trait:Sire!trait_y1:y1   27.20936  26.59364  1.0231532     P   0
trait:Sire!trait_y3:y1  -12.81265  41.71667 -0.3071350     P   0
trait:Sire!trait_y3:y3  124.88915 125.13874  0.9980055     P   0
units:trait!R             1.00000        NA         NA     F   0
units:trait!trait_y1:y1 132.36803  25.00112  5.2944842     P   0
units:trait!trait_y3:y1 -59.97696  39.91635 -1.5025662     P   0
units:trait!trait_y3:y3 647.80434 122.32450  5.2957858     P   0

注意,这里的方差组分是component一列,前三列为两个性状的遗传方差组分,后三列为环境方差组分。

遗传相关 = cov_sire_12/sqrt(V_sire_11*V_sire_22)

表型相关 = (cov_sire_12 + cov_units_12)/sqrt((V_sire_11 + V_units_11)*(V_sire_22 + V_units_22))

因此相关的计算方法:

> # 遗传相关
> vpredict(mod, rg ~ V2/sqrt(V1*V3))
     Estimate        SE
rg -0.2197948 0.6668593
> # 表型相关
> vpredict(mod, rp ~ (V2 + V6)/sqrt((V1+V5)*(V3+V7)))
     Estimate        SE
rp -0.2072908 0.1434625

如果只是单纯的使用统计的相关系数,计算如下:

> # 相关系数:统计的方法
> cor(harvey$y1,harvey$y3)
[1] -0.3208209
> cor.test(harvey$y1,harvey$y3)

  Pearson's product-moment correlation

data:  harvey$y1 and harvey$y3
t = -2.6886, df = 63, p-value = 0.009171
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.52373856 -0.08345172
sample estimates:
       cor 
-0.3208209

可以看到结果偏高,所以注意,如果你是分析的育种数据,还是要按照上面数量遗传的方法,通过计算方差组分和协方差组分,计算出的遗传相关和表型相关才是正确的。

本文分享自微信公众号 - 育种数据分析之放飞自我(R-breeding),作者:邓飞2013

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-03-29

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 如何检测遗传相关的显著性:LRT检验操作方法

    3. 软件实现遗传相关计算 3.1 数据格式 前三列是系谱,有3个性状: y1, y2, y3 , 想要计算y1和y3的遗传相关,并用LRT检验显著性

    邓飞
  • asreml 设定初始值 固定初始值

    一个朋友问我,如何固定asreml的初始值,现在分为单性状和多性状进行说明。 为何要固定初始值: 1,由于群体较小,估算的方差组分不准确,需要手动设定初始值,直...

    邓飞
  • 将多个性状画一张箱线图上R实战

    邓飞
  • 自定义 View 实战 04 - 圆形进度条

    code_horse
  • 自定义 View 实战 02 - 圆弧形进度条

    code_horse
  • 有关程序员的几个爆笑段子

    2、世界上最遥远的距离不是生与死,而是你亲手制造的BUG就在你眼前,你却怎么都找不到她。

    民工哥
  • 记录python一些技巧

    公众号---志学Python
  • 熊逸《唐诗50讲》怀古篇 - 学习笔记与感想

      在怀古篇中,熊逸老师选择了6首诗,分别cover三大怀古热门人物:贾谊、王昭君及项羽,以及一大网红怀古热门景点:秦淮河,最后是唐朝人自己的怀古—白居易的《长...

    Edison Zhou
  • Python元编程:控制你想控制的一切

    很多人不理解“元编程”是个什么东西,关于它也没有一个十分准确的定义。这篇文章要说的是Python里的元编程,实际上也不一定就真的符合“元编程”的定义。只不过我...

    小小科
  • Python元编程:控制你想控制的一切

    專 欄 ❈松直,Python中文社区专栏作者,计算机在读,Python拥趸,知乎专栏:从Python开始❈ 很多人不理解“元编程”是个什么东西,关于它也没有一...

    Python中文社区

扫码关注云+社区

领取腾讯云代金券