hapmap 格式 转化为Plink格式方法

1.需求说明

进行重测序或者GBS时,hapmap是比较常见的格式,生信中经常使用这种格式。但是在GWAS和GS中,数据筛选,质控,构建矩阵都是使用的plink的格式。本文介绍如何tassel和vcftools两个软件,将hapmap格式的数据转化为plink格式的数据。

环境:linux系统

2. 安装软件

首先,要安装anaconda或者miniconda,然后使用conda install进行软件安装, 安装conda方法:https://docs.anaconda.com/anaconda/install/linux/

2.1 安装tassel

tassel的按照方法,使用git将文件copy到本地,然后将里面的内容(可执行文件) copy到home下的bin文件中, 不用设置路径了。

git clone https://bitbucket.org/tasseladmin/tassel-5-standalone.git

2.2 安装vcftools

https://anaconda.org/bioconda/vcftools

conda install -c bioconda vcftools

2.3 安装R语言

https://anaconda.org/r/r

conda install -c r r

3. 文件格式

3.1  hapmap格式:genotype.hmp.txt

行头:

assembly# center protLSID assayLSID panelLSID QCcode sample1 sample2 ...

内容:

rs# alleles chrom pos strand assembly# center protLSID assayLSID panelLSID QCcode Sample_YCX334/12Sample_ya>1:1151 C 1 1151 + NA NA NA NA NA NA N N N C C C N C>1:1203 T/C 1 1203 + NA NA NA NA NA NA T N T N T T T T>1:1249 A/C 1 1249 + NA NA NA NA NA NA A N A N A A A A>1:1266 G/A 1 1266 + NA NA NA NA NA NA G N G G G G G N>1:1277 T/C 1 1277 + NA NA NA NA NA NA T T T T T T T N>1:1325 C/T 1 1325 + NA NA NA NA NA NA C N N N N N C N>1:1335 G/T 1 1335 + NA NA NA NA NA NA G G G G G G G G>1:1362 G/A 1 1362 + NA NA NA NA NA NA G G G G G G G G>

3.2 plink格式:ped和map

plink格式是基因组选择中经常用到的文件类型, plink软件功能强大,运行速度快。

3.2.1.map格式

格式说明链接: http://zzz.bwh.harvard.edu/plink/data.shtml#map

map格式的文件, 主要是图谱文件信息, 主要包括染色体名称, 所在的染色体和所在染色体的坐标.

1, map文件没有行头

2, map文件包括四列: 染色体, SNP名称, SNP位置,  碱基对坐标

染色体编号为数字, 未知为0

SNP名称为字符或数字, 如果不重要, 可以从1编号, 注意要和bed文件SNP列一一对应

染色体的摩尔未知(可选项, 可以用0)

SNP物理坐标

3, 如果只有SNP名称, 可以手动构建map文件, 第二列为SNP名称, 其它三列为0即可.

Example:

1snp1011snp2021snp303

这里有3个SNP, 分别名为snp1, snp3, snp3(第二列)

这三个SNP在第一个染色体上(第一列)

第三列为0

第四列为SNP所在染色体的坐标

3.2.2.ped格式

格式说明链接:http://zzz.bwh.harvard.edu/plink/data.shtml#ped

bed格式的文件, 主要包括SNP的信息, 包括个体ID, 系谱信息, 表型和SNP的分型信息.

1, 数据没有行头, 空格或者tab隔开的文件 2, 必须要有六列, 包括系谱信息, 表型信息

第一列: Family ID # 如果没有, 可以用个体ID代替

第二列: Individual ID # 个体ID编号

第三列:  Paternal ID # 父本编号

第四列:  Maternal ID # 母本编号

第五列:  Sex (1=male; 2=female; other=unknown) # 性别, 如果未知, 用0表示

第六列:  Phenotype # 表型数据, 如果未知, 用0表示

第七列以后: 为SNP分型数据, 可以是AT CG或11 12, 或者A T C G或1 1 2 2

3, 上面六列, 必须要有, 如果没有相关数据, 用0表示.

Example:

110010G G22C C120020A A00A C1312120012A C210010A A2200220022A A2200231212A A22A A

数据包括两个家系(第一列)

每个家系有三个个体(第二列)

第三列父本编号

第四列母本编号

第五列性别

第六列表型值

第七列, 第八列为一个基因型

第九列, 第十列为第二个基因型

第十一列, 第十二列为第三个基因型

4. 测试数据

将下面文件保存为:hmp.txt 注意, 知乎中是.,但数据应该是#,下面这个代码是正确的。

