Flagging

How a statistic becomes a Pass, Warning, Fail or Unknown.

testing

Turn a per-class statistic into a Pass, Warning, Fail or Unknown flag.

Every check asks one question: can a single feature tell the two classes apart? The feature's values are read as if they were a classifier's scores, and the AU-ROC is the answer - 0.5 for a feature that says nothing about the class, 1.0 for one that settles it by itself. A fixed boundary makes that a flag: Pass up to 0.6, Warning up to 0.7, Fail above that.

Two floors sit in front of the scoring, and a check has to clear both. Nothing is scored on fewer than MIN_SEQUENCES_PER_CLASS sequences per class - counted per position for the per-position checks, since a position is compared only on the sequences long enough to have it. And a per-position check additionally stops where fewer than DEFAULT_MIN_COVERAGE of a class reaches the position. A check that misses either floor reports Unknown, which says the comparison was not made, not that it came out clean.

flag_significant_differences runs every check on one pair of classes, and is what both commands call. Why the boundaries and the floors are where they are is on the How a flag is decided page.

flag_significant_differences

flag_significant_differences(stats1, stats2)

Run every check on one pair of classes.

Parameters:

Name	Type	Description	Default
stats1	SequenceStatistics	Computed statistics for the first class.	required
stats2	SequenceStatistics	Computed statistics for the second class.	required

Returns:

Type	Description
dict	Tuple of (summary_statuses, failed_by_feature). summary_statuses is
dict	every flag and metric in report order: the headline checks first, then
tuple[dict, dict]	their sub-checks. failed_by_feature is what the plots shade - only the
tuple[dict, dict]	sub-checks that came out Warning or Fail, keyed by feature and then by
tuple[dict, dict]	what was flagged: { 'Per sequence nucleotide content': {'A': 'Warning', 'G': 'Fail'}, 'Per sequence dinucleotide content': {'GG': 'Fail'}, 'Per position nucleotide content': {'A': {52: 'Warning'}, 'G': {66: 'Fail', 70: 'Fail'}}, 'Per position reversed nucleotide content': {}, }

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

def flag_significant_differences(stats1: 'SequenceStatistics',
                                 stats2: 'SequenceStatistics') -> tuple[dict, dict]:
    """Run every check on one pair of classes.

    Args:
        stats1: Computed statistics for the first class.
        stats2: Computed statistics for the second class.

    Returns:
        Tuple of `(summary_statuses, failed_by_feature)`. `summary_statuses` is
        every flag and metric in report order: the headline checks first, then
        their sub-checks. `failed_by_feature` is what the plots shade - only the
        sub-checks that came out Warning or Fail, keyed by feature and then by
        what was flagged:

            {
                'Per sequence nucleotide content': {'A': 'Warning', 'G': 'Fail'},
                'Per sequence dinucleotide content': {'GG': 'Fail'},
                'Per position nucleotide content':
                    {'A': {52: 'Warning'}, 'G': {66: 'Fail', 70: 'Fail'}},
                'Per position reversed nucleotide content': {},
            }
    """
    results = {}

    ordered_stats = [
        'Unique bases',
        'Sequence Duplications within Labels',
        'Duplicate Sequences between Labels',
        'Sequence lengths',
        'Per sequence GC content',
        'Per sequence nucleotide content',
        'Per sequence dinucleotide content',
        'Per position nucleotide content',
        'Per position reversed nucleotide content',
    ]

    all_results = {}

    all_results['Unique bases'] = {'Flag': _flag_unique_bases(stats1, stats2)}
    all_results['Sequence Duplications within Labels'] = _flag_duplicate_sequences(stats1, stats2)
    all_results['Duplicate Sequences between Labels'] = {
        'Flag': _flag_duplication_between_datasets(stats1.sequences, stats2.sequences)
    }

    model_results = direct_feature_model(stats1, stats2)
    all_results.update(model_results)

    _warn_about_unscored_checks(stats1, stats2, all_results, ordered_stats)

    # Order: aggregates first, then details
    for stat_name in ordered_stats:
        if stat_name in all_results:
            results[stat_name] = all_results[stat_name]

    for stat_name in ordered_stats:
        for key in all_results:
            if key.startswith(f"{stat_name} - ") and key != stat_name:
                results[key] = all_results[key]

    # Build failed_by_feature dict for visualization
    failed_by_feature = _extract_failed_features(all_results)

    return results, failed_by_feature

direct_feature_model

direct_feature_model(stats1, stats2)

Score every feature of one class against the same feature of another.

Parameters:

