GBS hapmap 格式 转化为Plink格式方法
GBS hapmap 格式 转化为Plink格式方法
邓飞
发布于 2019-10-24 19:34:00
发布于 2019-10-24 19:34:00
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 vcftools2.3 安装R语言
https://anaconda.org/r/r
conda install -c r r3. 文件格式
3.1 hapmap格式：genotype.hmp.txt
行头：
rs# alleles chrom pos strand 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:
1 snp1 0 1
1 snp2 0 2
1 snp3 0 3- 这里有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:
1 1 0 0 1 0 G G 2 2 C C
1 2 0 0 2 0 A A 0 0 A C
1 3 1 2 1 2 0 0 1 2 A C
2 1 0 0 1 0 A A 2 2 0 0
2 2 0 0 2 2 A A 2 2 0 0
2 3 1 2 1 2 A A 2 2 A A- 数据包括两个家系 (第一列)
- 每个家系有三个个体 (第二列)
- 第三列父本编号
- 第四列母本编号
- 第五列性别
- 第六列表型值
- 第七列, 第八列为一个基因型
- 第九列, 第十列为第二个基因型
- 第十一列, 第十二列为第三个基因型
4. 测试数据
将下面文件保存为：hmp.txt 注意， 知乎中是.，但数据应该是#，下面这个代码是正确的。
rs# alleles chrom pos strand assembly# center protLSID assayLSID panelLSID QCcode Box302.YX.2017.06.002 Box302.YX.2017.06.003 Box302.YX.2017.06.004 Box302.YX.2017.06.005 Box302.YX.2017.06.11 Box302.YX.2017.06.013 Box302.YX.2017.06.018 Box302.YX.2017.06.019 Box302.YX.2017.06.020 Box302.YX.2017.06.21
10000235 A/C 6 80388997 + NA NA NA NA NA NA AC AA AC AA AA AA AA AA AC AC
10000345 G/A 8 17865885 + NA NA NA NA NA NA GG GG AA GG GG GG GG AA AA GG
10004575 G/T 7 32792755 + NA NA NA NA NA NA GG GG GG GG GG GT GG GG GG GG
10006974 T/C 1 14418789 + NA NA NA NA NA NA CT CC TT CT CT CT CT TT TT CC
10006986 A/G 2 76416246 + NA NA NA NA NA NA AA AA AA AA AA AG AG GG GG AA
10007074 G/C 14 105043500 + NA NA NA NA NA NA GG GC GC GC GC GC GC GC GG GC
10007097 C/A 15 118863849 + NA NA NA NA NA NA CC CC CC CC CC CC AC CC CC AC
10007113 G/A 9 127852320 + NA NA NA NA NA NA GG GG GG AG GG GG AG AA AA GG
10007117 C/T 2 150034851 + NA NA NA NA NA NA CC CT CC CT CT CT TT CC CT CT
10007153 A/C 6 145256960 + NA NA NA NA NA NA AA AC AA AC AA AA AA AC AA AA
10007668 G/A 10 19391507 + NA NA NA NA NA NA GG GG GG AG AG AA GG GG GG GG
12784072 G/A 17 55229968 + NA NA NA NA NA NA AG AG AG GG AG GG GG AG AA GG
12784632 T/C 6 28132948 + NA NA NA NA NA NA TT TT TT TT TT TT TT CT TT TT
12784633 T/C 6 28132948 + NA NA NA NA NA NA TT TT TT TT TT TT TT CT TT TT
12784634 T/C 6 28132948 + NA NA NA NA NA NA TT TT TT TT TT TT TT CT TT TT
12784635 T/C 6 28132948 + NA NA NA NA NA NA TT TT TT TT TT TT TT CT TT TT
12784636 G/A 1 160773437 + NA NA NA NA NA NA AG GG GG GG AG AG GG AG AG GG
12784637 G/A 1 160773437 + NA NA NA NA NA NA AG GG GG GG AG AG GG AG AG GG
12784638 G/A 1 160773437 + NA NA NA NA NA NA AG GG GG GG AG AG GG AG AG GG
12784639 G/A 1 160773437 + NA NA NA NA NA NA AG GG GG GG AG AG GG AG AG GG5. 处理代码
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.ped
Box302.YX.2017.06.002 Box302.YX.2017.06.002 0 0 0 -9 C 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 G
Box302.YX.2017.06.003 Box302.YX.2017.06.003 0 0 0 -9 C 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 G
Box302.YX.2017.06.004 Box302.YX.2017.06.004 0 0 0 -9 T 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 G
Box302.YX.2017.06.005 Box302.YX.2017.06.005 0 0 0 -9 C 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 G
Box302.YX.2017.06.11 Box302.YX.2017.06.11 0 0 0 -9 C 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 G
Box302.YX.2017.06.013 Box302.YX.2017.06.013 0 0 0 -9 C 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 G
Box302.YX.2017.06.018 Box302.YX.2017.06.018 0 0 0 -9 C 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 G
Box302.YX.2017.06.019 Box302.YX.2017.06.019 0 0 0 -9 T 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 G
Box302.YX.2017.06.020 Box302.YX.2017.06.020 0 0 0 -9 T 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 A
Box302.YX.2017.06.21 Box302.YX.2017.06.21 0 0 0 -9 C 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 G G
(base) [dengfei@ny01 hmp_2_plink]$ cat result.map
1 10006974 0 14418789
1 12784636 0 160773437
1 12784637 0 160773437
1 12784638 0 160773437
1 12784639 0 160773437
2 10006986 0 76416246
2 10007117 0 150034851
6 12784632 0 28132948
6 12784633 0 28132948
6 12784634 0 28132948
6 12784635 0 28132948
6 10000235 0 80388997
6 10007153 0 145256960
7 10004575 0 32792755
8 10000345 0 17865885
9 10007113 0 127852320
10 10007668 0 19391507另外，可以编写R代码，提取map文件，将代码转化，然后转置，但是效率较低。
6. 参考资料
https://zhuanlan.zhihu.com/p/38590403
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2019-10-23，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读
目录