作者： 劉峻宇
1、BLUPF90
首先講一下BLUPF90家族，在這個家族中包括將近20種小程序，構(gòu)建H矩陣需要用到renumf90.exe和preGSf90.exe兩個程序。Windows系統(tǒng)下按照正常的軟件是不能打開的，需要在cmd控制臺上后臺運行（開始→搜索cmd→回車），你會看到這個界面：
百度搜索會有很多命令，這里只需要一些簡單的：
cd + 空格 +下一層子目錄（cd E:\software\blupf90）→表示進入下一層的子目錄
cd + 空格 + ..（cd ..）→表示回上一層目錄
cd加空格加/ （cd /）→回到根目錄
之所以要學這些指令是因為我們需要把控制臺的路徑設(shè)置與文件所在的路徑一致，這樣軟件才能讀取文件
3、準備文件
一共需要4個文件：表型，系譜，snp，參數(shù)文件（txt文件）
表型和系譜的數(shù)據(jù)不再多說，格式上有一點，不需要表頭（列名），第一行就是數(shù)據(jù)
snp文件：格式為0,1,2,5的形式，NA設(shè)置為5，形式如下：
構(gòu)建參數(shù)文件，表型，系譜和snp數(shù)據(jù)都在參數(shù)文件中呈現(xiàn)，格式如下：
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、構(gòu)建H矩陣
（1）在cmd中隨時提取blupf90中的程序，需要設(shè)置一下環(huán)境變量，將blupf90中的程序所在的路徑添加在電腦環(huán)境變量中；或者直接將blupf90程序與所保存的準備文件放在同一個文件夾，但是這樣比較麻煩。
環(huán)境變量自己可以百度設(shè)置，不再多說；假設(shè)文件所在的路徑為：E:\ssGBLUP\BLUPF90，在cmd中的路徑設(shè)置：
E：
cd ssGBLUP\BLUPF90
OK
在已經(jīng)設(shè)置好環(huán)境變量和將準備文件所在的路徑設(shè)置好后,運行renumf90.exe程序
輸入"renumf90.exe"，回車，會讓你輸入?yún)?shù)文件，輸入剛才保存的參數(shù)文件
（2）待程序跑完后，你會發(fā)現(xiàn)原文件夾中有一個renf90.par的文件，打開如下：
紅色方框中的是需要手動添加的，保存H逆矩陣
再回到cmd中，運行pregsf90.exe，回車后，輸入renf90.par的文件
（3）最后原文件夾中會有一個Hinv.txt文件和一個renaddxx.ped的文件 。這兩個文件需要再次修改一下
renadd03.ped在Excel打開后， 按照第一列排序后，只留個體編號一列，另存為sortedAnimalID.txt
Hinv.txt 中是上三角矩陣，需要將它改為下三角矩陣，才能正確的在asreml中運行。如果有需要可以繼續(xù)求逆，成為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
