【GS课堂】如何利用系谱计算近交系数和亲缘关系系数

来源: 发表日期:2019-04-10 浏览量:5856

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 查看结果

     

    结果如下:


  • ID1ID2r
    111.0000
    120.0000
    130.5000
    140.5000
    150.4714
    160.2357
    210.0000
    221.0000
    230.5000
    240.0000
    250.2357
    260.5893
    310.5000
    320.5000
    331.0000
    340.2500
    350.5893
    360.5303
    410.5000
    420.0000
    430.2500
    441.0000
    450.5893
    460.2946
    510.4714
    520.2357
    530.5893
    540.5893
    551.0000
    560.6111
    610.2357
    620.5893
    630.5303
    640.2946
    650.6111
    661.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



 



 





分享: