Skip to content

association

Each allele's frequency regressed on a phenotype you measured per pool, weighted by how precisely each pool's frequency was measured, with a p-value read off a permutation rather than a table.

PoolSeqFlow analysis association

It reads the depth tables and your metadata.csv. The phenotype is a pt_ column, declared in analysis.metadata.phenotypes and chosen for this module in analysis.modules.association.phenotypes.

What it assumes, and why that sentence comes first

This is a model-based test. It is not an assumption-light one, and there is no assumption-light one available: any test that weights pools by their precision has already conceded that they differ, and a null that then treats them as interchangeable contradicts its own statistic. So the assumptions are stated rather than hidden, and there are three.

  1. How precisely a frequency was measured. For one allele at one site in one unit, the variance is p(1−p) × (dispersion + 1/n_eff) — read sampling and pool sampling through n_eff, plus whatever excess uneven pooling left behind. Only the ratio between units is ever a claim: multiplying every weight by a constant leaves the result unchanged, so a pool size wrong by the same factor everywhere costs nothing, and one wrong differently per unit costs calibration.
  2. Units are independent. A unit is one independent biological entity; two pools are two units unless analysis.design.technicalRep says they are one material measured twice.
  3. A site's alleles are several views of one comparison, which is what the next section is about.

Where those do not describe your experiment, the numbers will still be produced and will still be wrong. permutations.tsv is what tells you so, and if your design is one this does not fit, please come and say so — the roster grows by someone describing a problem this cannot answer, and fitting your data to the nearest available analysis serves nobody.

A site with more than two alleles

Every allele is fitted, the reference included, and the site's answer is the strongest of them. Nothing is collapsed, nothing is dropped, and there is no correction for how many there were.

That is worth unpacking, because it is where this module departs from what most tools do.

One regression per allele. A triallelic site gives three fits — reference against the phenotype, first alternate against it, second alternate against it. Each allele has its own frequency and therefore its own variance under assumption 1 above; the weights are the same for all of them, because a weight describes how deeply that unit was read at that site and depth is a property of the site.

The site statistic is the largest |t| among them. Written S in association.tsv.

Why the reference row is not skipped. A site's frequencies sum to 1, so the reference carries the negated sum of every alternate. When a phenotype pushes two alternates up together, neither is impressive alone and the reference falls by both — so the reference row is where that signal is, and an analysis that only looks at alternates cannot see it. Measured against a planted signal of exactly that shape, reading the alternates alone recovers 40% fewer of the sites than reading every allele.

MajorAlleleToRef.py has already made the reference the cohort's major allele, so "reference" here is not a privileged biological state — it is simply the commonest allele, and skipping it would privilege "not the commonest" for no reason.

Why more alleles is not an unfair advantage. A site of four alleles offers four chances at a large |t| where a biallelic site offers two, which sounds like it should favor multiallelic sites. It does not, because the p-value is not read off a table: it is read against a null built by rearranging that same site, which has the same alleles, the same constraint that they sum to 1, and therefore the same number of chances. The comparison is like for like at every arity, so nothing needs correcting — no Bonferroni, no union bound.

That is not an argument from first principles alone. Measured on null sites of two, three and four alleles, the false-positive rate is flat across arity, and the sites the module selects have the same allele counts as the genome it selected them from — which is what arity_mean and arity_selected in permutations.tsv let you check on your own data. The union-bound alternative — take the smallest allele p and multiply by k − 1 — measured 37% above its nominal rate at triallelic sites, which is why it is not used.

What it costs. A site of k alleles carries k − 1 free tests, not k: at a biallelic site the reference and alternate rows are the same test with the sign flipped, and their p-values are identical. So one site is one finding, whatever its arity, and the allele rows tell you which allele carries it rather than giving you several findings.

The smallest p your design can reach

Before any of the output: a design of n units cannot report a p below 2 / n!, however strong the signal. Reversing a phenotype negates every slope and leaves the statistic alone, so the reversed arrangement always ties with the observed one and two of the rearrangements are always counted.

units smallest p can a site be significant at 0.05?
4 0.083 no
5 0.017 yes
6 0.0028 yes
8 0.00005 yes

At four units the run says so in a sentence and the table is a ranking of effect sizes, not a test. That is a true statement about the experiment rather than a fault of the method, and it is the honest answer to a well-known failure of this kind of study: a modestly powered pool-seq scan read as though it were well powered produces scattered genome-wide significance from a single causal variant, which looks like a polygenic architecture and is noise Long et al. 2026. A floor is a structural refusal to draw that picture.

The floor binds hardest when the pools are read to similar depths. Where depths differ a great deal the design carries more information than a bare count of units suggests, and the run can reach further down; design_floor and floor in permutations.tsv are those two numbers side by side.

association.tsv — one row per site

Always complete: every site the module could read has a row, whether or not it produced a statistic.

