GATK4简明用法

GATK4是最新的GATK版本，它在算法上进行了优化，运行速率得到提高，而且整合了picard。GATK4依然是用java 语言开发的，但使用方式上更加人性化，比如所有命令都是gatk cmd方式，这里的cmd是任何可以用的cmd。GATK4 的最佳实践给出了5套pipeline: Germline SNP/Indel, Somatic SNV/Indel, RNAseq SNP/Indel, Germline CNV, Somatic CNV。本文是前段时间参与Broad和Intel中国在北京的培训班上的精简记录，供自己参考用，主要是我所关注的SNV/Indel。

背景知识

■ 建库library


注意几个名字的理解：

  • fragment修补: 末端补平、3‘加A、加接头(插有barcode)、PCR扩增(引入P5/P7)
  • adaptor: 测序primer结合(Read 1 primer，Read 2 primer)，其中反向adaptor也是
  • index barcode:
  • P5/P7: flowcell表面结合(桥式扩增时两个都用，测序时之前切断P7)

■ BAM/SAM file

■ VCF file

■ 不同类型的variant

下面将要介绍的内容有：

Module 1: Data Pre-processing

Prepare mapped, cleaned and sorted BAM
Step 1: MAPPING
BWA for DNA, STAR for RNA-seq
Step 2: Sort

1
gatk SortSam  -I  -O

Step 3: Mark duplicates

1
gatk MarkDuplicates -I  -O

Step 4: BQSR
BaseRecalibrator

1
2
3
gatk BaseRecalibrator -R ref.fasta -I sample.bam -O recal.table
          -knownSites snps.vcf.gz -knownSites indels.vcf.gz
OR        -O recal_2nd.table  -bqsr recal.table #generate the 2nd recal table

ApplyBQSR

1
2
3
4
gatk ApplyBQSR -R ref.fasta -I sample.bam -O sample_bqsr.bam
               -bqsr recal.table
               -SQQ 10 -SQQ 20 -SQQ 30 -SQQ40 # bin quals
               --emit_original_quals  # emit original quals(OQ) to OQ tag

Tool:AnalyzeCovariates

1
gatk AnalyzeCovariates -before recal.table -after recal_2nd.table -plots plots.pdf

Notes:
alt contigs in GRCh38; the unmapped bam workflow; RNA-seq mapping

Module 2: Best Practice for Germiline SNP & InDel



■ Step 1: Call variants with HaplotypeCaller (consolidate/joint-call with cohort)
Note: VCF和GVCF格式略不同。GVCF称为genomic VCF，里面含有Non-variant site信息，不仅仅是该样本中的SNV。相当于提前把坐标对齐了，方便后面joint-calling。包括四步操作:

  • identify ActiveRegions
  • assemble plausible haplotypes (local realignment, collect likely haplotypes, Smith-Waterman align)
  • score haplotypes using PairHMM
  • genotype each sample at each potential variant site

  • Default mode

    1
    2
    gatk HaplotypeCaller -R ref.fasta -I sample.bam -O sample.vcf
    -L 20:10000-11000

  • Realignment mode and ensemble haplotypes:

    1
    2
    gatk HaplotypeCaller -R ref.fasta -I sample.bam -O sample.vcf
    #可产生多个ArtificialHaplotype

  • GVCF mode

    1
    2
    3
    gatk HaplotypeCaller -R ref.fasta -I sample.bam -O sample.g.vcf
    -ERC GVCF|BP_RESOLUTION -L 20:10000-11000
    #GVCF for block, BR_RESOLUTION for site

GenomicsDBImport Consolidate GVCFs

1
2
3
4
5
gatk GenomicsDBImport -R ref.fasta -V sample1.g.vcf -V sample2.g.vcf -V sample3.g.vcf
--genomicsdb-workspace-path output/trio
--intervals 20:10000000-10200000 | chr20
#run on a single interval at a time 每次只能一个interval，如果多个用WDL scripts
or gatk CombineGVCFs -R ref.fasta -V sample1.g.vcf -V sample2.g.gvcf -O combined.g.vcf

GenotypeGVCFs Jointly genotype

