表型相关系数与标准误

表型相关系数
邓飞    2020-08-31    468

  原创 邓飞2013  育种数据分析之放飞自我

  来自专辑

  R语言灌水

  一个同学发信给我:

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

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

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

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

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

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

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

  首先我们生成两个数据:

  > 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

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

本文为专栏作者授权科易网发表,版权归原作者所有。文章系作者个人观点,不代表科易网立场,转载请联系原作者。如有任何疑问,请联系ky@1633.com。

熟练使用统计软件和遗传评估软件,负责开发了具有数据管理、常规遗传评估、MAS、GWAS和GS的育种云平台,具有丰富的全基因组选择分析及平台建设经验。 [查看更多]

13163阅读数 9237访问数 924订阅数
热门观点

向邓飞提问