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.
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.
- How precisely a frequency was measured. For one allele at one site in one unit, the variance is
p(1−p) x (dispersion + 1/n_eff)(read sampling and pool sampling throughn_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. - Units are independent. A unit is one independent biological entity; two pools are two units unless
analysis.design.technicalRepsays they are one material measured twice. - 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_<phenotype>.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 but 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. There is no Bonferroni and 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 1 / n!, however strong the signal. The observed arrangement is always among the rearrangements counted, so the smallest p is the one where it beats every other outright.
| units | smallest p | under 0.05? |
|---|---|---|
| 3 | 0.17 | no |
| 4 | 0.042 | yes |
| 5 | 0.0083 | yes |
| 6 | 0.0014 | yes |
| 8 | 0.000025 | yes |
At three 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.
A site read in fewer units than the design holds is rearranged among the units it was read in, so it has the floor of a design that size: a site read in three of six cannot reach 0.05 either. n_observed in association_<phenotype>.tsv is that count, and the run says in a sentence how many tested sites it applies to.
Ties keep a site above the floor. It is reached only by a site whose observed arrangement beats every other outright. Where every unit carries the same weight, as equal depths and pool sizes give, some arrangements always tie with the observed one: the reversal of a phenotype symmetric about its mean, and every reshuffle inside a level of a categorical one, with the swap of the two levels when they are the same size. Three cases against three controls read to equal depths cannot go below 0.1. design_floor in permutations.tsv is 1 / n!, and floor beside it is the smallest p this run could report, which differs only when the rearrangements were sampled.
association_<phenotype>.tsv: one row per site¶
One file per phenotype, named after it: fitting pt_wingspan and pt_resistant publishes association_pt_wingspan.tsv and association_pt_resistant.tsv. A character in a phenotype's name other than a letter, a digit, ., _ or - becomes _ in its file name, and two phenotypes whose names would then be the same are refused. To compare phenotypes at a site, join their files on kind, chrom and pos.
Always complete: every site the module could read on the sequences it fitted (every sequence, unless chromosomes names some) has a row, whether or not it produced a statistic.
| Column | What |
|---|---|
kind, chrom, pos |
which table, and where |
k |
how many alleles the site holds, the reference included |
n_observed, n_units |
units read at this site, and units in the design. With fewer than three read, or none of them differing in the phenotype, the site has no statistic |
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 weight of the units read at this site: their effective sample size when dispersion is 0, and less once a dispersion is estimated |
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.
alleles_<phenotype>.tsv: one row per allele¶
The effect sizes, one file per phenotype beside its site table, for the sites reportBelow and reportTop select from that phenotype's: every site at a perm_p of reportBelow or under, and the first reportTop by perm_p, ties going to the larger S with a flagged site first. That is the order the report ranks them in, so with reportTop at 20 or more every site it shows has its alleles. 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. The published association.R names the selection in its header.
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 binary scale, in the units at the level phenotype.tsv names as coded 1, and for an ordinal one, further along its order, which is read as each level's place in it.
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. NA where no site was read in two units, which leaves nothing to test either |
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 |
the median of S² over the tested sites, against a χ² on one degree of freedom. It is not 1 on a clean table; what it reads there is below |
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 |
lambda_gc reads above 1 when nothing is wrong. S is not a one-degree χ² even at a null site: read in m units it is a t on m − 2 degrees of freedom, and at a site of several alleles it is the largest of several. On a perfectly calibrated null it reads:
- 1.21 at six units, 1.29 at five, 1.47 at four and 2.20 at three
- about 3.8 at three alleles and 6 at four, read in six units
A table mixes these, so read lambda_gc against the mix of units and alleles in yours, not against 1. With vcffilter.keepLowDepthAsZero on, sites move to fewer units and it rises with nothing wrong. None of this touches perm_p, which is read against rearrangements of each site's own units, not against a χ².
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 for each sequence you name in analysis.modules.association.chromosomes, which also restricts the fit to those sequences. A reference holding scaffolds, unplaced contigs or an organelle would otherwise put all of them in the test; naming the chromosomes keeps them out. A run that names none fits every sequence and draws no plot.
Restricting changes the numbers, not only the plots. The dispersion is estimated from the named sequences' sites and fdr_p adjusts across them, so a site's fdr_p from a run over two chromosomes is not the one a run over the whole genome would give it.
A flat ceiling across one of these is your design's floor and not a feature of the data.
The report: the results laid out¶
The report reads the tables above and nothing else, in the order this page asks you to read them. The phenotype as it was read comes first, each pool's cell beside the number the fit used. Then what the design can show: for each phenotype, SNPs and indels apart, the units, the rearrangements, both floors, the sites read, tested and selected, and the diagnostics from permutations.tsv. When the design's floor is above 0.05, the report says so in bold before any site is shown. Every phenotype of a run shares that floor, because the units alone set it.
Then the strongest sites of each phenotype: the twenty with the smallest permutation p, or every tested site when there are fewer, with ties going to the larger S. Beside each is the allele with the largest |t| and its slope, read from the phenotype's alleles_<phenotype>.tsv, so a site that table did not select shows NA there. The caption counts how many tested sites reach the smallest p the run could report. qq.png comes last, and after it any Manhattan plots in the order of the sequences they draw.
A flagged site prints flagged where its S would be, and comes first among the sites it ties with. Its t is not comparable, so the allele named for it is the one with the steepest slope. A biallelic site's two alleles carry one |t| with opposite slopes, and the report names the one whose slope is positive: the allele that rises with the phenotype, more frequent at the level coded 1 of a binary scale or further along an ordinal one. An allele longer than 20 bases is shown by its first 12 and its length.
What it can be set to¶
| Setting | Default | What |
|---|---|---|
phenotypes |
no default | which pt_ column to fit. Required |
permutations |
10000 | the budget. The set is enumerated whole while it fits inside this and the design has eight units or fewer, and sampled otherwise |
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 of each phenotype reach its allele table |
chromosomes |
none | the sequences to fit, each drawn as a Manhattan plot; none fits every sequence and draws none |
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.