【原】如何利用系譜計算近交系數(shù)和親緣關(guān)系系數(shù)

育種數(shù)據(jù)分析 2021-11-18

展開全文

《線性模型在動物育種值預(yù)測中的應(yīng)用》第二章：親屬間的遺傳協(xié)方差，P19

1, 概念定義

近交系數(shù)： 近交系數(shù)（inbreeding coefficient）是指根據(jù)近親交配的世代數(shù)，將基因的純化程度用百分?jǐn)?shù)來表示即為近交系數(shù)，也指個體由于近交而造成異質(zhì)基因減少時，同質(zhì)基因或純合子所占的百分比也叫近交系數(shù)，個體中兩個親本的共祖系數(shù)。

親緣系數(shù): 是指兩個個體間加性基因效應(yīng)間的相關(guān).

兩者的區(qū)別和聯(lián)系:

近交系數(shù)是個體的值
親緣系數(shù)是兩個個體之間的值

兩者的計算方法:

可以使用通徑分析的方法進(jìn)行計算
也可以采用由系譜構(gòu)建親緣關(guān)系A(chǔ)矩陣的形式進(jìn)行計算, 這種方法在數(shù)據(jù)比較大時更為方便

2, 系譜數(shù)據(jù)

這里我們模擬了四個個體的系譜關(guān)系, 想要計算一下每個個體的近交系數(shù), 以及個體間的親緣系數(shù), 使用R語言實現(xiàn).

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函數(shù).

library(nadiv)
pped = prepPed(ped)
pped

完整的系譜如下, NA表示親本未知.

3, 計算親緣關(guān)系A(chǔ)矩陣

as.matrix(makeA(pped))

這里我們使用makeA函數(shù), 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, 計算近交系數(shù)

用親緣關(guān)系矩陣A的對角線-1，即是個體的近交系數(shù)

diag(A)-1

可以看出, 1,2,4,3的近交系數(shù)為0. 個體5和6的近交系數(shù)為0.125.

5, 計算親緣系數(shù)

根據(jù)計算的親緣關(guān)系A(chǔ)矩陣，這個矩陣時個體間的方差協(xié)方差矩陣, 對角線為每個個體的方差, 非對角線為個體間的協(xié)方差. 公式為:

rij = cov(i,j)/sqrt(var(i)*var(j))

因此如果我們要計算個體1和2的親緣系數(shù), 可以用A12/(sqrt(A11*A22)) = 0/sqrt(1*1) = 0, 因此個體1和2 的親緣系數(shù)為0.

因為共有6個個體, 1和2的親緣系數(shù) = 2和1的親緣系數(shù), 因此他們之間的親緣系數(shù)一共有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 計算出矩陣的行, 確定循環(huán)數(shù)

#2 計算出個體的ID名在矩陣中的順序, 因為有些ID可能是字符或者沒有順序, 主要用于后面的個體編號的確定

#3 為了計算更快, 我們生成一個6*6的矩陣

#4 寫一個循環(huán), 因為矩陣時對稱的, 所以我再第二個for循環(huán)時從i開始, 而不是從1開始, 后面mat[j,i] = mat[i,j]再賦值, 這樣更快.

#5 生成的mat矩陣查看

#6 根據(jù)ID生成兩列, 表示他們之間的親緣系數(shù), 這個和矩陣變?yōu)橄蛄亢笠灰粚?yīng).

#7 查看結(jié)果

結(jié)果如下:

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函數(shù), 更方便.

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

結(jié)果和上面是一致的.

所有代碼匯總?cè)缦?

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)