GATK练习题

~ admin

因为本教程年代有些久远，无法保证其正确性，只是用于交流讨论。

教程适用于GATK v4， release 2.0， 源自http://www.broadinstitute.org/gatk/guide/topic?name=best-practices

流程图：

GATK SNP call流程图
GATK SNP call流程图

在SNP分析中，最基本的流程分为连续的多个步骤，初始mapping，精练初始短序，多样品插入删除及SNP分析，较准事件评分。这个过程也适用于再测序，全外显子组分析，深度基因组测序分析，及低通全基因组分析。

要完成这个练习，需要准备：
1.环境：基于UNIX/Linux计算机（包括MAC OS），了解和掌握基本的Unix命令，熟悉字符操作界面，会翻译运行二进制程序，了解什么是PATH，安装有C/C++编译库（mac机需要安装Xcode, Linux/Unix会预装编译环境）。
2.：

3.材料：已经mapping好的NA12878 BAM文件及其indexed文件。其它相关材料，比如基因组，dbSNP等等。

阶段I 准备文件
第一步：使用BWA完成fastq原文件至基因组的映射（mapping）。（略）

第二步：使用Picard的MarkDuplicate工具去除重复，使用GATK的realignment工具重新完成local alignment工作，使用GATK的recalibration工具较准评分。这一步需要运行两次。第一次，对每个测序lane的数据先做一次。而后第二次对每个样品的数据集合在一起再做一遍。

去除重复的作用是希望去除在测序过程中产生的一个片断被多次测序的影响。

java -jar MarkDuplicates.jar \
     INPUT=NA12878_2x250_PCRfree_20.bam \
     OUTPUT=original.bam \
     METRICS_FILE=picard_info.txt \
     REMOVE_DUPLICATES=true \
     ASSUME_SORTED=true \
     VALIDATION_STRINGENCY=LENIENT

那为什么要做local realignment工作呢？假设以下几种情况：
a.发生了insertion或者deletion。
b.如果有连续出现的均聚合物（比如说连续的多个T，或者A等），在测序结果中7个T有可能被记为6个或者8个。这时就会影响最终的计算。如下图：
GATKwh0-BP-2-Realignment_Page_05
在经过local relaignment之后，通过插入空隔，将这种影响消除。
GATKwh0-BP-2-Realignment_Page_06
那么问题就产生了，如何判断这是一个正常的insertion或者deletion，而不是测序错误呢？GATK通过评分的办法来判定，评分的基准有三：
i）已知的SNP库，比如dbSNP, 1K Genomes等。ii）上述情况中的a，以初始比对的过程中就被标定为insertion或者deletion的（记录在SAM文件中CIGAR标记中）。iii）测序结果显示极有可能是一个indel。
从上述基准中我们可以看出，local realignment需要准备好已知的SNP库。在操作上，GATK的local realignment被分为两步走：
A）使用RealignerTargetCreator标记哪些区域需要local realignment。

## find out the contigs shared in bam file and reference. If the reference is same for mapping, skip this step and -XL parameter in next step 
## -L and -XL are opposite arguments GATK, the later one is to exclude any contig not exist in reference. 
samtools view -H original.bam | grep "^@SQ"
cat human_g1k_v37.fasta | grep "^>"
## remove NC_007605 and hs37d5 from original.bam
samtools view -h original.bam | grep -v "NC_007605" | grep -v "hs37d5" | samtools view -bS - > new.bam
## sort the bam file by reference
## index the bam file
samtools index new.bam
 
java -jar GenomeAnalysisTK.jar -T RealignerTargetCreator \
                               -R human_g1K_v37.fasta \
                               -I new.bam \
                               --known dbsnp_137.b37.vcf \
                               --known 1000G_omni2.5.b37.vcf \
                               --known 1000G_phase1.indels.b37.vcf \
                               --known 1000G_phase1.snps.high_confidence.b37.vcf \
                               -o realigner.intervals

如果只对已知的indel感兴趣，这一步时可以不输入BAM文件。在这一步时，可以使用已知的SNP库，这样会大大提高运算效率以及准确性。输出文件使用-o参数来确定。更多参数。
B) 使用IndelRealigner实现local realignment。

