多性状选择，很有必要。
第一种多性状选择：
分别计算出单性状的育种值，然后根据权重进行选择。这种方法有一定效果，但是模型中没有考虑到性状间的协方差，误差较大。
第二种多性状选择：
在计算时，用多性状模型，直接得到性状的育种值，然后再根据权重进行计算综合育种值进行选择。准确性高。
肯定推荐第二种多性状模型，因为准确性高，一步到位。
1. 多性状模型为何更优
在育种过程中，经常对多个性状进行选择，这些性状可能有遗传相关。
多性状模型可以充分利用性状的表型相关和遗传相关等信息，针对两性状或多性状对个体进行评估。
多性状模型最佳线性无偏预测（MBLUP）的优点之一是增加了评估的精确性。精确度的增加取决于遗传相关和残差相关大小。这些相关的差异越大，精确性增加越多。
低遗传力和高遗传力同时进行多性状分析时，低遗传力性状收益更多。而且多性状模型分析时，因为考虑了性状间的残差协方差，所以增加了评定的准确性。
另外，多性状模型可以避免「淘汰偏差（culling bias）」，比如断奶体重和周岁体重这两个性状，很多个体都会根据断奶体重进行淘汰，因此很多个体都没有周岁体重。如果只针对周岁体重的个体进行遗传评估（选择后的个体），结果会产生偏差，因为它没有利用断奶体重的信息。所以，对断奶体重和周岁体重进行双性状模型评估，可以消除这种偏差。
2. 遗传相关的原因
遗传相关的分子原因，可以分为两类：
- 第一大类是基因的一因多效和基因间的连锁造成的性状间的遗传相关。它分为两种：
- 不同基因的连锁导致的遗传相关会随着基因间的重组而改变，这种遗传相关会随着世代的增加，基因连锁逐渐消失，由此导致遗传相关也逐渐减小。除非是高紧密连锁，一般而言，由基因连锁导致的遗传相关是不稳定遗传的
- 基因的一因多效导致的遗传相关是能够稳定遗传的。
- 第二大类是由于两性状受个体所在的相同环境导致的相关，称为环境相关。其中，由于等位基因的显性和非等位基因的上位性造成的相关也不能遗传，他们都并入环境相关。
3. 遗传相关的用途
3.1 间接选择
遗传相关可用于确定间接选择的依据和预测间接选择反应大小。间接选择是指当一个性状（如X）不能直接选择或者选择效果很差时，借助与之相关的另一个性状（如Y）的选择来达到对性状X的选择目的
间接选择在育种实践中具有很重要的意义，比如有些性状只有屠宰后才能度量，有些性状只能在个体的生命晚期才能度量。有些性状受性别限制无法度量，还有些性状直接选择效果不理想，在这些情况下都可以考虑采用间接选择
3.2 不同环境下的选择
遗传相关可用于比较不同环境条件下的选择效果。实际上，不但不同性状可以来估计遗传相关，而且同一性状在不同环境下的表现也可以作为不同的性状来估计遗传相关。
不同环境下的遗传相关，为解决育种工作中的一个重要实际问题提供了理论依据，即在条件优良的种畜场选育的优良品种，推广到条件较差的其它条件生产厂是否能保持其优良特性。
实际上，遗传相关可以推断同一性状在不同环境下的反应是否一致。
3.3 多性状选择
一般而言，只要涉及两个性状以上的选择问题，无不需要用到遗传相关这一参数，用于指定相关形状的选择指数，这也是遗传相关主要的用途之一。
4. 示例演示
「数据来源：」
数据来源我编写的learnasreml包，原数据是一篇论文中的示例数据，我将其放到包中，方便调用。
代码语言：javascript复制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，随机因子为加性效应。
代码语言：javascript复制# 对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)
查看模型迭代情况和方差组分：
代码语言：javascript复制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
可以看到，模型收敛，根据方差组分，计算遗传相关和表型相关：
代码语言：javascript复制> 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)
结果：可以看出遗传相关和表型相关都达到了极显著。
代码语言：javascript复制> 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
「输出多性状育种值」
代码语言：javascript复制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值是这个样子的：
代码语言：javascript复制> 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格式。
代码语言：javascript复制> 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. 完整代码
代码语言：javascript复制# 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")
