basicstats

What is in your published tables, per pool. It is the module to run first: everything it reports is a description of the data rather than a test of anything, and the numbers it prints are the ones every other module weights by.

PoolSeqFlow analysis basicstats

It reads the frequency and depth tables, and it needs your metadata.csv — the pool sizes, the ploidy and the experimental design all come from the project as the pipeline recorded it. A project whose metadata file was never copied across is refused rather than described as having no design.

design.tsv — one row per pool

Column	What
pool	the RG_Sample a column of every published table is named after
libraries	the SampleIDs merged into it, ;-separated — one for most pools, more where reads were merged
n_libraries	how many, so a merged pool is visible without reading the list
one per exp_ column	that pool's value for each of your experimental variables, in the order metadata.csv gives them; blank where the cell was
pool_size	individuals in the pool: its param_poolSize, or the project's poolSize where that cell was blank
ploidy	copies of the genome per individual, from parameters.config
n_chrom	ploidy × pool_size, the number of chromosomes sampled. This is the only place ploidy enters any statistic here
detection_limit	1 / (2 × ploidy × pool_size) — the frequency below which step 7's filter took a call to be error. Not a property of this analysis: it is what already happened to the tables, repeated here so the numbers beside it can be read against it

Nothing in this file is estimated. It is the design the results were produced under, printed where an analysis can be read against it — which is the point, because a frequency means nothing without the pool size it was read from.

A pool of one chromosome is refused, and it is the only combination that is. A single haploid genome has no segregating sites, and the correction every diversity estimate here applies divides by n_eff - 1, which is zero there. Correct ploidy, or that pool's param_poolSize.

sites.tsv — what the pipeline called

One row per sequence per kind of variant. It carries no pool, because it is a property of the tables rather than of any one column.

Column	What
chrom	the sequence, in the order the tables give them — never sorted, so chr10 sits where your reference puts it rather than between chr1 and chr2
kind	snp or indel. Step 7 publishes the two as separate tables and nothing here pools them
sites	rows of the depth table: records, not alleles. A record carrying three alleles is one site
alleles	rows of the frequency table: one per allele including the reference, so a triallelic site contributes three

A sequence with no surviving variants has no row here at all, and is indistinguishable from one your reference does not have. That is a limit of reading the published tables: they hold what was called, and nothing records what was looked at and found invariant.

depth.tsv — depth at the called sites

One row per pool per sequence, over the SNP table.

Column	What
pool, chrom	which pool, which sequence
sites	SNP sites on that sequence
depth_mean, depth_median	the ordinary summaries of that pool's depth
depth_harmonic	the harmonic mean, which is the one every effective sample size below is computed from

A pool's depth at a site is the sum of its cell in the depth table — the reads supporting any allele there, after step 7's depth, quality and false-positive filters. It is not coverage, and it is not what Output/Reports/Depth measured: the sites here are the ones that survived calling, every one of them carries at least vcffilter.minDP reads in every sample by construction, and mapping and base quality minima applied to the pileup that they did not. The two numbers are not one quantity measured twice, and this one is always the larger.

diversity.tsv — gene diversity and effective sample size

One row per pool, over the called SNP sites of the whole project.

Column	What
pool	the RG_Sample
n_chrom	ploidy × pool_size, as in design.tsv
sites	the called SNP sites
segregating	how many of them are segregating for this pool, by the rule below
depth_harmonic	the harmonic mean of this pool's depth over those sites
n_eff_harmonic	the pool's effective sample size over them, from depth_harmonic
h_sum	the numerator: the sum over sites of the corrected gene diversity
pi_per_called_site	h_sum divided by the sites that contributed a value to it. That is sites, unless a cell of the depth table was missing — which a published table cannot hold, since vcffilter.minDP removes a site where any sample falls short

The diversity is Nei's, over every allele at a site, corrected for the pool's effective sample size at that site:

\[\hat\pi = \frac{1}{\text{sites}}\sum_{\text{sites}} \frac{n_{\text{eff}}}{n_{\text{eff}} - 1}\left(1 - \sum_j p_j^2\right), \qquad n_{\text{eff}} = \frac{n_{\text{chrom}} \cdot d}{n_{\text{chrom}} + d - 1}\]

Not 2p(1-p): that form assumes two alleles and a privileged reference, so a triallelic site has to be collapsed or dropped before it can be used. \(1 - \sum_j p_j^2\) needs no such choice, is defined for any number of alleles, and reduces to 2p(1-p) where there are two. See Nei 1973 for the statistic, Ferretti et al. 2013 for why a pooled estimate needs correcting at all, and Hivert et al. 2018 for the effective sample size — including which of the two forms in circulation this is.

h_sum and sites are published, and the ratio is published beside them, because the ratio is the part that is contestable. Per-site diversity summed over the sites that were called is not the same as per-base diversity over the sites that were callable, and our tables hold only the former. A project's real π is h_sum over the number of positions that could have produced a call, which needs a per-position pass over the ready BAMs that no module does yet. pi_per_called_site is therefore an upper bound on π and will typically be far above it — invariant sites are missing from the denominator entirely. Quote h_sum and the count, not the ratio, unless you mean per called site and say so.

Our number will not match PoPoolation's or poolfstat's on a project with triallelic sites, and the reason is deliberate. Kofler et al. 2011a implements estimators defined over two alleles, and poolfstat's reader skips multiallelic records outright; the diversity here sums over every allele the site has. On a strictly biallelic project the two agree. On a viral quasispecies, where the three- and four-allele sites are the signal, discarding them discards the result.

