How to Run a Somatic Variant Calling Pipeline on Your Own Sequencing Data

By the end of this you will have a filtered, annotated VCF of candidate somatic SNVs and indels from a pair of aligned read sets (or one, if you have no matched normal), a contamination estimate, a read-orientation artifact model, and a concrete sense of the variant allele fraction floor your data supports. You need: a reference FASTA matching the one used for alignment (GRCh38 with alt contigs, or the no-alt analysis set if your aligner ignored them), coordinate-sorted, duplicate-marked BAM or CRAM files with a valid @RG line, GATK 4.4+, bcftools, samtools, roughly 8 GB RAM per Mutect2 shard, and 2–3× your BAM size in scratch. If you are working from a Personal Molecular Profile, your input is blood or saliva WGS, which changes what “somatic” means in practice. More on that in step 1.

1. Decide what you are calling

A somatic variant is a change present in some of your cells and absent from the zygote you started as. It arises after fertilization: during development (mosaicism), during ordinary tissue turnover (clonal expansions in blood, skin, esophagus), or in a tumor. A germline variant is in every nucleated cell and gets transmitted. The BRCA question that shows up constantly has both answers: a pathogenic BRCA1 allele inherited from a parent is germline and sits at ~50% allele fraction in every tissue, while a BRCA1 mutation confined to tumor tissue is somatic and sits at whatever fraction the tumor clone occupies. The sequencing is identical. Only the comparison distinguishes them.

That comparison is the whole design of the pipeline. HaplotypeCaller assumes a diploid genome and tests genotypes {0/0, 0/1, 1/1}, so it wants allele fractions near 0, 0.5, or 1.0 and penalizes everything else. Mutect2 drops the ploidy prior entirely. It models the tumor as a mixture at an unknown allele fraction f, computes the likelihood of the read pileup under “variant present at some f > 0” versus “reference plus sequencing error,” and emits a log-odds score (TLOD). That is the only structural difference that matters, and it is why you cannot get 5% allele fraction calls out of HaplotypeCaller no matter how you tune it.

If your samples are tumor and matched normal, this is the standard workflow. If your only sample is blood WGS from a healthy person, the realistic target is clonal hematopoiesis and post-zygotic mosaicism, which live mostly below 10% allele fraction, and 30× WGS will not see them. Skip to step 8 before you spend compute.

2. Get the alignments right before you call anything

Somatic callers punish bad alignment harder than germline callers do, because a 5% allele fraction call has 2–3 supporting reads and a single misalignment cluster looks exactly like a real subclone. Minimum bar:

bwa-mem2 mem -t 32 -K 100000000 -Y \
  -R '@RG\tID:tum1\tSM:TUMOR\tPL:ILLUMINA\tLB:lib1' \
  GRCh38.fa r1.fq.gz r2.fq.gz \
  | samtools sort -@ 8 -m 2G -o tumor.bam -
samtools index tumor.bam
gatk MarkDuplicates -I tumor.bam -O tumor.md.bam -M tumor.dup.txt \
  --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500
gatk BaseRecalibrator -R GRCh38.fa -I tumor.md.bam \
  --known-sites Homo_sapiens_assembly38.dbsnp138.vcf \
  --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
  -O tumor.recal.table
gatk ApplyBQSR -R GRCh38.fa -I tumor.md.bam -bqsr tumor.recal.table -O tumor.final.bam

The -K 100000000 fixes the chunk size so alignment is deterministic across thread counts, which you want if you ever need to reproduce a call. -Y uses soft clipping for supplementary alignments, which keeps structural-variant-adjacent reads recoverable.

Coverage targets: 30× for the normal, 60–100× for a tumor at typical purity, and much deeper for anything low-fraction. Clinical sequencing practice pushes higher and requires per-base depth reporting rather than mean depth, because mean depth hides the exons sitting at 8×.¹ If you would rather not assemble this yourself, nf-core/sarek wraps the whole thing (alignment, BQSR, multiple somatic callers, annotation) in a Nextflow pipeline that runs identically on a laptop, a cluster, or cloud, and it is what we reach for when we want a reproducible end-to-end run instead of a hand-built one.²

3. Build a panel of normals

