用到的软件包：asreml

运行GetASremlInfor函数

将此代码在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
  ##############################求近交系数###################################
  inbreeding_coefficient<-diag(biao_zhun_xi_pu_ju_zhen)-1#求出矩阵的对角线，-1就是近交系数
  inbreeding_coefficient[inbreeding_coefficient < 0.000000001] <- 0
  matrix_number<-(1:n)###生成矩阵位置编号
  animal_code<-ped[,1]
  data<-data.frame(matrix_number,animal_code,inbreeding_coefficient)###生成一个数据框，包括序号和近交系数
  write.csv(data,"inbreeding_coefficient_ data.csv")
  ################################下面是求相关系数###############################
  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)######生成一个数据框，包括行、列、相关系数
  real_xishu<-en[which(en$xiang_guan_xi_shu >= 0.0000001),]##把近亲系数为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

调用GetASremlInfor函数，会在工作路径下输出两个Excel

近交系数：inbreeding_coefficient_ data.csv

亲缘关系系数：coancestry_file.csv

打印近交系数前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>

打印亲缘关系系数前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>