How many reads do I need for variant calling?

It depends on the variant class and the target. For SNPs the packet guidance is about 20x per base, with 15x described as cost-effective for whole-genome SNVs. For indels, plan on 50x or more across targets. Check the fraction of bases at or above 20x, not just the mean.

Is more depth or more replicates better for variant calling?

Past saturation, replicates win. Merging replicates gave better accuracy than one high-depth replicate, with Ti/Tv of 2.00 versus 1.39 for a single high-depth library. More depth from the same library mostly re-samples the same fragments, so it does not fix uneven capture.

How do I make a saturation or rarefaction curve for variant calling?

Downsample the BAM at several fractions with samtools view -s, call variants at each with identical filters, and plot variant counts against mean depth. If the last step adds few new variants, you are on the plateau. Do it separately for SNVs and indels, and for your genes of interest.

My mean coverage is 30x but I am missing variants. Why?

Mean depth hides low-coverage exons. Roughly 7 to 11% of genes show problematic coverage on common exome platforms, linked to GC content, repeats and segmental duplications. Compute per-target depth over your BED and look at which exons fall below 20x.

Sanity check · Variant Calling

How to Tell If You Sequenced Deep Enough in Variant Calling

Mean coverage looks fine, the VCF looks clean, and you still cannot say whether the missing variants are absent or just unsampled.

By Ming "Tommy" Tang, Director of Bioinformatics in Big Pharma · Reviewed October 2026 · 5 min read

You got a VCF back from an exome or genome run, the report says mean 30x, and someone asks whether the sample needs to be re-sequenced. Or you are planning the next batch and have to decide between more reads per sample and more samples. The mean answers neither question.

What is at stake: shallow or uneven coverage fails quietly. The failure is false negatives. The variants you did call stay mostly reliable, so the VCF looks fine while a causal variant in a poorly covered exon simply never appears. A "no variant found" result then gets read as "no variant present".

This page gives you a one-hour workflow. You will compute coverage per target instead of per sample, build a downsampling curve to see which side of saturation you are on, and decide between more depth, more replicates, or neither.

What it looks like when it's happening

The sequencing report says mean 30x, but a gene you care about has zero or near-zero calls across its exons.
samtools depth over your target BED shows a long tail of positions at very low depth next to a few positions at several times the average.
Adding more reads to a re-run barely changes the number of variants called, but the same few exons stay empty.
Indel calls are far sparser than SNP calls in the same sample, and indel sites look ragged in IGV with few supporting reads.
Ti/Tv ratio of a single low-depth sample sits well below what you expect from a good call set, much lower than a merged call set from the same individual.
Heterozygous sites appear as homozygous reference or homozygous alt in one replicate and heterozygous in another.
A somatic caller reports more low-VAF calls in your deepest sample than in shallower ones, and many are C to T or G to A changes.

Why it happens

Variant calling is sampling. At each position the caller sees a handful of reads drawn from two chromosomes (germline) or from a mixture of cells (somatic). With few reads, a heterozygous site can be sampled from one allele only and called homozygous, or missed entirely. That is why the GIAB depth benchmark on HG002 finds precision stable across depths while recall climbs with coverage, and why genotype errors such as heterozygous allele dropout account for most false positives, not spurious sites.

Returns diminish fast for SNVs. Published numbers: SNV concordance exceeds 95% at 17.6x and 99% at 13.7x in ultra-deep whole-genome benchmarks, and about 15x is described as cost-effective for SNVs. Indels are harder. They need roughly 50x or more across targets, and indel performance stays below SNP performance at every depth tested. Going from 40x to 60x costs 2.5x the sequencing for about the same absolute gain as going from 2x to 10x.

The mean hides where reads land. In exome data, capture is uneven: 7 to 11% of genes show problematic coverage across NimbleGen, Agilent and Illumina TruSeq platforms, tied to high GC content, repeats and segmental duplications. Longer exons are more uneven. A sample can hit mean 30x while a GC-rich first exon sits at a depth where nothing can be called. More reads from the same library re-sample the same favored fragments, so this tail barely moves.

