Sequence statistics

Where the per-sequence and per-position statistics are computed, and where the constants that decide how far the per-position checks reach are defined.

seq_stats

Per-class sequence statistics: the features the class comparison compares.

One SequenceStatistics holds the sequences of a single class and the features computed from them - lengths, GC content, nucleotide and dinucleotide composition, per-position composition, duplication levels. The comparison in genomic_benchmarks_qc.utils.testing reads these, and so do the plots and the HTML report.

Two windows bound the per-position features, because a position is compared only on the sequences long enough to reach it. scored_end_position is as far as a position may be flagged: it needs a cohort of at least MIN_SEQUENCES_PER_CLASS sequences and at least DEFAULT_MIN_COVERAGE of the class. It is also the window the figures draw. end_position reaches further, as far as MIN_SEQUENCES_PER_REPORTED_POSITION sequences: the positions in between are named in the report as Unknown, so a dataset whose sequences simply end is not read as one whose tail was too thin to compare.

What the two windows mean for reading a report is on the per-position checks page.

SequenceStatistics

SequenceStatistics(
    sequences,
    filename,
    filepath,
    label,
    seq_column=None,
    end_position=None,
    slug=None,
    min_coverage=DEFAULT_MIN_COVERAGE,
)

The sequences of one class, and the statistics computed from them.

Statistics are computed on demand by compute, not on construction, and cached in stats.

Hold one class's sequences together with how to identify it.

Parameters:

Name	Type	Description	Default
sequences	list[str]	The sequences of this class, uppercased by the reader.	required
filename	str	Name of the file they came from, shown in the report.	required
filepath	str	Full path they came from, shown in the report.	required
label	str	The class name, shown verbatim in reports and plots.	required
seq_column	str | None	Sequence column they came from, or None for FASTA.	None
end_position	int | None	Last position the per-position checks reach, 1-based and inclusive. Defaults to the last position at least MIN_SEQUENCES_PER_REPORTED_POSITION of these sequences reach. It cannot widen what gets flagged - the scored window decides that - so an explicit value only ever trims.	None
slug	str | None	Path form of label; derived from it when not given, but normally passed in by the caller, which is the only place that can tell whether it collides with another class.	None
min_coverage	float	Fraction of these sequences that must reach a position before it may set a flag, on top of the MIN_SEQUENCES_PER_CLASS sequences every scored position needs. 0 leaves only that count. Default: 0.25 (DEFAULT_MIN_COVERAGE).	DEFAULT_MIN_COVERAGE

Source code in src/genomic_benchmarks_qc/utils/seq_stats.py

def __init__(self, sequences: list[str], filename: str, filepath: str, label: str,
             seq_column: str | None = None, end_position: int | None = None,
             slug: str | None = None, min_coverage: float = DEFAULT_MIN_COVERAGE):
    """Hold one class's sequences together with how to identify it.

    Args:
        sequences: The sequences of this class, uppercased by the reader.
        filename: Name of the file they came from, shown in the report.
        filepath: Full path they came from, shown in the report.
        label: The class name, shown verbatim in reports and plots.
        seq_column: Sequence column they came from, or None for FASTA.
        end_position: Last position the per-position checks reach, 1-based and
            inclusive. Defaults to the last position at least
            [MIN_SEQUENCES_PER_REPORTED_POSITION][genomic_benchmarks_qc.utils.seq_stats.MIN_SEQUENCES_PER_REPORTED_POSITION]
            of these sequences reach. It cannot widen what gets flagged - the
            scored window decides that - so an explicit value only ever trims.
        slug: Path form of `label`; derived from it when not given, but normally
            passed in by the caller, which is the only place that can tell whether it
            collides with another class.
        min_coverage: Fraction of these sequences that must reach a position before it
            may set a flag, on top of the
            [MIN_SEQUENCES_PER_CLASS][genomic_benchmarks_qc.utils.testing.MIN_SEQUENCES_PER_CLASS]
            sequences every scored position needs. 0 leaves only that count.
            Default: `0.25`
            ([DEFAULT_MIN_COVERAGE][genomic_benchmarks_qc.utils.seq_stats.DEFAULT_MIN_COVERAGE]).
    """
    self.filename = filename
    self.filepath = filepath
    self.label = label
    # `label` is shown verbatim in reports and plots; `slug` is the
    # filesystem-safe, collision-free form used to build report paths.
    self.slug = slug if slug is not None else slugify(label)
    self.seq_column = seq_column
    self.sequences = sequences
    self.min_coverage = min_coverage

    self.end_position = end_position
    """Last position the per-position checks reach, resolved by `compute`."""

    self.scored_end_position = None
    """Last position that may be flagged, and the last one the figures draw.

    Resolved by `compute`, alongside `end_position`.
    """

    self.stats = {}
    """The statistics `compute` produces, empty until it has run."""

end_position instance-attribute

end_position = end_position

Last position the per-position checks reach, resolved by compute.

scored_end_position instance-attribute

scored_end_position = None

Last position that may be flagged, and the last one the figures draw.

Resolved by compute, alongside end_position.

stats instance-attribute

stats = {}

The statistics compute produces, empty until it has run.

compute

compute()

Compute this class's statistics and resolve its per-position windows.

Results are cached on self.stats. The statistics dictionary holds:

