利用SHAPEIT将vcf文件进行基因型（genotype）定相（phasing）：查看两个突变是否来源于同一条链（染色体或父本或母本）

Posted 2020-11-19 橙子牛奶糖

首先，下载SHAPEIT.

按照里面的步骤安装完后，将vcf文件进行基因型定相，分四步走。

 

第一步，将vcf文件转化为plink二进制文件（.bed, .bim, .fam）。

这一步需要用到GATK里的GenomeAnalysisTK工具，见如下命令：

java -Xmx8g -jar GenomeAnalysisTK.jar -T VariantsToBinaryPed  -R GRCh37.fa -V file.vcf --metaData sampleID.fam -mgq 0 -bed file.bed -bim file.bim  -fam file.fam　

　注：这里的metaData，输入的fam文件跟plink的fam格式一样，如果你没有fam文件，则需要自己手动生成，第一列和第一列可以都为样本的ID;

FAM files

The fields in a FAM file are

• Family ID
• Sample ID
• Paternal ID
• Maternal ID
• Sex (1=male; 2=female; other=unknown)
• Affection (0=unknown; 1=unaffected; 2=affected)

Example of a FAM file of the binary PLINK format:
FAM1	NA06985	0	0	1	1
FAM1	NA06991	0	0	1	1
0	NA06993	0	0	1	1
0	NA06994	0	0	1	1
0	NA07000	0	0	2	1
0	NA07019	0	0	1	1
0	NA07022	0	0	2	1
0	NA07029	0	0	1	1
FAM2	NA07056	0	0	0	2
FAM2	NA07345	0	0	1	1

 

######号外，在这一步的时候，如果出现报错信息：ERROR MESSAGE: No metadata provided for sample XXXXXX，很有可能是fam文件在作怪，需要对fam文件进行样本排序，使其的顺序跟vcf文件保持一致；

第二步，去除高缺失率的基因型和样本

plink --noweb --bfile file --keep-allele-order --me 1 1 --set-me-missing --make-bed --out file.nomendel

 

plink --noweb --bfile file.nomendel --keep-allele-order --geno 0.05 --make-bed --out file.nomendel.filter

　　

　####对数据进行高缺失率去除是很有必要的，如果跳过这一步，直接定相的话，会出现：“ERROR: XXX fully missing individuals (=100%)” 报错　

第三步，对基因型进行定位，需要下载遗传图谱（genetic map）

shapeit --duohmm -W 5 --thread 1 --input-bed file.nomendel.filter.bed file.nomendel.filter.bim file.nomendel.filter.fam --input-map genetic_map_b37.txt -O hapData

第四步，将定相的文件生成vcf格式：

shapeit -convert --input-haps hapData --output-vcf hapData.vcf

　　至此，vcf中没有定相的变异会由“0/0 0/1 1/0 1/1 ”变为“0|0 0|1 1|0 1|1”