Past saturation, extra depth can also hurt. Somatic callers (VarScan, SomaticSniper, Strelka, MuTect2) show rising false positive discovery rates at higher coverage, so you need tighter filtering above 50x. Regions with more than 2x the average depth are enriched for artifacts, and GATK HaplotypeCaller can call nothing at positions with tens of thousands of reads. FFPE adds damage on top: deamination produces C to T and G to A artifacts that look like low-VAF variants, and fewer than 90% of the exome reaches 20x.

The checks

Run them in order. Each one tells you what healthy looks like and what the problem looks like.

0/8 checked · saved in this browser

01Read the depth report as a distribution, not a mean
Take the coverage summary from your pipeline or compute it yourself, and look for the fraction of target bases at or above 20x, not the mean. For SNP calling the rule of thumb in the packet is 20x per base; for indels, 50x or more across targets. Write down the fraction of target bases that clear 20x.
Healthy
A large majority of target bases above 20x, with the fraction close to the mean-based expectation, and a mean depth that sits near the median depth.
Red flag
Mean looks healthy but the fraction of bases at or above 20x is low, or the median is far below the mean. A few extremely deep regions are pulling the mean up.
02Compute per-target coverage over your BED
Run samtools depth with -a over the capture or gene-panel BED so zero-depth positions are included, then summarize per position. Count how many target bases sit below 20x and which exons they belong to. Zero-depth positions must be in the denominator or you will overestimate coverage.
bash
```
samtools depth -a -b targets.bed sample.bam \
  | awk '{n++; s+=$3; if($3>=20) ok++} END{printf "mean=%.1f frac>=20x=%.3f\n", s/n, ok/n}'
```
Healthy
The fraction of target bases at or above 20x is high and the low-coverage positions are scattered, not clustered in your genes of interest.
Red flag
Low-depth positions cluster in specific exons or genes, especially ones you were asked to evaluate. Those variants cannot be called, whatever the mean says.
03List the genes and exons that fail, and check why
Join the positions below 20x back to the BED and aggregate by exon and gene. Then check whether failures line up with GC content, repeats or segmental duplications, which are the known drivers of problematic exome coverage. Also check exon length, since longer exons are more uneven.
Healthy
A short list of failing exons that are explained by high GC, repeats or segmental duplications, and that are not on your gene list of interest.
Red flag
Failures concentrated in clinically relevant genes, or failures that repeat in every sample from the same capture kit. That is a library or panel property, not a one-sample accident.
04Look at depth outliers on the high side
Flag regions with depth above 2x the regional average. The bcftools documentation treats these as enriched for artifacts. Inspect a few in IGV for pileups of reads with mapping quality problems or clipped alignments. If you use bcftools mpileup, remember the default --max-depth is 250 reads per file and may truncate high-coverage projects.
Healthy
Few regions with depth more than 2x the local average, and calls in those regions look like ordinary heterozygous or homozygous sites in IGV.
Red flag
Hotspots of very high depth with ragged pileups, or a caller that returns no variants where the pileup clearly shows an alternate allele. For HaplotypeCaller this appears at extreme depth and is a caller limit, not absence of signal.
05Build a downsampling curve of calls versus reads
Subsample the BAM at several fractions with samtools view -s (seed.fraction, for example 42.25 for 25%), re-run your caller on each, and plot the number of high-confidence variants, and separately the number in your target genes, against mean depth. Use the same filters at every point. This is the rarefaction curve for variant calling.
bash
```
for f in 25 50 75; do
  samtools view -b -s 42.${f} sample.bam > sub_${f}.bam
  samtools index sub_${f}.bam
done
```
Healthy
The curve flattens: the last step from 75% to 100% of reads adds very few new SNVs. You are on the plateau and more depth from this library buys little.
Red flag
The curve is still climbing at 100%, which means you are on the steep side and more depth would help. Or it flattens but a specific gene stays empty at every fraction, which means the problem is capture, not depth.
06Separate SNPs from indels in the curve
Split your downsampling results by variant type. Indels need more depth than SNVs, so the two curves saturate at different points. Compare the point where each curve stops rising to the 20x (SNP) and 50x or more (indel) guidance.
Healthy
SNV counts plateau early and indel counts plateau later or lower, consistent with indels being the harder class.
Red flag
You are reporting indels from a run that was only deep enough for SNVs. Indel recall at low depth in the GIAB benchmark is a small fraction of what it is at 60x.
07Check Ti/Tv and genotype consistency
Compute the transition to transversion ratio of your SNVs and compare across replicates or against an earlier call set. In the replicate-merging study, merged replicates reached a Ti/Tv of 2.00 plus or minus 0.13, while single high-depth and low-depth replicates were 1.39 and 1.22. If you have two libraries from the same individual, compare genotypes at shared sites.
Healthy
Ti/Tv near that of a good call set for your data type, and genotype agreement between replicates that is high at sites where both have adequate depth.
Red flag
Ti/Tv well below a merged reference, or heterozygous sites flipping to homozygous between libraries. That is allele dropout from under-sampling, not biology.
08For FFPE or somatic work, check the artifact signature
Tabulate the substitution types of your low-VAF calls. Deamination in FFPE gives an excess of C to T and G to A changes, mostly below 10% VAF. Compare their share with an unaffected control or a fresh-frozen sample if you have one. Also check whether the number of low-VAF calls rises with depth.
Healthy
Low-VAF calls spread across substitution classes, and call counts that stabilize as depth rises once filters are applied.
Red flag
A spike in C to T and G to A at VAF below 10%, and a call count that rises with depth. Extra depth is amplifying damage and noise, not finding mutations.

