GS课堂第五期
【GS课堂】如何利用系谱计算近交系数和亲缘关系系数
育种过程中选种选配的原则是: 优中选优, 控制近交。那么如何控制近交呢?这就需要计算近交系数和亲缘关系系数。具体怎样进行计算,就请跟着我一起往下看吧!
1.概念定义
近交系数: 近交系数是指根据近亲交配的世代数,将基因的纯化程度用百分数来表示;也指个体由于近交而造成异质基因减少时,同质基因或纯合子所占的百分比,个体中两个亲本的共祖系数。
亲缘系数: 是指两个个体间加性基因效应间的相关。
两者的区别和联系:
近交系数是个体的值
亲缘系数是两个个体之间的值
两者的计算方法:
可以使用通径分析的方法进行计算
也可以采用由系谱构建亲缘关系A矩阵的形式进行计算, 这种方法适用于数据量较大时使用
2.系谱数据
这里我们模拟了四个个体的系谱关系, 想要计算一下每个个体的近交系数, 以及个体间的亲缘系数, 使用R语言实现。
ped <- data.frame(ID=c(3,4,5,6),Sire=c(1,1,4,5),Dam=c(2,NA,3,2))ped
首先, 我们看到个体1和2的亲本未知, 所以我们先将系谱补充完整, 使用nadiv的prepPed函数。
library(nadiv)pped = prepPed(ped)pped
完整的系谱如下, NA表示亲本未知。
as.matrix(makeA(pped))
这里我们使用makeA函数, A矩阵如下:
4.计算近交系数
用亲缘关系矩阵A的对角线-1,即是个体的近交系数。
diag(A)-1
可以看出, 1,2,4,3的近交系数为0。个体5和6的近交系数为0.125。
5.计算亲缘系数
根据计算的亲缘关系A矩阵,这个矩阵是个体间的方差协方差矩阵, 对角线为每个个体的方差, 非对角线为个体间的协方差。公式为:
rij = cov(i,j)/sqrt(var(i)*var(j))
因此,如果我们要计算个体1和2的亲缘系数, 可以用:
A12/(sqrt(A11*A22)) = 0/sqrt(1*1) = 0
因此个体1和2 的亲缘系数为0。
因为共有6个个体, 1和2的亲缘系数 = 2和1的亲缘系数, 因此他们之间的亲缘系数一共有6*5/2 = 15个。这里我们都计算,共有36行。
第一种方案:
n = dim(A)[1] #1id = row.names(A) #2mat = matrix(0,n,n) #3for(i in 1:n){ #4 for(j in i:n){mat[i,j] = A[i,j]/(sqrt(A[i,i]*A[j,j]))mat[j,i] = mat[i,j]}}mat #5re = data.frame(ID1= rep(id,n),ID2=rep(id,each = n),y = round(as.vector(mat))) #6re#7
这里我们的编程思路如下:
#1 计算出矩阵的行, 确定循环数
#2 计算出个体的ID名在矩阵中的顺序, 因为有些ID可能是字符或者没有顺序, 主要用于后面的个体编号的确定
#3 为了计算更快, 我们生成一个6*6的矩阵
#4 写一个循环, 因为矩阵是对称的, 所以我在第二个for循环时从 i 开始, 而不是从1开始, 后面mat[j,i] = mat[i,j]再赋值, 这样更快
#5 生成的mat矩阵查看
#6 根据ID生成两列, 表示他们之间的亲缘系数, 这个和矩阵变为向量后一一对应
#7 查看结果
结果如下:
ID1 ID2 r 1 1 1.0000 1 2 0.0000 1 3 0.5000 1 4 0.5000 1 5 0.4714 1 6 0.2357 2 1 0.0000 2 2 1.0000 2 3 0.5000 2 4 0.0000 2 5 0.2357 2 6 0.5893 3 1 0.5000 3 2 0.5000 3 3 1.0000 3 4 0.2500 3 5 0.5893 3 6 0.5303 4 1 0.5000 4 2 0.0000 4 3 0.2500 4 4 1.0000 4 5 0.5893 4 6 0.2946 5 1 0.4714 5 2 0.2357 5 3 0.5893 5 4 0.5893 5 5 1.0000 5 6 0.6111 6 1 0.2357 6 2 0.5893 6 3 0.5303 6 4 0.2946 6 5 0.6111 6 6 1.0000 第二种方案:
这里也可以用我写的learnasreml包的: mat_2_coefficient函数, 更方便。
library(learnasreml)re2 = mat_2_coefficient(mat)head(re2)
结果和上面是一致的。
6.代码汇总
模拟系谱数据:
ped <- data.frame(ID=c(3,4,5,6),Sire=c(1,1,4,5),Dam=c(2,NA,3,2))ped
library(nadiv)
pped = prepPed(ped)pped
计算近交系数:
A = as.matrix(makeA(pped))round((diag(A) -1),3
计算亲缘系数:
n = dim(A)[1]id = row.names(A)mat = matrix(0,n,n)for(i in 1:n){ for(j in i:n){mat[i,j] = A[i,j]/(sqrt(A[i,i]*A[j,j]))mat[j,i] = mat[i,j]}}matre = data.frame(row = rep(id,n),col=rep(id,each = n),y = round(as.vector(mat)))re
用learnasreml包计算亲缘系数:
library(learnasreml)re2 = mat_2_coefficient(mat)head(re2)
参考文献:
《线性模型在动物育种值预测中的应用》 第二章:亲属间的遗传协方差,P19