Column What
phenotype, kind, chrom, pos which fit, which table, and where
k how many alleles the site holds, the reference included
n_observed, n_units units carrying data here, and units in the design
S the site statistic: the largest |t| over every allele
perm_p the published p, from the permutation
fdr_p perm_p corrected across sites, by default Benjamini–Hochberg Benjamini & Hochberg 1995
mean_weight the average effective sample size at this site
max_leverage how much of the slope one unit is carrying, which at six of them is the question a reader has
zero_variance 1 where an allele had no residual left; see below

A flagged site's S is not comparable between implementations and must not be quoted. Two alleles of a perfectly separated site are algebraically one test, and whether each lands on exactly zero residual or on 1e-32 decides between an infinite t and a merely enormous one — we have seen two builds of the same language disagree. perm_p is unaffected and remains the answer.

association_alleles.tsv — one row per allele

The effect sizes, for the sites reportBelow and reportTop select. This table is a selection and the site table is not: at a hundred million sites the allele rows are a twenty-gigabyte file, which is a data dump rather than a result.

b1 is the slope, se its standard error, and t and p the closed-form diagnostics. b1 > 0 means the allele is more frequent at the higher phenotype value — for a categorical scale, at the level phenotype.tsv names as coded 1.

The parametric p is a diagnostic and not the answer. It is well calibrated on an ordinary site and makes impossible claims on a degenerate one — six reads of evidence will produce a p of 1e-9 that no design of six pools can support.

permutations.tsv — what the run assumed and what it absorbed

A guard that passes tells you nothing. These are published so that someone reading your results, without rerunning them, can see what the model was asked to do on your behalf.

Column What it means, and when to worry
permutations, exhaustive how many rearrangements, and whether that was all of them
floor, design_floor the smallest p this run could reach, and the smallest the design can
dispersion the excess variance uneven pooling left behind, in units of p(1−p). Large means the weights absorbed a lot; it is estimated from your data unless you set it. A published 0 is a floor, not necessarily a measurement — the method of moments can return a negative excess when the units scatter less than their own sampling variance predicts, which is noise and not a unit measured better than its depth allows, and it is reported as 0
depth_phenotype_cor depth lined up with the phenotype. This breaks label-based tests badly and is the reason this module does not use one
lambda_gc genomic inflation. If it is 3, nothing in the table is a p-value
arity_mean, arity_selected, selected the allele counts of the sites selected against all sites. They should agree; multiallelic sites at the top of a table are a biological claim, and this is how you check it is one

phenotype.tsv — the phenotype as it was read

One row per pool: the cell as you wrote it, the number the fit used, the unit it belongs to, and for a categorical scale which level became 1.

Read this first. A reversed case/control coding has no more of a check than dd/MM against MM/dd does, and printing the resolved answer is the only place either is caught.

qq.png — the p-values against the uniform they should follow

Always drawn. A permutation p is discrete, so the points sit on a ladder rather than a line, and the ladder's ceiling is your design's floor. A curve bending above the diagonal is inflation; the diagonal itself is what a null genome looks like.

manhattan_<sequence>.png — association along a sequence

Drawn only for the sequences you name in analysis.modules.association.chromosomes, on the same rule basicstats follows: a genome has more sequences than anyone wants plots for.

A flat ceiling across one of these is your design's floor and not a feature of the data.

What it can be set to

Setting Default What
phenotypes which pt_ column to fit. Required
permutations 10000 the budget. The set is enumerated whole while it fits inside this, and sampled above it
fdr BH anything p.adjust takes
dispersion estimated the excess variance term. 0 weights by n_eff alone
reportBelow, reportTop 0.05, 1000 which sites reach the allele table
chromosomes none which sequences get a plot
binSize, workers, usecpp 100000, task cpus, true how the work is divided and which implementation does it

Enumerating is a correctness requirement, not a speed one. A sampled p is (1 + reached) / (1 + draws) Phipson & Smyth 2010, and sampling a set small enough to enumerate can return a p below the smallest value your design supports. Leave permutations above n! for small designs.

The compiled path

The per-site parse runs through allele_frequencies.cpp, compiled on your machine, and that is the default:

PoolSeqFlow analysis association          # compiled, the default
PoolSeqFlow analysis association nocpp    # plain R, for one run

Both produce the same numbers, which the test suite asserts by comparing the published tables byte for byte. allele_frequencies.cpp is published beside your results whether or not the run used it, and the header of association.R names the path that produced the numbers.

What to cite

Every folder association publishes carries a CITATIONS.md and a references.bib with these already filled in, and PoolSeqFlow analysis cite association prints them without running anything. Cite the statistics, not only the software.

Hivert et al. 2018 the effective sample size each unit is weighted by, and which of the two forms in circulation this is
Phipson & Smyth 2010 why a sampled permutation p is (1 + reached) / (1 + draws) and never the raw share
Benjamini & Hochberg 1995 the correction across sites, applied to the permutation p and never to the per-allele one
Long et al. 2026 what a modestly powered pool-seq scan produces when it is read as though it were well powered, which is why this module prints a floor rather than assuming one

Full entries, with what each is cited for, are in the Bibliography.