Filtering & Frequency Calculations

Step 7 takes the joint VCF that step 6 produced and turns it into the frequency tables you actually read. Between the two sit the filters that decide which of the called sites survive.

The order is fixed (stages 5 to 9 of The Filter Chain): normalize to the major allele, drop alleles too few samples support, drop sites too thin or too poorly called, split SNPs from indels, then convert allelic depths to frequencies. Only two of those stages have parameters, and both are analysis-affecting.

Parameter	What it does
poolSize, ploidy	Set the smallest allele frequency worth believing, which is what the false-positive filter tests against. Both can be set per sample or per run
filterFalsePositives.sampleThreshold	Removes alternate alleles without support across enough samples. The filter that makes the pipeline pool-aware, and the one most worth understanding before changing anything
vcffilter.minDP, vcffilter.minQUAL	Remove whole sites that are too poorly called, or too shallow in any one sample unless keepLowDepthAsZero is on
vcffilter.keepLowDepthAsZero	Off by default. On, a cell below minDP no longer removes its site: it is written as unread, and the site is kept if enough of its cells reach the floor
vcffilter.minSamples	With keepLowDepthAsZero on, how many of a site's cells must reach minDP for it to be kept. 2 by default, from 1 to the number of pools with reads
vcf.fileName	Removes nothing; it names the files this step writes, after outputPrefix unless you give it a name of its own

filterFalsePositives.sensitivity is computed from poolSize and ploidy rather than set, and can be overridden per sample in metadata.csv; see Metadata for the per-sample form.

For the stage-by-stage account with counts, read The Filter Chain.

poolSize and ploidy

poolSize = 50    // individuals in one pool
ploidy   = 2     // copies of the genome per individual

These two are the pipeline's most consequential settings. They are not used to model anything; they set the minimum credible allele frequency:

\[S = \frac{1}{2 \times \text{ploidy} \times \text{poolSize}}\]

An alternate allele must reach S in a sample for that sample to count as supporting it.

poolSize is per pool, not per run. With eight pools of 50 individuals each, it is 50, not 400.

The factor of two is deliberate. A pool of 50 diploids holds 100 chromosomes, so one chromosome is frequency 0.01. S comes out at 0.005 (half of that), so a genuine singleton clears the threshold with margin rather than sitting exactly on it. The observed fraction of a singleton is itself noisy, and a threshold placed exactly at the expected value would reject half of them.

poolSize	ploidy	One chromosome is	S
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
50	1	0.0200	0.0100
50	4	0.0050	0.0025

If your pools differ in size, give each one its own param_poolSize in metadata.csv and every column is judged at its own threshold. The global poolSize here stays as the default for pools that do not state one. See Pool size belongs to the pool.

ploidy applies to a whole run. It describes the sequence rather than the sample, so a diploid animal's haploid mitochondrial genome belongs in a separate run with its own value. See Sequences with a different ploidy.

sampleThreshold

filterFalsePositives {
    sampleThreshold = 0.2
}

The fraction of samples that must independently support an allele at frequency S or above for it to be kept. This is what makes the filter a corroboration test rather than a frequency cutoff, and it is what allows S to be set so low without drowning in sequencing error: errors do not recur at the same position across independent libraries.

The comparison is count >= n_samples × sampleThreshold, so the effective requirement is the next whole number up:

Samples	M at 0.2	Must support
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

The default removes population-private alleles

With eight samples, an allele found in one pool is discarded regardless of how frequent it is there: one does not reach two. If private or population-specific variation is the subject of your study, this default is wrong for you.

Value	Behavior	Suits
0.2 (default)	Needs roughly a fifth of samples	Shared, evolving variation across comparable pools
Low enough to require 1 sample	Any single pool can carry an allele	Private variants, small sample counts, discovery runs
Higher, e.g. 0.5	Needs half the samples	Conservative core-variant sets; high-noise data

Remember that the denominator is the number of VCF columns, which is the number of distinct RG_Sample values and not the number of FASTQ pairs. Merging replicates by sharing an RG_Sample changes this threshold as a side effect (see Metadata).

Depth and quality

vcffilter {
    minDP              = 20
    minQUAL            = 30
    keepLowDepthAsZero = false
    minSamples         = 2
}

By default both are site-level filters, applied as two commands in sequence:

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

minDP → bcftools view -e "FMT/DP<N || FMT/DP==0"

Removes a site if any sample falls below the depth or has no reads at all. Read the negation carefully: the site survives only when every sample meets the floor.

minQUAL → vcftools --minQ

Removes sites whose QUAL falls below the value.

keepLowDepthAsZero, minSamples → bin/mask_depth.awk

With keepLowDepthAsZero = true the depth test moves from the site to the cell, as the walkthrough describes. A cell below minDP, or with no reads, is written as unread, and its site is kept when at least minSamples of its cells reach minDP. minSamples runs from 1 to the number of pools with reads, since a pool the metadata names but no reads reach has no column in the VCF, and a run asking for anything else is refused before it starts. While keepLowDepthAsZero is false, minSamples has no effect, though changing it still counts as a change to a project's parameters.

The weakest library sets the threshold for every site

Because the default test is "any sample below minDP", one under-sequenced pool removes sites for all of them. Three pools at depths 50/15/55, 60/12/28 and 10/12/20, with minDP = 20:

minDP = 20  →  0 of 3 sites kept   (poolB is 15/12/12, so it fails everywhere)
minDP = 12  →  2 of 3 sites kept
minDP = 20, keepLowDepthAsZero = true, minSamples = 2  →  2 of 3 sites kept, poolB unread at both

Check Output/Reports/Coverage/ for your weakest sample before choosing a value, and after the first run compare the site count in <name>.vcf against the distinct positions that reached the frequency tables. A near-total wipeout is this filter, not a broken pipeline.

If "every sample" is too strict for your design, turn on keepLowDepthAsZero rather than editing the expression in the step. A site filter alone, such as -i "COUNT(FMT/DP>=20)>=2", keeps the same sites but leaves the shallow cell's reads in the depth table that every analysis module reads, where they count as a measured frequency; the mask writes that cell as unread.

Output naming

outputPrefix = 'Test'

vcf {
    fileName = "${params.outputPrefix}"
}

outputPrefix starts the names your results carry. The VCFs and frequency tables take it through vcf.fileName: with the default, a finished run leaves Test.vcf, Test_snp_freq.tsv, Test_indel_freq.tsv, and Test_annotated.vcf if annotation ran. Every table, figure and report an analysis publishes starts with it too (File names). Reads, BAMs and the per-sample reports keep their samples' names. Worth setting to something descriptive, because it is the name your results carry from here on. It may hold letters, digits, ., _ and -, starting with a letter or a digit, and a run naming anything else is refused before it starts.

vcf.fileName can name the VCF on its own. As the template writes it, it follows outputPrefix, and in a run table it follows each run's own. Write a name there instead, such as fileName = 'calls_v2', and the VCF and its tables carry that while every analysis file keeps outputPrefix.

A run table that varies outputPrefix varies the VCF name with it, and runs whose VCFs are named differently call their variants separately, from the same BAMs. To label each run's analyses and still call once, write a fixed fileName in parameters.config: the runs then share the calling, and each prefix gets its own filtering and tables.

Both are recorded with the rest of your parameters, so changing either after a run stops the next run before it starts, the way any recorded change does, and nothing is renamed. Restore the name the results were produced under, or start again with PoolSeqFlow reset.