用到的軟件包：asreml
運(yùn)行GetASremlInfor函數(shù)
將此代碼在R語言里面運(yùn)行即可
GetASremlInfor<-function(ped){
require(asreml)
ainv <- asreml.Ainverse(ped)$ginv
ani <- ainv
head(ani)
n<-max(ani$Row,ani$Column)###查看逆矩陣中最大的號(hào)
mat=matrix(0,n,n)###生成N階值為的0矩陣
mat_xiang_guan=matrix(0,n,n)###生成N階值為的0矩陣
mat[cbind(ani$Row,ani$Column)]<-ani$Ainverse####把行號(hào)和列號(hào)對(duì)應(yīng)的值賦給空矩陣
mat[upper.tri(mat)]=t(mat)[upper.tri(t(mat))]####生成對(duì)稱矩陣生成A-1（asreml生成的逆矩陣）
nnn_biao_ni<-mat
biao_zhun_xi_pu_ju_zhen<-solve(nnn_biao_ni)###把原來的逆矩陣變換成系譜相關(guān)矩陣A
##############################求近交系數(shù)###################################
inbreeding_coefficient<-diag(biao_zhun_xi_pu_ju_zhen)-1#求出矩陣的對(duì)角線怎栽，-1就是近交系數(shù)
inbreeding_coefficient[inbreeding_coefficient < 0.000000001] <- 0
matrix_number<-(1:n)###生成矩陣位置編號(hào)
animal_code<-ped[,1]
data<-data.frame(matrix_number,animal_code,inbreeding_coefficient)###生成一個(gè)數(shù)據(jù)框晓避，包括序號(hào)和近交系數(shù)
write.csv(data,"inbreeding_coefficient_ data.csv")
################################下面是求相關(guān)系數(shù)###############################
nnn<-biao_zhun_xi_pu_ju_zhen
for (a in (1:n)){
for (b in (1:n)){
x<-nnn[a,b]/sqrt(nnn[a,a]*nnn[b,b])
mat_xiang_guan[a,b]<-x
}
}
xiang_guan_xi_shu<-mat_xiang_guan
dim(xiang_guan_xi_shu)<-NULL###把矩陣編程變量尘盼，一維的
hang<-rep(animal_code,each=n)####生成行
lie<-rep(animal_code,n)######生成列
length(lie);length(hang);length(xiang_guan_xi_shu)
xiang_guan_xi_shu
en<-data.frame(hang,lie,xiang_guan_xi_shu)######生成一個(gè)數(shù)據(jù)框，包括行、列、相關(guān)系數(shù)
real_xishu<-en[which(en$xiang_guan_xi_shu >= 0.0000001),]##把近親系數(shù)為0的列去掉
colnames(real_xishu)<-c("animal_code_1","animal_code_2","coefficient_of_coancestry")
real_xishu
write.csv(real_xishu,"coancestry_file.csv")
}
構(gòu)建系譜信息
ped <- data.frame(ID = c(1,2,3,4,5,6), Sire = c("NA","NA",1,1,4,5), Dam =c("NA","NA",2,"NA",3,2))
ped
<table>
<thead><tr><th scope=col>ID</th><th scope=col>Sire</th><th scope=col>Dam</th></tr></thead>
<tbody>
<tr><td>1 </td><td>NA</td><td>NA</td></tr>
<tr><td>2 </td><td>NA</td><td>NA</td></tr>
<tr><td>3 </td><td>1 </td><td>2 </td></tr>
<tr><td>4 </td><td>1 </td><td>NA</td></tr>
<tr><td>5 </td><td>4 </td><td>3 </td></tr>
<tr><td>6 </td><td>5 </td><td>2 </td></tr>
</tbody>
</table>
調(diào)用GetASremlInfor函數(shù)，會(huì)在工作路徑下輸出兩個(gè)Excel
近交系數(shù)：inbreeding_coefficient_ data.csv
親緣關(guān)系系數(shù)：coancestry_file.csv
打印近交系數(shù)前6行
inbreeding_value <- read.csv("inbreeding_coefficient_ data.csv")
head(inbreeding_value)
<table>
<thead><tr><th scope=col>X</th><th scope=col>matrix_number</th><th scope=col>animal_code</th><th scope=col>inbreeding_coefficient</th></tr></thead>
<tbody>
<tr><td>1 </td><td>1 </td><td>1 </td><td>0.000</td></tr>
<tr><td>2 </td><td>2 </td><td>2 </td><td>0.000</td></tr>
<tr><td>3 </td><td>3 </td><td>3 </td><td>0.000</td></tr>
<tr><td>4 </td><td>4 </td><td>4 </td><td>0.000</td></tr>
<tr><td>5 </td><td>5 </td><td>5 </td><td>0.125</td></tr>
<tr><td>6 </td><td>6 </td><td>6 </td><td>0.125</td></tr>
</tbody>
</table>
打印親緣關(guān)系系數(shù)前6行
jinjiaoxishu <- read.csv("coancestry_file.csv")
head(jinjiaoxishu)
<table>
<thead><tr><th scope=col>X</th><th scope=col>animal_code_1</th><th scope=col>animal_code_2</th><th scope=col>coefficient_of_coancestry</th></tr></thead>
<tbody>
<tr><td>1 </td><td>1 </td><td>1 </td><td>1.0000000</td></tr>
<tr><td>3 </td><td>1 </td><td>3 </td><td>0.5000000</td></tr>
<tr><td>4 </td><td>1 </td><td>4 </td><td>0.5000000</td></tr>
<tr><td>5 </td><td>1 </td><td>5 </td><td>0.4714045</td></tr>
<tr><td>6 </td><td>1 </td><td>6 </td><td>0.2357023</td></tr>
<tr><td>8 </td><td>2 </td><td>2 </td><td>1.0000000</td></tr>
</tbody>
</table>
