The Filter Chain

A variant that reaches a frequency table has come through nine stages spread across four steps — six that remove something and three that only rewrite or split. A read is subject to the first three; after that the chain works on sites and alleles. This page walks the whole chain in order: what each stage removes, which parameter controls it, and what you are trading when you move that parameter.

If you are trying to work out why a variant you expected is missing, read this page top to bottom — the answer is usually earlier in the chain than people look.

The chain at a glance

#	Stage	Operates on	Removes	Parameter
1	Alignment filter	Reads	Unmapped, non-paired, duplicate, secondary, supplementary, low-MAPQ reads	cleanBAM.filter, cleanBAM.required, cleanBAM.mapq
2	Depth capping	Reads at a position	Reads above a ceiling measured from the sample's own coverage	capBAM.maxDepth or param_capMaxDepth
3	Pileup filter	Reads at a position	Low mapping quality, low base quality; optionally caps depth again	variantCall.baseQualMin, variantCall.varQualMin, variantCall.maxDepth
4	Variant calling	Sites	Non-variant sites	variantCall.callOptions
5	Major-allele normalization	Allele order	Nothing — it rewrites	—
6	False-positive filter	Alternate alleles	Alleles without cross-sample support	poolSize or param_poolSize, ploidy, filterFalsePositives.sampleThreshold
7	Depth & quality filter	Sites	Sites where any sample is under-covered, and low-QUAL sites	vcffilter.minDP, vcffilter.minQUAL
8	SNP/INDEL split	Sites	Splits into two files; nothing is lost	—
9	Frequency conversion	—	Nothing	—

Stages 1–4 happen in steps 4, 5 and 6; stages 5–9 are the five sub-steps of step 7.

---

1. Alignment filter (step 4)

The last stage of BAM cleanup is a samtools view with three conditions:

samtools view -F 0xF0C -f 0x2 -q 30 -b

Flag	Source	Effect
-F 0xF0C	cleanBAM.filter	Excludes unmapped, mate-unmapped, secondary, QC-fail, duplicate and supplementary reads
-f 0x2	cleanBAM.required	Requires the read to be properly paired
-q 30	cleanBAM.mapq	Excludes reads with mapping quality below 30

The MAPQ floor is easy to miss because it is not part of either flag word. At 30 it is a strict filter — it discards reads that map ambiguously, which in a repetitive genome can be a substantial fraction. That is usually the right call for Pool-seq, because an ambiguously placed read contributes a read count to the wrong position and frequencies are read counts. But if your coverage reports from step 5 show much less depth than you sequenced for, this is the first place to look.

Duplicate removal happens just upstream (samtools markdup -r) and matters more here than in individual sequencing: a PCR duplicate is a second vote from a molecule that should only vote once, and in a frequency estimate every vote counts directly.

Full flag reference: Alignment & Cleaning.

2. Depth capping (steps 5 and 6)

Step 5 measures a per-position depth histogram for every sample and publishes it to Reports/Depth/. From that histogram it decides a ceiling for that sample, and step 6 applies the ceiling to the BAM before calling. Nothing is deleted: a position deeper than the ceiling is truncated to it, and the sample's ready BAM is untouched on disk.

Why the ceiling is measured rather than set. Pooled coverage is uneven by design, so there is no depth that is "too deep" in the abstract. What is worth truncating is a second population of positions at high depth — a collapsed repeat, or a PCR hill, where reads belonging to several places in the genome stack up on one. Those positions carry read counts that no single locus produced, and a frequency is a read count. A single hand-set number cannot find them: it cuts legitimate coverage in a deep sample and lets the pile-up through in a shallow one, and nothing in the output says which happened.

What the detector looks for. One rise and fall is a library. A rise, a fall to nothing, and a second rise is two populations, and the ceiling goes where coverage ran out:

Blue is the sample before capping, green after. Everything above the ceiling collapses onto it, which is the spike at 447. Both axes are logarithmic and both have to be — the second population here holds two thousandths of the covered genome, so on a linear count axis it is invisible.

The same shape at a different scale is the same decision. Organelle contamination is a small number of positions two orders of magnitude above the library:

A long tail is not a second population. This is the distinction the whole stage turns on, and it is why nothing here uses a mean, a median or a maximum — the maximum is precisely the thing being cut. An overdispersed library reaching 3000× has one population, and is left alone:

