【GS課堂】如何利用系譜計算近交系數和親緣關系系數

來源: 發表日期:2019-04-10 瀏覽量:290

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



 



 





分享:
富盈门财富