Precomputed GPN-Star scores#
This workflow annotates the real songlab/omim_traitgym Mendelian regulatory-variant benchmark from the GPN-Star collection with the public songlab/gpn-star-scores Parquet shards, then computes pooled global area under the precision-recall curve (AUPRC). It does not download a model or whole-genome alignment.
The join keys are chrom, one-based pos, ref, and alt. They must use the same reference assembly and chromosome spelling as the selected score set.
TraitGym benchmark#
This matched benchmark contains 338 OMIM regulatory variants and 3,042 controls. We pin the exact public dataset revision and evaluate all chromosomes together.
import polars as pl
from sklearn.metrics import average_precision_score
KEYS = ["chrom", "pos", "ref", "alt"]
TRAITGYM_REVISION = "9317562efb8c61f31bb5fc62a19f731b2f8b4384"
TRAITGYM_ROOT = (
"hf://datasets/songlab/omim_traitgym@"
f"{TRAITGYM_REVISION}"
)
variants = pl.read_parquet(f"{TRAITGYM_ROOT}/test.parquet")
assert variants.height == 3_380
assert variants.get_column("label").sum() == 338
variants.select(KEYS + ["label", "match_group"]).head()
| chrom | pos | ref | alt | label | match_group |
|---|---|---|---|---|---|
| str | i64 | str | str | bool | str |
| "1" | 1425822 | "C" | "G" | false | "PLS_4" |
| "1" | 1615869 | "C" | "T" | false | "PLS_0" |
| "1" | 1659060 | "G" | "A" | false | "PLS_4" |
| "1" | 1659114 | "A" | "G" | false | "PLS_5" |
| "1" | 2050958 | "T" | "C" | false | "5_prime_UTR_variant_7" |
Query the required positions remotely#
The released tables are partitioned by chromosome. Each lazy scan filters positions before collection, so the remote query reads only the needed columns and Parquet row groups without downloading chromosome shards. This 19-chromosome benchmark can still take several minutes because its variants are dispersed across the genome.
DATASET_REVISION = "5c799b2ec6aa089f0caa8294ae72adb4510f81ae"
SCORE_SET = "gpn-star-hg38-m447-200m"
SCORE_ROOT = (
"hf://datasets/songlab/gpn-star-scores@"
f"{DATASET_REVISION}/data/{SCORE_SET}/llr"
)
results = []
for chrom in variants.get_column("chrom").unique(maintain_order=True):
chrom_variants = variants.filter(pl.col("chrom") == chrom)
positions = chrom_variants.get_column("pos").unique().to_list()
scores = (
pl.scan_parquet(f"{SCORE_ROOT}/llr_chr{chrom}.parquet")
.filter(pl.col("pos").is_in(positions))
)
results.append(
chrom_variants.lazy()
.join(scores, on=KEYS, how="left")
.collect(engine="streaming")
)
annotated = pl.concat(results)
missing = annotated.get_column("llr_calibrated").null_count()
assert missing == 0, f"{missing} variants were absent from the score table"
annotated
| chrom | pos | ref | alt | OMIM | consequence | label | tss_dist | match_group | llr_calibrated | abs_llr_calibrated |
|---|---|---|---|---|---|---|---|---|---|---|
| str | i64 | str | str | str | str | bool | i64 | str | f32 | f32 |
| "1" | 1425822 | "C" | "G" | null | "PLS" | false | 48 | "PLS_4" | 4.75 | 0.213 |
| "1" | 1615869 | "C" | "T" | null | "PLS" | false | 35 | "PLS_0" | -1.236 | 1.236 |
| "1" | 1659114 | "A" | "G" | null | "PLS" | false | 101 | "PLS_5" | 5.75 | 1.559 |
| "1" | 1659060 | "G" | "A" | null | "PLS" | false | 47 | "PLS_4" | 0.626 | -0.626 |
| "1" | 2074688 | "G" | "A" | null | "5_prime_UTR_variant" | false | 0 | "5_prime_UTR_variant_2" | -1.869 | 1.869 |
| … | … | … | … | … | … | … | … | … | … | … |
| "X" | 155613005 | "C" | "T" | null | "PLS" | false | 52 | "PLS_52" | -0.78 | 0.78 |
| "X" | 155719093 | "C" | "A" | null | "5_prime_UTR_variant" | false | 4 | "5_prime_UTR_variant_101" | -0.644 | 0.644 |
| "X" | 155881342 | "A" | "C" | null | "PLS" | false | 2 | "PLS_57" | -1.853 | 1.853 |
| "X" | 155881414 | "C" | "T" | null | "5_prime_UTR_variant" | false | 35 | "5_prime_UTR_variant_110" | 0.198 | -0.198 |
| "X" | 156020793 | "A" | "G" | null | "upstream_gene_variant" | false | 167 | "upstream_gene_variant_9" | 0.523 | -0.523 |
llr_calibrated is the mutation-rate-calibrated alternate-versus-reference log-likelihood ratio; more-negative values indicate greater constraint or predicted effect. Therefore the classifier ranking score is -llr_calibrated. abs_llr_calibrated is an independently calibrated score and must not be recomputed as abs(llr_calibrated).
annotated = annotated.with_columns(
(-pl.col("llr_calibrated")).alias("effect_score")
)
auprc = average_precision_score(
annotated.get_column("label").to_numpy(),
annotated.get_column("effect_score").to_numpy(),
)
print(f"Genome-wide join global AUPRC: {auprc:.4f}")
Genome-wide join global AUPRC: 0.7644
annotated.write_parquet("gpn_star_scored_variants.parquet")