So is a library on a reference that most of the genome does not map to. Two thirds of the covered positions here sit at depth 1–5. That is not coverage and it is not an anomaly either, and the ceiling has to stay above the real lobe at 60× rather than below it:

The three outcomes, and the third is the one to read. Under the default capBAM.maxDepth = -1 a sample is capped where the detector finds a second population, and left uncapped where it does not — including when the histogram is one the detector cannot read. Every sample's decision is published as a sentence beside its histogram, whether it was capped or not, because "nothing was done" is the one outcome you cannot see in the output.

The detector declines deliberately when the deep population is too large to call an artefact:

Four fifths of the covered genome is at 20000× and one fifth at 200×. Which of those is the artefact cannot be read off a depth histogram, and capping at 500 would truncate most of a real genome on a guess. PoolSeqFlow reports it and does nothing. If you know which population is real, param_capMaxDepth in metadata.csv sets that sample's ceiling by hand — see Metadata.

The knob.

capBAM.maxDepth	What happens
-1 (default)	Measure a ceiling per sample from its own histogram; leave the sample uncapped where there is nothing to cut
a positive number	Cap every sample at that depth
0	Do not cap at all

param_capMaxDepth overrides it for one sample and takes the same three values. The histogram is published on every run whichever you choose — it costs little and it is what tells you whether the sample needed capping.

3. Pileup filter (step 6)

bcftools mpileup -B -C 50 -q 30 -Q 30 -d 0 -a AD,DP,SP,INFO/AD -Ou

Flag	Parameter	Effect
-q 30	variantCall.varQualMin	Minimum mapping quality for a read to be counted
-Q 30	variantCall.baseQualMin	Minimum base quality for a base to be counted
-C 50	variantCall.scaleMapQ	Downgrades mapping quality for reads with excessive mismatches
-d 0	variantCall.maxDepth	A second, flat cap on reads per file per position. 0 is no limit
-B	fixed	Disables BAQ (base alignment quality) recalculation
-a AD,DP,SP,INFO/AD	fixed	Emits the allelic-depth fields everything downstream depends on

variantCall.maxDepth is not the same setting as capBAM.maxDepth, and the two zeros do not mean the same thing. capBAM.maxDepth = 0 means do not cap; variantCall.maxDepth = 0 means impose no ceiling, because that is what -d 0 means to bcftools mpileup. Both ship as a state where nothing is truncated flatly, and the measured per-sample ceiling from stage 2 does the work instead.

It ships as 0 because a flat number applied to every sample is what stage 2 exists to replace. Setting it to a positive value restores a backstop under the measured ceiling — one number for every sample of every run, applied after capping. If you set it, compare it against your step 5 coverage reports first: a cap that bites truncates the read counts frequencies are computed from, and nothing downstream flags it. See maxDepth.

-B is a deliberate choice for pooled data. BAQ downweights bases near indels to suppress false positives that arise from misalignment in a single diploid genome. In a pool, the same signal may be a genuine low-frequency indel, and BAQ's correction assumes a genotype model that does not apply. Disabling it keeps the raw evidence and leaves the decision to the cross-sample filter at stage 6.

4. Variant calling (step 6)

bcftools call -m -A -v -Ov

Flag	Effect	Why it matters here
-m	Multiallelic caller	Handles sites with more than one alternate allele, which biallelic-only calling would collapse
-A	Keep all alternate alleles from the pileup	Without this, bcftools discards alternates it considers unlikely under a genotype model — exactly the low-frequency alleles a pool is meant to detect
-v	Variant sites only	Invariant sites are dropped

-A is the flag that makes this a Pool-seq caller rather than a general one. The default behavior prunes alternate alleles that no plausible genotype supports, which is sound for an individual and wrong for a pool, where a true allele at frequency 0.01 supports no genotype at all.

5. Major-allele normalization (step 7)

bin/MajorAlleleToRef.py re-encodes the VCF so the most-read allele is the reference. It does not remove anything; it rewrites.

For each site it sorts the alleles by INFO/AD — the read count summed across every sample — and reorders REF, ALT, INFO/AD, INFO/DP, FORMAT/AD and FORMAT/DP to match. DP is recomputed as the sum of the reordered AD, so depth and allelic depth cannot disagree.

