GATK变异检测

生信喵

1.变异检测

1. #生成gvcf  
   
2. time gatk HaplotypeCaller --emit-ref-confidence GVCF -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -I merge.sorted.markdup.BQSR.bam -O merge.HC.g.vcf.gz
   
3. 
   
4. #合并gvcf  
   
5. time gatk GenotypeGVCFs -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.g.vcf.gz -O merge.HC.vcf.gz

复制代码

2.结果过滤
2.1 VQSR
       准备的已知变异集作为训练集，可以是 Hapmap、OMNI，1000G，dbsnp，瓶中基因组计划等这些国际性项目的数据，然后利用训练集对每一个位点进行过滤。利用 VariantRecalibrator工具进行机器学习，ApplyVQSR 工具进行处理。 VQSR 过滤 SNP 和 InDel 分别进行，首先处理 SNP，得到结果后，在进行 InDel 处理。

1. #处理SNP  
   
2. # --max-gaussians默认值为8，报错提示需要降低
   
3. gatk VariantRecalibrator --max-gaussians 6 -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.vcf.gz --resource:hapmap,known=false,training=true,truth=true,prior=15.0 /share/home/xiehs/data/GATK/hg38/hapmap_3.3.hg38.vcf.gz --resource:omni,known=false,training=true,truth=false,prior=12.0 /share/home/xiehs/data/GATK/hg38/1000G_omni2.5.hg38.vcf.gz --resource:1000G,known=false,training=true,truth=false,prior=10.0 /share/home/xiehs/data/GATK/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /share/home/xiehs/data/GATK/hg38/dbsnp_146.hg38.vcf.gz -an QD -an MQ -an MQRankSum -an ReadPosRankSum -an FS -an SOR -mode SNP -O merge.HC.snps.recal --tranches-file output.tranches --rscript-file output.plots.R
   
4. 
   
5. 
   
6. gatk ApplyVQSR -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.vcf.gz -O merge.HC.snps.VQSR.vcf.gz --recal-file merge.HC.snps.recal --tranches-file merge.HC.snps.tranches -mode SNP

复制代码

      1、HapMap，它来自国际人类单倍体型图计划，数据集包含了大量家系数据，并且有非常严格的质控和严密的实验验证，因此它的准确性是目前公认最高的。
       2、Omni，这个数据源自 Illumina 的 Omni 基因型芯片，它的验证结果常常作为基因型的金标准。
       3、1000G 千人基因组计划（1000 genomes project）质控后的变异数据，质控后，它包含的绝大部分都是真实的变异，但由于没办法做全面的实验验证，并不能排除含有少部分假阳的结果。
       4、dbSNP。dbSNP 收集的数据，实际都是研究者们发表了相关文章提交上来的变异，这些变异很多是没做过严格验证的
处理 InDel

1. #处理InDel  
   
2. gatk VariantRecalibrator -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.snps.VQSR.vcf.gz --max-gaussians 4 --resource:mills,known=false,training=true,truth=true,prior=12.0 /share/home/xiehs/data/GATK/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz --resource:dbsnp,known=true,training=false,truth=false,prior=2.0 /share/home/xiehs/data/GATK/hg38/dbsnp_146.hg38.vcf.gz -an QD -an DP -an FS -an SOR -an ReadPosRankSum -an MQRankSum -mode INDEL -O merge.HC.snps.indel.recal --tranches-file merge.HC.snps.indel.tranches --rscript-file merge.HC.snps.indel.plots.R
   
3. 
   
4. gatk ApplyVQSR -R /share/home/xiehs/data/GATK/hg38/Homo_sapiens_assembly38.fasta -V merge.HC.snps.VQSR.vcf.gz -O merge.HC.snps.indel.VQSR.vcf.gz --truth-sensitivity-filter-level 99.0 --tranches-file merge.HC.snps.indel.tranches --recal-file merge.HC.snps.indel.recal -mode INDEL

复制代码

2.2 Hard-filter
       Hard-filter 硬过滤，可以根据以下标准来进行过滤，gatk 过滤的时候，snp 与 indel 是分别进行的。也可以选择 bcftools 进行简单过滤。
       QualByDepth（QD）
       FisherStrand (FS)
       StrandOddsRatio (SOR)
       RMSMappingQuality (MQ)
       MappingQualityRankSumTest (MQRankSum)
       ReadPosRankSumTest (ReadPosRankSum)
       此脚本为硬过滤(hard-filter)的方法，主要用于不能进行 VQSR 的情况，例如非人物种，或外显子，芯片数据等。

1. # 使用SelectVariants，选出SNP  
   
2. time gatk SelectVariants -select-type SNP -V merge.HC.vcf.gz -O merge.HC.vcf.snp.gz
   
3.   
   
4. # 为SNP作硬过滤  
   
5. time gatk VariantFiltration -V merge.HC.vcf.snp.gz --filter-expression "QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O merge.HC.vcf.snp.filter.gz
   
6.   
   
7. # 使用SelectVariants，选出Indel  
   
8. time gatk SelectVariants -select-type INDEL -V merge.HC.vcf.gz -O merge.HC.indel.vcf.gz
   
9.   
   
10. # 为Indel作过滤  
    
11. time gatk VariantFiltration -V merge.HC.indel.vcf.gz --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name "PASS" -O merge.HC.indel.filter.vcf.gz
    
12.   
    
13. # 重新合并过滤后的SNP和Indel  
    
14. time gatk MergeVcfs -I merge.HC.snp.filter.vcf.gz -I merge.HC.indel.filter.vcf.gz -O merge.HC.filter.vcf.gz

复制代码