用到的軟件包：asreml

運行GetASremlInfor函數(shù)

將此代碼在R語言里面運行即可

GetASremlInfor<-function(ped){
  require(asreml)
  ainv <- asreml.Ainverse(ped)$ginv
  ani <- ainv
  head(ani)
  n<-max(ani$Row,ani$Column)###查看逆矩陣中最大的號
  mat=matrix(0,n,n)###生成N階值為的0矩陣
  mat_xiang_guan=matrix(0,n,n)###生成N階值為的0矩陣
  mat[cbind(ani$Row,ani$Column)]<-ani$Ainverse####把行號和列號對應的值賦給空矩陣
  mat[upper.tri(mat)]=t(mat)[upper.tri(t(mat))]####生成對稱矩陣生成A-1（asreml生成的逆矩陣）
  nnn_biao_ni<-mat
  biao_zhun_xi_pu_ju_zhen<-solve(nnn_biao_ni)###把原來的逆矩陣變換成系譜相關矩陣A
  ##############################求近交系數(shù)###################################
  inbreeding_coefficient<-diag(biao_zhun_xi_pu_ju_zhen)-1#求出矩陣的對角線，-1就是近交系數(shù)
  inbreeding_coefficient[inbreeding_coefficient < 0.000000001] <- 0
  matrix_number<-(1:n)###生成矩陣位置編號
  animal_code<-ped[,1]
  data<-data.frame(matrix_number,animal_code,inbreeding_coefficient)###生成一個數(shù)據(jù)框，包括序號和近交系數(shù)
  write.csv(data,"inbreeding_coefficient_ data.csv")
  ################################下面是求相關系數(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)######生成一個數(shù)據(jù)框，包括行、列、相關系數(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")
}

構建系譜信息

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ù)，會在工作路徑下輸出兩個Excel

近交系數(shù)：inbreeding_coefficient_ data.csv

親緣關系系數(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>

打印親緣關系系數(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>