When to use PoolSeqFlow¶
This page is about fit. It describes what pooled sequencing changes about an analysis, what PoolSeqFlow assumes your design looks like, and the cases where a different tool is the better answer. If you already know Pool-seq is what you want, skip to Getting Started.
What pooling buys and what it costs¶
In a pooled library, DNA from many individuals is combined before sequencing. Every read still comes from exactly one chromosome in one individual, but nothing in the data records which. What survives is the proportion of reads carrying each allele, which is an estimate of the allele frequency in the pool.
| You gain | You give up |
|---|---|
| Frequency estimates for many individuals at a fraction of the per-individual cost | Individual genotypes — and therefore heterozygosity, relatedness, and anything phased |
| Depth concentrated where it matters: on the frequency estimate, not on calling a genotype confidently | The ability to separate sampling noise from real low-frequency variation without care |
| A design that scales to large populations and many time points | Straightforward variant filtering — the usual genotype-based heuristics do not apply |
The second row is the one that shapes this pipeline. In individual sequencing, a variant seen on two reads out of a hundred is almost certainly an error. In a pool of 50 diploids, one chromosome out of 100 is a real allele at frequency 0.01, and it looks exactly the same. Distinguishing the two is not something a fixed cutoff can do, which is why PoolSeqFlow derives its threshold from your pool size and ploidy rather than hard-coding one. See The Filter Chain.
What PoolSeqFlow assumes¶
These are structural assumptions. If your data does not match, the pipeline will either refuse to start or produce something that is not what you meant.
| Assumption | Where it comes from | If you do not match |
|---|---|---|
| Paired-end Illumina reads | readPattern matches an R1/R2 pair; step 4 requires properly-paired alignments (0x2) |
Single-end data will not survive the pairing filter |
| One reference genome per run | Variant calling is a single joint bcftools mpileup over all BAMs |
Samples on different references cannot be called together. Running one set of reads against several references is a different thing, and supported — that is what a run table is for |
| One ploidy per run | ploidy sets the detection threshold for every pool in the run |
Split the work into runs, each with its own ploidy — see below |
| Several comparable pools per run | The false-positive filter keeps a variant only if a fraction of samples support it | A single-sample run needs sampleThreshold reconsidered — see below |
| A reference FASTA, and a GFF if you annotate | referenceFile, gffFile |
Either may be gzipped or plain; the pipeline takes both and unpacks what it needs |
| Compute and storage reachable from one process | mainDir and storageDir can be on different filesystems, and both must be mounted |
Cloud object storage without a filesystem mount is not supported |
Where PoolSeqFlow ends¶
PoolSeqFlow produces allele frequency tables. It deliberately stops there. It does not compute FST, run CMH or other tests for allele frequency change, generate sync/mpileup formats for PoPoolation, call structural variants or copy number, or perform any population-genetic modeling.
That is a scope decision, not an omission: the statistics you want depend entirely on your design, and a per-site frequency table is the input nearly all of them take. Annotation via SnpEff is included because it operates per-variant and needs the same reference build the pipeline already indexed.
Cases that need a second look¶
Pools of different sizes¶
Pool size feeds the sensitivity threshold:
A pool holds \(\text{ploidy} \times \text{poolSize}\) chromosomes, so a single one of them represents a frequency of \(1/(\text{ploidy} \times \text{poolSize})\); the threshold is set at half that. A pool of 10 and a pool of 500 therefore have very different limits, and judging both at one threshold is either too permissive for the small pool or too strict for the large one.
So give each pool its own. param_poolSize in metadata.csv is a column of pool sizes, and each pool's threshold is computed from its own. Because the size describes the pool rather than the row, rows sharing an RG_Sample must agree on it, and a blank cell means the global poolSize rather than agreement with anything.
The global poolSize in parameters.config stays as the default for any pool that does not state one, so a run where every pool really is the same size needs nothing extra.
Sequences with a different ploidy¶
ploidy applies to a whole run, and it belongs to the sequence rather than to the sample: a diploid animal carries a haploid mitochondrial genome, and one threshold cannot be right for both.
Separate them into runs. A run table can vary any parameter, so give the nuclear chromosomes and the organellar sequences their own references and their own ploidy:
RunID,referenceFile,gffFile,ploidy
nuclear,chromosomes.fasta.gz,chromosomes.gff.gz,2
mitochondrial,mito.fasta.gz,mito.gff.gz,1
Both runs read the same FASTQ files, and the work they share — trimming and quality control — is done once rather than twice. Changing ploidy re-derives that run's sensitivity threshold automatically, so there is nothing else to keep in step.
The same argument applies to anything else whose copy number differs from the nuclear genome: chloroplast sequences, and unplaced scaffolds whose ploidy you are not confident about. Splitting them out costs one row in the table and makes each threshold defensible.
Single-sample runs¶
The false-positive filter is a cross-sample consistency check: an allele is kept only if at least sampleThreshold (default 0.2) of the samples show it at or above \(f_{\min}\). With one sample, that fraction rounds to a requirement that the one sample support it, so the filter degrades to a plain per-site frequency cutoff. It still works, but it is doing much less than it does on a multi-sample run, and the cross-sample corroboration that justifies a permissive \(f_{\min}\) is gone.
Genuinely private variants¶
At the default sampleThreshold = 0.2, an allele present in only one pool out of eight (12.5% of samples) is removed, even at high frequency in that pool. This is the correct default for detecting shared, evolving variation; it is the wrong default if population-private alleles are the point of your study. Lower sampleThreshold accordingly, and see Variant Calling for what that costs you.
Very high depth¶
Deep is not a problem in itself: no flat depth ceiling ships any more. Each sample gets one measured from its own coverage, and a library that is uniformly deep is left alone — see Depth capping. What a very deep run does change is what an anomaly looks like, so read Output/Reports/Depth/ rather than assuming: a sample reported uncapped had nothing separable to cut, which on a deeply and unevenly sequenced library is worth knowing.
Choosing between run layouts¶
| If your samples are… | Do this |
|---|---|
| Independent pools you want compared | One run, one row per pool in metadata.csv, distinct RG_Sample values |
| One pool split across lanes or runs to reach depth | One run, one row per FASTQ pair, sharing an RG_Sample so the reads are combined |
| Technical replicates you want treated as one observation | Share an RG_Sample |
| Technical replicates you want to compare against each other | Distinct RG_Sample values |
| Pools of different sizes | One run — give each pool its own param_poolSize in metadata.csv |
| To be compared against different reference genomes | One project, one row per reference in a run table |
| Nuclear and organellar sequences together | One project, one row per ploidy in a run table |
Only the last two need a run table; everything above them is a single run.
The RG_Sample decision is the one most often got wrong by accident, because two FASTQ pairs from one pool look exactly like two ordinary samples. It is covered in full in Metadata.