Why do it at all. The reference genome is one individual's assembly. There is no reason its allele should be the common one in your population, and when it is not, every frequency in that row is reported against a rare baseline. Two studies on the same species then report mirror-image frequencies for the same site. Normalizing to the major allele makes rows comparable across samples, across runs and across projects.

Two consequences worth knowing.

The ordering is cohort-wide, not per-sample. REF is the allele most read across all samples combined. In an individual sample where the cohort-minor allele is locally dominant, that sample's frequency for REF will be below 0.5 — which is meaningful, not an error.

FORMAT/GT is set to ./. on every genotype. A pool has no genotype, and leaving bcftools' diploid call in place would invite downstream tools to read it as one. Making it explicitly missing is the honest encoding. Keep it in mind if you point a genotype-based tool at these VCFs: it will find nothing, by design.

The script runs twice — once before the false-positive filter and again after it, because splitting and rejoining multiallelic records can change allele order.

6. False-positive filter (step 7)

This is the filter that makes the pipeline pool-aware, and the one most worth understanding before you change anything.

bcftools norm -m -  vcf                    # 1. split multiallelic sites into one line per ALT
| bcftools view -i "INFO/AD[1]>0"          # 2. drop alleles with no supporting read at all
| awk  …                                   # 3. count samples that clear their own threshold
| sed  '*' → 'X'                           # 4. mask the spanning-deletion allele
| bcftools norm -m+                        # 5. rejoin into multiallelic records
| sed  'X' → '*'

What the filter keeps

An alternate allele survives only if both hold:

INFO/AD[1] > 0

The allele has at least one supporting read somewhere in the cohort.

At least M samples show it at or above their own threshold S

Each sample is judged against the threshold for the pool it belongs to, not against a single number for the run.

where, for a given pool,

\[S = \frac{1}{2 \times \text{ploidy} \times \text{poolSize}} \qquad M = n_{\text{samples}} \times \text{sampleThreshold}\]

This is a cross-sample corroboration filter, not a frequency cutoff. A site is not kept because it is frequent; it is kept because several independent pools saw it. That distinction is what allows S to be set very low without drowning in sequencing error: errors are random and do not recur in the same place across libraries, whereas real low-frequency variants do.

Why this is not one bcftools expression

Earlier releases did the counting inside bcftools view -i, which can only apply one threshold across every sample. Per-pool sensitivity does not fit in that expression: a per-sample threshold would have to become a per-sample FORMAT field, annotated onto every record and stripped off again afterwards.

Doing the counting in awk instead buys a property worth more than the tidiness it costs. Each threshold is bound to the sample's name, not to its column position. A column that changes place — because you reordered metadata.csv, which you are free to do — still gets its own pool's threshold. Binding by position would have applied the wrong pool's threshold after a reorder and produced a perfectly ordinary-looking table.

It also fails loudly rather than guessing. If the VCF contains a sample column that your pool sizes do not mention, the filter stops rather than judging that column at some other pool's threshold.

What was applied, recorded in the VCF

Each filtered VCF carries a header line per pool, stating the size used and the threshold derived from it:

##PoolSeqFlowPool=<ID=Sample1T1,PoolSize=50,Sensitivity=0.005>

These are provenance, written beside the data they were applied to. Note that the filtered VCFs are intermediates and most do not survive a completed run — the durable record of your pool sizes is .poolseqflow_metadata, kept beside the results.

Where S comes from

A pool of poolSize individuals, each carrying ploidy copies of the genome, contains \(\text{ploidy} \times \text{poolSize}\) chromosomes, so a single chromosome carries a frequency of \(1 / (\text{ploidy} \times \text{poolSize})\). The extra factor of two puts S at half that — a true singleton clears the threshold with margin rather than sitting exactly on it, which matters because the observed fraction of a singleton is itself noisy.

poolSize	ploidy	One chromosome is	S (threshold)
10	2	0.0500	0.0250
25	2	0.0200	0.0100
50	2	0.0100	0.0050
100	2	0.0050	0.0025
200	2	0.0025	0.00125

Larger pools give a smaller S, so the filter admits rarer alleles — appropriately, since a larger pool really can contain rarer ones.