rs#  alleleschrom pos  strandassembly#centerprotLSIDassayLSIDpanelLSIDQCcodeBox302.YX.2017.06.002Box302.YX.2017.06.003Box302.YX.2017.06.004Box302.YX.2017.06.005Box302.YX.2017.06.11Box302.YX.2017.06.013Box302.YX.2017.06.018Box302.YX.2017.06.019Box302.YX.2017.06.020Box302.YX.2017.06.2110000235A/C680388997+NANANANANANAACAAACAAAAAAAAAAACAC10000345G/A817865885+NANANANANANAGGGGAAGGGGGGGGAAAAGG10004575G/T732792755+NANANANANANAGGGGGGGGGGGTGGGGGGGG10006974T/C114418789+NANANANANANACTCCTTCTCTCTCTTTTTCC10006986A/G276416246+NANANANANANAAAAAAAAAAAAGAGGGGGAA10007074G/C14105043500+NANANANANANAGGGCGCGCGCGCGCGCGGGC10007097C/A15118863849+NANANANANANACCCCCCCCCCCCACCCCCAC10007113G/A9127852320+NANANANANANAGGGGGGAGGGGGAGAAAAGG10007117C/T2150034851+NANANANANANACCCTCCCTCTCTTTCCCTCT10007153A/C6145256960+NANANANANANAAAACAAACAAAAAAACAAAA10007668G/A1019391507+NANANANANANAGGGGGGAGAGAAGGGGGGGG12784072G/A1755229968+NANANANANANAAGAGAGGGAGGGGGAGAAGG12784632T/C628132948+NANANANANANATTTTTTTTTTTTTTCTTTTT12784633T/C628132948+NANANANANANATTTTTTTTTTTTTTCTTTTT12784634T/C628132948+NANANANANANATTTTTTTTTTTTTTCTTTTT12784635T/C628132948+NANANANANANATTTTTTTTTTTTTTCTTTTT12784636G/A1160773437+NANANANANANAAGGGGGGGAGAGGGAGAGGG12784637G/A1160773437+NANANANANANAAGGGGGGGAGAGGGAGAGGG12784638G/A1160773437+NANANANANANAAGGGGGGGAGAGGGAGAGGG12784639G/A1160773437+NANANANANANAAGGGGGGGAGAGGGAGAGGG

5. 处理代码

5.1 排序

run_pipeline.pl-SortGenotypeFilePlugin-inputFile hmp.txt-outputFile test.sort.hmp.txt-fileType Hapmap

生成一个test.sort.hmp.txt的

5.2 生成 vcf4.0 格式的文件

run_pipeline.pl-fork1-h test.sort.hmp.txt-export-exportType VCF

生成一个test.vcf文件

5.3 使用vcftools生成plink文件

vcftools--vcf test.vcf--plink--out  tassel.test.vcf2plink

生成tassel.test.vcf2plink文件

5.4 使用plink将vcf文件, 变为bed文件

plink--file tassel.test.vcf2plink--make-bed--out tassel.test.vcf2plink

5.5 使用plink将bed文件转化为map和ped文件

plink--bfile tassel.test.vcf2plink--recode--out result

结果生成:result.ped和result.map文件:

(base)[dengfei@ny01 hmp_2_plink]$ cat result.pedBox302.YX.2017.06.002Box302.YX.2017.06.002000-9C T A G A G A G A G A A C C T T T T T T T T C A A A G G G G G G G G G G C C A GBox302.YX.2017.06.003Box302.YX.2017.06.003000-9C C G G G G G G G G A A T C T T T T T T T T A A C A G G G G G G G G C G C C A GBox302.YX.2017.06.004Box302.YX.2017.06.004000-9T T G G G G G G G G A A C C T T T T T T T T C A A A G G A A G G G G C G C C A GBox302.YX.2017.06.005Box302.YX.2017.06.005000-9C T G G G G G G G G A A T C T T T T T T T T A A C A G G G G A G A G C G C C G GBox302.YX.2017.06.11Box302.YX.2017.06.11000-9C T A G A G A G A G A A T C T T T T T T T T A A A A G G G G G G A G C G C C A GBox302.YX.2017.06.013Box302.YX.2017.06.013000-9C T A G A G A G A G G A T C T T T T T T T T A A A A T G G G G G A A C G C C G GBox302.YX.2017.06.018Box302.YX.2017.06.018000-9C T G G G G G G G G G A T T T T T T T T T T A A A A G G G G A G G G C G A C G GBox302.YX.2017.06.019Box302.YX.2017.06.019000-9T T A G A G A G A G G G C C C T C T C T C T A A C A G G A A A A G G C G C C A GBox302.YX.2017.06.020Box302.YX.2017.06.020000-9T T A G A G A G A G G G T C T T T T T T T T C A A A G G A A A A G G G G C C A ABox302.YX.2017.06.21Box302.YX.2017.06.21000-9C C G G G G G G G G A A T C T T T T T T T T C A A A G G G G G G G G C G A C GG(base)[dengfei@ny01 hmp_2_plink]$ cat result.map11000697401441878911278463601607734371127846370160773437112784638016077343711278463901607734372100069860764162462100071170150034851612784632028132948612784633028132948612784634028132948612784635028132948610000235080388997610007153014525696071000457503279275581000034501786588591000711301278523201010007668019391507

另外,可以编写R代码,提取map文件,将代码转化,然后转置,但是效率较低。

特别注意:在执行上述步骤之前要对数据进行质控,确保原始数据内不包含错误行,数据行的分隔符不要错误。

6. 参考资料

https://zhuanlan.zhihu.com/p/38590403

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

推荐阅读更多精彩内容