• Filename: str
• Filepath: str
• Label: str, or 'N/A'
• Sequence column: str, or 'N/A'
• Number of sequences: int
• Number of bases: int
• Unique bases: list of str
• %GC content: float
• Number of sequences left after deduplication: int
• Empty sequences: int
• Per sequence nucleotide content: pd.DataFrame (index: sequence_id, columns: nucleotides, values: frequency)
• Per sequence dinucleotide content: pd.DataFrame (index: sequence_id, columns: dinucleotides, values: frequency)
• Per position nucleotide content: pd.DataFrame (index: position, columns: nucleotides, values: frequency)
• Per position reversed nucleotide content: pd.DataFrame (index: position, columns: nucleotides, values: frequency)
• Per sequence GC content: dict pd.DataFrame (index: sequence_id, columns: GC content (%), values: GC content)
• Sequence lengths: pd.DataFrame (index: sequence_id, columns: Length, values: length of the sequence)
• Sequence duplication levels: dict {sequence: extra copies}, holding only the sequences that occur more than once

Returns:

Type	Description
dict	That statistics dictionary, and end_position - the last position
int	the per-position checks reach.

Source code in src/genomic_benchmarks_qc/utils/seq_stats.py

def compute(self) -> tuple[dict, int]:
    """Compute this class's statistics and resolve its per-position windows.

    Results are cached on `self.stats`. The statistics dictionary holds:

    - Filename: str
    - Filepath: str
    - Label: str, or 'N/A'
    - Sequence column: str, or 'N/A'
    - Number of sequences: int
    - Number of bases: int
    - Unique bases: list of str
    - %GC content: float
    - Number of sequences left after deduplication: int
    - Empty sequences: int
    - Per sequence nucleotide content: pd.DataFrame
      (index: sequence_id, columns: nucleotides, values: frequency)
    - Per sequence dinucleotide content: pd.DataFrame
      (index: sequence_id, columns: dinucleotides, values: frequency)
    - Per position nucleotide content: pd.DataFrame
      (index: position, columns: nucleotides, values: frequency)
    - Per position reversed nucleotide content: pd.DataFrame
      (index: position, columns: nucleotides, values: frequency)
    - Per sequence GC content: dict pd.DataFrame
      (index: sequence_id, columns: GC content (%), values: GC content)
    - Sequence lengths: pd.DataFrame
      (index: sequence_id, columns: Length, values: length of the sequence)
    - Sequence duplication levels: dict {sequence: extra copies},
      holding only the sequences that occur more than once

    Returns:
        That statistics dictionary, and `end_position` - the last position
        the per-position checks reach.
    """
    message = f"Computing statistics for {self.filename}"
    if self.label is not None:
        message += f", label {self.label}"
    if self.seq_column is not None:
        message += f", sequence column: {self.seq_column}"
    logger.info(message)

    self._compute_basic_statistics()
    self._compute_per_sequence_statistics()
    self._compute_sequence_duplication_levels()

    self._resolve_position_windows()

    return self.stats, self.end_position

coverage_curve

coverage_curve(end_position)

Fraction reaching each 1-based position up to and including end_position.

The denominator behind every per-position statistic, as a curve. The report needs it three times over - the figure draws it, the interactive viewer carries it, and the prose quotes single points of it - and it used to be worked out three ways, one of them a Python loop over the class per position. One answer, so the figure and the viewer cannot come to draw different curves.

Source code in src/genomic_benchmarks_qc/utils/seq_stats.py

def coverage_curve(self, end_position: int) -> np.ndarray:
    """Fraction reaching each 1-based position up to and including `end_position`.

    The denominator behind every per-position statistic, as a curve. The
    report needs it three times over - the figure draws it, the interactive
    viewer carries it, and the prose quotes single points of it - and it used
    to be worked out three ways, one of them a Python loop over the class per
    position. One answer, so the figure and the viewer cannot come to draw
    different curves.
    """
    return _coverage_from_lengths(
        self.stats['Sequence lengths'].values.flatten(),
        np.arange(1, end_position + 1),
    )

coverage_at

coverage_at(position)

Fraction of this class's sequences that reach position (1-based).

This is the denominator behind every per-position statistic at that position, and it is what the report shows so a reader can tell how much data stands behind the far end of the per-position plots.

Source code in src/genomic_benchmarks_qc/utils/seq_stats.py

def coverage_at(self, position: int) -> float:
    """Fraction of this class's sequences that reach `position` (1-based).

    This is the denominator behind every per-position statistic at that
    position, and it is what the report shows so a reader can tell how much
    data stands behind the far end of the per-position plots.
    """
    lengths = self.stats['Sequence lengths'].values.flatten()
    return float(_coverage_from_lengths(lengths, np.array([position]))[0])

cohort_floor

cohort_floor(stats1, stats2)

The share of a class a position's cohort has to reach to be flagged.

Each class requires its own number of sequences behind a position - the larger of MIN_SEQUENCES_PER_CLASS and min_coverage of the class - so as a share of a class the floor differs between the two, and the binding one is the larger share: a position has to clear the floor in both classes.

Returns:

Type	Description
float	Fraction of a class, or 0.0 when neither class has a floor.

Source code in src/genomic_benchmarks_qc/utils/seq_stats.py

def cohort_floor(stats1, stats2) -> float:
    """The share of a class a position's cohort has to reach to be flagged.

    Each class requires its own number of sequences behind a position - the
    larger of `MIN_SEQUENCES_PER_CLASS` and `min_coverage` of the class - so as a
    share of a class the floor differs between the two, and the binding one is
    the larger share: a position has to clear the floor in both classes.

    Returns:
        Fraction of a class, or 0.0 when neither class has a floor.
    """
    binding = 0.0
    for stats in (stats1, stats2):
        count = stats.stats['Number of sequences']
        if not count:
            continue
        needed = stats._required_cohort(count)
        binding = max(binding, needed / count)
    return binding

DEFAULT_MIN_COVERAGE module-attribute

DEFAULT_MIN_COVERAGE = 0.25

MIN_SEQUENCES_PER_REPORTED_POSITION module-attribute

MIN_SEQUENCES_PER_REPORTED_POSITION = 50