Skip to content

mecfs_bio.build_system.task.genome_reference_harmonization.palindromes

Strand resolution of palindromic SNVs in untrusted tables: gwaslab's rule plus a distance check.

Rows arrive oriented so that NEA is the reference base, so a panel record with REF = NEA and ALT = EA describes the same variant. A variant is resolvable only when:

  • its EAF is at most palindrome_maf_threshold, or at least 1 - palindrome_maf_threshold;
  • the panel record exists;
  • the panel MAF is at most panel_maf_threshold;
  • after the decision, its EAF is within palindrome_max_af_distance of the panel AF.

It is kept when EAF and panel AF lie on the same side of 0.5, and strand-flipped (statistics flipped, alleles unchanged) otherwise. The distance check is not part of gwaslab's rule. On DecodeME, panels contradict each other's strand call for about one in ten palindromes beyond a distance of 0.1, against about one in five thousand within 0.02; see experiments/claude/gnomad_af_reference/palindrome_distance_vs_panel_disagreement.md.

Functions:

Attributes:

FREQUENCY_EPSILON module-attribute 

FREQUENCY_EPSILON = 1e-06

decide_palindrome_strands 

decide_palindrome_strands(
    rows: DataFrame,
    panel: DataFrame,
    options: GenomeReferenceHarmonizationOptions,
) -> pl.Series

One PalindromeDecision per row of rows (POS, EA, NEA, EAF), in row order.

Source code in mecfs_bio/build_system/task/genome_reference_harmonization/palindromes.py 
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
def decide_palindrome_strands(
    rows: pl.DataFrame,
    panel: pl.DataFrame,
    options: GenomeReferenceHarmonizationOptions,
) -> pl.Series:
    """One PalindromeDecision per row of rows (POS, EA, NEA, EAF), in row order."""
    records = panel.select(
        GWASLAB_POS_COL,
        pl.col(PANEL_REF_COL).alias(GWASLAB_NON_EFFECT_ALLELE_COL),
        pl.col(PANEL_ALT_COL).alias(GWASLAB_EFFECT_ALLELE_COL),
        pl.col(PANEL_AF_COL).cast(pl.Float64).alias(_PANEL_AF),
    )
    joined = rows.join(records, on=_KEYS, how="left", maintain_order="left")
    assert joined.height == rows.height, "duplicate panel records for a palindromic SNV"
    eaf = pl.col(GWASLAB_EFFECT_ALLELE_FREQ_COL).cast(pl.Float64)
    af = pl.col(_PANEL_AF)
    threshold = options.palindrome_maf_threshold
    eaf_informative = (eaf <= threshold + FREQUENCY_EPSILON) | (
        eaf >= 1 - threshold - FREQUENCY_EPSILON
    )
    panel_informative = (
        pl.min_horizontal(af, 1 - af) <= options.panel_maf_threshold + FREQUENCY_EPSILON
    )
    same_side = ((af < 0.5) & (eaf < 0.5)) | ((af > 0.5) & (eaf > 0.5))
    # The EAF the variant would have after the decision: unchanged if kept, 1 - EAF if flipped.
    decided_eaf = pl.when(same_side).then(eaf).otherwise(1 - eaf)
    close_to_panel = (decided_eaf - af).abs() <= (
        options.palindrome_max_af_distance + FREQUENCY_EPSILON
    )
    resolvable = (eaf_informative & panel_informative & close_to_panel).fill_null(False)
    return joined.select(
        pl.when(~resolvable)
        .then(pl.lit(PALINDROME_UNRESOLVED))
        .when(same_side)
        .then(pl.lit(PALINDROME_KEEP))
        .otherwise(pl.lit(PALINDROME_STRAND_FLIP))
        .alias(_DECISION)
    )[_DECISION]