【原】GBS hapmap 格式轉化為Plink格式方法

育種數(shù)據(jù)分析 2021-11-18

展開全文

1.需求說明

進行重測序或者GBS時，hapmap 是比較常見的格式，生信中經(jīng)常使用這種格式。但是在GWAS和GS中，數(shù)據(jù)篩選，質控，構建矩陣都是使用的plink的格式。本文介紹如何tassel 和vcftools兩個軟件，將hapmap格式的數(shù)據(jù)轉化為plink格式的數(shù)據(jù)。

環(huán)境：linux系統(tǒng)

2. 安裝軟件

首先，要安裝anaconda或者miniconda，然后使用conda install進行軟件安裝, 安裝conda方法：https://docs./anaconda/install/linux/

2.1 安裝tassel

tassel的按照方法，使用git將文件copy到本地，然后將里面的內容（可執(zhí)行文件） copy到home下的bin文件中，不用設置路徑了。

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

2.2 安裝vcftools

https:///bioconda/vcftools

conda install -c bioconda vcftools

2.3 安裝R語言

https:///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格式：`ped`和`map`

plink格式是基因組選擇中經(jīng)常用到的文件類型， plink軟件功能強大，運行速度快。

3.2.1 `.map格式`

格式說明鏈接: http://zzz.bwh./plink/data.shtml#map

map格式的文件, 主要是圖譜文件信息, 主要包括染色體名稱, 所在的染色體和所在染色體的坐標.

1, map文件沒有行頭

2, map文件包括四列: 染色體, SNP名稱, SNP位置, 堿基對坐標

染色體編號為數(shù)字, 未知為0
SNP名稱為字符或數(shù)字, 如果不重要, 可以從1編號, 注意要和bed文件SNP列一一對應
染色體的摩爾未知(可選項, 可以用0)
SNP物理坐標

3, 如果只有SNP名稱, 可以手動構建map文件, 第二列為SNP名稱, 其它三列為0即可.

Example:

1 snp1 0 11 snp2 0 21 snp3 0 3

這里有3個SNP, 分別名為snp1, snp3, snp3 (第二列)
這三個SNP在第一個染色體上 (第一列)
第三列為0
第四列為SNP所在染色體的坐標

3.2.2 `.ped格式`

格式說明鏈接:http://zzz.bwh./plink/data.shtml#ped

bed格式的文件, 主要包括SNP的信息, 包括個體ID, 系譜信息, 表型和SNP的分型信息.

1, 數(shù)據(jù)沒有行頭, 空格或者tab隔開的文件 2, 必須要有六列, 包括系譜信息, 表型信息

第一列: Family ID # 如果沒有, 可以用個體ID代替
第二列: Individual ID # 個體ID編號
第三列: Paternal ID # 父本編號
第四列: Maternal ID # 母本編號
第五列: Sex (1=male; 2=female; other=unknown) # 性別, 如果未知, 用0表示
第六列: Phenotype # 表型數(shù)據(jù), 如果未知, 用0表示
第七列以后: 為SNP分型數(shù)據(jù), 可以是AT CG或11 12, 或者A T C G或1 1 2 2

3, 上面六列, 必須要有, 如果沒有相關數(shù)據(jù), 用0表示.

Example:

1 1 0 0 1 0 G G 2 2 C C1 2 0 0 2 0 A A 0 0 A C1 3 1 2 1 2 0 0 1 2 A C2 1 0 0 1 0 A A 2 2 0 02 2 0 0 2 2 A A 2 2 0 02 3 1 2 1 2 A A 2 2 A A

數(shù)據(jù)包括兩個家系 (第一列)
每個家系有三個個體 (第二列)
第三列父本編號
第四列母本編號
第五列性別
第六列表型值
第七列, 第八列為一個基因型
第九列, 第十列為第二個基因型
第十一列, 第十二列為第三個基因型

4. 測試數(shù)據(jù)

將下面文件保存為：hmp.txt 注意，知乎中是.，但數(shù)據(jù)應該是#，下面這個代碼是正確的。

rs#alleleschromposstrandassembly#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.21

10000235A/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文件，變?yōu)閎ed文件

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 11000697401441878911278463601607734371127846370160773437112784638016077343711278463901607734372100069860764162462100071170150034851612784632028132948612784633028132948612784634028132948612784635028132948610000235080388997610007153014525696071000457503279275581000034501786588591000711301278523201010007668019391507

另外，可以編寫R代碼，提取map文件，將代碼轉化，然后轉置，但是效率較低。