育种数据分析之放飞自我分享 http://blog.sciencenet.cn/u/yijiaobai 关注:生物统计,数量遗传,混合线性模型,生物信息,R,Perl,Python,GWAS,GS相关方法,文章及代码

博文

表型相关与遗传相关及其显著性检验(数据及代码分享)

已有 7861 次阅读 2016-7-31 17:13 |个人分类:农学统计|系统分类:科研笔记

2018-10-09 修改备注:

1, 将代码修改, AAfun包不太好装, 使用VSNR包

2, VSNR包中有model.comp和asreml.lrt两个函数, 可以进行模型的比较和LRT检验

3, 代码correlation.R

correlation.R

附件重新上传.

4, 贴上了我的最新的微信公众号


表型相关与遗传相关及其显著性检验

数据及代码下载:

correlation.R

df.Rdata(注意数据下载后放在D盘根目录下才可运行该程序。)


# 导入数据
load("d:/df.Rdata")

# 查看数据前六行
head(df)

# 查看数据结构
str(df)

# 载入asreml软件
library(asreml)

# 载入VSNR软件, 如果没有, 使用devtools进行安装
# devtools::install_github("VSNC/VSNR") # 使用这行代码进行安装VSNR包
library(VSNR)

mod <- asreml(cbind(h1,h2) ~ trait + trait:Rep, # Rep为固定因子
               random=~ us(trait):Fam, # Fam为随机因子
               rcov=~ units:us(trait), # 残差
               data=df)
summary(mod)$varcomp

####表型相关##############
# 随机因子us变为diag,残差us变为diag
VSNR::pin(mod,pCORR ~ (V2+V6)/sqrt((V1+V5)*(V3+V7)))

mod1.1 <- asreml(cbind(h1,h3) ~ trait + trait:Rep, # Rep为固定因子
              random=~ diag(trait):Fam, # Fam为随机因子
              rcov=~ units:diag(trait), # 残差
              data=df)
VSNR::model.comp(mod,mod1.1)
VSNR::asreml.lrt(mod,mod1.1)
############环境相关###############
# 残差us变为diag
mod1.2 <- asreml(cbind(h1,h3) ~ trait + trait:Rep, # Rep为固定因子
                 random=~ us(trait):Fam, # Fam为随机因子
                 rcov=~ units:diag(trait), # 残差
                 data=df)
VSNR::model.comp(mod,mod1.2)
VSNR::asreml.lrt(mod,mod1.2)

##############基因型相关#############
# 随机因子us变为diag
mod1.3 <- asreml(cbind(h1,h3) ~ trait + trait:Rep, # Rep为固定因子
                 random=~ diag(trait):Fam, # Fam为随机因子
                 rcov=~ units:us(trait), # 残差
                 data=df)
VSNR::model.comp(mod,mod1.2)
VSNR::asreml.lrt(mod,mod1.2)




结果:

> load("d:/df.Rdata")

> head(df)

 TreeID Spacing Rep   Fam Plot    dj    dm    wd h1  h2  h3  h4  h5

1  80001       3   1 70048    1 0.334 0.405 0.358 29 130 239 420 630

2  80002       3   1 70048    2 0.348 0.393 0.365 24 107 242 410 600

3  80004       3   1 70048    4 0.354 0.429 0.379 19  82 180 300 500

4  80005       3   1 70017    1 0.335 0.408 0.363 46 168 301 510 700

5  80008       3   1 70017    4 0.322 0.372 0.332 33 135 271 470 670

6  80026       3   1 70002    2 0.359 0.450 0.392 30 132 258 390 570

> str(df)

'data.frame':827 obs. of  13 variables:

$ TreeID : Factor w/ 827 levels "80001","80002",..: 1 2 3 4 5 6 7 8 9 10 ...

$ Spacing: Factor w/ 2 levels "2","3": 2 2 2 2 2 2 2 2 2 2 ...

$ Rep    : Factor w/ 5 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...

$ Fam    : Factor w/ 55 levels "70001","70002",..: 44 44 44 15 15 2 2 10 10 10 ...

$ Plot   : Factor w/ 4 levels "1","2","3","4": 1 2 4 1 4 2 4 1 2 3 ...

$ dj     : num  0.334 0.348 0.354 0.335 0.322 0.359 0.368 0.358 0.323 0.298 ...

$ dm     : num  0.405 0.393 0.429 0.408 0.372 0.45 0.509 0.381 0.393 0.361 ...

$ wd     : num  0.358 0.365 0.379 0.363 0.332 0.392 0.388 0.369 0.347 0.324 ...