The panel of normals (PoN) is the single highest-yield filter in the pipeline. It captures recurrent artifacts: mapping errors in segmental duplications, systematic miscalls at homopolymers, library-prep damage signatures. Any site that shows a non-reference allele in multiple unrelated normals is an artifact, regardless of what your sample says.

Use 25–40 normals sequenced on the same platform, same library chemistry, same read length as your sample. A PoN built from a different chemistry is worse than no PoN, because it filters the wrong things and passes your actual artifacts.

# per normal, in tumor-only mode
for n in normals/*.bam; do
  gatk Mutect2 -R GRCh38.fa -I "$n" -max-mnp-distance 0 \
    -O "pon/$(basename "$n" .bam).vcf.gz"
done

gatk GenomicsDBImport -R GRCh38.fa -L wgs_calling_regions.hg38.interval_list \
  --genomicsdb-workspace-path pon_db \
  --merge-input-intervals true \
  $(for v in pon/*.vcf.gz; do echo -n "-V $v "; done)

gatk CreateSomaticPanelOfNormals -R GRCh38.fa \
  -V gendb://pon_db \
  --germline-resource af-only-gnomad.hg38.vcf.gz \
  --min-sample-count 2 \
  -O pon.vcf.gz

-max-mnp-distance 0 is mandatory here. Without it, adjacent SNVs get merged into MNPs that will not match the single-nucleotide records in your sample VCF, and the filter silently does nothing. GATK’s public 1000g_pon.hg38.vcf.gz is a reasonable fallback if you cannot build your own, but it is 1000 Genomes chemistry and it will miss artifacts specific to your library prep.

4. Call with Mutect2

Scatter by interval. A whole genome on one core takes days. Twenty-four shards on one chromosome each takes a few hours.

gatk --java-options "-Xmx8g" Mutect2 \
  -R GRCh38.fa \
  -I tumor.final.bam \
  -I normal.final.bam \
  -normal NORMAL_SM \
  --germline-resource af-only-gnomad.hg38.vcf.gz \
  -pon pon.vcf.gz \
  --f1r2-tar-gz f1r2.chr1.tar.gz \
  --native-pair-hmm-threads 4 \
  -L chr1 \
  -O unfiltered.chr1.vcf.gz

-normal NORMAL_SM must be the SM value from the normal’s read group, not the filename. Getting this wrong is the most common silent failure in the whole pipeline: Mutect2 treats both inputs as tumors, and you get a VCF full of germline heterozygotes.

The germline resource is an allele-frequency-only VCF. Mutect2 uses the population AF as a prior to compute GERMQ, the Phred-scaled probability that a call is germline rather than somatic. It does not hard-filter on it.

Then merge the shards:

gatk MergeVcfs $(for v in unfiltered.chr*.vcf.gz; do echo -n "-I $v "; done) -O unfiltered.vcf.gz
gatk MergeMutectStats $(for s in unfiltered.chr*.vcf.gz.stats; do echo -n "-stats $s "; done) -stats merged.stats

MergeMutectStats is easy to forget and FilterMutectCalls will refuse to run without the merged stats file, since it needs the total callable site count to calibrate its false-discovery threshold.

5. Model the artifacts, then filter

Three inputs feed the filter: read orientation bias, cross-sample contamination, and the raw call stats.

gatk LearnReadOrientationModel \
  $(for f in f1r2.chr*.tar.gz; do echo -n "-I $f "; done) \
  -O read-orientation-model.tar.gz

gatk GetPileupSummaries -I tumor.final.bam \
  -V small_exac_common_3.hg38.vcf.gz -L small_exac_common_3.hg38.vcf.gz \
  -O tumor.pileups.table
gatk GetPileupSummaries -I normal.final.bam \
  -V small_exac_common_3.hg38.vcf.gz -L small_exac_common_3.hg38.vcf.gz \
  -O normal.pileups.table

gatk CalculateContamination \
  -I tumor.pileups.table -matched normal.pileups.table \
  --tumor-segmentation segments.table \
  -O contamination.table

gatk FilterMutectCalls -R GRCh38.fa \
  -V unfiltered.vcf.gz \
  --contamination-table contamination.table \
  --tumor-segmentation segments.table \
  --ob-priors read-orientation-model.tar.gz \
  --stats merged.stats \
  -O filtered.vcf.gz

The orientation model matters more than people expect. FFPE fixation causes C>T deamination and Illumina’s acoustic shearing plus oxidation causes the OxoG G>T artifact, and both show a strong strand imbalance: the alt allele appears almost entirely on F1R2 reads or almost entirely on F2R1 reads. LearnReadOrientationModel fits the prior from your own data and FilterMutectCalls uses it to assign the orientation filter. Without it, an FFPE sample can produce thousands of spurious low-fraction C>T calls.

CalculateContamination estimates the fraction of reads from another individual by looking at allele fractions at common homozygous SNP sites. Above about 0.05 you have a real problem and should re-examine sample handling. Below 0.01 it is noise.

6. Read the output

Mutect2 emits VCF 4.2, bgzipped and tabix-indexed. Per record you care about a handful of fields:

• FORMAT/AD — reference and alt read counts, per sample.
• FORMAT/AF — Mutect2’s allele fraction estimate for the tumor.
• FORMAT/F1R2 and FORMAT/F2R1 — strand-specific alt counts, the raw material for the orientation model.
• INFO/TLOD — log-odds that the variant is present in the tumor. Default call threshold is 3.0.
• INFO/NLOD — log-odds that the site is diploid-het-free in the normal.
• INFO/GERMQ, INFO/SEQQ, INFO/STRANDQ — Phred-scaled confidence that the call is not germline, not sequencing error, not strand artifact.

FILTER is where the interpretation lives. PASS means it survived everything. Common non-PASS values: weak_evidence (TLOD below threshold), germline, panel_of_normals, normal_artifact (alt reads in the normal above what contamination explains), orientation, contamination, slippage (short tandem repeat), haplotype (linked to a filtered variant on the same haplotype), clustered_events, multiallelic.

Extract PASS calls into something tabular:

bcftools view -f PASS -Oz -o pass.vcf.gz filtered.vcf.gz
bcftools index -t pass.vcf.gz
bcftools query -f '%CHROM\t%POS\t%REF\t%ALT\t%FILTER\t%INFO/TLOD\t[%AD\t%AF\t]\n' \
  -s TUMOR pass.vcf.gz > pass.tsv

Do not delete the non-PASS records. When you go back to ask why a known variant is missing, the filter tag is the answer, and it is gone if you filtered at the bcftools view stage and kept only the output.

7. Verify before you believe

Single-caller somatic results are not trustworthy on their own. A large multi-center benchmark using whole-genome tumor-normal pairs with validated truth sets found wide variation in both sensitivity and precision across pipelines, with indels substantially harder than SNVs and a significant share of disagreement traceable to alignment and filtering choices rather than the caller’s core model.³ Run at least one structurally different second caller (Strelka2 or MuSE for SNVs, which use different statistical machinery from Mutect2’s local reassembly) and treat the intersection as high confidence and the symmetric difference as the set that needs eyes on it.

If you want a fast second opinion, Sentieon’s TNhaplotyper reimplements the Mutect algorithm with the same mathematics and produces near-identical calls at a large speed advantage, which makes it useful as a determinism check and for reprocessing at scale.⁴

For specific variants you care about, alignment-free validation is cheap and independent of every choice you made above. Build k-mers that span the mutation and count their exact occurrences in the raw FASTQ, bypassing the aligner entirely. This works for substitutions, insertions, and deletions and cleanly separates real support from mapping artifacts.⁵ Then open the site in IGV with soft clips shown. Real variants scatter across read positions and both strands. Artifacts cluster at read ends, appear on one strand, or sit inside a homopolymer run.

8. Know your allele fraction floor

This is where most self-directed somatic analysis goes wrong. The detection limit is set by depth and by the per-base error rate of your chemistry, not by the caller.

At 30× WGS, a variant at 5% allele fraction has an expected 1.5 supporting reads. You will miss most of them and the few you call will be indistinguishable from error. At 100×, 5% gives you 5 reads, which is callable but sits close to the Illumina substitution error floor of roughly 0.1–1% per base. To go below 1% you need unique molecular identifiers and consensus calling across duplicate families, which means targeted panels at several thousand × raw depth, not WGS.

Allele fraction also varies enormously across real samples for biological reasons rather than technical ones. Analysis of somatic TP53 mutations across tumor datasets showed mutant allele fraction spanning a wide range driven by tumor purity, copy number at the locus, and subclonal structure, so a low fraction does not mean a low-confidence call and a high one does not mean a clonal event.⁶ You cannot read clonality off allele fraction without matched copy number and purity estimates.

If you have no matched normal, every common germline variant becomes a candidate somatic call and you are left leaning entirely on the population AF prior plus your PoN. That works reasonably for variants absent from gnomAD and badly for everything else. Long reads change this calculation, because phasing lets you separate haplotypes and identify variants present on only a fraction of molecules without a normal comparison.⁷ For short-read tumor-only work, expect a false positive rate high enough that any specific call needs orthogonal confirmation.

Common problems

Both inputs treated as tumor. The -normal argument takes the read group SM value. Verify with samtools view -H normal.bam | grep '^@RG'. Symptom: tens of thousands of PASS calls at exactly 50% allele fraction.

PoN does nothing. You built it without -max-mnp-distance 0, so adjacent SNVs merged into MNP records that never match your sample’s single-nucleotide calls. Rebuild.

FilterMutectCalls errors on missing stats. You scattered the calling step but forgot MergeMutectStats. The .stats files sit next to each shard VCF.

Reference mismatch. Mutect2 checks sequence dictionaries and will fail loudly if contig names differ (chr1 vs 1), but it will run happily on a reference with matching names and different content, such as analysis-set GRCh38 versus full GRCh38 with alts. Compare md5 values in the .dict file against the BAM header.

Thousands of low-fraction C>T calls. FFPE deamination or old DNA. Confirm by checking whether the alt reads pile onto one of F1R2/F2R1, and make sure --ob-priors was passed to the filter.

Calls clustered in one region. Segmental duplication or a collapsed repeat. Intersect with the UCSC genomicSuperDups track and the ENCODE blacklist and drop the overlaps.

Recurrent low-fraction calls in blood at genes like DNMT3A, TET2, ASXL1. These are the canonical clonal hematopoiesis loci, and at sufficient depth they are common in blood from adults. They are a real biological finding about your cells, not a pipeline artifact. What they mean for you individually is a question for a clinician, specifically a hematologist, and not something to resolve from a VCF.

Everything here is measurement. A filtered VCF is evidence about which molecules were in a tube, and turning that into any statement about health requires a clinical laboratory, orthogonal confirmation, and a physician who can see the rest of you.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Daniel C. Koboldt. Best practices for variant calling in clinical sequencing. Genome Medicine, 2020. https://doi.org/10.1186/s13073-020-00791-w ↩

2. M. Garcia, Szilveszter Juhos, Malin Larsson, et al. Sarek: A portable workflow for whole-genome sequencing analysis of germline and somatic variants. F1000Research, 2020. https://doi.org/10.12688/f1000research.16665.2 ↩

3. Tyler Alioto, Ivo Buchhalter, Sophia Derdak, et al. A comprehensive assessment of somatic mutation detection in cancer using whole-genome sequencing. Nature Communications, 2015. https://doi.org/10.1038/ncomms10001 ↩

4. Donald Freed, Rafael Aldana, Jessica A. Weber, et al. The Sentieon Genomics Tools - A fast and accurate solution to variant calling from next-generation sequence data. bioRxiv (Cold Spring Harbor Laboratory), 2017. https://doi.org/10.1101/115717 ↩

5. HoJoon Lee, Ahmed Shuaibi, John Bell, et al. Unique k -mer sequences for validating cancer-related substitution, insertion and deletion mutations. NAR Cancer, 2020. https://doi.org/10.1093/narcan/zcaa034 ↩

6. Alexander J. Cole, Ying Zhu, Trisha Dwight, et al. Comprehensive analyses of somatic TP53 mutation in tumors with variable mutant allele frequency. Scientific Data, 2017. https://doi.org/10.1038/sdata.2017.120 ↩

7. Jared T. Simpson. Detecting Somatic Mutations Without Matched Normal Samples Using Long Reads. bioRxiv (Cold Spring Harbor Laboratory), 2024. https://doi.org/10.1101/2024.02.26.582089 ↩