How to Call somatic SNVs and indels using MuTect2

IMPORTANT: This is the legacy GATK documentation. This information is only valid until Dec 31st 2019. For latest documentation and forum click here

created by shlee

on 2017-03-15

This tutorial covers the GATK3 version of the MuTect2 tool which is in a BETAstate and will not be developed further. A newer and much-improved version of this tool is included in GATK4; we recommend you use that tool rather than this one. Accordingly, this document is provided as is and with no guarantees that everything will work correctly.

Tutorial example data are just that--examples we use to illustrate tool features. Data are inappropriate for deriving biological significance. Although we have aligned and post-alt processed reads to GRCh38, in this tutorial we focus only on results from the primary assembly.

1. Create a sites-only panel of normals (PoN) with MuTect2 and CombineVariants

The PoN allows additional filtering of calls, e.g. those that arise from technical artifacts. Therefore, it is important that the PoN consist of samples that are technically similar to and representative of the tumor samples, e.g. that were sequenced on the same platform using the same chemistry and analyzed using the same toolchain. For the tutorial, we have already created a PoN using forty 1000 Genomes Project exome samples aligned to GRCh38. We chose these randomly, downloaded CRAM files from the 1000 Genomes Project GRCh38 FTP site, converted to BAMs and indexed, and used the data directly with MuTect2 without further modification.

Within the tutorial bundle, the 1kg40m2ponsubset50k.vcf.gz_ and 1kg40m2ponsitesonlysubset50k.vcf.gz files represent the 40-sample PoN that we will use in Sections 5.

We emulate the PoN creation steps below but with small substitute data.

1.1. First, per sample BAM, invoke MuTect2’s --artifact_detection_mode to generate a per-normal-sample callset.

java -jar $GATK \ -T MuTect2 \ -I:tumor hcc1143_N_subset50k.bam \ --dbsnp dbSNP142_GRCh38_subset50k.vcf.gz \ --artifact_detection_mode \ -o 1_normalforpon.vcf.gz \ -L chr6:33,413,000-118,315,000 \ -R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta

We include the dbSNP resource so as to annotate sites already present in the database. Here we make the motions of this step using our matched normal. In the real world, you would perform this step on other germline samples that are NOT your tumor's matched normal.

1.2. Second, we use CombineVariants on all the per-normal-sample VCFs at once.

java -jar $GATK \ -T CombineVariants \ --arg_file inputs.list \ -minN 2 \ --setKey "null" \ --filteredAreUncalled \ --filteredrecordsmergetype KEEP_IF_ANY_UNFILTERED \ -o 2_pon_combinevariants.vcf.gz \ -R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta \ --genotypemergeoption UNIQUIFY

We keep only the variant sites that are in at least two samples, using the -minN 2 option. We omit filtered calls with the --filteredAreUncalled option but ensure a site is kept if at least two samples pass the site (KEEP_IF_ANY_UNFILTERED, -minN 2). Also, instead of supplying each of the forty VCFs with -V, we pass in an argument file ending in .list with the --arg_file option. This file lists all the input arguments in the same line, e.g. “-V file1.vcf -V file2.vcf …”.

In the real world, you would NOT use the last option--genotypemergeoption UNIQUIFY. The --genotypemergeoption UNIQUIFY allows us to use the same normal VCF twice, as different samples, to illustrate this step with the data at hand. Without it, the command will error because the read group sample RGSM tags are identical for our two VCFs. This generates a PoN with 13 records.

1.3. Finally, generate a sites-only VCF with Picard's MakeSitesOnlyVcf.

In the command in section 1.2, we could have used the --minimalVCF option to generate a simplified sites-only VCF. However, in this tutorial we deviate from the standard practice because we are in learning-mode, and so we care to retain information when possible. Using the --minimalVCF option with CombineVariants is faster but produces a qualitatively different VCF that removes TLOD scores. Here we generate a sites-only VCF that retains TLOD scores by removing the sample columns with Picard's MakeSitesOnlyVcf.

java -jar $PICARD MakeSitesOnlyVcf \ I=2_pon_combinevariants.vcf.gz \ O=3_pon_siteonly.vcf.gz

