专栏首页育种数据分析之放飞自我GBS hapmap 格式 转化为Plink格式方法

GBS hapmap 格式 转化为Plink格式方法

1.需求说明

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

行头:

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格式:pedmap

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	GG

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.pedresult.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

本文分享自微信公众号 - 育种数据分析之放飞自我(R-breeding),作者:邓飞2013

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-10-23

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • asreml 设定初始值 固定初始值

    一个朋友问我,如何固定asreml的初始值,现在分为单性状和多性状进行说明。 为何要固定初始值: 1,由于群体较小,估算的方差组分不准确,需要手动设定初始值,直...

    邓飞
  • 笔记 | GWAS 操作流程3:plink关联分析--完结篇

    注意,这里我使用的是ped和map格式,如果ped文件中有表型数据(第六列),如果想指定表型数据,用--pheno,包括三列:家系,个体,表型值。

    邓飞
  • 如何利用系谱进行家系划分并可视化?

    概念定义共祖系数:共祖系数为概率fAB,表示一个来自个体A,另一个来自个体B的两个同源基因(或等位基因)在系谱上是一致或相同的概率,也就是说来自同一祖先基因的概...

    邓飞
  • R语言naniar包(新名词:阴影矩阵;Shadow matrices)

    因为ggplot2不能处理缺失值,所以我们得到了一个warning message ,我们可以使用geom_miss_point() 去展示缺失数据。

    用户1359560
  • R语言中回归模型预测的不同类型置信区间应用比较分析

    我们正在这里做出一个预测。正如在R课堂上(以及在预测模型的过程中)所回顾的,当我们要为预测提供一个置信区间时,建议您为预测器确定置信区间(这将取决于预测误差)参...

    用户5031023
  • 数据处理第一节:选取列的基本到高级方法选取列列名

    博客原文:https://suzan.rbind.io/2018/01/dplyr-tutorial-1/ 作者:Suzan Baert

    用户1359560
  • 学徒讨论-在数据框里面使用每列的平均值替换NA

    他认为替换不干净,应该是循环有问题。希望我们帮忙检查,我通常是懒得看其他人写的代码,所以让群里的小伙伴们有空的都尝试写一下。

    生信技能树
  • R语言第二章数据处理⑥dplyr包(1)列选取目录选取列

    =========================================

    用户1359560
  • TidyFriday 每天 5 分钟,轻轻松松上手 R 语言(五)

    今天我们依旧利用 msleep 数据集来探讨 dplyr 的列筛选,并在最后补充几个行筛选的例子。

    王诗翔呀
  • RNA-seq(10):KEGG通路可视化:gage和pathview

    开始用gage包进行富集分析,gage()函数需要fold change 和Entrez gene IDs

    Y大宽

扫码关注云+社区

领取腾讯云代金券