Skip to content

mecfs_bio.build_system.task.protein_coloc.trait_loci

Genome-wide significant loci of a trait, as windows around lead variants.

Greedy distance clumping: take the most significant remaining variant below the p threshold as a lead, open the window lead +- half_window_bp, remove every significant variant inside it, and repeat. Windows on one chromosome that then overlap are merged; a merged locus keeps the most significant of its leads. This needs no LD.

A locus whose window overlaps any excluded interval (by default the extended MHC, whose long range LD breaks the single-causal-variant assumption) gets status excluded_mhc and is reported but not tested.

p is the two-sided normal p-value of BETA / SE, kept in log space.

Functions:

  • find_trait_loci –

    One row per merged locus window, in LOCI_SCHEMA, sorted by chromosome and start.

Attributes:

LOCI_SCHEMA module-attribute 

LOCI_SCHEMA = pl.Schema(
    {
        LOCUS_ID_COL: pl.String,
        LOCUS_CHR_COL: pl.Int64,
        LOCUS_START_COL: pl.Int64,
        LOCUS_END_COL: pl.Int64,
        LOCUS_LEAD_POS_COL: pl.Int64,
        LOCUS_LEAD_EA_COL: pl.String,
        LOCUS_LEAD_NEA_COL: pl.String,
        LOCUS_LEAD_MLOG10P_COL: pl.Float64,
        STATUS_COL: pl.String,
    }
)

find_trait_loci 

find_trait_loci(
    trait: DataFrame,
    p_threshold: float,
    half_window_bp: int,
    exclude: Sequence[GenomicInterval],
) -> pl.DataFrame

One row per merged locus window, in LOCI_SCHEMA, sorted by chromosome and start.

Source code in mecfs_bio/build_system/task/protein_coloc/trait_loci.py 
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
def find_trait_loci(
    trait: pl.DataFrame,
    p_threshold: float,
    half_window_bp: int,
    exclude: Sequence[GenomicInterval],
) -> pl.DataFrame:
    """One row per merged locus window, in LOCI_SCHEMA, sorted by chromosome and start."""
    loci = _merge_overlapping(
        _greedy_leads(_significant(trait, p_threshold), half_window_bp)
    )
    return pl.DataFrame(
        [
            {
                LOCUS_ID_COL: f"chr{lead.chrom}:{lead.start}-{lead.end}",
                LOCUS_CHR_COL: lead.chrom,
                LOCUS_START_COL: lead.start,
                LOCUS_END_COL: lead.end,
                LOCUS_LEAD_POS_COL: lead.pos,
                LOCUS_LEAD_EA_COL: lead.ea,
                LOCUS_LEAD_NEA_COL: lead.nea,
                LOCUS_LEAD_MLOG10P_COL: lead.mlog10p,
                STATUS_COL: (
                    STATUS_EXCLUDED_MHC
                    if _overlaps_any(lead, exclude)
                    else STATUS_CANDIDATE
                ),
            }
            for lead in loci
        ],
        schema=LOCI_SCHEMA,
    )