如何根据系谱计算近交系数并保存

有朋友问, 如何根据系谱计算近交系数, 可以输出excel的那种操作, 之前写过一篇是根据asreml的, 这里使用免费的nadiv包, 将结果输出到D盘. 同时不仅仅有输出近交系数的结果, 还有输出亲缘关系系数的结果到D盘.

• 近交系数结果:inbreeding.csv
• 亲缘关系系数结果: coeff.csv

为了证明程序的有效性, 以教科书的示例作为演示.

《动物遗传原理与育种方法》 陈国宏 张勤 主编 中国农业出版社

1, 根据系谱计算近交系数程序

第十四章 第四节：近交系数与亲缘系数计算

近交系数的公式：

书中的系谱关系

转化为系谱：ID,Sire,Dam的写法

ID <- c("X","B","A","C","D","E","H","G","I")
Sire <- c("B","E","C","E","E","H",0,0,0)
Dam <- c("A",0,"D",0,"G","G","I","I",0)
ped <- data.frame(ID,Sire,Dam)
ped

> ped
  ID Sire Dam
1  X    B   A
2  B    E   0
3  A    C   D
4  C    E   0
5  D    E   G
6  E    H   G
7  H    0   I
8  G    0   I
9  I    0   0

计算每个个体的近交系数，根据亲缘关系矩阵A，它的对角线减去1即为所在个体的近交系数

library(nadiv)
    
pped =prepPed(ped)

A = as.matrix(makeA(pped));A

re = data.frame(ID=row.names(A),inbreeding = diag(A)-1)

write.csv(re,"d:/inbreeding.csv",row.names=F)

将结果保存到D盘根目录

可以看到X的近交系数为0.1797，和书中的结果一致：

2, 亲缘关系系数

书中的系谱

转化为三列形式：ID,Sire,Dam

# 亲缘关系系数
id <- c(289,135,181,90,188,49)
sire <- c(135,90,49,16,0,16)
dam <- c(181,188,188,0,0,0)
ped <- data.frame(id,sire,dam)
ped

> ped
   id sire dam
1 289  135 181
2 135   90 188
3 181   49 188
4  90   16   0
5 188    0   0
6  49   16   0

计算方法

最方便的计算方法，计算亲缘关系矩阵A，则a和b的亲缘系数为：cov(a,b)/sqrt(var(a)*var(b))，两者的协方差除以两者方差积的开方

library(nadiv)
pped = prepPed(ped)
mat = as.matrix(makeA(pped))
# 135 VS 181
mat[5,6]/sqrt(mat[5,5]*mat[6,6])
# 289 VS 16
mat[1,7]/sqrt(mat[1,1]*mat[7,7])

0.3125
0.232495277487639
135和181的亲缘关系系数为:0.3125, 和书中结果一致.

将结果保存到D盘

id = row.names(mat)

id1 = rep(id,length(id))
id2 = rep(id,each=length(id))
coeff = matrix(0,dim(mat)[1],dim(mat)[1])
for(i in 1:dim(mat)[1]){
    for(j in 1:dim(mat)[1]){
        coeff[i,j] = mat[i,j]/sqrt(mat[i,i]*mat[j,j])
    }
}

coeff_value = as.vector(coeff)

re = data.frame(id1,id2,coeff_value)

write.csv(re,"d:/coeff.csv",row.names=F)

结论

使用通径连可以很好地理解计算方法，但是更为简单的方法是通过程序计算亲缘关系矩阵，计算近交系数和亲缘关系系数，通过nadiv计算更方便，也更快。

我的微信公众号