gain.genomic_resources.statistics package

Submodules

gain.genomic_resources.statistics.alleles module

Allele-content statistics for allele-score resources.

Vocabulary per GLOSSARY.md and ADR 0020. Where the coverage statistic answers where a score holds data, this one answers what its rows are: how many alleles a resource carries, how those alleles distribute over the five allele classes, and what the classes with structure look like inside – the 4x4 ref->alt matrix for substitutions, the exact length maps for the two anchored classes, and the complex grid for the rest.

An allele score counts no covered positions at all (gain#1118). Its rows collapse to points, so the span union gain.genomic_resources.statistics.coverage computes never applied to the kind – which is why it is deliberately absent from that statistic’s kinds – and the DISTINCT-position count this module used to keep in its place answered a question the page never asked. What an allele score’s rows ARE is the whole of what this statistic says.

Raw counts only. Anything needing a denominator is computed at render time, as the coverage statistic’s fractions are.

Laid out like its coverage twin: the per-region accumulator, the resource-wide statistic, the fold that merges a scan’s regions into one, and the write. The scan wiring that feeds it is in implementations/genomic_scores_impl/scan.py.

ALLELE_FORMAT_VERSION = 1

The format_version serialize stamps.

ALLELE_INSERTION_LENGTHS_IMAGE_FILE = 'statistics/allele_insertion_lengths.png'

The global images the statistics build renders beside the file. One each, never per chromosome (ADR 0020): the per-chromosome numbers are stored as data and rolled up for the picture.

ALLELE_STATISTIC = ('statistics/alleles.json', 1, <function <lambda>>)

The statistic’s one declaration – file, version, and the kind it is written for – read by the scan’s gate and by the repair check.

class AlleleChromosomeRow(chrom: str, allele_count: int, shares: dict[str, ClassShare] | None)

Bases: NamedTuple

One chromosome’s allele counts, as the info page renders them.

The chromosome is a FIELD rather than a mapping key, so a row is self-contained exactly as the coverage and fragment rows are: the template reads a row, never a pair it has to keep together.

allele_count: int

Alias for field number 1

chrom: str

Alias for field number 0

shares: dict[str, ClassShare] | None

What this chromosome’s alleles are MADE OF, one entry per class.

Keyed by CLASS_NAMES, and None for a chromosome with no alleles to take a share of – the row then renders empty cells carrying no sort key, exactly as the Coverage table’s row with no resolvable denominator does. percentages_over() owns that rule; the difference is only where it is asked. Coverage resolves a denominator per ROW and so does this, while the classes total that used to sit below the table resolved one for the whole table.

class AlleleCounts(allele_count: int, class_counts: dict[str, int], substitution_matrix: dict[tuple[str, str], int] | None, insertion_lengths: ExactLengths | None = None, deletion_lengths: ExactLengths | None = None, complex_grid: dict[tuple[int, int], int] | None = None)

Bases: NamedTuple

One chromosome’s – or a whole resource’s – allele counts.

class_counts is keyed by the class names in CLASS_NAMES and sums to allele_count: every row classifies, other absorbing what does not parse as alleles.

substitution_matrix is keyed by the ref/alt cells in MATRIX_CELLS, upper-cased, and sums to the substitution class count: a cell holds the rows whose pair classifies as a substitution of that cell, and nothing else lands in any cell. It is None – data unknown, distinct from a matrix of zeros – when restored from a file written before the matrix existed.

allele_count: int

Alias for field number 0

class_counts: dict[str, int]

Alias for field number 1

complex_grid: dict[tuple[int, int], int] | None

Alias for field number 5

deletion_lengths: ExactLengths | None

Alias for field number 4

display() → AlleleDisplay | None

This entry’s render payload, None if EVERY group is unknown.

One seam for the whole Alleles section, with the groups independently optional inside it: collapsing the payload the moment any single group is missing would hide the groups a file does carry, and the rollout guarantees no particular combination. A genuinely empty group still renders – zeros for the matrix, an empty grid for the complex cells – which is why “unknown” and “empty” are different answers here.

Everything derived is computed here, at render time, and never stored: transitions are the four A<->G / C<->T cells and transversions the eight remaining OFF-DIAGONAL cells – never “substitutions minus transitions”, which would silently count the diagonal’s identity rows as transversions.

insertion_lengths: ExactLengths | None

Alias for field number 3

substitution_matrix: dict[tuple[str, str], int] | None

Alias for field number 2

class AlleleDisplay(substitution_matrix: dict[tuple[str, str], int] | None, transitions: int | None, transversions: int | None, ts_tv: float | None, substitution_percentages: dict[tuple[str, str], str] | None, insertion_lengths: ExactLengths | None = None, deletion_lengths: ExactLengths | None = None, complex_grid: dict[tuple[int, int], int] | None = None)

Bases: NamedTuple

The Alleles section’s render payload, one field per stored group.

Raw cells come from the stored statistic; the derived numbers are computed by AlleleCounts.display() – which builds this – and never stored, as the coverage display derives its fractions.

Every group is INDEPENDENTLY optional, because the statistics roll out lazily and a resource may have been rebuilt under any one of them: a file written between gain#778 and gain#779 carries a matrix and no lengths, and must render its matrix rather than losing the whole section. Each of the page’s sections therefore asks for its own group, and this payload exists at all whenever ANY group is known.

complex_grid: dict[tuple[int, int], int] | None

Alias for field number 7

property complex_grid_renders_as_table: bool

Whether the complex cells render as a table, not a heatmap.

False without a grid, which the page never asks: it gates on the grid itself first, as it does for the matrix.

complex_rows() → list[tuple[str, str, int, str]]

The occupied complex cells as table rows, most alleles first.

A row is reference length, alternative length, alleles and the share of the complex class – the lengths labelled as the heatmap’s axes label them, so a clamped cell reads ≥64 in both. Empty without a grid; the page gates on the grid itself rather than reading this to find out.

The shares come from percentages_over(), the one rule the Alleles section writes a share by, so a rare cell reads <0.01% and a cell that is all but the whole class reads >99.99%, here exactly as they do in the classes column. Its denominator is the grid’s own total, which the TOTAL clamp makes exactly the complex class count: every complex row lands in one cell, so these SHARES sum to 100% – the rendered strings, as in the matrix, round off it.

That denominator is zero only when no cell is occupied, and then there are no rows to carry a share anyway – so the helper’s “no percentage at all” answer and this method’s empty result are the same answer, and it is returned as one.

deletion_lengths: ExactLengths | None

Alias for field number 6

indel_rows() → list[LengthStatisticsRow]

The indel statistics table: one row per known group.

Global only. Per-chromosome indel statistics would add eight columns to the Alleles table to answer a question nobody asks; the comparison these exist for is between the two GROUPS – is the deletion tail longer than the insertion tail – and that wants them side by side.

A group the file does not carry is left out rather than shown empty; the page says “not computed” for it, as it does for every other unknown group.

insertion_lengths: ExactLengths | None

The three gain#779 groups, None when the file predates them. An empty grid is KNOWN and empty, which is not the same thing.

matrix_rows() → list[tuple[str, list[MatrixCell]]]

The matrix as table rows: reference base, then a cell per alt.

Empty without a matrix; the page gates on the matrix itself rather than reading this to find out. This only pairs each stored count with the share AlleleCounts.display() already computed for it – nothing is derived here.

property nucleotides: tuple[str, ...]

The matrix’s axis labels, in stored order.

substitution_matrix: dict[tuple[str, str], int] | None

None when the file predates the matrix (gain#778).

substitution_percentages: dict[tuple[str, str], str] | None

Each cell’s share of the substitution class, formatted by percentages_over(). The denominator is read off the matrix rather than taken from class_counts: the two are equal by AlleleCounts’s invariant, and dividing the cells by their own total is what makes the sixteen SHARES come to exactly 100%. The sixteen rendered strings do not: two decimals round independently, and the floor and ceiling round further still. None with no matrix, and with no substitutions to take a share of.

transitions: int | None

The four A<->G / C<->T cells; None with no matrix.

transversions: int | None

The eight off-diagonal cells that are not transitions; the diagonal’s identity pairs are neither. None with no matrix.

ts_tv: float | None

None when there are no transversions – the template renders “not applicable” rather than dividing – and with no matrix.

class AlleleSectionDisplay(rows: list[AlleleChromosomeRow], class_counts: dict[str, int], class_percentages: dict[str, str] | None, detail: AlleleDisplay | None)

Bases: NamedTuple

The Alleles section’s render payload: the table and its totals.

Built in the implementation layer, as CoverageDisplay and FragmentDisplay are, so the template renders fields off an inert record rather than calling methods on the statistic itself.

property allele_count: int

The table’s total, summed off the rows it shows.

Derived rather than stored, as CoverageDisplay. global_covered and FragmentDisplay.global_fragments are: the counts have one source, so a total cannot drift from the rows under it.

class_counts: dict[str, int]

Alias for field number 1

property class_names: tuple[str, ...]

The class columns, in the order ADR 0020 states them.

Reached for THROUGH the payload rather than as a template global, because the template layer deliberately does not import gain.genomic_resources – which is why natural_chromosome_key lives in gain.utils instead. The constant itself is what this returns; the property is the seam that keeps the import out of the template environment.

class_percentages: dict[str, str] | None

Each class name’s share of allele_count, or None.

Computed HERE rather than on AlleleDisplay because the two numbers a share needs – the class counts and the allele total – are stored fields that every file carries, while that payload collapses when every OPTIONAL group is unknown. Computing them there dropped this column for a file written before the matrix (gain#777) whose shares were perfectly resolvable.

None when there are no alleles to take a share of, which drops the column; percentages_over() owns that rule.

detail: AlleleDisplay | None

matrix, indel lengths, complex grid.

None when the file carries none of them; each group inside it stays independently optional, as AlleleDisplay explains.

Type:

The optional groups

rows: list[AlleleChromosomeRow]

Alias for field number 0

class AlleleStatistics

Bases: RegionFoldedStatistic[RegionAlleles]

A resource’s allele content, per chromosome and global.

Folds RegionAlleles the way the base class does, and serializes to the resource’s ALLELE_STATISTICS_FILE as raw counts.

by_chromosome() → dict[str, AlleleCounts]

The per-chromosome counts, in natural chromosome order.

Ordered here rather than downstream because this order REACHES the info page: build_allele_section_display() turns these entries into the Alleles table’s rows as they come, and does not re-sort them the way its coverage and fragment siblings sort theirs. What it replaces is not arrival order but the plain string sort regions_in_genomic_order() applies before the fold – exactly the order iossifovlab/gain#983 calls wrong.

serialize() reads it too, so the order reaches statistics/alleles.json. Nothing downstream reads that file positionally: deserialize() and _total() are both order-blind, and calc_statistics_hash hashes the config and the source files, never the statistics content.

static deserialize(content: str) → AlleleStatistics

Create a statistic from serialized data.

global_counts() → AlleleCounts

The roll-up over every chromosome.

Off the regions directly rather than through by_chromosome(): _total is order-blind, so the ordering that accessor does would be paid and thrown away.

repeated_allele() → RepeatedAllele | None

The resource’s first repeated allele key, in genomic order.

serialize() → str

Return a serialized version of this statistic.

CLASS_NAMES: tuple[str, ...] = ('substitution', 'insertion', 'deletion', 'complex', 'other')

The five class names, in the order ADR 0020 states them – which is the order AlleleClass declares them in, so this reads that order rather than restating it. The serialized class map is written in it, so two builds of one resource produce byte-identical JSON however the rows arrived.

COMPLEX_GRID_TABLE_MAX_CELLS = 32

How many occupied cells the complex grid may hold and still render as a table of those cells rather than as the heatmap (gain#989).

A judgement call, not a derivation: 32 rows is about where a table stops being scannable, and a 64x64 grid starts having enough lit cells to show shape. Unlike COMPLEX_LENGTH_CLAMP this is a RENDERING choice and no part of the stored format.

RAISING it is free: a resource built under the old value keeps a PNG its page no longer references, which is the leftover image save_allele_statistics() already documents. LOWERING it needs the resources rebuilt with --force: a grid between the two values would start asking for a PNG that was never written, and a plain repo-stats will not notice – the statistics hash covers the table config, the score definitions and the data files, none of which a constant here moves.

COMPLEX_LENGTH_CLAMP = 64

The longest allele length the complex grid resolves exactly. A length at or above it lands in the grid’s top row or column, which therefore reads “this many bases or more”. The clamp is TOTAL – every complex row lands in exactly one cell, so the grid’s total is the complex class count and no overflow counter is needed.

Exact lengths rather than the shared log2 ladder (gain#779): that ladder’s first bin is exactly length 1, which no complex pair can have – a 1->1 pair is a substitution – and its second bin is {2, 3}, which would put a 2->3 complex in the same cell as a 3bp MNV and empty the diagonal of its meaning. Part of the stored format: it must not change once resources carry grids built from it.

class ClassShare(alleles: int, percentage: str, fraction: float)

Bases: NamedTuple

One class’s slice of one chromosome’s alleles.

Three renderings of one count, because the cell shows one, sorts on another and titles itself with the third: percentage is the text, fraction is the sort key – the NUMBER, so that a column sorts by size rather than by the string <0.01% – and alleles is the exact count the hover title carries.

The count is what the global classes table used to show as a column of its own (gain#1118 removed it). Keeping it here is what lets a reader recover the exact number the share rounds off, per chromosome rather than only for the resource.

alleles: int

Alias for field number 0

fraction: float

Alias for field number 2

percentage: str

Alias for field number 1

MATRIX_CELLS: tuple[tuple[str, str], ...] = (('A', 'A'), ('A', 'C'), ('A', 'G'), ('A', 'T'), ('C', 'A'), ('C', 'C'), ('C', 'G'), ('C', 'T'), ('G', 'A'), ('G', 'C'), ('G', 'G'), ('G', 'T'), ('T', 'A'), ('T', 'C'), ('T', 'G'), ('T', 'T'))

The sixteen ref->alt cells of the substitution matrix, row-major in NUCLEOTIDES order. The identity pairs are cells like any other: ADR 0020 classifies A>A as a substitution.

class MatrixCell(alleles: int, percentage: str | None)

Bases: NamedTuple

One substitution-matrix cell as the page renders it.

alleles is the stored count – spelled as the column it sits under rather than count, which a tuple already means something else by. The percentage is its share of the substitution class, None when no denominator resolves, which the page renders as no second line rather than as 0.00%.

alleles: int

Alias for field number 0

percentage: str | None

Alias for field number 1

NUCLEOTIDES: tuple[str, ...] = ('A', 'C', 'G', 'T')

The nucleotides an allele may be written with, in the order the classifier’s alphabet states them. The substitution matrix’s cells are keyed and serialized in this one order, so two builds of one resource produce byte-identical JSON however the rows arrived.

class RegionAlleles(chrom: str, start: int | None, end: int | None, *, track_repeats: bool = False)

Bases: object

The allele content of one scanned region, accumulated row by row.

Counts each ROW as an allele – duplicate (chrom, pos, ref, alt) rows are legitimate per-transcript data and each is one allele. Built with track_repeats – for a resource declaring allele_multiplicity: one – it also notes the first such key among the rows it owns, which repeated_allele() reports.

A region owns the rows whose point falls inside it, which is what makes the statistic chunk-invariant: the regions of a contig partition it, so a position carries rows in exactly one of them and no merge can double-count it. A row’s optional pos_end takes no part – an allele’s value stands for its ref/alt pair, not for the bases such a column may reach over – so ownership is the shared clip_span() asked about the point, and gain#636’s edge is answered there rather than again here.

add_allele(pos: int, ref: str | None, alt: str | None) → None

Fold one row, read as the point at pos.

add_allele_batch(pos_begin: ndarray, reference: ndarray, alternative: ndarray) → None

Fold a batch of column arrays, the counting vectorized.

The same rule add_allele() applies row by row. Ownership is vectorized outright; the classification cannot be – a class is a property of one ref/alt PAIR, and an array statement of it would be a second spelling of classify_allele() – so instead each DISTINCT pair in the batch is classified once and its multiplicity added. Same function, same answer, called once per pair rather than once per row: a real allele score is overwhelmingly substitutions, so a 100,000-row batch usually holds a handful of distinct pairs, and this is ~7x the row-by-row fold over whole-genome data. (A batch of entirely distinct pairs pays a small tally overhead instead.)

The nucleotide arrays are RAW, as AlleleScore.fetch_region_allele_arrays() yields them, so this path hands the classifier exactly the strings the per-record path does.

add_record(record: tuple[Any, ...]) → None

Fold one raw record.

The nucleotides come off the RAW record: the region-values transform an allele score applies collapses a row to the point it sits at and drops them (values_from_records).

counts() → AlleleCounts

This region’s counts, class map keyed by class name.

classmethod frozen(chrom: str, allele_count: int, class_counts: dict[str, int], *, substitution_matrix: dict[tuple[str, str], int] | None = None, insertion_lengths: ExactLengths | None = None, deletion_lengths: ExactLengths | None = None, complex_grid: dict[tuple[int, int], int] | None = None) → RegionAlleles

A region restored from serialized counts, with no scan state.

substitution_matrix is None for a file written before the matrix existed – data unknown, not a matrix of zeros.

merge(other: RegionAlleles) → None

Fold the adjacent region to the right into this one.

It is the adjacency – asserted by refuse_unmergeable – that lets the counts simply add: a row belongs to exactly one of two adjacent regions, so none is counted twice. It is also why the earlier region’s repeated key, when it has one, is the one kept: all rows of a key sit at one position, so each region’s count of its own key is complete.

repeated_allele() → RepeatedAllele | None

The first allele key this region’s rows repeat, if any.

Only a region built with track_repeats looks; a region restored from a file carries none.

class RepeatedAllele(chrom: str, pos: int, ref: str | None, alt: str | None, rows: int)

Bases: NamedTuple

An allele key held by several rows, and how many rows hold it.

alt: str | None

Alias for field number 3

chrom: str

Alias for field number 0

error(resource_id: str) → AlleleMultiplicityError

The refusal of a resource declaring allele_multiplicity: one.

pos: int

Alias for field number 1

ref: str | None

Alias for field number 2

rows: int

Alias for field number 4

allele_arrays_folded_into(score: AlleleScore, chrom: str, start: int, end: int, score_ids: list[str], *, batch_size: int, alleles: RegionAlleles) → Generator[tuple[ndarray, ndarray, dict[str, ndarray]], None, None]

The bulk read, nucleotides folded off it, validated as usual.

The array twin of records_folded_into(), and the same shape: a transducer that folds each batch and yields it onward. What it yields is the batch’s [:3] slice – a plain RecordArrays – because the scan’s array door (validate_record_arrays(), ADR 0008) unpacks three names and raises on the five an AlleleRecordArrays carries.

Folding on the way IN is what lets the nucleotides reach this statistic without the door having to carry them: nothing downstream ever sees the widened batch, so nothing has to pair the two back up.

That the fold precedes the door’s verdict is unobservable: a region the door refuses raises out of the scan, and its accumulator is discarded with the failed task rather than merged.

build_allele_section_display(statistics: AlleleStatistics) → AlleleSectionDisplay

Turn the stored allele counts into the section’s payload.

The rows are built over AlleleStatistics.by_chromosome() rather than replacing it: that accessor already orders the chromosomes naturally (gain#983) and serialize reads it too, so the ordering has one owner.

One walk of the regions serves both the rows and the roll-up, as AlleleStatistics.serialize() does it: _total is order-blind, so folding the ordered entries it already has costs nothing over AlleleStatistics.global_counts(), which would rebuild every entry a second time.

merge_region_alleles(resource_id: str, regions: Iterable[RegionAlleles | None]) → AlleleStatistics | None

Fold the regions’ counts, or None for a kind that has none.

plot_complex_grid(outfile: IO, grid: dict[tuple[int, int], int]) → None

Render the complex (len_ref, len_alt) cells as a heatmap.

Drawn over the FULL clamped square rather than only the occupied cells, so the diagonal – where the MNVs sit – is visible as a diagonal and an empty region reads as empty rather than as a missing axis. The axes are exact lengths, which is the whole point of the cell scheme (gain#779): a 2->3 complex sits one cell off the diagonal from a 3bp MNV, and a binned axis would have hidden that.

Counts span orders of magnitude on real scores, so the COLOUR is log-scaled – the same choice the length histograms make on their count axis – with empty cells left as the background rather than coloured as a genuine zero.

records_folded_into(records: Iterator[tuple[Any, ...]], alleles: RegionAlleles) → Generator[tuple[Any, ...], None, None]

Yield a record stream through, folding each row into alleles.

A transducer over the very stream the per-record scan is already pulling: this statistic rides that one read rather than re-reading the region, exactly as the coverage statistic rides it.

region_alleles_for(score: GenomicScore, chrom: str, start: int | None, end: int | None) → RegionAlleles | None

A region accumulator for an allele score, None for other kinds.

Gated on the BUILT SCORE CLASS rather than the resource type string.

Until 2026.8.5 that was load-bearing: allele_score and the deprecated np_score both built an AlleleScore while equivalent_resource_types aliased neither to the other, so a gate written on type strings skipped np_score silently (gain#777). np_score is gone (gain#920) and only one spelling reaches here now, but the class gate stays: it is the property this statistic actually depends on – that the score reads alleles – and a type-string gate would have to be revisited by the next spelling that builds an AlleleScore.

save_allele_statistics(resource: GenomicResource, statistics: AlleleStatistics | None) → None

Write the statistics into the resource, with their global images.

Laid out like save_and_plot_coverage: the file first, then one image per group that has something to draw. A group the resource publishes nothing for writes no image – the info page’s section is what says whether that is “not computed” or “genuinely none”.

An EMPTY group is skipped just as an unknown one is, in all three twins (write_length_chart states the rule once). Every group here applies to every allele score, so plotting the empty ones would put an all-zero deletion histogram on each of the many scores that carry only substitutions – and a logarithmic count axis cannot draw one at all.

serves_allele_arrays(score: GenomicScore, score_ids: list[str]) → bool

Whether the bulk read can hand this score’s rows their nucleotides.

Asked BEFORE a path is chosen rather than caught after: a region whose backend will not serve the ref/alt arrays must fall back to the per-record scan, which reads the nucleotides off the record, rather than produce a statistic with no class data.

Ask it of an OPEN score. On an unopened one the answer is merely conservative – a table naming its key columns nowhere but inside its own data file cannot be known to have them until that header is read – and a spurious False costs the whole region the bulk scan for no gain in correctness.

gain.genomic_resources.statistics.base_statistic module

class MergeableRegion(*args, **kwargs)

Bases: ScannedRegion, Protocol

A scanned region that can absorb its neighbour to the right.

What RegionFoldedStatistic needs of the accumulators it holds, and no more: where the region was scanned (from ScannedRegion) and how to fold the next one in. Whether that fold adds counts, unions spans or sums histograms is the region’s own business.

merge(other: Self) → None

Fold the adjacent region to the right into this one.

NON_NUMERIC_ERRORS = (<class 'TypeError'>, <class 'ValueError'>)

What a numeric accumulator’s nan skip raises for a value it cannot read as a number – the complaints non_numeric_error() re-words, stated once so the twins that catch them cannot drift apart. TypeError is np.isnan on text or Decimal; ValueError is np.isnan on a SEQUENCE, which answers with an array whose coercion to the bool the skip needs is what raises. Catching only the first let the second escape every statistics pass’s per-score containment and abort the whole build (gain#1337). A one-element sequence is the exception: its array coerces, and as_python_number() refuses it below the skip instead.

PYTHON_NUMBER_TYPES = (<class 'float'>, <class 'int'>)

What one value folded into a numeric accumulator may be, once numpy’s own scalars have been normalized to the Python value they hold. bool rides in through int, which it subclasses – a bool score under a number histogram really does produce a 0/1 histogram.

Shared by NumberHistogram and MinMaxValue for the reason non_numeric_error() below is: the twins fold the same values, so the allow-list is stated once (gain#1338, gain#1358). It lives here rather than with the histogram because histogram.py imports min_max.py and not the reverse.

Named for Python’s types to keep it distinct from NUMBER_HISTOGRAM_VALUE_TYPES: that one is declared score value_type STRINGS, this one is the type of a single folded value.

Public rather than folded into as_python_number(), because each twin tests it INLINE before calling that: add_value runs per value of every record, a Python float is what the scan hands it, and a function call on that path costs measurably more than one isinstance against a hoisted tuple (gain#1358 measured +8% against +3%). Hoisted rather than written as a tuple literal at each check, because a tuple of names is rebuilt on every call: 1M isinstance calls, best of 5 (gain#1338), 112 ns/call for a 3-member literal against 52 ns hoisted. float first because isinstance tests a tuple in order.

class RegionFoldedStatistic(statistic_id: str, description: str)

Bases: Statistic, Generic

A statistic accumulated one region per scanned window.

The coverage, allele and fragment statistics are all built this way: the scan produces one region per window, and the statistic keeps one region per chromosome, folding each new window into the region already held for its contig. Distinct chromosomes accumulate side by side; same-chromosome windows merge, which is where the regions’ adjacency rule applies.

Subclasses supply their id and description, their readers over the folded regions, and their own serialization. Everything about the FOLD itself lives here.

fold_region(region: R) → None

Fold one region in, keyed by its chromosome.

class ScannedRegion(*args, **kwargs)

Bases: Protocol

What every per-region accumulator shows the merge step.

A statistic accumulated per scanned REGION – coverage, allele content and fragments alike – carries the extent it was scanned over, because the merge is only sound over a contig’s non-overlapping windows. What it accumulates INSIDE that extent is the statistic’s own business and is deliberately absent here.

chrom: str

end: int | None

start: int | None

class Statistic(statistic_id: str, description: str)

Bases: object

Base class genomic resource statistics.

Statistics are generated using task graphs and aggregate values from a large amount of data. Each statistic should have a clearly defined single unit of data to process (for example, a nucleotide in a reference genome).

description: str

abstractmethod static deserialize(content: str) → Statistic

Create a statistic from serialized data.

finish() → None

Perform final calculations for the statistic.

This step is optional.

This is called when resource iteration is complete.

Can also be used when creating more complex resources via deserialization.

abstractmethod serialize() → str

Return a serialized version of this statistic.

statistic_id: str

class StoredStatistic(file: str, format_version: int, writes_for: Callable[[GenomicScore], bool])

Bases: NamedTuple

One versioned statistics file, and which scores a build writes it for.

The statistic’s one declaration: the scan asks writes_for whether to accumulate it, its writer stamps format_version into the file, and the repair flow reads both back to tell a resource whose file is missing or older than the schema – and report it, never rebuild it. Not part of calc_statistics_hash.

file: str

Alias for field number 0

format_version: int

Alias for field number 1

stored_version(content: bytes) → int | None

The version a stored file carries; 0 if none, None if unreadable.

Anything but an integer format_version – absent, null, a boolean, a non-object document – is 0, older than every schema. A file that is not a JSON document at all, or not text, is None: whatever reads it reports that, not this check. Raises nothing.

writes_for: Callable[[GenomicScore], bool]

Alias for field number 2

as_python_number(value: Any, what: str) → float | int

The Python number a value that is NOT already one folds as.

The slow path of the numeric contract NumberHistogram and MinMaxValue share, reached only by a value np.isnan accepted and the caller’s inline PYTHON_NUMBER_TYPES test refused – text and Decimal never get here, and neither does a Python number.

numpy’s own scalars fold as the Python value they hold, which is what an enumerated allow-list checked BEFORE normalizing got wrong: np.float32 is not a float (only np.float64 subclasses it) and np.bool_ is not an np.integer, so all of them were refused as non-numeric even though add_batch folds them and a gene score’s column really does arrive as one (gain#1338). item() rather than a wider allow-list: widening alone would not make the histogram’s two arms agree, because numpy 2 keeps np.float32 - <python float> in float32 and that picks a different bin at the edges – the witness is in test_add_batch_matches_add_value_loop_float32_at_bin_edges.

Anything else is refused: np.isnan accepting a value does not make it one a reducer can fold. A complex or a 0-d array passes it, and min() would order either without complaint and leave the extremum holding it (gain#1358). The refusal names what the CALLER handed over, not what it was normalized to, so a nullified score’s reason does not report a np.complex128 as a plain complex.

non_numeric_error(value: Any, what: str) → TypeError

The one wording of a numeric accumulator’s refusal.

what names the accumulator in the message; the rule is one – the same shape as refuse_unmergeable(), and here for the same reason: NumberHistogram and MinMaxValue are twins that fold the same values under the same numeric contract, and a reader of a nullified score’s reason should not be able to tell which of them produced it.

Returned rather than raised, so a caller can raise it from the numpy error it is re-wording – one of NON_NUMERIC_ERRORS.

refuse_unmergeable(what: str, left: ScannedRegion, right: ScannedRegion) → None

Refuse a pair that is not adjacent-and-in-order on one contig.

Region statistics are only ever produced over a contig’s non-overlapping windows, so anything else reaching a merge is a wiring error – and it is exactly the adjacency that lets the per-region counts simply add. Refusing it loudly is the difference between a failed build and a silently wrong statistic.

what names the statistic in the message; the rule is one.

regions_in_genomic_order(regions: Iterable[R | None]) → list[R]

The scan’s regions in genomic order, kinds carrying none dropped.

Ordered here rather than trusting the order the task arguments arrived in: the merge asserts adjacency, and a correct set of regions handed over shuffled would fail that assertion.

gain.genomic_resources.statistics.chart_style module

Where the style shared by the statistics charts is set.

histogram draws the number and categorical histograms; in this package, alleles draws the complex-allele grid and length_histogram the length ladder. The charts agree on their label font and figure size so they read as one set on the resource page, and this module is where that agreement lives – a leaf that imports nothing from gain, so both tiers reach it with a module-level import. test_architecture keeps it one and says why (gain#1486).

CHART_FIGSIZE = (15, 10)

Figure size, in inches, shared by the statistics charts – all but the complex-allele grid, which is square and keeps its own.

CHART_LABEL_FONT_SIZE = 20

Label font size shared by the statistics charts.

gain.genomic_resources.statistics.coverage module

Covered-position statistics for tabular genomic scores.

Vocabulary per GLOSSARY.md and ADR 0020: a covered position is a position spanned by at least one table row — value-blind, union semantics. A segment is a maximal run of touching-or-overlapping rows carrying equal values (the whole scanned score tuple, NA equal to NA, floats exact).

The whole of this statistic lives here: the per-region accumulator and the resource-wide statistic, the fold that merges a scan’s regions into one, the write, and the render payload the info page reads. Its allele twin is laid out the same way in gain.genomic_resources.statistics.alleles; the scan wiring that feeds either is in implementations/genomic_scores_impl/scan.py.

COVERAGE_FORMAT_VERSION = 2

The format_version serialize stamps.

COVERAGE_STATISTIC = ('statistics/coverage.json', 2, <function <lambda>>)

The statistic’s one declaration – file, version, and the kind it is written for – read by the scan’s gate and by the repair check.

class CoverageDisplay(rows: list[CoverageRow], global_length: int | None, uncovered: UncoveredContigs | None, segment_lengths: ExactLengths | None)

Bases: NamedTuple

The Coverage section’s render payload, shares resolved.

Raw counts come from the stored statistic; shares are computed at render time and never stored. global_length is the denominator the section answers what part of the reference genome has values against: the WHOLE resolved reference, including contigs the score never touched (gain#1041). It is None unless every covered chromosome resolved a length – a covered contig the reference does not list is proof the reference is the wrong one, and a global percent over a partial denominator would be misleading.

property global_covered: int

property global_fraction: float | None

The whole score’s share of the reference, as a number.

global_length: int | None

Alias for field number 1

property global_percent: str | None

The whole score’s share, as the page writes it.

The same rule the rows are written through, over the same two integers: a reference all but entirely covered reads >99.99% here exactly as one of its chromosomes does above (gain#1057).

property global_segments: int | None

The segment total, or None when any row lacks segments.

All-or-nothing like the stored statistic: a global over a partial set would silently understate.

property has_fractions: bool

Whether the section renders a Covered % column at all.

The summary rows carry percentages too, so this cannot be read off rows alone: a score with no values ANYWHERE has every contig rolled up and no rows left, and a resolved global fraction of 0.0 would be computed and then dropped for want of a column to print it in.

property has_segments: bool

rows: list[CoverageRow]

Alias for field number 0

segment_lengths: ExactLengths | None

The global segment-length record, or None if unknown.

The section’s table is read off this record and its image is drawn from the ladder derived from it, so the page decides whether to show either from the same record the build wrote or declined to draw – a proxy such as the segment total could disagree with what was actually written.

property segment_row: LengthStatisticsRow | None

The Segment lengths table’s one row, None when unknown.

Global only, as the indel rows are: per chromosome this would be four more columns on the Coverage table above, to answer a question nobody asks.

uncovered: UncoveredContigs | None

The untouched part of the reference, or None when unknowable.

None – rather than a zero roll-up – whenever global_length is: “these contigs have no values” is a claim about the resolved reference being the right one, and it is not made under a denominator already known to be wrong.

class CoverageRow(chrom: str, covered: int, length: int | None, segments: int | None)

Bases: NamedTuple

One chromosome’s rendered coverage: raw counts, share derived.

length is the denominator resolved for this chromosome, or None when none was – the row then renders its raw count only. segments is None when the stored statistic carries no segment data for the resource (an old file, or a kind that publishes none).

The share is held as the two INTEGERS it is a share of rather than as either rendering of it, because the page needs both and they must not disagree: the cell sorts on fraction and displays percent, and the boundaries percent respects – covered none of it, covered all of it – are decided on the counts (gain#1057).

chrom: str

Alias for field number 0

covered: int

Alias for field number 1

property fraction: float | None

The share as a number, for the cell’s sort key.

length: int | None

Alias for field number 2

property percent: str | None

The share as the page writes it, None without a length.

segments: int | None

Alias for field number 3

class CoverageStatistics

Bases: RegionFoldedStatistic[RegionCoverage]

A resource’s covered positions, per chromosome and global.

Folds RegionCoverage the way the base class does, and serializes to the resource’s statistics/coverage.json as raw counts. Fractions are deliberately not computed here: they need chromosome lengths, which belong to a reference genome resolved at render time.

covered_by_chromosome() → dict[str, int]

covered_global() → int

static deserialize(content: str) → CoverageStatistics

Create a statistic from serialized data.

segment_lengths_by_chromosome() → dict[str, ExactLengths]

Per-chromosome exact length records – the read API for the per-chromosome data the statistics file stores (rendered consumers use the global roll-up; gain#776 reads these). Empty when any chromosome’s lengths are unknown.

segment_lengths_global() → ExactLengths | None

The fold of the per-chromosome records, unknown if any is.

segments_by_chromosome() → dict[str, int]

segments_global() → int | None

serialize() → str

Return a serialized version of this statistic.

class RegionCoverage(chrom: str, start: int | None, end: int | None, *, publishes_segments: bool = True)

Bases: object

Coverage of one scanned region, accumulated row by row.

Consumes [begin, end] spans in non-decreasing begin order – the order the scan validators guarantee – and counts each position once. The rows it is fed are pairwise disjoint: since gain#1127 the only coverage-scanned kinds are position scores, whose validators refuse a row beginning at or before its predecessor’s end (adjacent rows are legal and common, and the segment algebra depends on that). So the scan hands over each row at its FULL extent, unclipped – disjoint spans cannot double-count a position, the union stays additive across parallel regions, and the segment runs are measured at their true length rather than the region’s (gain#1175 retired the clip that overlapping rows once needed).

add_interval(begin: int, end: int, values: tuple) → None

Fold one row span into the coverage and its run bookkeeping.

The union first – a running maximum over the right edge, so a row is counted once whatever it overlaps – then the runs: the row joins the open run while it touches or overlaps it and carries equal values, and closes it otherwise. Refuses a region that publishes no segments – see _refuse_without_segments().

add_interval_batch(left: ndarray, right: ndarray, cells: list[ndarray]) → None

Fold a batch of row spans, collapsed into runs.

The vectorized statement of the rule add_interval() applies row by row — it lives HERE, beside that rule, so the equality algebra has one home: rows collapse into a run while they touch or overlap the positions covered so far and every column compares equal, nan equal to nan (ADR 0020). Only the batch’s first run then goes through add_interval(); the rest close together in _close_through().

The touching test reads the running maximum end, which is exact for a position score (whose validators refuse overlap, so the previous row IS the running maximum).

left/right are the rows’ own full extents (see the class docstring) and cells is one kept column per scanned score, all equally long.

classmethod frozen(chrom: str, covered: int, segments: SegmentSummary | None) → RegionCoverage

A region restored from serialized counts, with no scan state.

segments of None marks the segments wholly unknown – the file predates them – and is the one way a region comes to publish no segments.

merge(other: RegionCoverage) → None

Fold the adjacent region to the right into this one.

Refuses a pair that is not adjacent-and-in-order on one chromosome – see refuse_unmergeable, which states that rule for this statistic and its two twins alike.

property publishes_segments: bool

Whether this region has segment numbers to answer with.

The one predicate behind both segment_summary()’s None and the accessors’ refusal, so the two gates cannot drift apart. False only for a region frozen() from a statistics file that carried no segment data.

property segment_count: int

How many segments the region holds.

Refuses a region that publishes none – see _refuse_without_segments().

segment_lengths() → ExactLengths | None

The exact record of the region’s segment lengths.

Finalizes the still-open bookkeeping: the first and the open run are folded in on top of the interior tally, so the record’s total is exactly segment_count. None for a region read from a file that stored the count without the record. Refuses a region that publishes no segments at all – see _refuse_without_segments().

segment_summary() → SegmentSummary | None

Segment count and exact lengths, or None if unknown.

Unknown means the region was deserialized from a statistics file that predates segments altogether. This is the ASKING form of the gate the count and lengths accessors refuse through – None here, an exception there, because a caller that asks may not know and one that reaches straight for a number has asserted it does.

class SegmentSummary(total: int, lengths: ExactLengths | None)

Bases: NamedTuple

A region’s segments as its statistics file stores them.

lengths is None for a region read from a file that stored the count without the exact record (format version 1): the count is still true, the length distribution is unknown. The count is spelled total for the reason ExactLengths gives – a tuple already means something else by count.

lengths: ExactLengths | None

Alias for field number 1

total: int

Alias for field number 0

class UncoveredContigs(contigs: int, length: int)

Bases: NamedTuple

The contigs of the reference that carry no values at all.

One roll-up rather than a row each: a reference genome routinely carries hundreds of contigs a score never touches (alts, decoys, an unplaced scaffold), and per-contig zero rows would bury the contigs that do have values. The count and the base pairs are what the global fraction is measured against but has nothing to show for.

Membership is zero covered positions, not absence from the stored statistic. The two differ by backend and by nothing else – a bigWig scan visits every header contig and stores a 0 for the empty ones, a tabix scan visits only contigs the index lists – so rolling up by absence would render the same data two ways.

contigs: int

Alias for field number 0

length: int

Alias for field number 1

property percent: str

none of it.

Written through the same rule as every other cell in the column rather than as a formatted constant, so the row cannot drift from its neighbours (gain#1057). Membership is zero covered positions, so this is the rule’s exact zero by construction.

Unguarded, and returning str rather than str | None, where CoverageRow.percent and CoverageDisplay.global_percent both check their denominator first: those hold an int | None, because a denominator is what may fail to resolve, while a roll-up exists only where one DID – and _plausible_lengths() has already dropped every contig that could contribute a zero to length.

Type:

The share of these contigs that is covered

accumulate_coverage(arrays: tuple[ndarray, ndarray, dict[str, ndarray]], coverage: RegionCoverage, region: tuple[str, int | None, int | None]) → None

Fold one batch of column arrays into the region’s coverage.

Rides owned_records_mask(), the record partition every other statistic reads: a region owns the rows whose pos_begin falls inside it and measures them whole. That is exact for coverage because the rows are pairwise disjoint (see RegionCoverage): they cannot double-count a position, so the union is additive across parallel regions at full span, and the segment runs the same feed builds are measured at their true length rather than the region’s. A record beginning past the region’s end is not owned and covers nothing – the gain#636 verdict, reached here by the partition rather than by a clip.

The spans reach RegionCoverage.add_interval_batch(), which owns the run-collapse algebra; nothing here knows what “equal values” means. The batches the backends return rarely carry a row outside the queried region, so the all-kept batch skips the mask copies entirely.

build_coverage_display(resource_id: str, statistics: CoverageStatistics, lengths: dict[str, int]) → CoverageDisplay

Turn stored counts plus a resolved denominator into the payload.

What this resolves is each share’s DENOMINATOR; the shares themselves are derived on the payload, from that denominator and the count it bounds (see CoverageRow). Either way nothing is stored: the statistic stays raw counts (see CoverageStatistics). A denominator that cannot bound what it must is withheld, degrading that row to a raw count rather than rendering a zero-division or a >100%.

lengths is the whole reference the score is measured against, not merely the contigs it touched, so the global share answers what part of the reference genome has values and the untouched remainder is reported as one roll-up (gain#1041).

merge_region_coverage(resource_id: str, regions: Iterable[RegionCoverage | None]) → CoverageStatistics | None

Fold the regions’ coverage, or None for an uncovered kind.

normalize_values(values: Iterable[Any]) → tuple

A row’s score values as the tuple segment equality compares.

Every spelling of “no value” – None on the per-record path, nan in a bulk float column – becomes None, so NA equals NA whichever path produced the row (ADR 0020). Everything else is compared as stored: floats exactly, no tolerance.

region_coverage_for(score: GenomicScore, chrom: str, start: int | None, end: int | None) → RegionCoverage | None

A region accumulator for a position score, None for other kinds.

Gated on the built score’s class rather than on the resource type string, for the reason region_alleles_for() gives.

A position score is the one kind whose rows are pairwise disjoint – the rule it is registered under in record_validation refuses a row beginning at or before its predecessor’s end – and that is the property this statistic depends on: only then does the union of the spans answer “is there data at this position at all?” exactly, so that the count is a genuine measure of what the resource covers and the fraction a genuine completeness figure. What disjointness buys the accumulator – full-span rows, an additive union, published segments – is RegionCoverage’s own docstring.

An allele score’s rows are points, so there is no span to union. A fragment score’s rows overlap by design; fragments and ADR 0020 say why that keeps the kind out of both coverage and segments.

resolve_chrom_lengths(resource: GenomicResource, ref_genome: ReferenceGenome | None, score_lengths: Callable[[], Mapping[str, ChromLength]], chroms: Iterable[str]) → dict[str, int]

Resolve chromosome lengths for the render-time denominator.

The ladder: the reference_genome the caller resolved from the resource’s label; else the score’s own lengths, kept where their source is exact (a bigWig header’s contig sizes); else raw counts (an empty mapping).

What comes back is the whole universe the fraction is measured against, not only the contigs the score touched (gain#1041): every contig of the resolved genome, or every contig the score rung has an exact length for. chroms – the covered contigs – is still passed in so that a covered contig the resolved source does NOT list is visible to the caller by its absence, which is what degrades the fraction.

The two rungs are not interchangeable: only the genome rung answers “what part of the reference genome has values”, while the score rung answers “what part of what this file declares”. docs/adr/0020 carries the worked example and is the record to amend if this changes; the user-facing half is in docs/source/grr.rst.

Both rungs are the CALLER’s to supply: the genome needs a repository, which only exists during a page build, and the score’s records come from the implementation’s ladder (GenomicScoreImplementation.get_chrom_lengths()), which opens the score to ask – so they are asked for, through score_lengths, only once the genome rung has nothing. Each record’s best answer is the candidate, and whether its source may serve as a denominator is ChromLengthSource.is_exact’s call, not this function’s (gain#1414): a tabix probe’s upper bound and an in-memory table’s extent are dropped here, and a contig with no length at all (proven empty, or one the probe could not bracket) with them.

save_and_plot_coverage(resource: GenomicResource, statistics: CoverageStatistics | None) → None

Write the coverage statistics and their histogram images.

Does nothing for a kind that has no coverage, and skips a group’s image when there is nothing to draw – whether the group is unknown or known and empty. The same rule as save_and_plot_alleles.

gain.genomic_resources.statistics.exact_lengths module

A length statistic’s exact record: stored, folded, read and charted.

ADR 0020 gave segments, fragments and indels one log2 binning for their length histograms. The ladder lumps {2, 3} into one bin and {4, 5, 6, 7} into the next, so no exact minimum, maximum, mean or median survives it – and those four are what the statistics table on an info page exists to show. So each statistic stores an ExactLengths record instead: an exact {length: count} map clamped at LENGTH_MAP_CLAMP, beside four scalars accumulated on the unclamped length.

The ladder is still what the CHART is drawn on; it is derived from the map at render time by length_ladder() rather than stored beside it, so the picture and the numbers beneath it cannot drift.

Nothing here knows what KIND of thing has a length. The indel groups of an allele score were the first users (gain#1118), a position score’s segments the second (gain#1543) and a fragment score’s fragments the third (gain#1544). The record’s stored key for the count is count whatever the kind, and the row formatter takes the group label from its caller, so “alleles”, “segments” and “fragments” all fit.

class ExactLengths(lengths: dict[int, int], total: int, sum: int, min: int | None, max: int | None)

Bases: NamedTuple

One group’s lengths: an exact map and four scalars.

lengths maps a length in base pairs to how many items have it, clamped at LENGTH_MAP_CLAMP – the key AT the clamp means “that long or longer”. The map is what the chart’s bins and the median are derived from.

The scalars are what keep the clamp from becoming a lie. total, sum, min and max are all accumulated on the UNCLAMPED length, so min, max and the mean stay exact however far the tail runs; only a median landing in the overflow bucket degrades, and it says so. A clamped map alone would understate the mean and cap the max – the one statistic that exists to describe the tail.

They are stored rather than derived even where the map could give them, because that is the whole point: sum cannot be recovered from a clamped map at all, and a max recovered from one would read 8192 for a 40,000 bp deletion. total does equal the map’s total and is kept beside the other three so the four merge as one thing and the file says outright what the mean is over.

min and max are None exactly when total is zero: a group with no items has no shortest and no longest, and 0 is not a length anything here can have.

classmethod from_stored(stored: dict[str, Any]) → ExactLengths

The record a statistics file holds under one group’s key.

min and max are read tolerantly and the other three are not, which is deliberate rather than sloppy: those two are legitimately null for a group that was scanned and found nothing, so absent and null must read alike. A group missing its map, count or sum is a malformed file, and raising names it.

property has_counts_to_plot: bool

Whether this group has anything to draw a chart of.

Asked of a KNOWN record; an unknown group is the callers’ None and gets the same answer, because the counts axis is logarithmic and can render neither, and a chart of nothing under a lengths heading states nothing either.

One spelling, because the statistics build and the page must agree exactly – a build that skips the image while the page links it leaves a dangling thumbnail, and the reverse leaves a file nothing references.

lengths: dict[int, int]

Alias for field number 0

max: int | None

Alias for field number 4

property mean: float | None

The mean length, exact past the clamp. None if empty.

property median: float | None

The middle length, or the mean of the middle two if even.

The standard convention, stated on the ITEMS rather than on the distinct lengths: {2, 3} is 2.5, and a group of one 2 and nine 3s has a median of 3, not 2.5.

Read off the map, so it degrades where the map does. When median_is_clamped this is a LOWER BOUND rather than the median, and the page renders it as one.

property median_is_clamped: bool

Whether either middle item fell in the overflow bucket.

The one statistic here the clamp can blunt, so it is asked outright rather than left to a reader to notice that a suspicious number is suspicious.

Asked of the two middle LENGTHS, never of their average. Every key in the map is at most the clamp, so an average reaching it means BOTH middles were in the overflow bucket; a group with one middle below the clamp and one above averages to something under it and would otherwise be published as though exact. On a group of one 1 bp and one 40,000 bp deletion that reads 4096.5, where the truth is 20,000.5 – a fabricated number beside three exact ones.

Since the clamped side is a floor, the average is a floor too: the true median is at least what median computes, whichever of the two middles was clamped.

middle_lengths() → tuple[int, int] | None

The one or two lengths the median is taken over.

Both indices land on the same item when the count is odd, so the pair collapses and the median is that length. Kept as the PAIR rather than folded straight into an average, because whether the clamp blunted the answer is a property of these two lengths and cannot be read back off their mean – see median_is_clamped.

The rank arithmetic is the same convention MedianAggregator applies to a weighted value list, restated here rather than shared: that one is a private method on a stateful annotation-layer aggregator keyed on its own accumulated values, and reaching it would cost more than the six lines it saves. If a third caller ever wants it, the extraction is a free weighted_median over (value, count) pairs in gain.utils.

min: int | None

Alias for field number 3

stored() → dict[str, Any]

This record as the statistics file writes it.

Keys SORTED and written as strings: the map is sparse and two chunkings of one resource meet its lengths in different orders, so sorting is what makes the file byte-identical however the rows arrived. The count is written under count.

sum: int

Alias for field number 2

total: int

How many items the map counts. Spelled total rather than count, which a tuple already means something else by – the reason MatrixCell spells its count alleles; a record that does not know what it counts needs a neutral word. The STORED key is count, which is what it is in the file.

LENGTH_MAP_CLAMP = 8192

The longest length the stored map resolves exactly. A length at or above it folds into one overflow bucket keyed by the clamp, which therefore reads “this many bases or more”. The clamp is TOTAL in the sense COMPLEX_LENGTH_CLAMP is: every length lands in exactly one bucket, so the map’s values sum to the record’s count.

Part of the stored format: it must not change once resources carry maps built from it.

It must never fall BELOW LENGTH_HISTOGRAM_DISPLAY_CAP, because the chart’s bins are derived from this map and a bin between the two would be drawn from lengths the map had already folded away. They are EQUAL, which is the tightest the rule allows and is what makes the derived chart identical to one drawn from a stored ladder: the plot sums every bin at or above the cap into one overflow bar anyway, so folding at the same length loses nothing it would have drawn separately. The consequence to know is the other direction – the display cap is documented as free to change, and it is not free to be RAISED. Raising it means raising this first, which is a stored-format change and needs every resource rebuilt.

class LengthArrayTally

Bases: object

A mutable group of lengths, accumulated a whole batch at a time.

The other scan-side counterpart to ExactLengths, for a statistic that meets its lengths as column arrays and in numbers the dict tally cannot afford – a position score’s segments run to billions. One int64 array of LENGTH_MAP_CLAMP + 1 counters stands in for the map, so a batch folds in with one bincount and the cost is bounded by the clamp, not by the run count. The map appears only when the group is frozen, sparse and with plain Python ints, so a record from here is indistinguishable from one the dict tally froze.

total, sum, min and max are the same four scalars the dict tally keeps, accumulated on the UNCLAMPED batch.

add(length: int) → None

Fold one length in, with no numpy vector call.

The scalar twin of add_batch(), for the lengths a caller meets one at a time – the per-record scan’s runs, and the runs a region can only measure once its neighbours are known. It runs once per segment, and a position score’s segments run to billions: a vector call here, even a one-element bincount, would cost microseconds where a counter increment costs a few hundred nanoseconds. The extremes are updated longhand for the reason LengthTally.add() gives.

No check that length is positive, again like the dict tally: a run’s length is end - begin + 1 over a row span, at least 1 by construction, and the check belongs where lengths arrive from outside – add_batch().

add_batch(lengths: ndarray) → None

Fold a whole batch of lengths in.

lengths is an integer array of any width; the fold widens to int64 itself. Every length must be at least 1: bincount would accept a 0 silently and put it in counter 0, which is not a length anything can have, so the batch is checked before it is folded – and a refused batch leaves the tally as it was.

frozen() → ExactLengths

This group as the inert record a region hands out.

The map is SPARSE – a counter that stayed at zero is no key – because that is the map the dict tally builds, and the file must not say which tally wrote it.

max: int | None

merge(other: LengthArrayTally) → None

Fold another group of the same kind into this one.

LengthTally.merge() over arrays: the counters add elementwise, total and sum add, and the extremes take the extreme through the same rule.

min: int | None

class LengthStatisticsRow(group: str, total: str, min: str, max: str, mean: str, median: str)

Bases: NamedTuple

One group as the statistics table renders it.

Every field is already TEXT. The table shows five different kinds of number – a count, two exact lengths and two averages, one of which can be a floor – and formatting each where it is built keeps the template free of a second rule about how a length is written.

group is whatever the caller labels the row: “insertions” on an allele page, “segments” or “fragments” elsewhere. The count column is total for the same reason the record’s is.

group: str

Alias for field number 0

max: str

Alias for field number 3

mean: str

Alias for field number 4

median: str

Alias for field number 5

min: str

Alias for field number 2

classmethod of(group: str, lengths: ExactLengths) → LengthStatisticsRow

One group’s row; the LENGTH cells empty when it holds none.

A group that was scanned and found nothing has a genuine count of zero, which renders as 0 – that is the answer, not a missing one. What it has no answer for is a shortest, longest or average length, so those four render empty rather than as zeros, which would read as items of length nothing.

total: str

Alias for field number 1

class LengthTally

Bases: object

A mutable group of lengths, accumulated one distinct length at a time.

The scan’s counterpart to ExactLengths, which is what a region hands out. Separate because the map is updated in place: a tuple rebuilt per row would copy the whole map each time, which on a resource whose lengths run to thousands of distinct values is quadratic in the lengths seen.

Dict-backed, so its cost is in the DISTINCT lengths seen, not in the clamp: right for a statistic that meets few distinct lengths per batch and folds each with its multiplicity. A statistic that meets its lengths as whole arrays wants LengthArrayTally.

add(length: int, multiplicity: int) → None

Fold multiplicity items of one exact length in.

Written out longhand rather than with min()/max(), which is not style: this runs once per DISTINCT length per scanned batch, and the two builtin calls measured 158ns of the 258ns this method had added over a fixed-bin increment. Comparing in place brings it to ~20ns over.

frozen() → ExactLengths

This group as the inert record a region hands out.

lengths: dict[int, int]

max: int | None

merge(other: LengthTally | ExactLengths) → None

Fold another group of the same kind into this one.

The ONE statement of how two groups come together, so the region merge and the global roll-up cannot drift: maps add per length, total and sum add, and the extremes take the extreme. No re-clamping – both maps are already keyed on clamped lengths, while min/max are exact on both sides and stay exact here. The other side may be a stored record as well as a tally: it is only read.

min: int | None

classmethod restored(lengths: ExactLengths) → LengthTally

A tally holding what a stored group already counted.

NO_LENGTHS = ({}, 0, 0, None, None)

A group that was scanned and holds nothing – distinct from a group that was never scanned, which is None.

Its map is shared, so nothing may mutate it. Nothing does: a roll-up accumulates in a fresh LengthTally and only ever READS a record, this one included. Add a path that mutates a group in place and it must start from LengthTally(), not from this.

folded_lengths(records: Iterable[ExactLengths | None]) → ExactLengths | None

The fold of every chromosome’s record, None if any is unknown.

What a statistic’s global entry IS for a length group: the merged_lengths() rule over all the chromosomes, so a partial roll-up can never silently understate. A file that stored the counts without the records (format version 1 of coverage.json and fragments.json) passes its count gate and fails this one.

Accumulated in one tally, so the fold costs each record’s own keys and nothing per record beyond them.

length_ladder(lengths: ExactLengths) → list[int]

The exact map binned onto the shared log2 ladder, for the chart.

Derived at render time rather than stored, so the chart and the statistics beside it cannot drift: there is one source of truth and the picture is a view of it.

The result is identical to a histogram stored on the ladder, not merely close. plot_length_histogram sums every bin at or above its display cap into a single overflow bar, and LENGTH_MAP_CLAMP equals that cap – so the lengths the map folded together are exactly the ones the chart was going to add up anyway. Below the cap the binning is the same function on the same lengths.

merged_lengths(left: ExactLengths | None, right: ExactLengths | None) → ExactLengths | None

The sum of two groups, unknown if either is.

The _merged_matrix() rule over the exact length maps: an unknown side makes the whole merge unknown rather than a smaller number. The arithmetic itself is LengthTally.merge(), so the roll-up here and the region merge in the scan are one rule.

merged_tallies(left: LengthTally | None, right: LengthTally | None) → LengthTally | None

The same rule between two SCANNED groups, folded left in place.

For a statistic whose regions hold a dict tally that may be unknown – the allele score’s indel groups.

stored_lengths(entry: Mapping[str, Any], key: str) → ExactLengths | None

A stored record under key, None when the entry carries none.

The record is the ONLY thing read. A file predating it – one written with the log2 ladder histogram the record replaced (insertion_length_histogram, segment_length_histogram, fragment_length_histogram) – carries no record, so the group reads as unknown and the page says “not computed” until the resource is rebuilt.

That is one reader rather than a compatibility branch, deliberately (ADR 0020 as amended by gain#1118). A branch that read the old histograms would have to publish them as an ExactLengths whose exact map, sum, min and max are all unrecoverable, so every statistic in the table would be a guess at bin resolution presented as a number.

write_length_chart(resource: GenomicResource, filename: str, lengths: ExactLengths | None, kind: str) → None

Draw the group’s chart into the resource, exactly when it has one.

A group that is unknown, or known and empty, writes no image – the info page’s section is what says which, and it asks ExactLengths.has_counts_to_plot too, so the build and the page cannot disagree about whether an image exists. What skipping costs is a previous build’s image left behind when a group empties out, and nothing links the leftover: the page reads the stored record, not the directory.

gain.genomic_resources.statistics.fragments module

Fragment count and exact fragment-length statistics for fragment scores.

Vocabulary per GLOSSARY.md: a fragment is a table row AS STORED. Overlapping, nested and duplicate rows each count once, at their own unclipped span — there is no run algebra and no stitching here, because a row is owned whole by exactly one scanned region.

This statistic used to ride inside gain.genomic_resources.statistics.coverage, as one more optional group in the coverage file. It does not any more (gain#1127): a fragment score’s rows deliberately overlap, so the union of their spans measures nothing a reader wants, and the kind is no longer coverage-scanned at all. While the two shared a carrier, dropping the union would have dropped the tally with it — which is why the tally moved out first.

Laid out like its two twins, coverage and alleles: the per-region accumulator, the resource-wide statistic, the fold that merges a scan’s regions into one, the write, and the render payload the info page reads. The scan wiring that feeds all three is in implementations/genomic_scores_impl/scan.py.

Lengths are kept EXACTLY, as the kind-neutral record in gain.genomic_resources.statistics.exact_lengths (ADR 0020 as amended by gain#1544, the way segments left the stored ladder in gain#1543 and the indel groups in gain#1118): the file stores the record, the table on the page is read off it, and the chart is drawn on the ladder derived from it at render time.

FRAGMENT_FORMAT_VERSION = 2

The format_version serialize stamps.

FRAGMENT_STATISTIC = ('statistics/fragments.json', 2, <function <lambda>>)

The statistic’s one declaration – file, version, and the kind it is written for – read by the scan’s gate and by the repair check.

class FragmentDisplay(rows: list[FragmentRow], fragment_lengths: ExactLengths | None)

Bases: NamedTuple

The Fragments section’s render payload.

Counts and the global length record – nothing to resolve: a fragment is a table row, and rows have no natural total to be a fraction of. The global count is the sum of the rows, exactly as the stored statistic’s global entry is the merge of its per-chromosome ones.

fragment_lengths: ExactLengths | None

The global fragment-length record, or None if unknown.

Unknown is a THIRD answer, distinct from a record that is known and empty: the counts and the lengths are read independently, so a file that stored the counts alone leaves the lengths unknown while the counts stay exact. The section renders all three apart – “not computed”, the table and chart, “no fragments” – because collapsing the first into the last would deny fragments the table beside it is counting.

property fragment_row: LengthStatisticsRow | None

The one row of the “Fragment lengths” table, if known.

property global_fragments: int

rows: list[FragmentRow]

Alias for field number 0

class FragmentRow(chrom: str, fragments: int)

Bases: NamedTuple

One chromosome’s fragment count, as the info page renders it.

chrom: str

Alias for field number 0

fragments: int

Alias for field number 1

class FragmentStatistics

Bases: RegionFoldedStatistic[RegionFragments]

A resource’s fragment counts and lengths, per chromosome and global.

Folds RegionFragments the way the base class does, and serializes to FRAGMENT_STATISTICS_FILE as raw counts beside each chromosome’s exact length record.

static deserialize(content: str) → FragmentStatistics

Create a statistic from serialized data.

fragment_lengths_by_chromosome() → dict[str, ExactLengths]

Per-chromosome fragment-length records, as stored.

A chromosome whose lengths are unknown is left OUT rather than given an empty record, which would read as “measured, and empty”.

fragment_lengths_global() → ExactLengths | None

The fold of every chromosome’s lengths – folded_lengths(), so the all-or-nothing rule is the coverage twin’s.

fragments_by_chromosome() → dict[str, int]

fragments_global() → int

serialize() → str

Return a serialized version of this statistic.

class RegionFragments(chrom: str, start: int | None, end: int | None, *, accumulates: bool = True)

Bases: object

The fragments of one scanned region, counted row by row.

Consumes row spans and counts each row once, at its own length. Unlike RegionCoverage this carries no opt-out flag: a region is built only for a kind whose rows ARE fragments, so every instance publishes a count.

add_fragment(length: int) → None

Count one fragment of that many base pairs.

The row’s OWN span, never clipped to the region: a region owns the rows beginning inside it and measures them whole, so a fragment is counted once at its true length however the contig was split.

add_fragment_batch(lengths: ndarray) → None

Count a whole batch of fragment lengths at once.

The vectorized statement of add_fragment(): one fold of the whole array into the tally, because a genome-scale fragment score has hundreds of thousands of rows, and this is the path ADR 0001 deleted the per-row object churn from. A length below 1 is refused by the tally itself.

fragment_lengths() → ExactLengths | None

The region’s fragment lengths as a record, None if unknown.

property fragments: int

How many rows this region counted.

classmethod frozen(chrom: str, fragments: int, lengths: ExactLengths | None) → RegionFragments

A region restored from serialized counts, with no scan state.

lengths of None marks the length record unknown – the file stored the count alone. The COUNT is unaffected and still reads, so a file like that renders its table and no image.

merge(other: RegionFragments) → None

Fold the adjacent region to the right into this one.

Refuses a pair that is not adjacent-and-in-order on one chromosome – see refuse_unmergeable, which states that rule for this statistic and its two twins alike.

No stitch is needed: a row is owned whole by exactly one region, so the merged count is a plain sum and the merged lengths are the tallies’ own merge. Only scanned regions get this far: a region restored from a file has no extents, and the adjacency rule refuses it.

accumulate_fragments(arrays: tuple[ndarray, ndarray, dict[str, ndarray]], fragments: RegionFragments, region: tuple[str, int | None, int | None]) → None

Fold one batch of column arrays into the region’s fragment tally.

Fragments partition RECORDS: the rows this region OWNS, each measured at its own unclipped span. So this rides owned_records_mask(), the record partition every statistic reads – a row is owned whole by exactly one region however the contig was split.

build_fragment_display(statistics: FragmentStatistics) → FragmentDisplay

The Fragments payload for a resource that has the statistic.

Always a payload: the file existing IS the statistic, where the group riding inside the coverage file used to have to answer “present but carrying no fragments” as well (gain#1127). Whether the section renders at all is decided by the file’s presence, in the implementation’s get_fragment_display.

merge_region_fragments(resource_id: str, regions: Iterable[RegionFragments | None]) → FragmentStatistics | None

Fold the regions’ counts, or None for a kind with no fragments.

region_fragments_for(score: GenomicScore, chrom: str, start: int | None, end: int | None) → RegionFragments | None

A region accumulator for a fragment score, None for other kinds.

Gated on the built score’s class rather than on the resource type string, for the reason region_alleles_for() gives: the builder has already resolved both spellings of the kind to FragmentScore, and asking the string again would restate its dispatch. The class is the property this statistic depends on – that the rows ARE fragments, each counted once at its own span.

save_and_plot_fragments(resource: GenomicResource, statistics: FragmentStatistics | None) → None

Write the fragment statistics and their histogram image.

Does nothing for a kind that has no fragments, and skips the image when there is nothing to draw. The same rule as save_and_plot_coverage and save_allele_statistics.

gain.genomic_resources.statistics.length_histogram module

The one fixed length ladder, and the chart drawn on it.

ADR 0020 gave every length histogram this binning as its STORED form, so that per-chromosome results merged into exact global ones and chunked scans merged exactly for the same reason. No kind stores it any more: indels left the stored ladder in gain#1118, segments in gain#1543 and fragments in gain#1544 (ADR 0020 as amended). All three keep an exact {length: count} map in gain.genomic_resources.statistics.exact_lengths and merge on that; what they use from here is the RENDERING – the map is projected onto these bins at draw time, so the chart keeps the shape the stored histograms drew. The ladder lives here rather than beside any one statistic because it is the shared contract every length chart is drawn on, not a detail of one kind.

What the ladder does NOT bin at all is the complex allele grid: its cells are exact lengths (ADR 0020 as amended by gain#779), for reasons that belong with that grid rather than here.

length_histogram_bin_index(length: int) → int

The fixed log2 bin a length of that many base pairs falls in.

plot_length_histogram(outfile: IO, histogram: list[int], item: str, display_cap: int = 8192) → None

Render a length histogram on the fixed log2 bins as PNG.

Styled to sit beside the per-score value histograms on the resource info page: the figure size and label font are the ones the statistics charts share, read from gain.genomic_resources.statistics.chart_style. item names what was measured – segments, fragments, insertions, deletions – and appears in both axis labels; the bins are the same ladder every time, which is what lets one renderer serve them all. Required, with no default: a fragment histogram silently labelled “segment” is the one mistake this parameter exists to prevent.

display_cap is the length the drawn axis stops at: every bin at or above it becomes one overflow bar. A parameter rather than a constant so a resource kind whose lengths genuinely run longer can raise its own axis, up to LENGTH_MAP_CLAMP. It is snapped down to its own bin on the ladder, so a cap between two edges caps at the lower one – pass a power of two to get the axis the number reads as.

gain.genomic_resources.statistics.min_max module

class MinMaxValue(score_id: str, min_value: float = nan, max_value: float = nan)

Bases: Statistic

Statistic that calculates Min and Max values in a genomic score.

add_value(value: float | generic | None) → None

Fold one value into the running extremum.

The same numeric contract as NumberHistogram.add_value, stated the same way and in the same order, so that a nullified score’s reason reads alike whichever twin refused (gain#1313, gain#1358). The agreement is pinned value by value in test_numeric_reducer_twins.

static deserialize(content: str) → MinMaxValue

Create a statistic from serialized data.

merge(other: MinMaxValue) → None

Widen this score’s range in place to cover other’s.

serialize() → str

Return a serialized version of this statistic.

class MinMaxValueStatisticMixin

Bases: object

static get_min_max_file(score_id: str) → str

class NullMinMaxValue(score_id: str, reason: str)

Bases: MinMaxValue

A score’s min/max, refused rather than measured.

The min/max twin of NullHistogram, and it exists for the same reason: a reducer that refuses a value has to leave something behind that the rest of the pass can keep feeding harmlessly, so one score costs one score rather than the resource’s whole statistics build (gain#1285, gain#1313).

add_value is a no-op, which is what makes the refusal cost ONE exception rather than one per record – the score is swapped out for this and every later record folds into nothing. Without that latch a caller catching per value would pay a raise per record, which on a genome-scale resource is multiples of the pass it is protecting.

It keeps the inherited nan min and max, so a caller that never learned about this type still reads it as a score with no values to bin. update_hist_confs does know about it, and nullifies it in a branch of its own – the same outcome as the nan branch beside it, reached without depending on the seed, and worded from reason: “the values could not be folded” and “there were no values” are different facts about a resource, and a curator reading a nullified score needs to know which.

It deliberately does NOT override the inherited MinMaxValue.serialize() / MinMaxValue.deserialize(), which would write the nan seed and read back a plain MinMaxValue – losing the refusal. Those two have no callers on any live path (ADR 0001), and a refused min/max never reaches a file; a refusal travels between the region tasks and the merge as an OBJECT. This is the one asymmetry with NullHistogram, which avoids the trap by being a sibling of its histogram rather than a subclass.

add_value(value: float | generic | None) → None

Fold one value into the running extremum.

The same numeric contract as NumberHistogram.add_value, stated the same way and in the same order, so that a nullified score’s reason reads alike whichever twin refused (gain#1313, gain#1358). The agreement is pinned value by value in test_numeric_reducer_twins.

merge(other: MinMaxValue) → None

Widen this score’s range in place to cover other’s.

gain.genomic_resources.statistics.moments module

The count, sum and sum of squares a number histogram folds beside its bars, and the mean and standard deviation read off them (gain#1589).

MOMENT_KEYS: Final = ('count', 'sum', 'sum_of_squares')

The keys the accumulators are stored under, beside a histogram’s min_value / max_value. Owned here with Moments.stored() and Moments.from_stored(), so the file format of the moments has one home.

class Moments(count: int = 0, sum: float = 0.0, sum_of_squares: float = 0.0)

Bases: object

Running count, sum and sum of squares of a weighted value stream.

Weighted exactly as a histogram’s bars are – by the count handed to add_value – so count is the bars’ total plus the out-of-range counts, and an out-of-range value still enters sum and sum_of_squares at its real value. Only the three accumulators are stored; mean and std are derived on read.

The two sums are exact: each per-term product rounds once, as a double, and the additions of those products do not round at all. sum and sum_of_squares are the exact totals rounded once, so they depend only on the multiset of folded terms, never on the order of add(), add_batch() and merge(), the batch boundaries or the numpy build. Two hosts that fold the same values store the same sums, to the bit.

A stream nobody has folded into is a Moments at zero; a histogram whose stored file predates these keys has NO moments, which its owner spells None rather than a zero here.

add(value: float, count: int) → None

Fold one value, count times.

The per-value twin of add_batch(): the products are value * count and value * value * count, in that order, so each term rounds exactly as the batch arm rounds it, and the exact sums make the two arms agree to the bit.

add_batch(values: ndarray, weights: ndarray) → None

Fold (value, weight) pairs, vectorized.

values are float64 and already free of nan. The per-term products are the ones add() folds – value * weight, then value * value * weight – and they are added exactly, so the result is the same, to the bit, as an add() per pair in any order.

One temporary, reused in place: a fresh 800 KB array per product of a 100k batch costs more than the reductions themselves.

count

static from_stored(data: dict[str, Any]) → Moments | None

Read the accumulators out of a histogram file’s mapping.

None when the file has none: one written before the accumulators existed has no count key, and nothing rebuilds it on its own (the statistics hash is unchanged), so its moments stay unknown until the resource’s statistics are next rebuilt.

property mean: float | None

The weighted mean; None when nothing was folded.

merge(other: Moments) → None

Add other’s accumulators to this one’s, exactly.

property std: float | None

The POPULATION standard deviation; None when nothing folded.

Divides by count, not count - 1: the histogram describes the whole resource, not a sample drawn from it. A variance that rounds below zero – possible only when the values barely vary around a large mean – is clamped to zero, with a warning.

static stored(moments: Moments | None) → dict[str, Any]

The accumulators under MOMENT_KEYS, for a histogram file.

As Python numbers whatever was folded (a numpy weight would otherwise reach json.dumps as an int64), and every key None for unknown moments – which from_stored() reads back as exactly that.

property sum: float

The weighted sum of the values, correctly rounded.

property sum_of_squares: float

The weighted sum of the squared values, correctly rounded.

summary() → tuple[tuple[str, str], ...] | None

n, mean and sd for a summary page, as (label, value) pairs; None when nothing folded.

n with thousands separators, the other two at the three significant digits a histogram’s values_domain uses. The page lays the pairs out one per line, so every info page formats them here and nowhere else.

MomentsSummary

The rendered moments, one (label, value) pair per line of a summary page’s cell: n, mean, sd in that order.

alias of tuple[tuple[str, str], …]

gain.genomic_resources.statistics.percentages module

How a share of a whole is written, for every table that writes one.

One rule in one place, so every info-page table that renders a count as a percentage of a total says the same thing the same way: the genomic-score page’s Coverage table (denominator: a chromosome’s length) and Alleles tables (denominator: an allele total), and the reference-genome page’s nucleotide and di-nucleotide tables (denominator: the bases, or the pairs, counted).

percentage_of() is the rule itself, and it is deliberately SCALAR – one (count, total) pair to one string. What varies between the callers is the contract around a MISSING denominator, and that is theirs to keep. Coverage resolves a denominator per row and degrades only that row, as do the Alleles table’s per-class share columns, whose denominator is each chromosome’s own allele count. percentages_over() is the whole-map answer – no total, no percentages at all – which the substitution matrix, the complex table and the reference-genome tables ask.

Every caller hands over the INTEGERS. The boundaries are decided on them, so a share stored already multiplied out cannot be written by this rule; the reference-genome page derives its shares from the stored counts for that reason.

percentage_of(count: int, total: int) → str

count as a percentage of total, floored and capped.

The result is text for an HTML page and can carry markup-significant characters – the floor and the ceiling begin with < and > – so a template rendering it must escape, which the .jinja HTML templates do and the Markdown ones deliberately do not.

total must be positive; a share of nothing is not a percentage, and what to show instead is the caller’s question rather than this one’s (see the module docstring).

Two answers a bare "%.2f%%" gets wrong:

• A nonzero count too small to survive two decimals renders <0.01%, never 0.00%. On a real score complex is 881 alleles out of 727,413,443 while other is genuinely empty – and telling those two apart is the whole reason a percentage is shown at all.

• A count that falls SHORT of the total but rounds up to it renders >99.99%, never 100.00% – the floor reflected (gain#990). On that same score the substitutions are all but 881 of the alleles, and a column reading substitution 100.00% beside complex <0.01% says the resource is entirely one class in the act of showing that it is not.

Both boundaries are decided on the INTEGERS, never on the rendered float: a count that IS the total renders 100.00% and a count of zero renders 0.00%, because only a share that is not the whole is written as short of it, and only a share that exists at all is written as too small to see.

percentages_over(counts: Mapping[K, int], total: int) → dict[K, str] | None

Each count as a percentage of total, None without one.

Every cell is written by percentage_of(). What this adds is the MAP contract around a missing denominator: a zero total has no percentage, and the answer is None for the WHOLE map rather than per cell, because the denominator is a property of the table. The caller decides what a table without one shows – the Alleles section drops the column, the reference-genome page leaves the cells empty.

gain.genomic_resources.statistics.png_bytes module

Make the bytes of a statistics image depend only on its pixels.

Matplotlib writes a tEXt chunk Software with its own version into every PNG it saves. Two hosts with different matplotlib versions then write different bytes from the same pixels, and the md5 of the image changes in .MANIFEST for no real change (gain#2048).

The chunk is removed after the image is drawn, not kept out with savefig(metadata=...): a resource plot_function calls savefig itself, and gain does not control that code.

A leaf that imports nothing from gain, as chart_style is: both histogram and this package write images, and the statistics package does not import histogram.

rendered_png(draw: Callable[[BinaryIO], object]) → bytes

Let draw write an image into memory, and return its bytes.

The bytes come back without the Software chunk, ready for the caller to write into the resource.

without_software_chunk(data: bytes) → bytes

Return data without its tEXt chunks keyed Software.

Every other chunk is copied byte for byte, so the image data and each CRC stay as they were. Bytes that are not a PNG come back unchanged: a plot_function can write anything.

gain.genomic_resources.statistics.record_validation module

What the statistics scan refuses, one rule per score kind (ADR 0027).

ADR 0008 gave the statistics scan sole ownership of validation – reads never validate – and gave each score kind its own rule rather than a shared one it did not choose. Both halves stand; this module is where the second one is written down.

The rules used to be @abstractmethod on GenomicScore with one body per kind beside the read code, while their only callers were the scan’s two doors (genomic_scores_impl/scan.py) and statistics/alleles.py. A read class carrying a rule only its consumer applies is what gain#1269 moved here.

Two functions, both functools.singledispatch(), both dispatching on the score’s class:

• validate_records() – the per-record door’s rule;

• validate_record_arrays() – the same rule over a batch’s columns.

Both are transducers: they hand back exactly what they were given, in order, and raise at the first record their kind cannot mean. A record the kind’s ORDERING rule refuses raises MalformedResourceError; a record whose own end precedes its own begin raises the plain OSError that inverted_span_error() builds, which is a claim about one record rather than about the resource’s order. Neither re-reads and neither materialises the region – the scan pays for one read, and these ride it.

Both read RAW records rather than the spans a kind yields, because a kind’s normalization destroys the evidence: an allele score collapses a record to the point it sits at, discarding its end entirely. Raw is also the only layer at which the per-record and the vectorized rule can say the same thing – clipping a record to the scanned region would tie the verdict to how the contig happened to be partitioned.

Three kinds, two rules. A position score promises one value per position, so it refuses records that overlap or merely touch. An allele score and a fragment score both allow several records at one position and refuse only a record that moves backwards; they differ solely in the noun their message names. So the backwards rule is written once per shape and registered twice, with the noun supplied at the registration. ADR 0008 accepted six copies to avoid a hidden shared statement – one class attribute two validators read and interpreted differently. A body called explicitly from two registrations hides nothing, so that reason does not reach this arrangement (ADR 0027).

A kind nobody wrote a rule for cannot be scanned. The undecorated body of each function – singledispatch’s default, which it keys on object rather than on GenomicScore – raises NotImplementedError naming the class it was handed. The three kinds are flat siblings, so a kind added later reaches that default rather than inheriting a rule chosen for something else.

That refusal is at run time – a missing registration has no static analogue the way a missing @abstractmethod override had – so test_every_buildable_kind_is_registered stands in for the check the type checker used to make. ADR 0001’s gain#1261 Amendment depends on this; see ADR 0027.

validate_record_arrays(score: GenomicScore, batches: Iterator[tuple[ndarray, ndarray, dict[str, ndarray]]], chrom: str) → Generator[tuple[ndarray, ndarray, dict[str, ndarray]], None, None]

validate_record_arrays(score: PositionScore, batches: Iterator[tuple[ndarray, ndarray, dict[str, ndarray]]], chrom: str) → Generator[tuple[ndarray, ndarray, dict[str, ndarray]], None, None]

validate_record_arrays(score: GenomicScore, batches: Iterator[RecordArrays], chrom: str, *, kind: str = "an allele score's") → Generator[RecordArrays, None, None]

validate_record_arrays(score: GenomicScore, batches: Iterator[RecordArrays], chrom: str, *, kind: str = "a fragment score's") → Generator[RecordArrays, None, None]

Yield a stream of raw column batches through, refusing a bad one.

The vectorized counterpart of validate_records(), over the batches the bulk scan is already pulling. It states the SAME ordering rule as its per-record twin, so a resource whose records are out of order is refused identically whichever path it was eligible for. Divergence between the two is what ADR 0008 records as the reason the shared class attribute was removed.

It also refuses a record whose end precedes its begin, as the per-record door does (see inverted_span_error()), although the backends the bulk path reads keep all but one such row out – a rule the backend happens to enforce is still the kind’s to state (ADR 0008; test_scan_array_door.py pins what each backend admits). The refusal is the per-record door’s, one suffix short of it for an allele record: the batches carry no ref/alt columns, so it cannot name them. A resource breaking both rules is refused for the same fault on both doors; _refuse_first_fault() says which.

chrom is what the batches were read for. A bulk scan reads one region, which lies within one contig, so the rules carry their ordering state across batches but never across contigs.

Refuses an unregistered kind, for the reason validate_records() gives.

validate_records(score: GenomicScore, records: Iterator[tuple[Any, ...]]) → Generator[tuple[Any, ...], None, None]

validate_records(score: PositionScore, records: Iterator[tuple[Any, ...]]) → Generator[tuple[Any, ...], None, None]

validate_records(score: GenomicScore, records: Iterator[Record], *, kind: str = "an allele score's") → Generator[Record, None, None]

validate_records(score: GenomicScore, records: Iterator[Record], *, kind: str = "a fragment score's") → Generator[Record, None, None]

Yield a raw record stream through, refusing a malformed one.

Dispatches on score’s class to the rule registered for that kind. A kind with no registration reaches this body and is refused: it would otherwise be validated by a rule nobody chose for it, which is the failure ADR 0008 exists to undo.

gain.genomic_resources.statistics.region_fold module

The scan-side fold that turns a kind’s regions into one statistic.

Not part of base_statistic because this fold reports through cli_errors, and the line that draws is a real one: base_statistic says what a statistic IS, and stays free of anything CLI-facing.

merge_regions(resource_id: str, regions: Iterable[MergeableRegion | None], build: Callable[[], S], what: str) → S | None

Fold one kind’s scanned regions into a single statistic.

None for a kind that contributed no regions at all – the caller writes no file for it, which is how “this resource has no such statistic” is said. A kind that produced nothing on SOME task is a None among the regions and is simply dropped.

what names the statistic in the failure line, and is the same string refuse_unmergeable() names it by – the two halves of what an operator reads on one failure. The failure is reported AND re-raised: reporting alone would leave the build writing a half-merged statistic.

Module contents