What to do about it

Re-capture or re-library instead of sequencing the same library deeper

When: Per-target coverage shows a stable set of exons below 20x across the downsampling curve, and the mean is already above what SNV calling needs.

Treat the failing exons as a library property. Prepare a new library, ideally on a platform or kit with better behavior in your GC-rich or repeat-rich targets, and sequence it to moderate depth. Then merge it with the original at the variant-calling stage so both contribute reads.

Caveat: Costs a new library and sample material. Some regions, such as segmental duplications, remain hard on any short-read capture.

Add replicates and call them together

When: Your curve is flat, genotypes disagree between libraries, or you care about rare variants.

Sequence a second independent library from the same individual, then run multisample or merged-replicate calling instead of picking one high-depth replicate. The replicate-merging study found merged calls beat a single high-depth replicate on accuracy, with Ti/Tv of 2.00 versus 1.39, and that joint calling reduced false positives for rare variants.

Caveat: Requires enough DNA and budget for a second library. Replicates share any systematic bias from the sample itself, such as FFPE damage.

Add depth only when the curve is still rising

When: The downsampling curve is climbing at 100% of reads and your fraction of bases at or above 20x is low, or you need indels and are well under 50x.

Plan the extra depth against the target: about 20x per base for SNPs, 50x or more across targets for indels. Re-run the downsampling curve after the top-up to confirm the plateau instead of assuming it.

Caveat: Cost per added variant rises quickly. Going from 40x to 60x needs 2.5x more sequencing for the same absolute recovery as 2x to 10x. Somatic callers need tighter filtering above 50x.

Raise caller limits where depth is extreme

When: Hotspots with very high read counts return no variants, or bcftools truncates at the default depth cap.

For bcftools mpileup, raise --max-depth above the default 250 reads per file to match your project. For GATK HaplotypeCaller, adjust max-reads-per-alignment-start as the GATK community thread suggests for extreme-depth failures, then compare calls before and after in IGV.

Caveat: These issues show up almost only at extreme depth. Raising limits costs runtime and memory, and does not fix regions that are deep because of mapping artifacts.

Filter deep regions and FFPE artifacts instead of resequencing

When: The problem is excess false positives at high depth or low-VAF deamination calls, not missing calls.

Flag and review regions above 2x the regional average depth. For FFPE, filter low-VAF C to T and G to A calls, and use a matched normal rather than relying on tumor-only calling, which labels germline variants somatic when databases are incomplete.

Caveat: Hard filtering at low VAF can remove real subclonal variants. Validate a sample of removed calls in IGV.

When not to "fix" it