1
2
gatk GenotypeGVCFs -R ref.fasta -V gendb://output/trio|combined.g.vcf -O final_varaint.vcf
[-L 20:10000000-10200000 | chr20]

■ Step 2: Filter (low quality)
to balance sensitivity and specifity.
two filtering approaches:

  • Hard-filters using binary thresholds: Applicable to all BUT requires expertise to define appropriately
  • Variant “recalibration” using machine learning: More powerful BUT requires well-curated known resources

VQSR: SNP and InDel must be separately handled

1
2
3
4
5
6
7
8
9
10
11
gatk VariantRecalibrator -R ref.fasta -V raw.SNPs.vcf -resource [tags]:filename.vcf
    -an DP -an QD -an FS -an MQRankSum -mode SNP 
    -recal-file raws.SNPs.recal 
    -tranches-file raw.SNPs.tranches 
    -rescript-file recal.plots.R
# tags can be hapmap, 1000G, omni, dbsnp
gatk ApplyVQSR -R ref.fasta -V raw.vcf -mode SNP 
    -recal-file raw.SNPs.recal -tranches-file raw.SNPs.tranches 
    -O recal.SNPs.vcf 
    -ts-filter-level 99.0

■ Step 3: Refinement

1
2
gatk CalculateGenotypePosteriors
gatk VariantFiltration

Subset callset (only SNP and only one sample)

1
2
gatk SelectVariants -R  ref.fasta -V trio.vcf.gz -sn NA12878 -select-type SNP 
--exclude-non-variants -O motherSNP.vcf.gz

■ Step 4: Annotate

1
2
3
gatk VariantAnnotator -R ref.fasta -V output/motherSNP.vcf.gz 
--resource giab:resources/motherGIABsnps.vcf.gz -E giab.callsets 
-O output/motherSNP.giab.vcf.gz

Tabulate annotation result and Visualize

1
2
3
4
5
6
gatk VariantsToTable -R ref/ref.fasta -V output/motherSNP.giab.vcf.gz
    -F CHROM -F POS -F QUAL -F BaseQRankSum -F MQRankSum
     -F ReadPosRandSum -F DP -F DP_Orig \
    -F FS -F MQ -F SOR -F giab.callsets -GF GO 
    -O output/motherSNP.giab.txt
Visualize in R

Filter variants with QUAL<30 (<10 by default)

1
gatk VariantFiltration -R ref/ref.fasta -V output/motherSNP.vcf.gz --filter-expression "QUAL <30" --fil ter-name "qual30" -O output/motherSNPqual30.vcf.gz



■ Step 5: Callset Evaluation

■ Step 6: check with IGV

Note: Use RNA-seq data (mapping with STAR)

Module 3: Best Practice for Somatic SNP & InDel


different from Germline variants, due to purity (contamination of N/T) and heterogeneity

■ Step 1: Call variants with Mutect2 (denoise with normal match and PoN)

Tumor-only mode

1
2
3
4
5
gatk Mutect2 -R ref/Homo_sapiens_assembly38.fasta
    -I bams/tumor.bam -tumor HCC1143_tumor
    --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter
    -L resources/chr17plus.interval_list
    -O output/1_somatic_m2.vcf.gz

Tumor with match mode (with realignment)

1
2
3
4
5
6
7
8
9
10
gatk Mutect2 -R ref/Homo_sapiens_assembly38.fasta
    -I bams/tumor.bam -I bams/normal.bam
    -tumor HCC1143_tumor -normal HCC1143_normal
    -pon resources/chr17_m2pon.vcf.gz
    --germline-resource resources/chr17_af-only-gnomad_grch38.vcf.gz
    --af-of-alleles-not-in-resource 0.0000025
    --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter
    -L resources/chr17plus.interval_list
    -O output/1_somatic_m2.vcf.gz
    -bamout output/2_tumor_normal_m2.bam

比对paired normal /a panel of normal /a germline population resource，更能发现真实的somatic SNV
■ Step 2: Filter(contamination)

1
2
3
4
5
6
7
8
gatk GetPileupSummaries -I bams/tumor.bam 
    -V resources/chr17_small_exac_common_3_grch38.vcf.gz 
    -O sandbo x/7_tumor_getpileupsummaries.table
