mecfs_bio.build_system.task.protein_coloc.candidate_selection

Candidate (trait locus, protein) pairs for colocalization, and their selection within a budget.

A protein is a candidate at a trait locus when one of its hit segments (from the protein hit table) overlaps the locus window on the same chromosome. Excluded loci get no candidates.

Each candidate is classified by where the protein's gene lies (Sun et al. 2023 ST3, GRCh38): - cis: the gene interval lies within cis_window_bp of the locus window; - trans: the gene is elsewhere; - unknown: the protein has no gene coordinates (for example a multi-gene complex).

Each pair gets a score, and pairs are taken in order of decreasing score: SCORE = SCORE_TRAIT + SCORE_PROTEIN - SCORE_PENALTY SCORE_TRAIT = min(trait lead -log10 p, clip) SCORE_PROTEIN = min(protein lead -log10 p, clip) SCORE_PENALTY = w * log2(1 + N_OTHER_PROTEINS) when the pair is not cis, else 0 N_OTHER_PROTEINS is the number of other proteins with a segment overlapping the same window. The penalty demotes trans pairs at pleiotropic hotspots, where many proteins share a pQTL; cis pairs are never penalized. Unknown is penalized like trans. Strong pairs often tie once both terms are clipped, so ties are broken by GENE_LEAD_DISTANCE_BP (distance from the trait lead to the protein's gene, nearest first, genes on other chromosomes or without coordinates last), then by the unclipped protein -log10 p, then by PROTEIN_ID and LOCUS_ID. A pair is selected while fewer than max_pairs are selected and its locus has fewer than max_per_locus selected pairs.

Classes:

• SelectionConfig –

Functions:

• build_candidates –

  One row per (candidate locus, protein with a segment overlapping its window).

• select_candidates –

  Score the candidates and mark each selected or not_selected_budget, in ranked order.

Attributes:

• CANDIDATE_SCHEMA –

CANDIDATE_SCHEMA module-attribute

CANDIDATE_SCHEMA = pl.Schema(
    {
        LOCUS_ID_COL: pl.String,
        PROTEIN_ID_COL: pl.String,
        HIT_GENE_COL: pl.String,
        PROTEIN_MLOG10P_COL: pl.Float64,
        CIS_STATUS_COL: pl.String,
        N_OTHER_PROTEINS_COL: pl.Int64,
        GENE_LEAD_DISTANCE_COL: pl.Int64,
    }
)

SelectionConfig

Attributes:

• cis_window_bp (int) –
• max_pairs (int) –
• max_per_locus (int) –
• mlog10p_clip (float) –
• trans_penalty_weight (float) –

cis_window_bp class-attribute instance-attribute

cis_window_bp: int = 1000000

max_pairs class-attribute instance-attribute

max_pairs: int = 30

max_per_locus class-attribute instance-attribute

max_per_locus: int = 5

mlog10p_clip class-attribute instance-attribute

mlog10p_clip: float = 50.0

trans_penalty_weight class-attribute instance-attribute

trans_penalty_weight: float = 2.0

build_candidates

build_candidates(
    loci: DataFrame,
    hit_segments: DataFrame,
    gene_coords: Mapping[ProteinId, GenomicInterval],
    cis_window_bp: int,
) -> pl.DataFrame

One row per (candidate locus, protein with a segment overlapping its window).

Source code in mecfs_bio/build_system/task/protein_coloc/candidate_selection.py

def build_candidates(
    loci: pl.DataFrame,
    hit_segments: pl.DataFrame,
    gene_coords: Mapping[ProteinId, GenomicInterval],
    cis_window_bp: int,
) -> pl.DataFrame:
    """One row per (candidate locus, protein with a segment overlapping its window)."""
    windows = loci.filter(pl.col(STATUS_COL) == STATUS_CANDIDATE).select(
        LOCUS_ID_COL,
        pl.col(LOCUS_CHR_COL).cast(pl.Int64),
        LOCUS_START_COL,
        LOCUS_END_COL,
        LOCUS_LEAD_POS_COL,
    )
    segments = hit_segments.select(
        PROTEIN_ID_COL,
        HIT_GENE_COL,
        pl.col(HIT_CHR_COL).cast(pl.Int64).alias(LOCUS_CHR_COL),
        HIT_SEG_START_COL,
        HIT_SEG_END_COL,
        HIT_LEAD_MLOG10P_COL,
    )
    overlapping = windows.join(segments, on=LOCUS_CHR_COL, how="inner").filter(
        pl.col(HIT_SEG_START_COL) <= pl.col(LOCUS_END_COL),
        pl.col(HIT_SEG_END_COL) >= pl.col(LOCUS_START_COL),
    )
    pairs = overlapping.group_by(
        LOCUS_ID_COL,
        LOCUS_CHR_COL,
        LOCUS_START_COL,
        LOCUS_END_COL,
        LOCUS_LEAD_POS_COL,
        PROTEIN_ID_COL,
    ).agg(
        pl.col(HIT_GENE_COL).first(),
        pl.col(HIT_LEAD_MLOG10P_COL).max().alias(PROTEIN_MLOG10P_COL),
    )
    return (
        pairs.with_columns(
            (pl.len().over(LOCUS_ID_COL) - 1).cast(pl.Int64).alias(N_OTHER_PROTEINS_COL)
        )
        .join(_gene_coords_frame(gene_coords), on=PROTEIN_ID_COL, how="left")
        .with_columns(_cis_status_expr(cis_window_bp), _gene_lead_distance_expr())
        .select(CANDIDATE_SCHEMA.names())
        .cast(CANDIDATE_SCHEMA)
        .sort(LOCUS_ID_COL, PROTEIN_ID_COL)
    )

select_candidates

select_candidates(
    candidates: DataFrame,
    loci: DataFrame,
    config: SelectionConfig,
) -> pl.DataFrame

Score the candidates and mark each selected or not_selected_budget, in ranked order.

Source code in mecfs_bio/build_system/task/protein_coloc/candidate_selection.py

def select_candidates(
    candidates: pl.DataFrame, loci: pl.DataFrame, config: SelectionConfig
) -> pl.DataFrame:
    """Score the candidates and mark each selected or not_selected_budget, in ranked order."""
    ranked = (
        candidates.join(
            loci.select(LOCUS_ID_COL, LOCUS_LEAD_MLOG10P_COL),
            on=LOCUS_ID_COL,
            how="left",
        )
        .with_columns(_score_columns(config))
        .with_columns(
            (
                pl.col(SCORE_TRAIT_COL)
                + pl.col(SCORE_PROTEIN_COL)
                - pl.col(SCORE_PENALTY_COL)
            ).alias(SCORE_COL)
        )
        .drop(LOCUS_LEAD_MLOG10P_COL)
        .sort(
            [
                SCORE_COL,
                GENE_LEAD_DISTANCE_COL,
                PROTEIN_MLOG10P_COL,
                PROTEIN_ID_COL,
                LOCUS_ID_COL,
            ],
            descending=[True, False, True, False, False],
            nulls_last=True,
        )
    )
    return ranked.with_columns(
        pl.Series(
            STATUS_COL,
            _budget_statuses(ranked[LOCUS_ID_COL].to_list(), config),
            dtype=pl.String,
        )
    )