作者: 刘峻宇
1、BLUPF90
首先讲一下BLUPF90家族,在这个家族中包括将近20种小程序,构建H矩阵需要用到renumf90.exe和preGSf90.exe两个程序。Windows系统下按照正常的软件是不能打开的,需要在cmd控制台上后台运行(开始→搜索cmd→回车),你会看到这个界面:
百度搜索会有很多命令,这里只需要一些简单的:
cd + 空格 +下一层子目录(cd E:\software\blupf90)→表示进入下一层的子目录
cd + 空格 + ..(cd ..)→表示回上一层目录
cd加空格加/ (cd /)→回到根目录
之所以要学这些指令是因为我们需要把控制台的路径设置与文件所在的路径一致,这样软件才能读取文件
3、准备文件
一共需要4个文件:表型,系谱,snp,参数文件(txt文件)
表型和系谱的数据不再多说,格式上有一点,不需要表头(列名),第一行就是数据
snp文件:格式为0,1,2,5的形式,NA设置为5,形式如下:
构建参数文件,表型,系谱和snp数据都在参数文件中呈现,格式如下:
snp格式的源代码如下:
rm(list=ls())
gc()
library(data.table)
snps <- fread("SNP_loci_single_v2.csv",header = TRUE,sep = ",",stringsAsFactors = FALSE)
snps <- snps[,-1]
snps_012 <- snps[,-c("AnimalID")]
setNAas5 <- function(x){
x[is.na(x)] <- 5
return(x)}
snps_0125 <- apply(snps_012,2,setNAas5)
setrowcollapse <- function(x){
return(paste(x,sep="",collapse = ""))
}
snps_0125_collapse <- apply(snps_0125,1,setrowcollapse)
snps_new <-data.table(snps[,c("AnimalID")],snps_0125_collapse)
colnames(snps_new) <- c("ID","SNP012")
fwrite(snps_new,file = "SNP0125forBLUPF90_2year_402_single_loci_v2.txt",quote = FALSE,sep=" ",row.names=FALSE,col.names = FALSE)
4、构建H矩阵
(1)在cmd中随时提取blupf90中的程序,需要设置一下环境变量,将blupf90中的程序所在的路径添加在电脑环境变量中;或者直接将blupf90程序与所保存的准备文件放在同一个文件夹,但是这样比较麻烦。
环境变量自己可以百度设置,不再多说;假设文件所在的路径为:E:\ssGBLUP\BLUPF90,在cmd中的路径设置:
E:
cd ssGBLUP\BLUPF90
OK
在已经设置好环境变量和将准备文件所在的路径设置好后,运行renumf90.exe程序
输入"renumf90.exe",回车,会让你输入参数文件,输入刚才保存的参数文件
(2)待程序跑完后,你会发现原文件夹中有一个renf90.par的文件,打开如下:
红色方框中的是需要手动添加的,保存H逆矩阵
再回到cmd中,运行pregsf90.exe,回车后,输入renf90.par的文件
(3)最后原文件夹中会有一个Hinv.txt文件和一个renaddxx.ped的文件 。这两个文件需要再次修改一下
renadd03.ped在Excel打开后, 按照第一列排序后,只留个体编号一列,另存为sortedAnimalID.txt
Hinv.txt 中是上三角矩阵,需要将它改为下三角矩阵,才能正确的在asreml中运行。如果有需要可以继续求逆,成为n*n的矩阵形式。
下三角矩阵源代码:
h_sort<- fread(input = "Hinv.txt",sep = " ",header = FALSE, stringsAsFactors = FALSE)
colnames(h_sort) <- c("id1","id2","value")
h_sort[,id1:=as.integer(id1)]
h_sort[,id2:=as.integer(id2)]
h <- h_sort[,c(2,1,3)]
colnames(h) <- c("id1","id2","value")
h_sort <- setorder(h,id1,id2)
fwrite(x=h_sort,file = "hinv.giv",sep = " ",row.names = FALSE,col.names = FALSE,quote = FALSE,na = "0")
最后asreml中用到的两个文件:hinv.giv 和 sortedAnimalID.txt