gatk CalculateContamination -I output/7_tumor_getpileupsummaries.table 
    -O output/8_tumor_calculationcontamination.table
gatk FilterMutectCalls -V output/1_somatic_m2.vcf.gz 
    --contamination-table output/8_tumor_calculationcontamination.table 
    -O output/9_somatic_oncefiltered.vcf.gz

■ Step3: check with IGV

Commands in GATK

Base Calling:

1
2
3
4
5
6
7
CheckIlluminaDirectory (P) 
CollectIlluminaBasecallingMetrics (P)
CollectIlluminaLaneMetrics (P)
ExtractIlluminaBarcodes (P)
IlluminaBasecallsToFastq (P)
IlluminaBasecallsToSam (P)
MarkIlluminaAdapters (P)

Copy Number Variant Discovery:

1
2
3
4
5
6
7
8
9
10
11
AnnotateIntervals
CallCopyRatioSegments
CombineSegmentBreakpoints
CreateReadCountPanelOfNormals
DenoiseReadCounts
DetermineGermlineContigPloidy
GermlineCNVCaller
ModelSegments
PlotDenoisedCopyRatios
PlotModeledSegments
PostprocessGermlineCNVCalls

Coverage Analysis:

1
2
3
4
5
6
7
8
9
10
ASEReadCounter
CollectAllelicCounts
CollectReadCounts
CountBases
CountBasesSpark
CountReads
CountReadsSpark
GetPileupSummaries
Pileup
PileupSpark

Diagnostics and Quality Control:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
AccumulateVariantCallingMetrics (P)
AnalyzeCovariates
BamIndexStats (P)
CalcMetadataSpark
CalculateContamination
CalculateReadGroupChecksum (P)
CheckFingerprint (P)
CheckPileup
CheckTerminatorBlock (P)
ClusterCrosscheckMetrics (P)
CollectAlignmentSummaryMetrics (P)
CollectBaseDistributionByCycle (P)
CollectBaseDistributionByCycleSpark
CollectGcBiasMetrics (P)
CollectHiSeqXPfFailMetrics (P)
CollectHsMetrics (P)
CollectIndependentReplicateMetrics (P)
CollectInsertSizeMetrics (P)
CollectInsertSizeMetricsSpark
CollectJumpingLibraryMetrics (P)
CollectMultipleMetrics (P)
CollectMultipleMetricsSpark
CollectOxoGMetrics (P)
CollectQualityYieldMetrics (P)
CollectQualityYieldMetricsSpark
CollectRawWgsMetrics (P)
CollectRnaSeqMetrics (P)
CollectRrbsMetrics (P)
CollectSequencingArtifactMetrics (P)
CollectTargetedPcrMetrics (P)
CollectVariantCallingMetrics (P)
CollectWgsMetrics (P)
CollectWgsMetricsWithNonZeroCoverage (P)
CompareBaseQualities
CompareDuplicatesSpark
CompareMetrics (P)
CompareSAMs (P)
ConvertSequencingArtifactToOxoG (P)
CrosscheckFingerprints (P)
CrosscheckReadGroupFingerprints (P)
EstimateLibraryComplexity (P)
EstimateLibraryComplexityGATK
FlagStat
FlagStatSpark
GetSampleName
MeanQualityByCycle (P)
MeanQualityByCycleSpark
QualityScoreDistribution (P)
QualityScoreDistributionSpark
ValidateSamFile (P)
ViewSam (P)

Intervals Manipulation:

1
2
3
4
5
6
BedToIntervalList (P)
IntervalListToBed (P)
IntervalListTools (P)
LiftOverIntervalList (P)
PreprocessIntervals
SplitIntervals

Metagenomics:

1
2
3
4
5
6
PathSeqBuildKmers
PathSeqBuildReferenceTaxonomy
PathSeqBwaSpark
PathSeqFilterSpark
PathSeqPipelineSpark
PathSeqScoreSpark

Other:

1
2
3
4
5
6
CreateHadoopBamSplittingIndex
FifoBuffer (P)
GatherBQSRReports
GatherTranches
IndexFeatureFile
ParallelCopyGCSDirectoryIntoHDFSSpark

