如何利用系譜計算近交係數和親緣關係係數

• 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)

希望可以幫到你.

轉發朋友圈是最大的支援.