The command also generates an index--3_pon_siteonly.vcf.gz.tbi. This produces records that have only eight columns. Compare the before and after using gzcat.

gzcat 3_pon_siteonly.vcf.gz | grep -v '##' | less

back to top

2. (Optional) Estimate cross-sample contamination with ContEst

ContEst detects cross-sample contamination and to some extent sample swaps. Even if our sample data are from cultured cell lines and we need not factor for tumor purity and tumor heterogeneity, we still need to consider contamination from other human sources and whether the extent of contamination is prohibitive or is in an acceptable range that allows us to proceed in analyzing a sample. We need to estimate contamination for both the tumor and the normal.

ContEst informs your downstream filtering. Consider tumor types in which you expect low purity. It would be good to know whether you have 0–2% or 2–5% contamination, as, depending on expected mutation rates, one result allows you to progress in analysis and the other requires planning manual review or even tossing the data.

For on-the-fly genotyping of the normal, ContEst will call any site with >80% bases showing an alternate allele with at least 50x coverage homozygous variant. If your normal sample has lower coverage overall, then you can alternately provide ContEst with a genotyped VCF using the --genotypes option.

java -jar $GATK -T ContEst \ -I:eval hcc1143_T_subset50k.bam \ -I:genotype hcc1143_N_subset50k.bam \ --popfile hapmap_3.3_grch38_pop_stratified_af_subset50k.vcf.gz \ --interval_set_rule INTERSECTION \ -o 4_T_contest.txt \ -R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta

For the normal, use a genotyped VCF of the normal. Our genotypes are from HaplotypeCaller.

java -jar $GATK -T ContEst \ -I:eval hcc1143_N_subset50k.bam \ --genotypes hcc1143_N_haplotypecaller50k.vcf.gz \ --popfile hapmap_3.3_grch38_pop_stratified_af_subset50k.vcf.gz \ --interval_set_rule INTERSECTION \ -o 5_N_contest.txt \ -R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta

If there is insufficient data to make a meaningful estimate, the resulting file will contain a warning. Take a look at each result using cat. Of the eight columns, the fourth column, labeled contamination, gives the percent contamination.

Converting 1.9% to a fraction gives us 0.019. So we know going forward that any calls with less than ~0.02 alternate allele fraction are suspect. For the tutorial, we omit the contamination estimate from the MuTect2 analysis.

Could this represent the fraction of a tumor subpopulation? If we estimate contamination for the normal sample using the tumor BAM to genotype on-the-fly, then ContEst estimates 42.2% contamination using 2,070 sites. Why does this happen?

☞ How does MuTect use the contamination estimate?

The handling differs for the different versions of MuTect. MuTect1 incorporates the contamination estimate with the --fraction_contamination option as a hard cutoff for calling in all samples. For the same parameter, v3.7 MuTect2 uses the fraction to downsample reads for each alternate allele. For example, if a site contains X coverage depth, then MuTect2 will remove X times the contamination fraction of reads for each alternate allele.

back to top

3. (Optional) Estimate extent of FFPE and OxoG artifacts with CollectSequencingArtifactMetrics

Picard’s CollectSequencingArtifactMetrics provides us with an estimate of the extent of FFPE and OxoG artifacts as well as other potential artifacts related to sequence context. These the tool categorizes as those that occur before hybrid selection (preadapter) and those that occur during hybrid selection (baitbias). The tool provides a global view of the genome that empowers decision making in ways that site-specific analyses cannot. For example, the metrics can help you decide whether you should consider downstream filtering.

java -jar $PICARD CollectSequencingArtifactMetrics \ I=hcc1143_T_subset50k.bam \ O=6_T_artifact \ R=~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta

java -jar $PICARD CollectSequencingArtifactMetrics \ I=hcc1143_N_subset50k.bam \ O=7_N_artifact \ R=~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta

Of the multiple metrics files (five each), examine the preadaptersummarymetrics_ using cat. Here are the metrics for the full callsets that show similar results to that of our small data.

☞ How do I interpret the TOTAL_QSCORE?

