Scientific Rationale

This page is about the tools this pipeline chose and why, not about pooled sequencing as a method. That method is established and well described elsewhere — see Futschik & Schlötterer 2010 for the estimation problem and Kofler et al. 2011a for the toolkit most people met it through.

What follows matters because a pooled allele frequency is a measurement, and every stage between the reads and that number can bias it. A choice that is harmless when you are calling genotypes in one individual is not necessarily harmless when the quantity of interest is a fraction.

Every run writes its own CITATIONS.md and references.bib beside the results, listing the tools it actually invoked with the versions that ran. Those are the citations to use in a methods section; this page is the argument for the choices.

Trimming: Trim Galore, then cutadapt

Trim Galore is an excellent adapter trimmer, and it is at its best when it is given the adapters rather than left to detect them. Its thread handling is the awkward part — --cores N runs rather more than N threads — which is why the pipeline works out the core allocation for you rather than passing a number straight through.

Running cutadapt afterwards costs time the pipeline would otherwise save. It is worth it because trimming can leave the retained portion of a read skewed in base composition, and the clipping step exists to keep the A/T and G/C ratios near 1 across the length of the read. For frequency estimation a compositional bias is not cosmetic — it moves the numbers, and it moves them in a way nothing downstream can detect.

Alignment: bwa, with a fixed batch size

bwa is the standard, and there is little to argue about. One detail is worth stating because it is easy to lose: the pipeline passes -K 10000000, which fixes how many input bases each batch processes regardless of the thread count. Without it, bwa's output depends on how many cores it was given, so the same reads aligned on two machines can differ. With it, the alignment is reproducible.

Cleaning and calling: samtools and bcftools, not GATK

This is a considered preference rather than a default. Working through GATK and Picard in earlier versions of this pipeline, the read-group handling in particular was a recurring source of trouble; the samtools equivalents are less forgiving to use but they produce a cleaner result. Everything Picard was doing in the earlier pipelines, here is now samtools.

The same goes for calling. Recent comparisons of germline callers put samtools/bcftools among the accurate ones and consistently the most efficient — short runtimes and low memory where the deep-learning callers need far more of both. That is the basis for using it here, and call quality matters more in this setting than it does when the output is a genotype: an allele frequency inherits every error in the underlying counts rather than rounding it away.

Those comparisons do not measure what this pipeline needs, and it is worth being clear about that. They benchmark genotype calls on a single diploid individual against a truth set. This pipeline wants allele frequencies from a pool, and precision on a genotype call says little directly about the bias in a frequency derived from allele depths. The choice here rests on efficiency, on the AD handling being what the rest of the pipeline reads, and on long familiarity with the tool's failure modes — not on a benchmark that measured a different quantity.

The cleaning step is deliberately strict: anything that is not a properly paired read is discarded. That costs a little depth, which is the right trade when the quantity being measured is a ratio.

Major-allele reassignment

After calling, each site is re-polarized so that the reference allele is the one most common across the whole cohort, rather than whichever allele happened to be in the reference genome.

The reason is that in pooled data the reference allele carries no special status. It is an accident of which individual was sequenced to build the assembly. Re-polarizing on cohort totals makes the frequency tables consistent with each other and stops downstream analyses from treating an arbitrary allele as a baseline. One consequence follows and is worth knowing: this is a property of the cohort, so a frequency from a six-pool run is not the same quantity as a frequency from a twelve-pool run.

Depth, and the false-positive filter

These two are simple to state and take a moment to sit with.

Each pool has a detection limit set by how many chromosomes went into it: 1 / (ploidy × poolSize) is the smallest allele fraction that pool can physically produce — one chromosome out of all of them. The limit is set to half of this value, (1 / (2 x ploidy × poolSize)) as a statistical threshold. Anything below it is noise rather than a rare allele, because the pool contains nothing that could generate it.

The filter therefore applies a threshold per pool, not one threshold for the whole file. A pool of ten individuals and a pool of a hundred have different resolutions, and judging both at the same cut-off either discards real variation in the large pool or keeps noise in the small one.

Clearing that limit in one pool is not by itself a reason to keep a site, so a second choice sits beside it: sampleThreshold, the proportion of pools in which an allele has to clear its own detection limit before the site survives. It defaults to 0.2, and it is converted to a count — with six pools, an allele must be detectable in at least two of them.

This one is deliberately yours to set, because the right value depends on what you are asking. A high threshold keeps variation that is shared across the cohort and discards anything private to one population — which is correct if you are estimating differentiation over standing variation, and exactly wrong if you are looking for a sweep or a locally adapted allele, since the sites carrying your signal are the ones a high threshold removes first. A low threshold keeps those, and admits more noise with them. The default is permissive on purpose: it is easier to notice that a table is noisy than to notice that the site you were looking for was discarded before you saw it.

The point that takes longer to absorb is that neither of these is a confidence filter. It is not saying an allele is probably wrong. It is saying that the experiment as designed could not have detected it, so a number reported there describes the sequencing rather than the population. Depths are recomputed alongside it because once the surviving allele set changes, the counts a frequency is derived from have to change with it.

Annotation: snpEff

Standard, optional, and unchanged in approach. It runs only when you ask for it, and a run that did not annotate does not cite it.

The analysis layer

The modules that read these tables — what each one estimates, what it assumes and what it cannot tell you — are documented with the modules themselves in Modules. This section will grow as they do.