Variant Calling

Step 6 is two things in order. First each ready BAM is capped at a depth ceiling; then the capped BAMs go into one joint pileup and bcftools calls across all of them at once. Both halves have parameters and they are easy to confuse, because both have one called maxDepth.

	Parameter	Applies to
First	capBAM.maxDepth	Each BAM separately, at a ceiling measured from that sample
Then	variantCall.maxDepth	The joint pileup, as one flat number for every sample

Everything here is analysis-affecting: change one and the next run stops rather than letting old and new results share a folder. What happens to the call set afterwards is Filtering & Frequency Calculations; for how the whole chain fits together, read The Filter Chain.

Capping each BAM

capBAM {
    maxDepth     = -1
    histogramMax = 100000    // how deep step 5's histogram looks
}

Value	What happens
-1 (default)	Step 5 measures a ceiling for each sample from its own depth histogram; step 6 truncates that sample's BAM to it. A sample with nothing worth cutting is left uncapped
positive N	Every sample is capped at N, measured or not
0	No capping at all — the BAMs reach the pileup as step 4 left them

param_capMaxDepth in metadata.csv overrides this for a single sample and takes the same three values, which is the way to handle one library that needs a different answer from the rest — see Metadata.

Nothing is deleted. A position deeper than the ceiling is truncated to it on the way into the pileup, and the sample's ready BAM is untouched on disk. Reads are dropped whole rather than trimmed, so a pair may lose one mate.

histogramMax is how deep step 5 looks, and a sample deeper than it stops the run. The histogram reaches 100000× by default, and everything above that would land in a single open bin — so a ceiling read from it would be read from a partial picture. Rather than choose on partial evidence, step 5 fails and names the sample and the depth it found. Raise histogramMax past that depth and run again. An organelle in a well-covered library is the case that reaches it in practice. Raising it costs nothing already produced: samtools stats reports only the depths that actually occur, so every value the run completes at gives the same histogram and the same ceiling, and step 0 does not compare it against previous runs.

Why the ceiling is measured rather than set. Pooled coverage is uneven by design, so no single depth is "too deep" in the abstract. What is worth truncating is a second population of positions — a collapsed repeat, or a PCR hill — carrying read counts no single locus produced. A hand-set number cannot find those: it cuts legitimate coverage in a deep sample and lets the pile-up through in a shallow one. The Filter Chain shows the histograms this is read off, including the two cases where the detector deliberately declines.

Every sample's decision is published, capped or not, beside its histogram in Output/Reports/Depth/ — "nothing was done" is the outcome you cannot otherwise see:

grep -H 'ceiling applied' Output/Reports/Depth/*_depth_report.txt

Pileup settings

variantCall {
    scaleMapQ   = 50     // -C  downgrade coefficient for mismatch-heavy reads
    varQualMin  = 30     // -q  minimum mapping quality
    baseQualMin = 30     // -Q  minimum base quality
    maxDepth    = 0      // -d  a flat depth cap on top of capBAM; 0 is no limit
}

-C 50 (scaleMapQ)

Downgrades mapping quality for reads carrying excessive mismatches. Reads that align poorly are more likely to be misplaced, and a misplaced read contributes its bases to the wrong position — which in a frequency estimate is a direct error, not just noise.

-B (fixed, not configurable)

Disables BAQ recalculation. BAQ downweights bases near indels to suppress misalignment artefacts under a single-genome model. In a pool, the same signal may be a genuine low-frequency indel, so the raw evidence is kept and the decision is left to the cross-sample filter.

maxDepth

-d caps the reads considered per file per position, and it ships as 0, which to bcftools mpileup means no limit. The BAMs reaching this point were already capped above, from each sample's own coverage — a second flat number on top of a measured one is what that stage exists to remove.

The two are not the same knob, and they do not even take the same values:

Value	capBAM.maxDepth	variantCall.maxDepth
-1	Measure a ceiling per sample	Not valid
0	Do not cap	No limit
positive N	Cap every sample at N	Cap the pileup at N

When to set it. Leave it at 0 unless you want a hard backstop under the measured ceiling — a memory limit on a shared machine is the usual reason, since a very deep pile-up costs mpileup memory whether or not it was capped. If you do set it, size it against your data rather than guessing: the depth reports named above say what each sample actually looked like.

A cap that bites truncates the read counts your frequencies are computed from, and nothing downstream flags it — which was the whole problem with the fixed 2000 this replaced.

Coming from 2.2.0 or earlier, your old bcftools.maxDepth is not carried across. It was the only depth control there was, and here it would sit as a flat ceiling under a measured one; migrate_config reports it under Format changed this release and explains the change. You get capBAM.maxDepth = -1 and variantCall.maxDepth = 0 — measured capping, no flat backstop. To reproduce older results exactly, put your old number back in variantCall.maxDepth and set capBAM.maxDepth = 0.

Calling

callOptions = "-m -A -v -Ov"

Flag	Effect
-m	Multiallelic caller — required for sites with more than one alternate allele
-A	Keep all alternate alleles from the pileup
-v	Output variant sites only
-Ov	Uncompressed VCF

-A is what makes this a Pool-seq caller. Without it, bcftools prunes alternate alleles that no plausible genotype supports — sound for an individual, and wrong for a pool, where a true allele at frequency 0.01 supports no genotype at all. If you edit callOptions, keep -A and -m.

Variant calling is a single joint task over all BAMs, not one task per sample. That is what produces a multi-sample VCF with comparable columns, and it is why every sample must share a reference.