The TOTAL_QSCORE is Phred-scaled such that lower scores equate to a higher probability the change is artifactual. E.g. 40 translates to 1 in 10,000 probability. For OxoG, a rough cutoff for concern is 30. If numbers are less than 30, then plan to filter your somatic calls further.

☞ What is FFPE? What is OxoG?

FFPE stands for formalin-fixed, paraffin-embedded. Formaldehyde deaminates cytosines and thereby results in C→T transition mutations.

For an explanation of OxoG see its Dictionary entry. Breifly, oxidation of guanine to 8-oxoguanine results in G→T transversion mutations during library preparation. Another Picard tool, CollectOxoGMetrics, similarly gives Phred-scaled scores for the 16 three-base extended sequence contexts. For even more details on this artifact, see the 2013 Costello et al, article in Nucleic Acids Research. For SNVs in our MuTect2 callset, the FOXOG annotation refers to fraction oxoG. Downstream tools such as Broad CGA’s D-ToxoG use this information. For more discussion, see Article#8771 and Article#8183.

back to top

4. Call somatic SNVs and indels with MuTect2

The above comparison is limited. How would you compare calling on GRCh37 versus on GRCh38? Given what you know about alt-aware alignment and post-alt processing, is there an additional step that this analysis is missing? Below are the six GRCh37 passing indel calls. Do they differ from the six GRCh38 passing indel calls? If so, how?

Examining the records, we find the GRCh38 calls differ qualitatively from that of GRCh37. This simplistic comparison is meant to show how deriving a good somatic callset involves comparing callsets, e.g. from different callers, manually reviewing passing and filtered calls and, as alluded to, additional filtering.

back to top

6. Generate BAMOUTs with MuTect2 and manually review alignments with IGV

Manual review extends from deciphering call record annotations to the nitty-gritty of reviewing read alignments using a visualizer.

6.1. Generate the BAMOUT with MuTect2

To review read alignments, we must generate those that reflect MuTect2’s graph-assembly results. We call these bamouts for the parameter we invoke. To demonstrate, we use a small interval.

java -jar $GATK -T MuTect2 \ -I:tumor hcc1143_T_subset50k.bam \ -I:normal hcc1143_N_subset50k.bam \ --dbsnp dbSNP142_GRCh38_subset50k.vcf.gz \ --cosmic CosmicCodingMuts_grch38_subset50k.vcf.gz \ --normal_panel 1kg_40_m2pon_sitesonly_subset50k.vcf.gz \ --disableOptimizations --dontTrimActiveRegions --forceActive \ -L 8_mutect2.vcf.gz \ -ip 1000 \ -bamout 10_m2_bamout.bam \ -R ~/Documents/ref/hg38/Homo_sapiens_assembly38.fasta

Notice the special options --disableOptimizations, --dontTrimActiveRegions and --forceActive (many thanks to @Sheila for these). To easily generate bamouts for regions surrounding call sites, we use the four-record VCF from step 4 with the -L option and add padding around the sites, e.g. 1000 bases before and after each site, with the -ip option.

Similarly, we have pre-generated a larger bamout for you using a VCF of the 76 solo-filtered indel sites and 22 passing indel sites.

6.2. Setup on IGV for manual review

In the remaining exercise, we will use data subset around the 76-solo-filtered and 22 passing indel sites. We will review data in IGV. First, load GRCh38 as the reference in IGV. Then load these files in order:

i. 98_m2_indels.vcf.gz

ii. 98_m2_bamout.bam

iii. hcc1143_T_subset50k.bam

iv. hcc1143_N_subset50k.bam

Go to the SLC35F1 locus at chr6:118314029 by typing in the genomic coordinates into IGV. If these alignments seem hard to decipher, it is because we need to tweak some settings.

a. In the View>Preferences>Alignments panel, turn off downsampling and turn on the center line.

b. Center the track on the deletion you see by dragging the track so that the center line is at the leftmost edge of the deletion. You may have to scroll to find the deletion.

c. In the tracks-view, right-click on the bamout track and select the following:

- Group by sample

- Color by tag: HC

- Sort by base

What are the three grouped tracks for the bamout? How do you feel about this indel call?