Name	Type	Description	Default
stats1	SequenceStatistics	Computed statistics for the first class.	required
stats2	SequenceStatistics	Computed statistics for the second class.	required

Returns:

Type	Description
dict	One entry per check and per sub-check, keyed by name - `'Sequence
dict	lengths','Per sequence nucleotide content - A','Per position
dict	nucleotide content - A position 52'` - each holding its AU-ROC, AU-PR,
dict	Accuracy and Flag. A check made of sub-checks also gets a headline entry
dict	under its own name, holding the worst of them.

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

def direct_feature_model(stats1: 'SequenceStatistics',
                         stats2: 'SequenceStatistics') -> dict:
    """Score every feature of one class against the same feature of another.

    Args:
        stats1: Computed statistics for the first class.
        stats2: Computed statistics for the second class.

    Returns:
        One entry per check and per sub-check, keyed by name - `'Sequence
        lengths'`, `'Per sequence nucleotide content - A'`, `'Per position
        nucleotide content - A position 52'` - each holding its AU-ROC, AU-PR,
        Accuracy and Flag. A check made of sub-checks also gets a headline entry
        under its own name, holding the worst of them.
    """
    indices_1 = np.arange(len(stats1.sequences))
    indices_2 = np.arange(len(stats2.sequences))

    results = {}

    # Scalar features
    results['Sequence lengths'] = _score_scalar_feature(
        stats1.stats['Sequence lengths'],
        stats2.stats['Sequence lengths'],
        'Sequence lengths',
        indices_1, indices_2,
    )

    results['Per sequence GC content'] = _score_scalar_feature(
        stats1.stats['Per sequence GC content'],
        stats2.stats['Per sequence GC content'],
        'Per sequence GC content',
        indices_1, indices_2,
    )

    # DataFrame features
    results.update(_score_dataframe_features(
        stats1.stats['Per sequence nucleotide content'],
        stats2.stats['Per sequence nucleotide content'],
        'Per sequence nucleotide content',
        indices_1, indices_2,
    ))

    results.update(_score_dataframe_features(
        stats1.stats['Per sequence dinucleotide content'],
        stats2.stats['Per sequence dinucleotide content'],
        'Per sequence dinucleotide content',
        indices_1, indices_2,
    ))

    # Position features (forward)
    bases = sorted(set(stats1.stats['Unique bases']) | set(stats2.stats['Unique bases']))
    end_position, scored_end_position = position_windows(stats1, stats2)

    pos_results, per_base_agg = _score_position_features(
        stats1.sequences, stats2.sequences, bases,
        'Per position nucleotide content',
        end_position=end_position, scored_end_position=scored_end_position,
        reverse=False,
    )
    results.update(pos_results)
    if per_base_agg:
        results['Per position nucleotide content'] = _aggregate_worst_case_metrics(
            per_base_agg.values())

    # Position features (reverse)
    pos_results_rev, per_base_agg_rev = _score_position_features(
        stats1.sequences, stats2.sequences, bases,
        'Per position reversed nucleotide content',
        end_position=end_position, scored_end_position=scored_end_position,
        reverse=True,
    )
    results.update(pos_results_rev)
    if per_base_agg_rev:
        results['Per position reversed nucleotide content'] = (
            _aggregate_worst_case_metrics(per_base_agg_rev.values()))

    return results

position_windows

position_windows(stats1, stats2)

The per-position windows a comparison of two classes runs in.

Each class resolves its own windows from its own sequence lengths, and a position belongs to a window only where it belongs to it in both: a position that half of one class does not reach cannot be flagged on the strength of the other class reaching it.

Returns:

Type	Description
int	Tuple of (end_position, scored_end_position), 1-based and inclusive - the
int	last position reported on and the last position allowed to set a flag,
tuple[int, int]	which is also the last position drawn.

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

def position_windows(stats1, stats2) -> tuple[int, int]:
    """The per-position windows a comparison of two classes runs in.

    Each class resolves its own windows from its own sequence lengths, and a
    position belongs to a window only where it belongs to it in both: a position
    that half of one class does not reach cannot be flagged on the strength of
    the other class reaching it.

    Returns:
        Tuple of (end_position, scored_end_position), 1-based and inclusive - the
        last position reported on and the last position allowed to set a flag,
        which is also the last position drawn.
    """
    return (min(stats1.end_position, stats2.end_position),
            min(stats1.scored_end_position, stats2.scored_end_position))

MIN_SEQUENCES_PER_CLASS module-attribute

MIN_SEQUENCES_PER_CLASS = 250