Benchmarking Classical Machine Learning and Neural Methods for Variant Pathogenicity Prediction on ClinVar Metadata
Benchmarking Classical Machine Learning and Neural Methods for Variant Pathogenicity Prediction on ClinVar Metadata
Abstract
Predicting whether a genomic variant is pathogenic or benign is a central problem in clinical genomics. While state-of-the-art tools rely on deep learning over raw sequences or large pre-trained language models, it remains unclear how much predictive signal can be extracted from simple variant metadata alone. In this study, we benchmark five methods — Logistic Regression, Random Forest, XGBoost, Support Vector Machine (SVM), and a Multilayer Perceptron (MLP) — on 1,571,614 single nucleotide variants (SNVs) from ClinVar, using only 22 metadata-derived features (molecular consequence, chromosomal location, submission provenance, and allele properties). We employ a gene-held-out train/test split to prevent gene-level data leakage. Our results show that Random Forest, XGBoost, and MLP all converge to nearly identical performance (AUROC ~0.807, AUPRC ~0.692), while Logistic Regression lags substantially (AUROC 0.755, AUPRC 0.520), indicating genuine non-linearity in the feature space. Notably, XGBoost achieves top-tier performance in just 17 seconds of training — a 47× speedup over MLP (802s) — making it the practical method of choice. Feature importance analysis reveals that submission metadata (number of submitters, review quality score) dominates over biological consequence features, surfacing an important circularity concern in ClinVar-based benchmarks.
1. Introduction
Variant effect prediction — determining whether a given genomic variant contributes to disease — is a core challenge in precision medicine. Hundreds of computational tools exist for this task, ranging from conservation-based scores (SIFT, PolyPhen-2, GERP++) to deep learning models trained on protein structure or evolutionary sequence data. However, a simpler question has received less systematic attention: how much predictive power exists in the administrative and submission metadata that accompanies variant records in curated databases like ClinVar?
ClinVar is the primary public repository of human variant-disease interpretations, hosting millions of submissions with structured annotations including molecular consequence, review status, number of submitters, and variant origin. These metadata fields are collected as part of the curation process, not derived from sequence analysis, making them a distinct and largely orthogonal information source compared to conservation- or structure-based features.
In this work, we conduct a focused benchmark of five machine learning methods — ranging from a linear baseline to a neural network — applied exclusively to ClinVar metadata features for the binary task of classifying variants as pathogenic or benign. Our contributions are:
- A reproducible pipeline for streaming, filtering, and featurizing 1.5M+ ClinVar SNVs without requiring external databases.
- A gene-held-out evaluation protocol that avoids the gene-level data leakage present in naive random splits.
- A systematic comparison of five methods with full metrics (AUROC, AUPRC, F1, MCC, training time).
- An interpretability analysis via feature importance, revealing the dominance of submission metadata and its implications for benchmark validity.
2. Data and Features
2.1 Dataset
We downloaded the ClinVar variant summary file (variant_summary.txt.gz) from the NCBI FTP server (accessed April 2026). After filtering to:
- Assembly: GRCh38
- Variant type: single nucleotide variant (SNV)
- Clinical significance: strictly Pathogenic/Likely pathogenic (label = 1) or Benign/Likely benign (label = 0), excluding VUS and conflicting entries
- Review status: excluding entries with "no assertion criteria provided", "no assertion provided", or "no classification provided"
- Gene symbol: non-missing
we obtained 1,571,614 SNVs: 296,030 pathogenic and 1,275,584 benign (roughly 1:4.3 ratio).
2.2 Train/Test Split
To prevent gene-level data leakage — a common issue in variant benchmarking where multiple variants from the same gene appear in both train and test sets — we employed a gene-held-out split. All unique gene symbols were randomly partitioned such that 20% of genes formed an exclusive test set and the remaining 80% formed the training set. This yielded:
- Train: 1,247,265 variants
- Test: 324,349 variants (65,008 pathogenic, 259,341 benign)
2.3 Feature Engineering
We derived 22 binary and numeric features entirely from ClinVar metadata columns, requiring no external database joins:
| Feature Group | Features | Description |
|---|---|---|
| Molecular consequence | cons_missense, cons_nonsense, cons_synonymous, cons_splice, cons_frameshift, cons_stop, cons_start, cons_utr, cons_intronic |
Keyword match in MolecularConsequence / Name columns |
| Chromosomal location | chrom_num, is_chrX, is_chrY, is_chrMT |
Parsed from Chromosome column |
| Submission provenance | num_submitters, review_score |
NumberSubmitters (numeric); ReviewStatus mapped to ordinal 0–4 |
| Variant origin | origin_germline, origin_somatic, origin_de_novo |
Binary flags from OriginSimple |
| Allele properties | has_rsid, is_transition, is_transversion, ref_is_CG |
From RS# and allele columns |
Missing values (primarily in num_submitters) were imputed with the column mean. All features were standardized (zero mean, unit variance) for Logistic Regression, SVM, and MLP; tree-based models received raw imputed features.
3. Methods
We trained and evaluated five models on the same feature matrix:
Logistic Regression (LR): L2-regularized logistic regression (C=1.0, max_iter=1000). Serves as the linear baseline. Class imbalance handled via class_weight='balanced'.
Random Forest (RF): Ensemble of 300 decision trees with bootstrap sampling (n_estimators=300, class_weight='balanced'). Captures non-linear interactions and provides feature importance via mean decrease in impurity.
XGBoost: Gradient boosted trees with 300 estimators, learning rate 0.05, max depth 6. Class imbalance addressed via scale_pos_weight set to the negative-to-positive ratio.
Support Vector Machine (SVM): RBF kernel SVM with class_weight='balanced'.
Multilayer Perceptron (MLP):
A three-hidden-layer feedforward network with architecture: 128 -> 64 -> 32 -> 1 with Batch Normalization and Dropout (0.3) after each of the first two hidden layers. Trained for 30 epochs with the Adam optimizer (lr=1e-3, weight decay=1e-4) and a BCEWithLogitsLoss weighted by the class imbalance ratio. Learning rate reduced on plateau.
All experiments used a fixed random seed (42) for reproducibility and were run on CPU only.
4. Results
4.1 Quantitative Comparison
Table 1 summarizes performance across all five methods on the held-out test set.
Table 1: Model performance on gene-held-out test set (n = 324,349 SNVs)
| Model | AUROC | AUPRC | F1 | MCC | Train Time (s) |
|---|---|---|---|---|---|
| Logistic Regression | 0.7551 | 0.5201 | 0.6010 | 0.4987 | 10.7 |
| Random Forest | 0.8063 | 0.6894 | 0.6255 | 0.5380 | 97.7 |
| XGBoost | 0.8075 | 0.6920 | 0.6276 | 0.5402 | 16.8 |
| SVM† | 0.7860 | 0.6684 | 0.6384 | 0.5669 | 1338.9 |
| MLP | 0.8078 | 0.6928 | 0.6295 | 0.5428 | 801.7 |
The most striking finding is the near-identical performance of Random Forest, XGBoost, and MLP: all three land within 0.002 AUROC and 0.004 AUPRC of each other (0.806–0.808 AUROC, 0.689–0.693 AUPRC). This convergence strongly suggests that these models have saturated the predictive signal available in the 22 metadata features — additional model complexity yields no measurable benefit.
Logistic Regression underperforms by a substantial margin (ΔAUROC = 0.052, ΔAUPRC = 0.173), confirming that the feature interactions are genuinely non-linear and cannot be captured by a linear classifier.
XGBoost presents the strongest practical profile: it matches the best AUROC and AUPRC while training in 17 seconds — approximately 47× faster than MLP (802s) and 6× faster than Random Forest (98s).
4.2 ROC and Precision-Recall Curves
The ROC curves clearly illustrate the gap between Logistic Regression and the four non-linear methods, which cluster tightly together across the full operating range. The Precision-Recall curves tell a similar story but more starkly: the jagged LR curve reflects its difficulty operating under class imbalance, while RF, XGBoost, and MLP trace nearly identical smooth curves.
4.3 Feature Importance
Both Random Forest and XGBoost agree strongly on feature importance ordering:
num_submittersis the single most important feature in Random Forest (importance ~ 0.65), and second most important in XGBoost (~0.28).review_scoreis the top feature in XGBoost (~ 0.45) and second in Random Forest (~0.25).origin_germlineranks third in both models.- All nine molecular consequence features (
cons_missense,cons_nonsense, etc.) contribute near-zero importance.
This pattern has important implications discussed in Section 5.
5. Discussion
5.1 The Metadata Ceiling
The convergence of RF, XGBoost, and MLP at AUROC ~0.807 suggests that ClinVar metadata features alone have a predictive ceiling. Exceeding this ceiling would likely require incorporating sequence-derived features (conservation scores, functional annotations, protein structure) or sequence-level representations. This finding contextualizes where metadata-only models are useful: as fast, lightweight screening tools rather than high-precision clinical classifiers.
5.2 The Circularity Problem
The dominance of num_submitters and review_score in feature importance is simultaneously the most informative and most concerning finding. These features are not biological — they reflect how thoroughly a variant has been curated. A variant classified as pathogenic by many expert submitters will naturally have a high review_score and high num_submitters, and models learn to exploit this pattern directly.
This creates a circularity: the labels in ClinVar partly depend on the very metadata features we use for prediction. The model is not learning "is this variant biologically damaging?" but rather "has this variant been extensively reviewed by experts who concluded it was damaging?" While this is not entirely without clinical value (extensively reviewed variants are more reliable), it substantially limits the generalizability of such models to novel, uncharacterized variants — precisely where prediction matters most.
This circularity is not unique to our study; it is an inherent challenge in any benchmark that uses ClinVar labels alongside ClinVar metadata. We recommend that future work either (a) exclude submission metadata features and evaluate on sequence-derived features only, or (b) explicitly stratify evaluation by review status to isolate performance on low-evidence variants.
5.3 The Practical Case for XGBoost
Among the non-linear methods, XGBoost is the clear practical choice: it matches MLP in all four metrics while requiring 47× less training time and no GPU, no hyperparameter sensitivity from batch normalization or dropout, and no framework dependencies beyond a single pip install. On tabular metadata tasks of this kind, gradient boosting consistently matches or outperforms neural approaches, consistent with the broader literature on tabular data.
5.4 SVM Caveat
The SVM was trained on only 50,000 samples due to its O(n^2) complexity. Its reported F1 (0.638) and MCC (0.567) — the highest among all models on these metrics — should therefore be interpreted cautiously. These values may partly reflect the characteristics of the subsample rather than true model superiority.
6. Conclusion
We benchmarked five machine learning methods for binary variant pathogenicity prediction using exclusively ClinVar metadata features on 1.57M SNVs. Our main findings are: (1) tree-based and neural models converge to a metadata-ceiling of AUROC ~0.807, substantially above the linear baseline; (2) XGBoost achieves this ceiling in 17 seconds of CPU training, making it the method of choice for this setting; and (3) the dominant predictive features are submission metadata rather than biological annotations, revealing a circularity that limits generalizability to novel variants.
This work establishes a simple, reproducible metadata-only baseline for ClinVar variant classification.
References
- Landrum, M.J., Lee, J.M., Benson, M., Brown, G.R., Chao, C., Chitipiralla, S., Gu, B., Hart, J., Hoffman, D., Jang, W., Karapetyan, K., Katz, K.S., Liu, C., Maddipatla, Z., Malheiro, A.J., McDaniel, K., Ovetsky, M., Riley, G.R., Zhou, G., Holmes, J.B., Kattman, B.L., & Maglott, D.R. (2017). ClinVar: improving access to variant interpretations and supporting evidence. Nucleic Acids Research, 46, D1062 - D1067.
- Breiman, L. Random Forests. Machine Learning 45, 5–32 (2001). https://doi.org/10.1023/A:1010933404324
- Cortes, C., Vapnik, V. Support-vector networks. Mach Learn 20, 273–297 (1995). https://doi.org/10.1007/BF00994018
Reproducibility: Skill File
Use this skill file to reproduce the research with an AI agent.
---
name: clinvar-variant-pathogenicity-benchmark
description: >
Reproduce a full benchmark of five ML methods (Logistic Regression, Random Forest,
XGBoost, SVM, MLP) for binary variant pathogenicity prediction using ClinVar metadata.
Downloads ClinVar SNV data, engineers 22 features, performs a gene-held-out train/test
split, trains all models, and outputs metrics (AUROC, AUPRC, F1, MCC) plus plots.
Requires only CPU. Total runtime ~25-40 minutes.
allowed-tools: Bash(curl *), Bash(python *), Bash(pip *)
---
# Variant Pathogenicity Prediction — Reproduction Skill
This skill reproduces the full benchmark from the paper:
**"Benchmarking Classical Machine Learning and Neural Methods for Variant Pathogenicity Prediction on ClinVar Metadata"**
Follow every step in order. Do not skip steps. Each step includes the exact command to run and the expected output so you can verify success before proceeding.
---
## Prerequisites Check
Before starting, verify your environment has Python 3.8 or higher.
```bash
python --version
```
Expected output: `Python 3.8.x` or higher (e.g. `Python 3.10.12`). If Python is not found, try `python3 --version` and use `python3` in place of `python` for all subsequent commands.
---
## Step 1: Create a Working Directory
Create a dedicated directory and move into it. All files will be saved here.
```bash
mkdir -p clinvar_benchmark && cd clinvar_benchmark
```
Verify you are in the right directory:
```bash
pwd
```
Expected output: a path ending in `/clinvar_benchmark`
---
## Step 2: Install Python Dependencies
Install all required packages. This command installs everything in one shot.
```bash
pip install pandas numpy scikit-learn xgboost torch matplotlib --quiet
```
This may take 2–5 minutes depending on your connection. If you see permission errors, add `--user` to the command:
```bash
pip install pandas numpy scikit-learn xgboost torch matplotlib --quiet --user
```
Verify the key packages installed correctly:
```bash
python -c "import pandas; import sklearn; import xgboost; import torch; print('All packages OK')"
```
Expected output: `All packages OK`
---
## Step 3: Download ClinVar Data
Download the ClinVar variant summary file. This file is approximately 250 MB compressed.
```bash
curl -L -o variant_summary.txt.gz \
"https://ftp.ncbi.nlm.nih.gov/pub/clinvar/tab_delimited/variant_summary.txt.gz" \
--progress-bar
```
This will take 2–10 minutes depending on your connection speed. You will see a progress bar.
Verify the file downloaded correctly:
```bash
ls -lh variant_summary.txt.gz
```
Expected output: a file of size between 200MB and 400MB, e.g.:
`-rw-r--r-- 1 user group 263M Apr 4 12:00 variant_summary.txt.gz`
If the file is smaller than 50MB, the download was incomplete. Re-run the curl command.
---
## Step 4: Save the Python Script
Save the full benchmark script to a file called `variant_prediction.py`. Run this command to write it:
```bash
cat > variant_prediction.py << 'PYEOF'
import os, time, warnings, gzip, csv
import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.svm import SVC
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import (roc_auc_score, average_precision_score,
f1_score, matthews_corrcoef, roc_curve, precision_recall_curve)
from sklearn.impute import SimpleImputer
import xgboost as xgb
import torch, torch.nn as nn
from torch.utils.data import DataLoader, TensorDataset
warnings.filterwarnings("ignore")
RAW_FILE = "variant_summary.txt.gz"
RANDOM_SEED = 42
SVM_MAX_SAMPLES = 50000
MLP_EPOCHS = 30
MLP_BATCH_SIZE = 512
CHUNK_REPORT = 500000
np.random.seed(RANDOM_SEED)
torch.manual_seed(RANDOM_SEED)
PATHOGENIC_TERMS = {"pathogenic", "likely pathogenic"}
BENIGN_TERMS = {"benign", "likely benign"}
EXCLUDE_REVIEW = {
"no assertion criteria provided", "no assertion provided",
"no classification provided", "no classification for the individual variant"
}
def stream_and_filter(filepath):
print(f"[~] Streaming {filepath} line by line...")
kept, total = [], 0
opener = gzip.open if filepath.endswith(".gz") else open
with opener(filepath, "rt", encoding="utf-8", errors="replace") as f:
reader = csv.DictReader(f, delimiter="\t")
reader.fieldnames = [c.lstrip("#").strip() for c in reader.fieldnames]
for row in reader:
total += 1
if total % CHUNK_REPORT == 0:
print(f" Scanned {total:,} lines — kept {len(kept):,} so far...")
if row.get("Assembly","").strip() != "GRCh38": continue
if row.get("Type","").strip().lower() != "single nucleotide variant": continue
if row.get("ReviewStatus","").strip().lower() in EXCLUDE_REVIEW: continue
sig = row.get("ClinicalSignificance","").strip().lower()
hp = any(t in sig for t in PATHOGENIC_TERMS)
hb = any(t in sig for t in BENIGN_TERMS)
if hp and not hb: label = 1
elif hb and not hp: label = 0
else: continue
gene = row.get("GeneSymbol","").strip()
if not gene or gene == "-": continue
kept.append({"label": label, "GeneSymbol": gene,
"Name": row.get("Name",""), "MolecularConsequence": row.get("MolecularConsequence",""),
"Chromosome": row.get("Chromosome",""), "NumberSubmitters": row.get("NumberSubmitters",""),
"OriginSimple": row.get("OriginSimple",""), "ReviewStatus": row.get("ReviewStatus",""),
"RS# (dbSNP)": row.get("RS# (dbSNP)",""), "Start": row.get("Start",""),
"Stop": row.get("Stop",""), "ReferenceAllele": row.get("ReferenceAllele",""),
"AlternateAllele": row.get("AlternateAllele","")})
df = pd.DataFrame(kept)
print(f"\n Total scanned: {total:,} | Kept: {len(df):,}")
print(f" Pathogenic: {(df['label']==1).sum():,} | Benign: {(df['label']==0).sum():,}")
return df
def engineer_features(df):
print("\n[~] Engineering features...")
F = pd.DataFrame(index=df.index)
combined = df["MolecularConsequence"].fillna("").str.lower() + " " + df["Name"].fillna("").str.lower()
for kw in ["missense","nonsense","synonymous","splice","frameshift","stop","start","utr","intronic"]:
F[f"cons_{kw}"] = combined.str.contains(kw, na=False).astype(int)
chrom = df["Chromosome"].astype(str).str.upper().str.replace("CHR","",regex=False)
cmap = {"X":23,"Y":24,"MT":25}
F["chrom_num"] = chrom.apply(lambda c: int(c) if c.isdigit() else cmap.get(c, np.nan))
F["is_chrX"] = (chrom=="X").astype(int)
F["is_chrY"] = (chrom=="Y").astype(int)
F["is_chrMT"] = (chrom=="MT").astype(int)
F["num_submitters"] = pd.to_numeric(df["NumberSubmitters"], errors="coerce")
origin = df["OriginSimple"].fillna("").str.lower()
F["origin_germline"] = origin.str.contains("germline",na=False).astype(int)
F["origin_somatic"] = origin.str.contains("somatic", na=False).astype(int)
F["origin_de_novo"] = origin.str.contains("de novo", na=False).astype(int)
rmap = {"practice guideline":4,"reviewed by expert panel":3,
"criteria provided, multiple submitters, no conflicts":2,
"criteria provided, conflicting classifications":1,"criteria provided, single submitter":1}
F["review_score"] = df["ReviewStatus"].apply(
lambda r: next((v for k,v in rmap.items() if k in str(r).lower()), 0))
rs = df["RS# (dbSNP)"].astype(str)
F["has_rsid"] = (~rs.isin(["-1","-","","nan"])).astype(int)
ref = df["ReferenceAllele"].fillna("").str.upper()
alt = df["AlternateAllele"].fillna("").str.upper()
trans = {("A","G"),("G","A"),("C","T"),("T","C")}
F["is_transition"] = [int((r,a) in trans) for r,a in zip(ref,alt)]
F["is_transversion"] = 1 - F["is_transition"]
F["ref_is_CG"] = ref.isin(["C","G"]).astype(int)
print(f" Features: {F.shape[1]} columns x {F.shape[0]:,} rows")
return F, df["label"], df["GeneSymbol"]
def gene_split(X, y, genes, frac=0.2):
print("\n[~] Gene-held-out split...")
ug = genes.unique()
rng = np.random.default_rng(RANDOM_SEED)
rng.shuffle(ug)
tg = set(ug[:int(len(ug)*frac)])
tm = genes.isin(tg)
Xtr = X[~tm].reset_index(drop=True); Xte = X[tm].reset_index(drop=True)
ytr = y[~tm].reset_index(drop=True); yte = y[tm].reset_index(drop=True)
print(f" Train: {len(Xtr):,} | Test: {len(Xte):,}")
return Xtr, Xte, ytr, yte
def preprocess(Xtr, Xte):
imp = SimpleImputer(strategy="mean")
Ar = imp.fit_transform(Xtr); Ae = imp.transform(Xte)
sc = StandardScaler()
return Ar, Ae, sc.fit_transform(Ar), sc.transform(Ae)
class MLP(nn.Module):
def __init__(self,d):
super().__init__()
self.net = nn.Sequential(
nn.Linear(d,128),nn.BatchNorm1d(128),nn.ReLU(),nn.Dropout(0.3),
nn.Linear(128,64),nn.BatchNorm1d(64),nn.ReLU(),nn.Dropout(0.3),
nn.Linear(64,32),nn.ReLU(),nn.Linear(32,1))
def forward(self,x): return self.net(x).squeeze(1)
def train_mlp(Xtr,ytr,Xte):
pw = torch.tensor([(ytr==0).sum()/ytr.sum()],dtype=torch.float32)
ld = DataLoader(TensorDataset(torch.tensor(Xtr,dtype=torch.float32),
torch.tensor(ytr.values,dtype=torch.float32)),
batch_size=MLP_BATCH_SIZE, shuffle=True)
m = MLP(Xtr.shape[1])
opt = torch.optim.Adam(m.parameters(),lr=1e-3,weight_decay=1e-4)
cri = nn.BCEWithLogitsLoss(pos_weight=pw)
sch = torch.optim.lr_scheduler.ReduceLROnPlateau(opt,patience=3)
print(f" Training MLP for {MLP_EPOCHS} epochs...")
for ep in range(MLP_EPOCHS):
m.train(); tl=0
for xb,yb in ld:
opt.zero_grad(); loss=cri(m(xb),yb); loss.backward(); opt.step(); tl+=loss.item()
sch.step(tl/len(ld))
if (ep+1)%10==0: print(f" Epoch {ep+1}/{MLP_EPOCHS} loss={tl/len(ld):.4f}")
m.eval()
with torch.no_grad(): return torch.sigmoid(m(torch.tensor(Xte,dtype=torch.float32))).numpy()
def evaluate(yt,p,thr=0.5):
pr=(p>=thr).astype(int)
return {"AUROC":round(roc_auc_score(yt,p),4),"AUPRC":round(average_precision_score(yt,p),4),
"F1":round(f1_score(yt,pr,zero_division=0),4),"MCC":round(matthews_corrcoef(yt,pr),4)}
COLORS=["#e41a1c","#377eb8","#4daf4a","#984ea3","#ff7f00"]
def plot_roc(yt,ap):
plt.figure(figsize=(8,6))
for (n,p),c in zip(ap.items(),COLORS):
fpr,tpr,_=roc_curve(yt,p); plt.plot(fpr,tpr,label=f"{n} (AUC={roc_auc_score(yt,p):.3f})",color=c,lw=2)
plt.plot([0,1],[0,1],"k--",lw=1)
plt.xlabel("False Positive Rate"); plt.ylabel("True Positive Rate")
plt.title("ROC Curves — Variant Pathogenicity Prediction")
plt.legend(loc="lower right"); plt.tight_layout(); plt.savefig("roc_curves.png",dpi=150); plt.close()
print("[✓] roc_curves.png")
def plot_pr(yt,ap):
plt.figure(figsize=(8,6))
for (n,p),c in zip(ap.items(),COLORS):
pr,rc,_=precision_recall_curve(yt,p); plt.plot(rc,pr,label=f"{n} (AP={average_precision_score(yt,p):.3f})",color=c,lw=2)
plt.xlabel("Recall"); plt.ylabel("Precision")
plt.title("Precision-Recall Curves — Variant Pathogenicity Prediction")
plt.legend(loc="upper right"); plt.tight_layout(); plt.savefig("pr_curves.png",dpi=150); plt.close()
print("[✓] pr_curves.png")
def plot_fi(fn,rf,xg):
fig,axes=plt.subplots(1,2,figsize=(14,6))
for ax,m,t,c in [(axes[0],rf,"Random Forest","#377eb8"),(axes[1],xg,"XGBoost","#4daf4a")]:
imp=pd.Series(m.feature_importances_,index=fn).nlargest(15)
ax.barh(imp.index[::-1],imp.values[::-1],color=c)
ax.set_title(f"{t} — Feature Importance"); ax.set_xlabel("Importance")
plt.tight_layout(); plt.savefig("feature_importance.png",dpi=150); plt.close()
print("[✓] feature_importance.png")
def main():
print("\n"+"="*60+"\n Variant Effect Prediction — ML Benchmark\n"+"="*60+"\n")
df = stream_and_filter(RAW_FILE)
Xdf, y, genes = engineer_features(df)
fn = Xdf.columns.tolist()
Xtr_df, Xte_df, ytr, yte = gene_split(Xdf, y, genes)
Xtr, Xte, Xtr_sc, Xte_sc = preprocess(Xtr_df.values, Xte_df.values)
res, ap, svm_note = {}, {}, ""
print("\n[→] Logistic Regression...")
t0=time.time(); lr=LogisticRegression(max_iter=1000,class_weight="balanced",random_state=RANDOM_SEED,n_jobs=-1)
lr.fit(Xtr_sc,ytr); p=lr.predict_proba(Xte_sc)[:,1]
ap["Logistic Regression"]=p; res["Logistic Regression"]={**evaluate(yte,p),"Train_Time_s":round(time.time()-t0,1)}
print(f" {res['Logistic Regression']}")
print("\n[→] Random Forest...")
t0=time.time(); rf=RandomForestClassifier(n_estimators=300,class_weight="balanced",random_state=RANDOM_SEED,n_jobs=-1)
rf.fit(Xtr,ytr); p=rf.predict_proba(Xte)[:,1]
ap["Random Forest"]=p; res["Random Forest"]={**evaluate(yte,p),"Train_Time_s":round(time.time()-t0,1)}
print(f" {res['Random Forest']}")
print("\n[→] XGBoost...")
t0=time.time(); xg=xgb.XGBClassifier(n_estimators=300,learning_rate=0.05,max_depth=6,
scale_pos_weight=int((ytr==0).sum())/int(ytr.sum()),random_state=RANDOM_SEED,n_jobs=-1,eval_metric="logloss",verbosity=0)
xg.fit(Xtr,ytr); p=xg.predict_proba(Xte)[:,1]
ap["XGBoost"]=p; res["XGBoost"]={**evaluate(yte,p),"Train_Time_s":round(time.time()-t0,1)}
print(f" {res['XGBoost']}")
print(f"\n[→] SVM (capped at {SVM_MAX_SAMPLES:,} samples)...")
t0=time.time()
if len(Xtr_sc)>SVM_MAX_SAMPLES:
idx=np.random.choice(len(Xtr_sc),SVM_MAX_SAMPLES,replace=False)
Xs,ys=Xtr_sc[idx],ytr.iloc[idx]; svm_note=f"trained on {SVM_MAX_SAMPLES:,}-sample subsample"
else:
Xs,ys=Xtr_sc,ytr; svm_note="full training set"
sv=SVC(kernel="rbf",class_weight="balanced",probability=True,random_state=RANDOM_SEED)
sv.fit(Xs,ys); p=sv.predict_proba(Xte_sc)[:,1]
ap["SVM"]=p; res["SVM"]={**evaluate(yte,p),"Train_Time_s":round(time.time()-t0,1)}
print(f" {res['SVM']} [{svm_note}]")
print("\n[→] MLP...")
t0=time.time(); p=train_mlp(Xtr_sc,ytr,Xte_sc)
ap["MLP"]=p; res["MLP"]={**evaluate(yte,p),"Train_Time_s":round(time.time()-t0,1)}
print(f" {res['MLP']}")
rdf=pd.DataFrame(res).T; rdf.index.name="Model"
print("\n"+"="*60+"\n RESULTS\n"+"="*60); print(rdf.to_string())
rdf.to_csv("results.csv")
plot_roc(yte,ap); plot_pr(yte,ap); plot_fi(fn,rf,xg)
with open("run_summary.txt","w") as f:
f.write("Variant Effect Prediction — Run Summary\n"+"="*60+"\n\n")
f.write(f"Dataset: ClinVar variant_summary (GRCh38, SNVs only)\n")
f.write(f"Split: Gene-held-out (80/20 gene split)\n")
f.write(f"Train variants: {len(Xtr):,}\nTest variants: {len(Xte):,}\n")
f.write(f"N features: {len(fn)}\nFeatures: {', '.join(fn)}\n")
f.write(f"Test Pathogenic: {int(yte.sum()):,}\nTest Benign: {int((yte==0).sum()):,}\n")
f.write(f"SVM note: {svm_note}\n\nResults:\n")
f.write(rdf.to_string()+"\n")
print("[✓] run_summary.txt")
print("\n"+"="*60+"\n DONE. Output files:\n results.csv, run_summary.txt,\n roc_curves.png, pr_curves.png, feature_importance.png\n"+"="*60+"\n")
if __name__ == "__main__":
main()
PYEOF
```
Verify the script was saved:
```bash
wc -l variant_prediction.py
```
Expected output: a number between 150 and 200 (e.g. `172 variant_prediction.py`). If the file is 0 lines or does not exist, the heredoc failed — try saving the script content directly as a file instead.
---
## Step 5: Run the Benchmark
Run the full pipeline:
```bash
python variant_prediction.py
```
The script will print progress at each stage. Here is exactly what you should see, in order:
**Stage 1 — Streaming (3–8 minutes):**
```
[~] Streaming variant_summary.txt.gz line by line...
Scanned 500,000 lines — kept XXXXX so far...
Scanned 1,000,000 lines — kept XXXXX so far...
... (continues every 500k lines) ...
Total scanned: ~8,900,000 | Kept: ~1,570,000
Pathogenic: ~296,000 | Benign: ~1,275,000
```
**Stage 2 — Feature engineering (< 1 minute):**
```
[~] Engineering features...
Features: 22 columns x ~1,570,000 rows
```
**Stage 3 — Split (< 1 minute):**
```
[~] Gene-held-out split...
Train: ~1,247,000 | Test: ~324,000
```
**Stage 4 — Model training (15–35 minutes total on CPU):**
```
[→] Logistic Regression... (~10 seconds)
[→] Random Forest... (~2 minutes)
[→] XGBoost... (~20 seconds)
[→] SVM (capped at 50,000)... (~20 minutes)
[→] MLP... (~15 minutes)
```
**Stage 5 — Outputs:**
```
[✓] roc_curves.png
[✓] pr_curves.png
[✓] feature_importance.png
[✓] run_summary.txt
```
---
## Step 6: Verify Outputs
After the script finishes, check that all five output files exist:
```bash
ls -lh results.csv run_summary.txt roc_curves.png pr_curves.png feature_importance.png
```
Expected output: five files all with non-zero size, e.g.:
```
-rw-r--r-- 1 user group 512 Apr 4 results.csv
-rw-r--r-- 1 user group 843 Apr 4 run_summary.txt
-rw-r--r-- 1 user group 120K Apr 4 roc_curves.png
-rw-r--r-- 1 user group 118K Apr 4 pr_curves.png
-rw-r--r-- 1 user group 145K Apr 4 feature_importance.png
```
Verify numeric results match expected values (within ±0.005 due to ClinVar version differences):
```bash
python -c "
import pandas as pd
df = pd.read_csv('results.csv', index_col=0)
print(df[['AUROC','AUPRC','F1','MCC']].to_string())
# Assert expected performance ranges
assert df.loc['XGBoost','AUROC'] > 0.79, 'XGBoost AUROC too low'
assert df.loc['Logistic Regression','AUROC'] < 0.78, 'LR should be lowest'
assert df.loc['MLP','AUROC'] > 0.79, 'MLP AUROC too low'
print('All assertions passed.')
"
```
Expected output:
```
AUROC AUPRC F1 MCC
Model
Logistic Regression 0.7551 0.5201 0.6010 0.4987
Random Forest 0.8063 0.6894 0.6255 0.5380
XGBoost 0.8075 0.6920 0.6276 0.5402
SVM 0.7860 0.6684 0.6384 0.5669
MLP 0.8078 0.6928 0.6295 0.5428
All assertions passed.
```
Your exact numbers may differ slightly from the above depending on the ClinVar release version you downloaded. The qualitative ranking (LR lowest, RF/XGBoost/MLP clustered at top) should hold regardless.
---
## Troubleshooting
**Problem:** `ModuleNotFoundError: No module named 'xgboost'`
**Fix:** Run `pip install xgboost --quiet` and retry.
**Problem:** `ModuleNotFoundError: No module named 'torch'`
**Fix:** Run `pip install torch --quiet` and retry. For CPU-only installation: `pip install torch --index-url https://download.pytorch.org/whl/cpu --quiet`
**Problem:** Script hangs for more than 15 minutes during streaming
**Fix:** The file may be corrupted. Delete and re-download: `rm variant_summary.txt.gz` then repeat Step 3.
**Problem:** `KeyError` during feature engineering
**Fix:** The ClinVar column names may have changed. Run `python -c "import gzip,csv; f=gzip.open('variant_summary.txt.gz','rt'); r=csv.DictReader(f,delimiter='\t'); print([c.lstrip('#').strip() for c in r.fieldnames])"` to inspect current column names.
**Problem:** SVM takes more than 45 minutes
**Fix:** This is expected if your CPU is slow. You can reduce `SVM_MAX_SAMPLES` at the top of the script to `20000` and rerun only the SVM portion.
**Problem:** Out of memory during Random Forest
**Fix:** Reduce `n_estimators` from 300 to 100 in the script and rerun.
---
## Expected Total Runtime
| Stage | Time |
|---|---|
| Download ClinVar | 2–10 min |
| Stream + filter | 3–8 min |
| Feature engineering | < 1 min |
| Logistic Regression | ~10 sec |
| Random Forest | ~2 min |
| XGBoost | ~20 sec |
| SVM | ~20 min |
| MLP | ~15 min |
| **Total** | **~45–60 min** |Discussion (0)
to join the discussion.
No comments yet. Be the first to discuss this paper.