head 10字节_优秀了！10万系谱，计算近交系数，不到1秒！

今天我们看一下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 Dam
1  101 Sire_1   0
2  102 Sire_1   0
3  103 Sire_1   0
4  104 Sire_1   0
5  105 Sire_1   0
6  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)
       Inbreeding
Sire_1          0
Sire_2          0
Sire_3          0
Sire_4          0
Sire_5          0
Sire_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 101
Sire_1    1.0      0      0      0      0      0      0      0      0 0.5
Sire_2    0.0      1      0      0      0      0      0      0      0 0.0
Sire_3    0.0      0      1      0      0      0      0      0      0 0.0
Sire_4    0.0      0      0      1      0      0      0      0      0 0.0
Sire_5    0.0      0      0      0      1      0      0      0      0 0.0
Sire_6    0.0      0      0      0      0      1      0      0      0 0.0
Sire_7    0.0      0      0      0      0      0      1      0      0 0.0
Sire_8    0.0      0      0      0      0      0      0      1      0 0.0
Sire_9    0.0      0      0      0      0      0      0      0      1 0.0
101       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秒