HIBLUP教程｜第8期：重复记录性状模型拟合

在生产实践中，对于窝产仔数、精子活力与密度等性状，需要对不同胎次或在不同的采集时间进行多次记录，获得同一个体的多个观测数据，进而构建重复记录模型进行遗传参数估计。HIBLUP也可以拟合重复记录性状模型，其方法与拟合单性状模型几乎相同，拟合单性状模型的具体方法可见HIBLUP教程｜第6期：单性状模型遗传参数估计。唯一的区别在于，在拟合重复记录模型时需要将表型文件第一列的ID分配给环境随机效应参数，即“–rand 1”，来估计永久环境效应。

1、One Step 拟合重复记录性状模型

以GBLUP模型为例，HIBLUP拟合重复记录模型的命令行如下：

hiblup \
  --single-trait \
  --pheno demo.repeat.phe \
  --pheno-pos 8 \
  --dcovar 2,3 \
  --qcovar 4,5 \
  --rand 1 \
  --bfile demo \
  --out demo 

与单性状模型相同，HIBLUP在拟合重复记录模型时通过–single-trait识别单性状模型，当检测到命令行中仅包含表型文件(–pheno)和基因型文件(–bfile)时，HIBLUP将自动拟合GBLUP模型。重复记录的表型文件如示例文件demo.repeat.phe所示：

重复记录的表型文件中包含了性状多次测量的表型值以及记录的环境因素，并含表头（即列名）。其中，第一列必须为个体ID；其后的第二至第七列中记录了每次测量时的性别、季节、日龄、初生重、场、母本等因素；第八列为个体T1性状多次测量的表型值，同一个体可以具有多个表型值。

–pheno-pos：用于指定需要分析的性状在表型文件中所在列的位置；

–rand：在拟合重复记录模型时需要指定表型文件第一列的ID作为环境随机效应来剖分永久环境效应方差。

程序运行效果如下：

如上图所示，在打印的日志信息中，HIBLUP会简单列出当前分析的模型公式，用户可通过模型公式检查输入的参数是否和需要拟合的模型一致，如下：

T1 = 1 + sex(F) + season(F) + day(C) + bornweight(C) + id(R) + GA(GR) + e

其中，T1为表型的名称，性别和季节为固定效应（F: fixed effect），日龄和初生重为协变量（C: covariate），个体多次记录ID作为永久环境效应（R: random），基因型关系矩阵(GA: genome-based relationship matrix)为加性遗传随机效应（GR: genetic random）。程序运行后生成demo.vars、demo.beta、demo.anova、demo.rand和demo.log共5个文件。

所有随机效应的估计方差和遗传力储存在demo.vars文件中，格式如下：

其中第一列为各随机效应对应的名称，包括永久环境效应（id），加性遗传随机效应（GA）和残差效应（e）；第二列为估计的随机效应方差；第三列为随机效应方差的标准误；第四列为随机效应的遗传力；第五列为随机效应遗传力的标准误。

所有固定效应及协变量的回归系数估计值保存在demo.beta文件中，格式如下：

文件的第一列为均值、固定效应的各个水平以及协变量的名称；第二列为对应的估计值；第三列为估计值的标准误。

方差分析表保存在demo.anova文件中，格式如下：

方差分析表的第一列为固定效应和协变量名称；第二列为自由度；第三列为总平方和；第四列为均方；第五列为F检验的统计量F值；第六列为显著性检验得到的P值，可以通过P值判断该效应的影响是否显著。

所有个体的随机效应预测值保存在demo.rand文件中，格式如下：

其中，首列为个体ID，其后的列包括永久环境效应、加性遗传随机效应（即个体育种值）和残差效应的估计值。

同理，当检测到输入命令行中仅包含表型文件（–pheno）和系谱文件(–pedigree)时，HIBLUP将自动拟合重复记录性状的PBLUP模型，命令行如下：

hiblup \
  --single-trait \
  --pheno demo.repeat.phe \
  --pheno-pos 8 \
  --dcovar 2,3 \
  --qcovar 4,5 \
  --rand 1 \
  --pedigree demo.ped \
  --out demo