java -Xmx4g -jar GenomeAnalysisTK.jar -T IndelRealigner \
                               -R human_g1K_v37.fasta \
                               -I new.bam \
                               -known dbsnp_137.b37.vcf \
                               -known 1000G_omni2.5.b37.vcf \
                               -known 1000G_phase1.indels.b37.vcf \
                               -known 1000G_phase1.snps.high_confidence.b37.vcf \
                               --targetIntervals realigner.intervals \
                               -o realigned.bam

上一步的结果在这一步中通过–targetIntervals参数传入。这一步需要较多的内存，通过Xmx4g设定java运行环境可以使用4G的内存。在可选参数中，–consensusDeterminationModel较为重要，使用KNOWNS_ONLY模式运行速度快，但只计算已知的SNP；默认使用USE_READS,对在原始比对中已经标记为indel的也会包括进来；USE_SW模式现在看来还没有实验支持其可用性。更多参数。

因为每个碱基的测序得分（Quality scores）对于下游的分析是至关重要的，然后测序结果直接给出的得分却有倾向性，所以必须通过较准的步骤。如下图所示。
GATKwh0-BP-3-Base_recalibration_Page_05
在操作上GATK的recalibration被分为两步走：较正（BaseRecalibrator）以及输出（PrintReads）。之所以分为两步，是因为可以通过连续两次调用BaseRecalibrator来生成一个recalibration前后比对图，从而实现输出可视化效果及输出结果兼顾。
GATKwh0-BP-3-Base_recalibration_Page_14
首先建立较正模型，然后应用较正模型到数据上。更多参数。

java -Xmx4g -jar GenomeAnalysisTK.jar -T BaseRecalibrator \
                                      -R human_g1K_v37.fasta \
                                      -I realigned.bam \
                                      --knownSites dbsnp_137.b37.vcf \
                                      --knownSites 1000G_omni2.5.b37.vcf \
                                      --knownSites 1000G_phase1.indels.b37.vcf \
                                      --knownSites 1000G_phase1.snps.high_confidence.b37.vcf \
                                      -o recal.grp

而后使用PrintReads将其结果输出至新的BAM文件中。

java -Xmx2g -jar GenomeAnalysisTK.jar -T PrintReads \
                                      -R human_g1K_v37.fasta \
                                      -I realigned.bam \
                                      --BQSR recal.grp \
                                      -o recal.bam

这一步生成的BAM文件将会得到新的测序评分。旧的评分会保留在OQ标签下。更多参数。

这一步过程中，如果希望将多个bam文件合并成一个文件，可以使用picard的MergeSamFiles工具。如果想将一个bam文件分成多个bam文件，可以使用GATK的SplitSamFile工具。
因为MarkDuplicate只与测序本身有关，如果你的样品只在一个lane中，只需要运行一次即可。如果你的一个样品在多个lane中，后来需要合并成一个样品，则需要运行两次。
对于多样品的realignment和recalibration,建议分开进行，因为这样在时效上更有利。但是，如果是成对的样品比如tumor/normal这种的，则最好同时运行，以避免不同次运行的误差。

第三步：使用GATK的ReduceReads工具缩小文件。
因为测序的结果文件都很大，即使比对结果的BAM文件也会少则几G多则几十G，在传输上非常不便。为此，GATK提供了ReduceReads工具来对BAM文件进行有损压缩，压缩后的文件甚至可以小至压缩前的1/100大小。它对于GATK的下游分析是有利的（因为文件小，读写上会有优势），但是并不保证其可以被其它工具正确识别和使用。所以经压缩后的文件并不适合用于存档。压缩后的文件样子如下图所示。
GATKwh0-BP-4-Compression_Page_09

java -Xmx4g -jar GenomeAnalysisTK.jar -T ReduceReads \
                                      -R human_g1K_v37.fasta \
                                      -I recal.bam \
                                      -o reduced.bam

阶段II 突变分析及基因分型：SNP和indel分析
经过了阶段I，我们将得到一个或者多个BAM文件。每一个BAM文件都是代表着一个样品，一个样品不应该有多个文件。每一个BAM文件都是经过realignment和recalibration的，并且可能的话经过压缩的。对于样品，如果测序深度在10X以上的话，推荐的突变分析的域值为Q30，如果小于10X的话，样品多于100可以使用Q4，等于或者小于100使用Q10。
对于突变分析，基本的问题就是怎么判断一个突变是测序误差（random machine noise）还是真实的突变，这其实就是一个large scale Bayesian inference problem。在GATK中提供了两种方法，一是不基于模型的UnifiedGenotyper，它速度非常快，不考虑相邻碱基的影响。一是基于体外拼接模型的HaplotypeCaller。后者是GATK推荐的方式，它还支持2倍体模式。虽然它不支持多线程模式，但是可以使用Queue来执行HaplotypeCaller。

