今天我们看一下asreml如何根据系谱计算近交系数和亲缘关系系数，asreml计算的特点是：

1，运算速度特别快，对于10万的系谱，运算时间小于10s，是其它软件不可比拟的。

2，对于一些植物或者群体分组的系谱，有些个体(或者家系)可能做父本，可能做母本，asreml软件可以处理这种系谱。

3，asreml直接计算近交系数，所以数据量大时也可以直接计算出来。

这里，我们先用示例数据演示一下，然后我们使用模拟数据(10万系谱)计算一下结果，做一个测试，结果运行时间不到1秒。

1. 载入asreml4-R包

library(asreml)

2. 查看示例系谱数据

> data(harvey)> ped = harvey[,1:3]> head(ped)  Calf   Sire Dam1  101 Sire_1   02  102 Sire_1   03  103 Sire_1   04  104 Sire_1   05  105 Sire_1   06  106 Sire_1   0

3. 计算A逆矩阵

> ainv = ainverse(ped)> head(ainv)     Row Column Ainverse[1,]   1      1 3.666667[2,]   2      2 3.666667[3,]   3      3 2.666667[4,]   4      4 3.666667[5,]   5      5 3.333333[6,]   6      6 3.000000

4. 提取近交系数

注意，提取Inbreeding，通过attr(ainv,"inbreeding")进行提取。

> re = data.frame(Inbreeding = attr(ainv,"inbreeding"))> head(re)       InbreedingSire_1          0Sire_2          0Sire_3          0Sire_4          0Sire_5          0Sire_6          0

5. 计算A矩阵

因为asreml计算的是A逆矩阵的三元组形式，需要首先转化为矩阵形式，然后求逆。asreml4-R中有相关转化的函数，简单快捷。

> mat_ainv > A_mat = zapsmall(solve(mat_ainv))> A_mat[1:10,1:10]       Sire_1 Sire_2 Sire_3 Sire_4 Sire_5 Sire_6 Sire_7 Sire_8 Sire_9 101Sire_1    1.0      0      0      0      0      0      0      0      0 0.5Sire_2    0.0      1      0      0      0      0      0      0      0 0.0Sire_3    0.0      0      1      0      0      0      0      0      0 0.0Sire_4    0.0      0      0      1      0      0      0      0      0 0.0Sire_5    0.0      0      0      0      1      0      0      0      0 0.0Sire_6    0.0      0      0      0      0      1      0      0      0 0.0Sire_7    0.0      0      0      0      0      0      1      0      0 0.0Sire_8    0.0      0      0      0      0      0      0      1      0 0.0Sire_9    0.0      0      0      0      0      0      0      0      1 0.0101       0.5      0      0      0      0      0      0      0      0 1.0

6. 完整示例代码

# 近交系数library(asreml)data(harvey)ped = harvey[,1:3]head(ped)ainv = ainverse(ped)head(ainv)

re = data.frame(Inbreeding = attr(ainv,"inbreeding"))head(re)

mat_ainv A_mat = zapsmall(solve(mat_ainv))A_mat[1:10,1:10]

7. 十万系谱数据计算近交系数

> library(data.table)> library(asreml)> ped = fread("ped_test.ped")> dim(ped)[1] 96280     3> head(ped)   Progeny Sire Dam1:       1    0   02:       2    0   03:       3    0   04:       4    0   05:       5    0   06:       6    0   0> system.time({ainv = ainverse(ped)}) 用户  系统  流逝 0.809 0.037 0.846 > dim(ainv)[1] 293068      3> head(ainv)     Row Column Ainverse[1,]   1      1       61[2,]   2      2       61[3,]   3      3       61[4,]   4      4       61[5,]   5      5       61[6,]   6      6       61> inbr = data.frame(Inbreeding = attr(ainv,"inbreeding"))> head(inbr)  Inbreeding1          02          03          04          05          06          0> tail(inbr)      Inbreeding96275 0.0457145996276 0.0457145996277 0.0457145996278 0.0457145996279 0.0457145996280 0.04571459

结论：

系谱个数：96280

计算近交系数时间：0.864秒