$ h1     : int  29 24 19 46 33 30 37 32 34 28 ...

$ h2     : int  130 107 82 168 135 132 124 126 153 127 ...

$ h3     : int  239 242 180 301 271 258 238 290 251 243 ...

$ h4     : int  420 410 300 510 470 390 380 460 430 410 ...

$ h5     : int  630 600 500 700 670 570 530 660 600 630 ...

> library(asreml)

> library(AAfun)

> fm.2 <- asreml(cbind(dj,h5) ~ trait + trait:Rep,

+ random=~ us(trait):Fam,

+ rcov=~ units:us(trait),

+ subset=Spacing=="3",

+ data=df,

+ maxit=100)

ASReml: Sun Jul 31 16:04:07 2016


US matrix updates modified 2 times to remain positive definite.

    LogLik         S2      DF      wall     cpu

 -5159.9423      1.0000  1105  16:04:07     0.0 (6 restrained)

 -3826.6081      1.0000  1105  16:04:07     0.0 (6 restrained)

US matrix updates modified 1 times to remain positive definite.

 -6789.9048      1.0000  1105  16:04:07     0.0 (6 restrained)

US matrix updates modified 1 times to remain positive definite.

 -2842.6010      1.0000  1105  16:04:07     0.0 (6 restrained)

US matrix updates modified 1 times to remain positive definite.

 -2285.0717      1.0000  1105  16:04:07     0.0 (1 restrained)

US matrix updates modified 1 times to remain positive definite.

 -1740.3287      1.0000  1105  16:04:07     0.0 (6 restrained)

 -1288.1124      1.0000  1105  16:04:07     0.0 (1 restrained)

Logliklihood decreased to   -1463.44 - trying again with reduced updates

 -1269.4379      1.0000  1105  16:04:07     0.0 (3 restrained)

Notice: NonPosDef US matrix modified

 -1193.0692      1.0000  1105  16:04:07     0.0 (3 restrained)

Notice: NonPosDef US matrix modified

 -1085.8035      1.0000  1105  16:04:07     0.0 (3 restrained)

 -1006.3407      1.0000  1105  16:04:07     0.0 (1 restrained)

 -1002.5588      1.0000  1105  16:04:07     0.0 (2 restrained)

Notice: NonPosDef US matrix modified

  -985.8516      1.0000  1105  16:04:07     0.0 (1 restrained)

  -984.2175      1.0000  1105  16:04:07     0.0 (3 restrained)

Notice: NonPosDef US matrix modified

  -968.4023      1.0000  1105  16:04:07     0.0 (1 restrained)

  -967.1408      1.0000  1105  16:04:07     0.0 (3 restrained)

Notice: NonPosDef US matrix modified

  -951.9053      1.0000  1105  16:04:07     0.0 (1 restrained)

  -950.5182      1.0000  1105  16:04:07     0.0 (2 restrained)

  -885.9153      1.0000  1105  16:04:07     0.0 (1 restrained)

Notice: NonPosDef US matrix modified

  -875.3218      1.0000  1105  16:04:07     0.0 (1 restrained)

Notice: NonPosDef US matrix modified

  -868.8819      1.0000  1105  16:04:07     0.0

  -868.1653      1.0000  1105  16:04:07     0.0

  -867.0717      1.0000  1105  16:04:07     0.0

  -867.0247      1.0000  1105  16:04:07     0.0

  -867.0246      1.0000  1105  16:04:07     0.0

  -867.0246      1.0000  1105  16:04:07     0.0

US variance structures were modified in 6 instances to make them positive definite


Finished on: Sun Jul 31 16:04:07 2016


LogLikelihood Converged

> summary(fm.2)$varcomp

                         gamma component std.error z.ratio constraint

trait:Fam!trait.dj:dj  6.39e-05  6.39e-05  2.21e-05    2.89   Positive

trait:Fam!trait.h5:dj  7.55e-02  7.55e-02  4.65e-02    1.62   Positive

trait:Fam!trait.h5:h5  4.56e+02  4.56e+02  1.90e+02    2.40   Positive

R!variance             1.00e+00  1.00e+00        NA      NA      Fixed

R!trait.dj:dj          4.98e-04  4.98e-04  3.14e-05   15.84   Positive

R!trait.h5:dj         -2.13e-01 -2.13e-01  7.36e-02   -2.89   Positive

R!trait.h5:h5          5.33e+03  5.33e+03  3.36e+02   15.83   Positive

> ####表型相关##############

> pin(fm.2,pCORR ~ (V2+V6)/sqrt((V1+V5)*(V3+V7)),signif=T)


     Estimate     SE sig.level

