一些育种的基本概念，需要看教科书理解一下。这次我学习教科书，加深对相关概念的理解，最后用实际数据演示一下如何计算遗传相关及显著性。

「参考书：」

• Raphael A.Mrode编著，于向春 张豪 主译《线性模型在动物育种值预测中的应用》 第三版
• 陈国宏, 张勤. 动物遗传原理与育种方法[M]. 中国农业出版社, 2009.

1. 多性状模型为何更优

在育种过程中，经常对多个性状进行选择，这些性状可能有遗传相关。

多性状模型可以充分利用性状的表型相关和遗传相关等信息，针对两性状或多性状对个体进行评估。

多性状模型最佳线性无偏预测（MBLUP）的优点之一是增加了评估的精确性。精确度的增加取决于遗传相关和残差相关大小。这些相关的差异越大，精确性增加越多。

低遗传力和高遗传力同时进行多性状分析时，低遗传力性状收益更多。而且多性状模型分析时，因为考虑了性状间的残差协方差，所以增加了评定的准确性。

另外，多性状模型可以避免「淘汰偏差（culling bias）」，比如断奶体重和周岁体重这两个性状，很多个体都会根据断奶体重进行淘汰，因此很多个体都没有周岁体重。如果只针对周岁体重的个体进行遗传评估（选择后的个体），结果会产生偏差，因为它没有利用断奶体重的信息。所以，对断奶体重和周岁体重进行双性状模型评估，可以消除这种偏差。

2. 遗传相关的原因

遗传相关的分子原因，可以分为两类：

• 第一大类是基因的一因多效和基因间的连锁造成的性状间的遗传相关。它分为两种：
• 不同基因的连锁导致的遗传相关会随着基因间的重组而改变，这种遗传相关会随着世代的增加，基因连锁逐渐消失，由此导致遗传相关也逐渐减小。除非是高紧密连锁，一般而言，由基因连锁导致的遗传相关是不稳定遗传的
• 基因的一因多效导致的遗传相关是能够文稳定遗传的。
• 第二大类是由于两性状受个体所在的相同环境导致的相关，称为环境相关。其中，由于等违纪基因见显性和非等位基因的上位性造成的相关也不能遗传，他们都并入环境相关。

3. 遗传相关的用途

3.1 间接选择

遗传相关可用于确定间接选择的依据和预测间接选择反应大小。间接选择是指当一个性状（如X）不能直接选择或者选择效果很差时，借助与之相关的另一个性状（如Y）的选择来达到对性状X的选择目的

间接选择在育种实践中具有很重要的意义，比如有些性状只有屠宰后才能度量，有些性状只能在个体的生命晚期才能度量。有些性状受性别限制无法度量，还有些性状直接选择效果不理想，在这些情况下都可以考虑采用间接选择

3.2 不同环境下的选择

遗传相关可用于比较不同环境条件下的选择效果。实际上，不但不同性状可以来估计遗传相关，而且同一性状在不同环境下的表现也可以作为不同的性状来估计遗传相关。

不同环境下的遗传相关，为解决育种工作中的一个重要实际问题提供了理论依据，即在条件优良的种畜场选育的优良品种，推广到条件较差的其它条件生产厂是否能保持其优良特性。

实际上，遗传相关可以推断同一性状在不同环境下的反应是否一致。

3.3 多性状选择

一般而言，只要涉及两个性状以上的选择问题，无不需要用到遗传相关这一参数，用于指定相关形状的选择指数，这也是遗传相关主要的用途之一。

4. 示例演示

「数据来源：」

数据来源我编写的learnasreml包，原数据是一篇论文中的示例数据，我将其放到包中，方便调用。

library(learnasreml)
library(asreml)
data("animalmodel.dat")
data("animalmodel.ped")
head(animalmodel.dat)
head(animalmodel.ped)


dat = animalmodel.dat
ped = animalmodel.ped
dat[dat==0] = NA
str(dat)

> str(dat)
'data.frame':	1084 obs. of  6 variables:
 $ ANIMAL: Factor w/ 1084 levels "1","2","3","5",..: 864 1076 549 989 1030 751 987 490 906 591 ...
 $ MOTHER: Factor w/ 429 levels "1","2","3","8",..: 362 268 216 375 396 289 328 255 347 240 ...
 $ BYEAR : Factor w/ 34 levels "968","970","971",..: 1 1 2 2 2 2 3 3 3 3 ...
 $ SEX   : Factor w/ 2 levels "1","2": 1 1 2 1 2 1 2 1 1 1 ...
 $ BWT   : num  10.77 9.3 3.98 5.39 12.12 ...
 $ TARSUS: num  24.8 22.5 12.9 20.5 NA ...

「构建模型：」

这里使用个体动物模型，双性状模型：BWT和TARSUS两个性状。固定因子为SEX，随机因子为加性效应。

# 对BWT和TARSUS两个性状进行多性状模型分析
## 固定因子：SEX
## 随机因子：ANIMAL 

ainv = ainverse(ped)
mod = asreml(cbind(BWT,TARSUS) ~ trait + trait:SEX,
             random = ~ us(trait):vm(ANIMAL,ainv),
             residual = ~ units:us(trait),
             data=dat)

查看模型迭代情况和方差组分：

summary(mod)$varcomp

mod$converge

> summary(mod)$varcomp
                                           component std.error  z.ratio bound %ch