Do not add reads when the downsampling curve has already flattened and your fraction of bases at or above 20x is high. A remaining empty region is then a capture, mappability or biology question. Segmental duplications and GC-rich exons stay thin at any depth, and a gene with no calls may simply have no variant. Do not chase a variant count that rises with depth in a somatic or FFPE sample without checking the substitution spectrum first: the increase may be damage. Also leave a shallow sample alone if it is a screening run where only SNVs matter and the 15x or 20x guidance is met; spending 2.5x more sequencing for 40x to 60x buys little.

Five things experienced analysts do here

Ask for the fraction of target bases at or above 20x before you look at the mean. If a report gives only mean depth, compute the rest yourself with samtools depth -a over the BED.
Always name the variant class in a depth requirement. About 20x covers SNPs; indels need 50x or more across targets, so one number for both is wrong.
Make downsampling a routine step on pilot samples. A curve that has flattened is the evidence you need to say no to a request for more reads.
Treat a negative result as weaker than a positive one. Shallow data fails by missing variants, so check coverage at the exact position before you write that a variant is absent.
Compare per-target coverage across samples from the same kit. If the same exons fail every time, the fix is the capture design, and no per-sample re-sequencing will repair it.

Questions people ask

How many reads do I need for variant calling?: It depends on the variant class and the target. For SNPs the packet guidance is about 20x per base, with 15x described as cost-effective for whole-genome SNVs. For indels, plan on 50x or more across targets. Check the fraction of bases at or above 20x, not just the mean.
Is more depth or more replicates better for variant calling?: Past saturation, replicates win. Merging replicates gave better accuracy than one high-depth replicate, with Ti/Tv of 2.00 versus 1.39 for a single high-depth library. More depth from the same library mostly re-samples the same fragments, so it does not fix uneven capture.
How do I make a saturation or rarefaction curve for variant calling?: Downsample the BAM at several fractions with samtools view -s, call variants at each with identical filters, and plot variant counts against mean depth. If the last step adds few new variants, you are on the plateau. Do it separately for SNVs and indels, and for your genes of interest.
My mean coverage is 30x but I am missing variants. Why?: Mean depth hides low-coverage exons. Roughly 7 to 11% of genes show problematic coverage on common exome platforms, linked to GC content, repeats and segmental duplications. Compute per-target depth over your BED and look at which exons fall below 20x.

Guide · How to Catch a Sample Swap in Variant Calling
Guide · How to Catch a Genome Build Mismatch in Variant Calling
Guide · How to Detect Contamination in Variant Calling
Guide · How to Fix Chromosome Naming Mismatches in Variant Calling
Guide · How to Tell If You Sequenced Deep Enough in ATAC-seq
Glossary · Variant calling

Sources

GIAB Coverage Benchmark: How sequencing depth affects germline variant calling accuracy — Recall versus depth on HG002, stable precision, shallow depth fails by false negatives, genotype errors dominate false positives.
Empirical evaluation of variant calling accuracy using ultra-deep whole-genome sequencing data — SNV concordance at 17.6x and 13.7x; roughly 15x as cost-effective whole-genome depth.
Novel metrics to measure coverage in whole exome sequencing datasets reveal local and global non-uniformity — Problematic coverage in 7 to 11% of genes, linked to GC content, repeats and segmental duplications; exon length and unevenness.
Improved Variant Calling Accuracy by Merging Replicates in Whole-Exome Sequencing Studies — Merged replicates beat a single high-depth replicate; Ti/Tv values.
Comparing variant calling algorithms for target-exon sequencing in a large sample — 20x for SNP calling and 50x or more for indels.
In-depth comparison of somatic point mutation callers based on different tumor next-generation sequencing depth data — False positive rates rise with coverage in somatic callers above 50x.
bcftools variant-calling howto — Regions above 2x average depth, mpileup default --max-depth 250, samtools depth.
GATK HaplotypeCaller - High Depth Issues — HaplotypeCaller failures at extreme depth and max-reads-per-alignment-start.
DeepVariant GitHub Repository — Accuracy across depths and the 'twenty is the new thirty' remark.

Part of the Sequencing depth and saturation series.