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

发布于

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

1, 概念定义

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

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

两者的区别和联系:

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

两者的计算方法:

  • 可以使用通径分析的方法进行计算
  • 也可以采用由系谱构建亲缘关系A矩阵的形式进行计算, 这种方法在数据比较大时更为方便

2, 系谱数据

这里我们模拟了四个个体的系谱关系, 想要计算一下每个个体的近交系数, 以及个体间的亲缘系数, 使用R语言实现.

代码语言：javascript
复制
ped <- data.frame(ID=c(3,4,5,6),Sire=c(1,1,4,5),Dam=c(2,NA,3,2))
ped

ID

Sire

Dam

3

1

2

4

1

NA

5

4

3

6

5

2

首先, 我们看到个体1和2的亲本未知, 所以我们先将系谱补充完整, 使用nadiv的prepPed函数.

代码语言：javascript
复制
library(nadiv)
pped = prepPed(ped)
pped

完整的系谱如下, NA表示亲本未知.

3, 计算亲缘关系A矩阵

代码语言：javascript
复制
as.matrix(makeA(pped))

这里我们使用makeA函数, A矩阵如下:

1

2

3

4

5

6

1

1.00

0.000

0.5000

0.5000

0.5000

0.2500

2

0.00

1.000

0.5000

0.0000

0.2500

0.6250

3

0.50

0.500

1.0000

0.2500

0.6250

0.5625

4

0.50

0.000

0.2500

1.0000

0.6250

0.3125

5

0.50

0.250

0.6250

0.6250

1.1250

0.6875

6

0.25

0.625

0.5625

0.3125

0.6875

1.1250

4, 计算近交系数

用亲缘关系矩阵A的对角线-1，即是个体的近交系数
代码语言：javascript
复制
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行.

第一种方案:

代码语言：javascript
复制
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 查看结果

结果如下:

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函数, 更方便.

代码语言：javascript
复制
library(learnasreml)
re2 = mat_2_coefficient(mat)
head(re2)

结果和上面是一致的.

所有代码汇总如下:

代码语言：javascript
复制
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)

希望可以帮到你.

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

#na#var