java -Xmx4g -jar GenomeAnalysisTK.jar -T HaplotypeCaller \
                                      -R human_g1K_v37.fasta \
                                      -I reduced.bam \
                                      --dbsnp dbsnp_137.b37.vcf \
                                      -stand_call_conf 30.0 \
                                      -stand_emit_conf 30.0 \
                                      --minPruning 3 \
                                      -o output.vcf

更多参数。

之后，我们可以为SNP和INDEL分析结果加入注释。这里需要使用的是snpEff以及VariantAnnotator。
下载并安装VariantAnnotator支持的snpEff。详细说明。

java -jar snpEff.jar download GRCh37.64
java -Xmx4g -jar snpEff.jar eff -v \
     -i vcf -o gatk GRCh37.64 output.vcf > output.snpEff.vcf

之后使用VariantAnnotator.

java -jar GenomeAnalysisTK.jar -T VariantAnnotator \
     -R human_g1k_v37.fasta \
     -A SnpEff \
     --variant output.vcf \
     --snpEffFile output.snpEff.vcf \
     -o annotated.vcf

阶段III 报告结果
第一步：较正
对于突变来说，之前的实验结果表明，多个突变会聚集在一起（cluster），这些聚集区是随机分布（Gaussianly distributed）。所以对于每一个突变，都需要依着这一原则对其重新评分。
我们先来看一看VCF文件是如何记录一个突变的。
GATKwh0-BP-6-Variant_recalibration_Page_07
对突变的recalibration，在GATK中分为两步实现：构建高斯混合模型（VariantRecalibrator），使用模型重写VCF文件（ApplyRecalibration）。

java -Xmx4g -jar GenomeAnalysisTK.jar \
   -T VariantRecalibrator \
   -R human_g1K_v37.fasta \
   -input annotated.vcf \
   -resource:hapmap,known=false,training=true,truth=true,prior=15.0 hapmap_3.3.b37.vcf \
   -resource:omni,known=false,training=true,truth=false,prior=12.0 1000G_omni2.5.b37.vcf \
   -resource:dbsnp,known=true,training=false,truth=false,prior=6.0 dbsnp_137.b37.vcf \
   -an QD -an HaplotypeScore -an MQRankSum -an ReadPosRankSum -an FS -an MQ -an InbreedingCoeff \
   -mode SNP \
   -recalFile output.recal \
   -tranchesFile output.tranches \
   -rscriptFile output.plots.R

之后重建一个新的VCF文件

java -Xmx3g -jar GenomeAnalysisTK.jar \
   -T ApplyRecalibration \
   -R human_g1K_v37.fasta \
   -input annotated.vcf \
   --ts_filter_level 99.0 \
   -tranchesFile output.tranches \
   -recalFile output.recal \
   -mode SNP \
   -o output.recalibrated.filtered.vcf

对于输出文件，需要打开它确认新文件已经写好。
同时VQSR会生成一个图像来描述recalibration的效果。

第二步：精练
对基因分型来说，如果有亲代及子代的关系数据，需要使用PhaseByTransmission来精练突变结果。关系文件描述可从http://pngu.mgh.harvard.edu/~purcell/plink/data.shtml#ped找到。

 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -T PhaseByTransmission \
   -R human_g1K_v37.fasta \
   -V output.recalibrated.filtered.vcf \
   -ped input.ped \
   -o tranmission.phased.vcf

无论是否经过PhaseByTransmission的精练，结果都需要使用ReadBackedPahsing来refine。

 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -T ReadBackedPhasing \
   -R human_g1K_v37.fasta \
   -I reduced.bam \
   -V tranmission.phased.vcf \
   --phaseQualityThresh 20.0 \
   -respectPhasInput \
   -o phased.vcf

在此之后通过使用统计学手段，来对结果进一步精练，通常使用Beagle。为此，GATK提供了从VCF至Beagle及从beagle至VCF的文件转换工具ProduceBeagleInput及BeagleOutputToVCF

java -Xmx2g -jar dist/GenomeAnalysisTK.jar -T ProduceBeagleInput\
      -R human_g1K_v37.fasta \
      -V phased.vcf \
      -o path_to_beagle_output/beagle_output