A site is segregating for a pool when an allele other than that pool's own major one reaches max(detection_limit, minReads / depth). Two limbs, because either alone is wrong: detection_limit is 1 / (2 × ploidy × pool_size), which at any ordinary depth is below one read and so admits every sequencing error; minReads / depth alone stops discriminating once the pool is large enough that one chromosome is rarer than a couple of reads. The crossover is at depth = minReads / detection_limit. The pool's own major allele and not the cohort's: a pool fixed for whatever the cohort calls alternate is not segregating, and reading the majority off the reference column would report that it is.

neff.tsv — effective sample size, at two levels and from two sources

How many independent chromosomes a frequency read off this data is actually worth. Everything else the analysis layer will ever weight by is this number, so it is reported at every level the data supports rather than collapsed to one.

Column	What
level	library — one sequencing run, your SampleID; or pool — one column of the published tables, your RG_Sample
id	the library or the pool
pool	the pool it belongs to, so a library row can be read against its pool's
source	histogram — every position the library covered, from step 5; called — the sites that survived filtering, from the depth table
positions	how many positions that figure was measured over
depth_harmonic	the harmonic mean depth over them
n_chrom, n_eff	the chromosomes sampled, and what they are worth at that depth
estimate	exact, or lower_bound — see below

The two sources are not one quantity measured twice, and called is always the larger. Five things separate them: the histogram is measured before the depth ceiling and the calls after it; samtools stats counts every spanning read where the pileup applied a mapping and base quality minimum; a histogram counts positions where the depth table sums AD; every called site carries at least vcffilter.minDP reads in every sample by construction, which censors the low end away; and a site had to be variable to be called at all. Compare a pool's two rows to see how much the filter chain concentrated the data — do not treat their difference as error.

There is no library row from the called source, and there cannot be. The published tables carry one column per RG_Sample, so a merged pool's libraries are already summed inside them and nothing can separate them again. Per-library figures come from the histograms or from nowhere.

A merged pool's genome-wide figure is a lower bound

When two libraries were merged into one pool, its depth at a position is their depths added — and the harmonic mean of that sum cannot be recovered from the two histograms. Adding the histograms would assume the libraries are independent across positions, which is false and errs optimistic. So the module adds the parts' harmonic means instead, and marks the row lower_bound.

That is sound rather than arbitrary: the harmonic mean is concave and homogeneous of degree one, therefore superadditive, so \(H(a+b) \ge H(a) + H(b)\) always. Measured over 100,000 positions:

the libraries	H(a)	H(b)	H(a) + H(b)	true H(a+b)
exactly proportional	29.98	17.99	47.97	47.97
realistically correlated	29.88	17.37	47.25	47.94
independent (unrealistic)	30.14	17.93	48.07	55.62

Exact when the libraries cover positions in proportion to one another, about 1.4% conservative under realistic correlation, and never optimistic. It understates precision, which is the safe direction — a test weighted by it is conservative, not anti-conservative. It does assume the libraries covered the same positions, which is why each row carries the count it was measured over.

A pool of one library is marked exact: there is nothing to add.

A pool with no histogram gets no row from that source, and none is guessed. DepthProfile skips a sample whose ceiling is already decided, so a sound project can lack them; the run names those pools rather than averaging over the libraries that did have one.

depth_<sequence>.png — depth along a sequence

One file per sequence you name, one panel per pool, a point per called site.

params {
    analysis {
        modules {
            basicstats {
                chromosomes = ['chr2L', 'chr3R']
            }
        }
    }
}

Nothing is drawn until you name something. A genome has more sequences than anyone wants plots for, and which of them is worth looking at is your question rather than one this could answer for you. With the setting empty the run lists the candidates and their called-site counts, so a name can be copied straight into the file.

A named sequence with no called site is refused, naming what the tables do have. A sequence on which nothing survived filtering is invisible in the published tables and cannot be told apart from one your reference does not contain, so an empty plot would be a claim the data does not support.

It is depth, not coverage, and the axis says so. Every point is a site that was called; the gaps between them are sites the pipeline did not call, which is not the same as sites with no reads. A run of thin points is a stretch where calling was sparse — read it against depth.tsv for that sequence rather than as a coverage trace.

What it can be set to

params {
    analysis {
        modules {
            basicstats {
                minReads    = 2         // alternate reads before a site counts as segregating
                binSize     = 100000    // sites handed to one worker at a time
                workers     = 0         // 0 means the cores Nextflow gave the task
                usecpp      = true      // the compiled path, above
                chromosomes = []        // sequences to draw a depth plot for
            }
        }
    }
}

minReads is the only one that changes a number, and chromosomes is the only one that decides whether a file appears. binSize, workers and usecpp decide how the work is divided and which implementation does it: every combination of them produces the same output, which is what the test suite asserts by running the corpus through all of them and comparing the published tables byte for byte.

The compiled path

basicstats computes depth and gene diversity through site_diversity.cpp, compiled on your machine, and that is the default. The other implementation is vectorized R, which needs nothing at all:

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

For a whole project, analysis.modules.basicstats.usecpp = false in basicstats.config does the same thing permanently. Both produce the same numbers; what the choice costs is in How long the per-site work takes, along with what to do when the compiled path will not build.

site_diversity.cpp is published in the results folder whether or not the run used it, and the header of basicstats.R beside it names the path that produced the numbers, the bin size and the number of workers.

What to cite

Every folder basicstats publishes carries a CITATIONS.md and a references.bib with these already filled in, and PoolSeqFlow analysis cite basicstats prints them without running anything. Cite the statistics, not only the software — what a reader needs to check is the estimator, and the estimator is not the program that ran it.

Nei 1973	the gene diversity itself, summed over every allele at a site rather than over two
Ferretti et al. 2013	why a pooled estimate needs correcting at all: individuals are sampled into the pool, then reads are sampled from it
Hivert et al. 2018	the effective sample size everything is weighted by, including which of the two forms in circulation this is

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