当检测到输入命令行中同时包含表型文件（–pheno）、系谱文件（–pedigree）以及基因型文件(–bfile)时，HIBLUP将自动拟合重复记录性状的SSGBLUP模型，命令行如下：

hiblup \
  --single-trait \
  --pheno demo.repeat.phe \
  --pheno-pos 8 \
  --dcovar 2,3 \
  --qcovar 4,5 \
  --rand 1 \
  --bfile demo \
  --pedigree demo.ped \
  --out demo

2、Two Steps 拟合重复记录性状模型

由于HIBLUP采用基于方差协方差V矩阵的计算策略，其矩阵维度等于表型记录数，对于非重复记录性状，表型记录数小于总个体数，在内存消耗和计算效率上具有明显优势；但对于重复记录性状，表型记录数远大于总个体数，如几百个体具有数万甚至几十万的表型记录，此时按照上面的方式运行HIBLUP，可能导致内存溢出等问题。

因此，对于数据量庞大的重复记录性状（如：表型记录数是个体数10倍以上），除了使用上述一步法直接拟合重复记录模型以外，还可以使用两步法进行计算：首先在R中使用lme4包的lmer函数构建不包含遗传随机效应的模型（重复力模型），即表型文件中的个体ID作为随机效应，同时加入固定效应和协变量等环境因素，计算遗传参数；再将计算得到的随机效应值作为表型值通过HIBLUP中的–pheno参数加入到模型中，并添加遗传随机效应进行遗传参数估计。具体流程如下：

第一步，在R中使用lme4包的lmer函数构建混合效应模型，提取随机效应值（即BLUP值），代码如下所示：

install.packages("lme4")  
library(lme4)  
phe <- read.table("demo.repeat.phe",header = TRUE)  
#拟合BLUP模型，同时加入固定效应和协变量等环境因素：
mod <- lmer(T1~1+sex+season+day+bornweight+(1|id), data=phe)
summary(mod)

使用summary函数查看模型的汇总信息，如下图所示：

图中Random effects表单包含所有随机效应的方差，其中id对应为性状的遗传方差和永久环境效应方差的总和（Vt），Residual对应性状的残差方差（Ve），性状的总方差为所有随机效应方差的总和，即Vp=Vt+Ve，性状的重复力可通过Vt/Vp计算获得。

将个体ID与得到的随机效应值合并，写入mod.phe文件中，作为运行HIBLUP新的表型文件：

blup <- ranef(mod)$id  
mod_rand <- data.frame(id=rownames(blup), T1_blup=blup[,1]) 
head(mod_rand,4)
      id   T1_blup
1 IND1001 -6.151284
2 IND1002  5.409942
3 IND1003  7.911197
4 IND1004 10.495767
write.table(mod_rand,"mod.phe",quote = FALSE,row.names = FALSE,col.names = TRUE)  

第二步：使用HIBLUP将上一步得到的个体BLUP值作为表型，剖分永久环境效应和遗传随机效应的命令如下：

hiblup \
  --single-trait \
  --pheno mod.phe \
  --pheno-pos 2 \
  --bfile demo \
  --out demo_lmer

注意：此时无需在命令行中加入“–rand 1”。

程序运行效果如下：

需要注意的是，与一步法拟合不同，上图结果中GA对应为性状的加性遗传方差Va，e对应为永久环境效应方差Vp，并不是性状的残差方差，性状的残差方差Ve在第一步拟合lmer函数结果中获得，此时的h2也并非性状的遗传力，而是加性遗传方差和永久环境效应方差的占比情况，性状的遗传力的计算需要结合第一步获得的表型方差Vp，即h2=Va/Vp。

程序运行后得到的育种值储存在.rand文件的GA列中，可以在R中通过以下代码计算一步法与两步法获得的育种值之间的相关系数：

demo <- read.table("demo.rand",header = TRUE)  
demo_lmer <- read.table("demo_lmer.rand",header = TRUE)
index <- match(demo_lmer$ID, demo$ID)   
cor(demo$GA[index], demo_lmer$GA)

两种方案获得的育种值之间的比较：

如上图所示，两种方案获得的育种值相关系数几乎为1，说明结果可靠。

感谢大家的持续关注，有关软件的任何建议与疑问，欢迎大家加入QQ交流群，交流群号：935348004