利用BLUPF90构建H矩阵

作者: 刘峻宇

1、BLUPF90
首先讲一下BLUPF90家族,在这个家族中包括将近20种小程序,构建H矩阵需要用到renumf90.exe和preGSf90.exe两个程序。Windows系统下按照正常的软件是不能打开的,需要在cmd控制台上后台运行(开始→搜索cmd→回车),你会看到这个界面:

CMD控制台
2、控制台指令
百度搜索会有很多命令,这里只需要一些简单的:
cd + 空格 +下一层子目录(cd E:\software\blupf90)→表示进入下一层的子目录
cd + 空格 + ..(cd ..)→表示回上一层目录
cd加空格加/ (cd /)→回到根目录
之所以要学这些指令是因为我们需要把控制台的路径设置与文件所在的路径一致,这样软件才能读取文件
3、准备文件
一共需要4个文件:表型,系谱,snp,参数文件(txt文件)
表型和系谱的数据不再多说,格式上有一点,不需要表头(列名),第一行就是数据
snp文件:格式为0,1,2,5的形式,NA设置为5,形式如下:

snp文件格式

构建参数文件,表型,系谱和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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 215,923评论 6 498
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 92,154评论 3 392
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 161,775评论 0 351
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,960评论 1 290
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,976评论 6 388
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,972评论 1 295
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,893评论 3 416
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,709评论 0 271
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 45,159评论 1 308
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,400评论 2 331
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,552评论 1 346
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 35,265评论 5 341
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,876评论 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,528评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,701评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,552评论 2 368
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,451评论 2 352

推荐阅读更多精彩内容