用到的软件包：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
<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>
调用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>