Read Data Manipulation:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
AddCommentsToBam (P)
AddOrReplaceReadGroups (P)
ApplyBQSR
ApplyBQSRSpark
BQSRPipelineSpark
BamToBfq (P)
BaseRecalibrator
BaseRecalibratorSpark
BaseRecalibratorSparkSharded
BuildBamIndex (P)
BwaAndMarkDuplicatesPipelineSpark
BwaSpark
CleanSam (P)
ClipReads
ConvertHeaderlessHadoopBamShardToBam
DownsampleSam (P)
ExtractOriginalAlignmentRecordsByNameSpark
FastqToSam (P)
FilterSamReads (P)
FixMateInformation (P)
FixMisencodedBaseQualityReads
GatherBamFiles (P)
LeftAlignIndels
MarkDuplicates (P)
MarkDuplicatesGATK
MarkDuplicatesSpark
MarkDuplicatesWithMateCigar (P)
MergeBamAlignment (P)
MergeSamFiles (P)
PositionBasedDownsampleSam (P)
PrintReads
PrintReadsSpark
ReorderSam (P)
ReplaceSamHeader (P)
RevertBaseQualityScores
RevertOriginalBaseQualitiesAndAddMateCigar (P
RevertSam (P)
SamFormatConverter (P)
SamToFastq (P)
SetNmAndUqTags (P)
SetNmMdAndUqTags (P)
SimpleMarkDuplicatesWithMateCigar (P)
SortSam (P)
SortSamSpark
SplitNCigarReads
SplitReads
SplitSamByLibrary (P)
SplitSamByNumberOfReads (P)
UmiAwareMarkDuplicatesWithMateCigar (P)
UnmarkDuplicates

Reference:

1
2
3
4
5
6
7
8
BaitDesigner (P)
BwaMemIndexImageCreator
CreateSequenceDictionary (P)
ExtractSequences (P)
FindBadGenomicKmersSpark
NonNFastaSize (P)
NormalizeFasta (P)
ScatterIntervalsByNs (P)

Short Variant Discovery:

1
2
3
4
5
6
7
CombineGVCFs
GenomicsDBImport
GenotypeGVCFs
HaplotypeCaller
HaplotypeCallerSpark
Mutect2
ReadsPipelineSpark

Structural Variant Discovery:

1
2
3
4
5
DiscoverVariantsFromContigAlignmentsSAMSpark
ExtractSVEvidenceSpark
FindBreakpointEvidenceSpark
StructuralVariationDiscoveryPipelineSpark
SvDiscoverFromLocalAssemblyContigAlignmentsSp

Variant Evaluation and Refinement:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
AnnotatePairOrientation
AnnotateVcfWithBamDepth
AnnotateVcfWithExpectedAlleleFraction
CNNScoreVariants
CNNVariantTrain
CNNVariantWriteTensors
CalculateGenotypePosteriors
CalculateMixingFractions
Concordance
CountFalsePositives
CountVariants
CountVariantsSpark
FilterVariantTranches
FindMendelianViolations (P)
Funcotator
GenotypeConcordance (P)
ValidateBasicSomaticShortMutations
ValidateVariants
VariantsToTable

Variant Filtering:

1
2
3
4
5
6
7
ApplyVQSR
CreateSomaticPanelOfNormals
FilterByOrientationBias
FilterMutectCalls
FilterVcf (P)
VariantFiltration
VariantRecalibrator

Variant Manipulation:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
FixVcfHeader (P)
GatherVcfs (P)
GatherVcfsCloud
LiftoverVcf (P)
MakeSitesOnlyVcf (P)
MergeVcfs (P)
PrintVariantsSpark
RemoveNearbyIndels
RenameSampleInVcf (P)
SelectVariants
SortVcf (P)
SplitVcfs (P)
UpdateVCFSequenceDictionary
UpdateVcfSequenceDictionary (P)
VariantAnnotator
VcfFormatConverter (P)
VcfToIntervalList (P)

赏

谢谢你请我吃糖果

支付宝
微信

扫一扫，分享到微信

微信分享二维码
这不是科学的乐园<br>而是科学的森林<br>无处不充满着引人入胜的奥秘<br>Email:sunnysean@qq.com<br>