trait:vm(ANIMAL, ainv)!trait_BWT:BWT        2.986952 0.5218132 5.724180     P 0.0
trait:vm(ANIMAL, ainv)!trait_TARSUS:BWT     2.494674 0.9920503 2.514665     P 0.1
trait:vm(ANIMAL, ainv)!trait_TARSUS:TARSUS 12.397268 3.0675900 4.041370     P 0.0
units:trait!R                               1.000000        NA       NA     F 0.0
units:trait!trait_BWT:BWT                   2.995429 0.4183571 7.159981     P 0.0
units:trait!trait_TARSUS:BWT                3.596542 0.8368501 4.297714     P 0.1
units:trait!trait_TARSUS:TARSUS            17.767578 2.6651323 6.666678     P 0.0
> 
> mod$converge
[1] TRUE

可以看到，模型收敛，根据方差组分，计算遗传相关和表型相关：

> vpredict(mod,rg ~ V2/sqrt(V1*V3))
    Estimate        SE
rg 0.4099554 0.1192531
> 
> vpredict(mod,rp ~ (V2+V6)/sqrt((V1+V5)*(V3+V7)))
    Estimate         SE
rp 0.4534364 0.03175156

「遗传相关和表型相关的显著性」

这里使用似然比检验（LRT），另外构建两个模型：

moda：

• 加性矩阵diag
• 残差us modb：
• 加性矩阵diag
• 残差diag

moda = asreml(cbind(BWT,TARSUS) ~ trait + trait:SEX,
             random = ~ diag(trait):vm(ANIMAL,ainv),
             residual = ~ units:us(trait),
             data=dat)

modb = asreml(cbind(BWT,TARSUS) ~ trait + trait:SEX,
              random = ~ diag(trait):vm(ANIMAL,ainv),
              residual = ~ units:diag(trait),
              data=dat)
lrt.asreml(mod,moda)
lrt.asreml(mod,modb)

结果：可以看出遗传相关和表型相关都达到了极显著。

> lrt.asreml(mod,moda) # 遗传相关显著性
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

         df LR-statistic Pr(Chisq)   
mod/moda  1       7.3347  0.003382 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> lrt.asreml(mod,modb) # 表型相关显著性
Likelihood ratio test(s) assuming nested random models.
(See Self & Liang, 1987)

         df LR-statistic Pr(Chisq)    
mod/modb  2       156.22 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

「输出多性状育种值」

blup = coef(mod)$random
head(blup)

library(tidyverse)
blup_re = blup %>% mutate(ID = sub(".*)_","",rownames(blup)), 
                Trait = sub("trait_","",sub(":vm.*","",rownames(blup)))) %>% 
  select(ID,Trait,effect) %>% pivot_wider(ID,names_from = Trait,values_from = effect)

head(blup_re)

fwrite(blup_re,"blup_re.csv")

注意，原始的blup值是这个样子的：

> head(blup)
                                    effect
trait_BWT:vm(ANIMAL, ainv)_1306 -0.9098018
trait_BWT:vm(ANIMAL, ainv)_1304  0.4085106
trait_BWT:vm(ANIMAL, ainv)_1298  0.2025910
trait_BWT:vm(ANIMAL, ainv)_1293  0.5439144
trait_BWT:vm(ANIMAL, ainv)_1290  1.2819816
trait_BWT:vm(ANIMAL, ainv)_1288  0.4501628

我这里用tidyverse对其进行了清洗，将ID和性状提取出来，并且重新进行了整理，结果为三列：ID，trait1-blup, trait2-blup，然后将其输出到csv格式。

> head(blup_re)
# A tibble: 6 x 3
  ID       BWT TARSUS
  <chr>  <dbl>  <dbl>
1 1306  -0.910 -0.538
2 1304   0.409 -0.218
3 1298   0.203  2.94 
4 1293   0.544  1.04 
5 1290   1.28   1.12 
6 1288   0.450  1.08

5. 完整代码

# devtools::install_github("dengfei2013/learnasreml")
library(learnasreml)
library(data.table)
library(asreml)
data("animalmodel.dat")
data("animalmodel.ped")
head(animalmodel.dat)
head(animalmodel.ped)


dat = animalmodel.dat
ped = animalmodel.ped
dat[dat==0] = NA
str(dat)

# 对BWT和TARSUS两个性状进行多性状模型分析
## 固定因子：SEX
## 随机因子：ANIMAL 

ainv = ainverse(ped)
mod = asreml(cbind(BWT,TARSUS) ~ trait + trait:SEX,
             random = ~ us(trait):vm(ANIMAL,ainv),
             residual = ~ units:us(trait),
             data=dat)

summary(mod)$varcomp

mod$converge

vpredict(mod,rg ~ V2/sqrt(V1*V3))

vpredict(mod,rp ~ (V2+V6)/sqrt((V1+V5)*(V3+V7)))


moda = asreml(cbind(BWT,TARSUS) ~ trait + trait:SEX,
             random = ~ diag(trait):vm(ANIMAL,ainv),
             residual = ~ units:us(trait),
             data=dat)

modb = asreml(cbind(BWT,TARSUS) ~ trait + trait:SEX,
              random = ~ diag(trait):vm(ANIMAL,ainv),
              residual = ~ units:diag(trait),
              data=dat)
lrt.asreml(mod,moda) # 遗传相关显著性
lrt.asreml(mod,modb) # 表型相关显著性

blup = coef(mod)$random
head(blup)

library(tidyverse)
blup_re = blup %>% mutate(ID = sub(".*)_","",rownames(blup)), 
                Trait = sub("trait_","",sub(":vm.*","",rownames(blup)))) %>% 
  select(ID,Trait,effect) %>% pivot_wider(ID,names_from = Trait,values_from = effect)

head(blup_re)

fwrite(blup_re,"blup_re.csv")

❝

欢迎关注我的公众号：育种数据分析之放飞自我。主要分享R语言，Python，育种数据分析，生物统计，数量遗传学，混合线性模型，GWAS和GS相关的知识。

❞