有朋友问, 如何根据系谱计算近交系数, 可以输出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 181mat[5,6]/sqrt(mat[5,5]*mat[6,6])# 289 VS 16mat[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计算更方便，也更快。