如何利用系谱计算近交系数和亲缘关系系数

• 2019 年 12 月 6 日
• 筆記

《线性模型在动物育种值预测中的应用》 第二章：亲属间的遗传协方差，P19

1, 概念定义

近交系数： 近交系数（inbreeding coefficient）是指根据近亲交配的世代数，将基因的纯化程度用百分数来表示即为近交系数，也指个体由于近交而造成异质基因减少时，同质基因或纯合子所占的百分比也叫近交系数，个体中两个亲本的共祖系数。

亲缘系数: 是指两个个体间加性基因效应间的相关.

两者的区别和联系:

• 近交系数是个体的值
• 亲缘系数是两个个体之间的值

两者的计算方法:

• 可以使用通径分析的方法进行计算
• 也可以采用由系谱构建亲缘关系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表示亲本未知.

3, 计算亲缘关系A矩阵

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] #1  id = row.names(A) #2  mat = matrix(0,n,n) #3  for(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 #5  re = data.frame(ID1= rep(id,n),ID2=rep(id,each = n),y = round(as.vector(mat))) #6  re#7

这里我们的编程思路如下:

#1 计算出矩阵的行, 确定循环数

#2 计算出个体的ID名在矩阵中的顺序, 因为有些ID可能是字符或者没有顺序, 主要用于后面的个体编号的确定

#3 为了计算更快, 我们生成一个6*6的矩阵

#4 写一个循环, 因为矩阵时对称的, 所以我再第二个for循环时从i开始, 而不是从1开始, 后面mat[j,i] = mat[i,j]再赋值, 这样更快.

#5 生成的mat矩阵查看

#6 根据ID生成两列, 表示他们之间的亲缘系数, 这个和矩阵变为向量后一一对应.

#7 查看结果

结果如下:

第二种方案:

这里也可以用我写的learnasreml包的: mat_2_coefficient函数, 更方便.

library(learnasreml)  re2 = mat_2_coefficient(mat)  head(re2)

结果和上面是一致的.

所有代码汇总如下:

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)  A  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]    }  }  mat  re = data.frame(row = rep(id,n),col=rep(id,each = n),y = round(as.vector(mat)))  re    library(learnasreml)  re2 = mat_2_coefficient(mat)  head(re2)

希望可以帮到你.

转发朋友圈是最大的支持.