pCORR  -0.0763 0.0449         *

---------------

Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1


> ############环境相关###############

> pin(fm.2,eCORR ~ V6/sqrt(V5*V7),signif = T)


     Estimate     SE sig.level

eCORR   -0.131 0.0441        **

---------------

Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1


> ##############基因型相关#############

> pin(fm.2,gCORR ~ V2/sqrt(V1*V3),signif = T)


     Estimate    SE sig.level

gCORR    0.442 0.257         *

---------------

Sig.level: 0'***' 0.001 '**' 0.01 '*' 0.05 'Not signif' 1




数据及模型介绍

两个性状xy在同一环境或者同一性状在不同年份地点会表现出不同程度的表型、遗传和环境相关,这种相关分为早晚相关和性状相关。

要计算相关先要计算出相关的方差、协方差组分,现用ASReml为例讲解如何得到相关结果。

数据来源见参考文献,截图如下:

01


模型如下:

02


模型解释:

要分析的性状为djh5Rep为固定因子,Fam为随机因子,trait代表两性状的联合,选择Spacing3的数据进行数据分析,数据集为df,最大迭代次数为100

方差和协方差结果:

03


结果解释:

trait1的基因型方差分量:VGtrait1,这里为6.38e-05V1

trait1trait2的基因型协方差分量:COVGtrait1.trait2,这里为7.55e-02V2

trait2的基因型方差分量:VGtrait2,这里为4.56e+02V3

trait1的误差方差分量:VRtrait1,这里为4.98e-04V5

trait2的误差方差分量:VRtrait2,这里为-2.13e-01V6

trait1trait2的误差协方差分量:COVRtrait1.trait2,这里为5.33e+03V7

表型相关定义、计算方法及软件代码:

是指用两个性状表型值求得的相关系数rp成为表型相关系数,它反映的是性状表现型之间的关系。

04



转化为:


05


表型相关以及显著性测验:

06


境相关定义、计算方法及软件代码:

环境相关是指两性状在同一环境下所体现的相关。

计算方法:

07



08


遗传相关定义、计算方法及软件代码:

是指同一遗传材料的两个性状间由于遗传原因所体现的相关,也即这个性状基因型值间的相关。广义的遗传相关(rg),用两性状的基因型协方差(covgxy)与各该性状的基因型标准差(σgxσgy)乘积之比度量。由于基因型值又可分解为加性的和非加性的两个部分,两性状加性效应值间的加性遗传相关称为狭义遗传相关。植物育种中常用广义遗传相关,动物育种中则常用狭义遗传相关。

遗传相关以及显著性结果:

09



10


表型相关和遗传相关的关系:

遗传相关是表型相关(γp)的一个组成部分,若相关的两个性状均具有较高的遗传力,则表型相关系数(γp)主要由遗传相关系数(γg)所决定。反之,若两性状的遗传力均甚低,则γp主要由环境相关系数(γe)所决定。由于遗传相关已经剔除了环境影响,故能比表型相关更确切地反映两个性状间的相关程度。

遗传相关的应用及注意事项:


育种工作中可应用遗传相关确定与育种目标性状有联系的某些性状的相对重要性。倘若某性状的遗传与育种目标性状的遗传高度相关,就可利用它对育种目标性状进行相关选择(或称间接选择)。这种方法特别适用于育种目标性状难以准确测量或遗传力甚低的情况,因为在这种情况下可通过对遗传力高且与育种目标性状高度相关的性状的选择,有效地改进育种目标性状。例如,林木的早期选择或林木性状的早期预测必须了解成熟性状与幼龄材性状之间的遗传相关,对目标性状进行间接选择必须了解目标性状与选择性状之间的遗传相关。




注意:在做双性状分析前,一般先对各性状进行独自的单性状分析,获得各自的最佳模型,并得到各性状的加性遗传方差和误差方差,这对于后续双性状分析的模型很有意义,因为可以用来构建R矩阵和G结构矩阵的初始值设置。

参考文献:

林元震. RASReml-R统计分析教程[M]. 中国林业出版社, 2014. 212~215

王富德, 张世苹. 表型相关系数和遗传相关系数的计算和应用[J]. 辽宁农业科学, 1982(1).


对农业数据感兴趣可以关注这个公众号:


qrcode_for_gh_0e98b5fcf43c_344.jpg





https://wap.sciencenet.cn/blog-2577109-993718.html

上一篇:R语言汇总统计,包括分组数据的汇总,结果友好(代码分享)
下一篇:无人机大批量采集玉米表型数据(文献解析)

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2021-12-9 18:25

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部