RYR1 Variant Pathogenicity Predictor

ENS210 (Computational Biology) asked for a predictor of variant pathogenicity. Most tools like SIFT, PolyPhen-2, and REVEL are trained genome-wide, but RYR1 - the calcium-release channel gene behind malignant hyperthermia and central core disease - has pathogenic variants that cluster tightly in a few structural hotspots. I built a gene-specific model instead: a logistic regression over 245 curated ClinVar missense variants that leans on where in the protein a substitution falls, not just what it is.

Why not just compare FASTA sequences

RYR1 is ~5,000 amino acids long, and almost every variant differs from the reference by exactly one residue - aligning full sequences just confirms "one position changed," which HGVS notation already tells you for free. It also throws away everything that actually predicts pathogenicity: allele frequency, conservation, review status, disease association. So instead of sequence alignment, every variant is parsed straight from its HGVS protein notation (L13R, p.Arg4893Gln, ...) into a (ref, position, alt) triple and scored against ClinVar/gnomAD-derived features.

A position-aware log-odds score

The main gene-specific feature is a per-position, per-substitution log-odds score built from ClinVar counts, with a pseudo-count alpha to avoid divide-by-zero on unseen substitutions:

def build_log_odds(df: pd.DataFrame):
    benign = np.zeros((max_pos + 1, len(AA_LIST)), dtype=float)
    path = np.zeros_like(benign)
 
    for _, r in df.iterrows():
        pos = int(r["prot_pos"])
        j = aa_to_idx.get(str(r["alt_aa"]).upper())
        if j is None or pos < 1 or pos > max_pos:
            continue
        if int(r["label"]) == 1:
            path[pos, j] += 1.0
        else:
            benign[pos, j] += 1.0
 
    benign += alpha
    path += alpha
    return np.log(path / benign), benign, path

A positive score means that substitution has historically skewed pathogenic at that exact position; zero means no prior evidence either way, which just pushes the model to lean on the other features - BLOSUM62 and Grantham distance for how chemically disruptive the substitution is, plus REVEL/PolyPhen/SIFT/CADD/phyloP where gnomAD has them (with explicit missingness flags where it doesn't).

Keeping the log-odds table leakage-free

The log-odds table is fit from labels, so recomputing it on the full dataset before cross-validation would leak test-set pathogenicity straight into a training feature. Every fold rebuilds the table from its training split only, then scores both splits against that fold-specific table:

for fold, (tr_idx, te_idx) in enumerate(skf.split(missense, missense["label"]), 1):
    train_df, test_df = missense.iloc[tr_idx].copy(), missense.iloc[te_idx].copy()
 
    log_odds_tr, _, _ = build_log_odds(train_df)
    train_df = attach_base_features(train_df, log_odds_tr)
    test_df = attach_base_features(test_df, log_odds_tr)
    # ... scale, fit LogisticRegression, score on the held-out fold

Standardized coefficients from the final model confirm the log-odds score dominates the decision - by a wide margin over REVEL, SIFT, and the biochemical substitution metrics.

Beating genome-wide predictors on their home gene

Evaluated on the same 245-variant set, the gene-specific model reached a leakage-free 5-fold ROC AUC of 0.943 ± 0.029 (accuracy 0.882, precision 0.944, recall 0.911), clearing every genome-wide baseline evaluated on the same variants - REVEL came closest at 0.90, with PolyPhen-2, SIFT, and CADD trailing between 0.79 and 0.83.

Using it

The trained model ships with the repo, so scoring a variant against the canonical RYR1 reference is one command:

python scripts/predict.py L13R R614C p.Arg4893Gln

Each variant gets a p_pathogenic probability, a benign / not_benign / uncertain call, and an evidence level reflecting how much ClinVar support exists at that exact position and substitution - so a confident-looking probability backed by zero training examples still reads as low-evidence.

Result

A gene-specific logistic regression that outperforms genome-wide predictors on RYR1 missense variants (ROC AUC 0.943 ± 0.029 vs. 0.90 for the best baseline), driven mostly by a leakage-free positional log-odds feature built directly from ClinVar, with a CLI that scores arbitrary HGVS variants against the shipped model.