技術標籤：head 10位元組真彩色影象資料量 計算

今天我們看一下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秒