java -Xmx4g -jar dist/GenomeAnalysisTK.jar \
      -R human_g1K_v37.fasta -T BeagleOutputToVCF \
      -V phased.vcf \
      -beagleR2:BEAGLE /myrun.beagle_output.r2 \
      -beaglePhased:BEAGLE /myrun.beagle_output.phased \
      -beagleProbs:BEAGLE /myrun.beagle_output.gprobs \
      -o output_vcf.vcf

经过了以上几步精练，其结果变化为：
GATKwh0-BP-7A-Genotype_refinement_Page_16

第三步：评价及过滤
如果我有两个样品，那么如何得到它们之间的差异突变呢？测序对误差有何影响呢？这些突变会影响激素水平吗？那些突变是个体特异的，还是有家族性的？这些都是在这一步需要解决的问题。
合并两个VCF文件。说明

java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T CombineVariants \
   --variant input1.vcf \
   --variant input2.vcf \
   -o output.vcf \
   -genotypeMergeOptions UNIQUIFY
 
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T CombineVariants \
   --variant:foo input1.vcf \
   --variant:bar input2.vcf \
   -o output.vcf \
   -genotypeMergeOptions PRIORITIZE
   -priority foo,bar

从一个VCF文件中提取属于某个样品的数据。说明

Select two samples out of a VCF with many samples:
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant input.vcf \
   -o output.vcf \
   -sn SAMPLE_A_PARC \
   -sn SAMPLE_B_ACTG
 
 Select two samples and any sample that matches a regular expression:
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant input.vcf \
   -o output.vcf \
   -sn SAMPLE_1_PARC \
   -sn SAMPLE_1_ACTG \
   -se 'SAMPLE.+PARC'
 
 Select any sample that matches a regular expression and sites where the QD annotation is more than 10:
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant input.vcf \
   -o output.vcf \
   -se 'SAMPLE.+PARC'
   -select "QD > 10.0"
 
 Select a sample and exclude non-variant loci and filtered loci:
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant input.vcf \
   -o output.vcf \
   -sn SAMPLE_1_ACTG \
   -env \
   -ef
 
 Select a sample and restrict the output vcf to a set of intervals:
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant input.vcf \
   -o output.vcf \
   -L /path/to/my.interval_list \
   -sn SAMPLE_1_ACTG
 
 Select all calls missed in my vcf, but present in HapMap (useful to take a look at why these variants weren't called by this dataset):
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant hapmap.vcf \
   --discordance myCalls.vcf
   -o output.vcf \
   -sn mySample
 
 Select all calls made by both myCalls and hisCalls (useful to take a look at what is consistent between the two callers):
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant myCalls.vcf \
   --concordance hisCalls.vcf
   -o output.vcf \
   -sn mySample
 
 Generating a VCF of all the variants that are mendelian violations:
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant input.vcf \
   -bed family.ped \
   -mvq 50 \
   -o violations.vcf
 
 Creating a set with 50% of the total number of variants in the variant VCF:
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant input.vcf \
   -o output.vcf \
   -fraction 0.5
 
 Select only indels from a VCF:
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant input.vcf \
   -o output.vcf \
   -selectType INDEL
 
 Select only multi-allelic SNPs and MNPs from a VCF (i.e. SNPs with more than one allele listed in the ALT column):
 java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R human_g1K_v37.fasta \
   -T SelectVariants \
   --variant input.vcf \
   -o output.vcf \
   -selectType SNP -selectType MNP \
   -restrictAllelesTo MULTIALLELIC

合理组合上述两个工具，可以实现很多操作。比如说找出某个样品中特异的突变，就可以先使用CombineVariants,而后SelectVariants。

java -jar GenomeAnalysisTK.jar
     –R human_g1K_v37.fasta \
     -T CombineVariants \
     -V:ESP ESP6500.chr1.snps.vcf \
     -V:T2D private.chr1.seq.imp.sites.vcf \
     -V:1000G 1000G_ALL_Sites.chr1.vcf \
     -o private.chr1.seq.imp.sites.COMBINED.vcf
java –jar GenomeAnalysisTK.jar
     –R human_g1K_v37.fasta \
     -T SelectVariants \
     -V private.chr1.seq.imp.sites.COMBINED.vcf \
     -select “set == ‘T2D’” \
     -o private.chr1.seq.imp.sites.PRIVATE.vcf

进一步过滤可以使用VariantFiltration。

java -Xmx2g -jar GenomeAnalysisTK.jar \
   -R uhuman_g1K_v37.fasta \
   -T VariantFiltration \
   -o output.vcf \
   --variant input.vcf \
   --filterExpression "AB < 0.2 || MQ0 > 50" \
   --filterName "Nov09filters" \
   --mask mask.vcf \
   --maskName InDel

For SNPs:

QD < 2.0 MQ < 40.0 FS > 60.0
HaplotypeScore > 13.0
MQRankSum < -12.5 ReadPosRankSum < -8.0 For indels: QD < 2.0 ReadPosRankSum < -20.0 InbreedingCoeff < -0.8 FS > 200.0

评价结果可以使用VariantEval。在评价时，通常使用一个标准文件来同时进行评价，最后统计两者的分布的差异，如果它与标准分布曲线一致，客观说明结果的质量很高。

java -jar GenomeAnalysisTK.jar \
     -R human_g1K_v37.fasta –L 20 \
     -T VariantEval \
     -eval:myCalls testCallSet.chr20.vcf \
     -eval:1000G_ALL_SITES 1000G_ALL_Sites.vcf.gz \
     -eval:1000G_FIN 1000G_62_FIN_Genotypes.bcf \
     -comp:OMNI_POLY_SITES omni_poly.vcf \
     -comp:OMNI_MONO_SITES omni_mono.vcf \
     -comp:1000G_ALL_SITES 1000G_ALL_Sites.vcf.gz \
     -D dbsnp_137_129_sites.vcf \
     -o myCalls_alongside_1000G.eval

其中第一个eval是自己的结果，第二及第三个都是两个标准文件。1000 Genomes sites有利于评价mates for expected Ti/Tv and indel size distribution，62—FIN有利于评价or number of calls, novelty, etc。具体评价可以参见http://www.broadinstitute.org/gatk//events/2038/GATKwh0-BP-8-Variant_analysis.pdf。

第四步：注释
需要使用到两个工具：snpEff以及VariantAnnotator.这里不做介绍。

发布者：admin

14 thoughts on “GATK练习题”

  1. 请问我开始bwa用的是ucsc的hg19, 在snpeff过程中用ensembl的数据库，会不会不匹配？
    另外-an InbreedingCoeff 是根据注释校正，可以不用吗？

    回复

    1. 大多数应该还是匹配的。但还是有很多不匹配。同时你也要注意，ensembl的注释是不包括chr这三个字符的，这一点与ucsc不同。inbreedingCoeff最好还是做。当然不做也只是增加假阳性。

      回复

  2. 感谢坛主回答，您的这个帖子非常实用，对我很有帮助。目前我遇到一个问题，似乎snpEff的database出了些问题，我在本地下不到数据库。而且我前期都是用ucsc hg19 的基因组序列比对，所以也尝试了从ucsc上下载refgene等gtf去人工构建snpEff，但总是出错。不知坛主能否提供一个利用ucsc基因组和gtf进行annotation的攻略？不胜感激！

    回复

  4. LZ你好，文章中local realignment那部分出现的reads覆盖图很鲜明的展现出了reads的比对情况以及snp，想问下那个图是是用软件实现的还是，自己写脚本实现的？多谢！

    回复

  5. 我所有需要做的事情欧老师都有讲！实在是太好了！
    目前需要做一件事情，是通过VCF文件重建许多个个体的二倍体基因组序列，还需要过滤一下每个sample的DP和GQ。您做过类似工作么？看了VCFx，貌似不能过滤sample的质量。目前考虑一个一个个体的使用SelectVariants+FastaAlternativeReferenceMarker，您觉得这样可以么？FastaAlternativeReferenceMarker可以输出二倍体序列么？

    回复

  6. LZ你好，能不能解释一下为什么在local relaignment时要用dbSNP数据库来确定其是否为真实的indel？dbSNP数据库是记录正常人的SNP吗？

    回复

  7. 您好，在进行RealignerTargetCreator时候，必须要有KNOWN site？如果没有known site的VCFfile， 咋办？谢谢

    回复

    1. 因为GATK是针对人而设计的，如果没有，跑它的pipeline会困难重重。我个人还是不建议使用。如果一定要用，也许可以跳过这个步骤。

      回复

发表评论

邮箱地址不会被公开。 必填项已用*标注