Each pool can be set to get its own row of that table. param_poolSize in metadata.csv can set the size per pool, and any pool without one falls back to the global poolSize. ploidy can be set per run rather than per pool, so a design mixing ploidies is split into runs — see Sequences with a different ploidy.

Where M comes from

M is a count of samples, computed as a fraction of however many are in the VCF. It is compared with >= against a non-integer value, so the effective requirement is the next whole number up:

Samples in run	M at sampleThreshold = 0.2	Samples that must support the allele
1	0.2	1
4	0.8	1
5	1.0	1
6	1.2	2
8	1.6	2
10	2.0	2
12	2.4	3
20	4.0	4

This filter removes population-private alleles

At the default, an allele found in only one pool out of eight is discarded no matter how frequent it is in that pool — one sample does not reach the two the threshold requires. If private variation is what you are studying, lower sampleThreshold before your first run, not after. See sampleThreshold.

Why the splitting and rejoining

FORMAT/AD[:1] indexes the first alternate allele. On a multiallelic record, a rare third allele would never be tested — it would simply ride along on whatever the second allele did. Splitting to one line per alternate (norm -m -) makes each allele stand on its own evidence; norm -m+ puts the survivors back together.

The * → X substitution around the rejoin masks the spanning-deletion allele, which norm -m+ does not handle in this position. It is restored immediately afterwards.

7. Depth and quality filter (step 7)

Two commands, both operating on whole sites:

bcftools view -e "FMT/DP<20" -Ov -o <name>_dp.vcf <input>
vcftools --vcf <name>_dp.vcf --minQ 30 --recode --recode-INFO-all --out <name>_dq

bcftools view -e "FMT/DP<20" (vcffilter.minDP)

Removes a site if any sample falls below the depth. A site survives only when every sample meets the floor, so the weakest library sets the threshold for the whole cohort — one under-sequenced pool removes sites for all of them.

vcftools --minQ 30 (vcffilter.minQUAL)

Removes sites whose QUAL falls below 30.

The depth test has to be site-level rather than per-sample. Both vcftools and bcftools express a genotype-level verdict by rewriting FORMAT/GT and nothing else — AD and DP survive untouched — and stage 9 reads AD. Stage 5 has already set every GT to ./. besides, so a genotype filter would have nothing left to mark.

Alternative expressions, what each trades, and how to pick a value: Depth and quality.

8. SNP/INDEL split (step 7)

Two vcftools passes over the same input — --remove-indels and --keep-only-indels — produce a SNP VCF and an INDEL VCF. Nothing is discarded; every surviving site lands in exactly one of the two files, and each is converted to its own frequency table.

9. Frequency conversion (step 7)

No filtering. bin/createDepthFile.sh extracts CHROM, POS, REF, ALT, INFO/AD and per-sample FORMAT/AD, and bin/depth2freq.awk divides each allele's read count by the row's total to give a frequency. The output format is described in Interpreting Results.

---

Tuning the chain

Work from the outside in. A variant lost at stage 1 cannot be recovered by loosening stage 6.

Symptom	Most likely stage	Parameter to examine
Far less depth than sequenced	1	cleanBAM.mapq, then duplicate rate in the step 5 reports
Depth plateaus at one number in one sample	2	That sample's Reports/Depth/ report — a measured ceiling was applied
Depth plateaus at the same number in every sample	3	variantCall.maxDepth, or capBAM.maxDepth set to a fixed depth
A sample you expected to be capped was not	2	Its Reports/Depth/ report says why; param_capMaxDepth overrules it
Low-frequency alleles absent everywhere	6	poolSize, or param_poolSize for the pool in question, and ploidy
Low-frequency alleles absent from one pool only	6	That pool's param_poolSize — a size set too low raises its threshold alone
Alleles present in one pool only, absent from output	6	filterFalsePositives.sampleThreshold
Whole sites missing despite good depth	7	vcffilter.minQUAL
Almost every site gone after filtering	7	vcffilter.minDP — one under-covered sample removes sites for all of them
Multiallelic sites reduced to two alleles	4	variantCall.callOptions — confirm -A is still present

Changing any of these invalidates existing outputs, and step 0 will stop the next run rather than mix results. That is covered in Design Decisions.