Catalog variants tutorial (variant_aware.py — ClinVar/TCGA/1000G batches)

Catalog variants tutorial (variant_aware.py — ClinVar/TCGA/1000G batches)#

This notebook is a tutorialized, step-by-step walkthrough of the unified variant_aware.py entry point. It decomposes the run into annotated, inspectable steps (explicit I/O expectations + lightweight sanity checks such as printing tensor shapes), while keeping the underlying execution logic identical to the script (we call the same functions and helpers).

CLI arguments#

Common arguments#

Argument	Type / choices	Required	Meaning
`--variation_mode`	`before` / `after`	yes	Score reference window (“before”) or ALT-substituted window (“after”).
`--fasta_sequence_path`	`PATH`	yes	Input FASTA containing window sequences.
`--variant_out_file`	`PATH`	yes	Output file (appended).
`--Transformer_path`	`PATH`	yes	Transformer path used by `build_Transformer_embeddings`.
`--model_save_path`	`PATH`	yes	Directory with BRIDGE checkpoints (`*.pth`).
`--device`	`cuda`, `cuda:N`, or `cpu`	no	Torch device (default: cuda if available else cpu).

Pipeline selection flags#

Argument	Type / choices	Required	Meaning
`--gwas`	`-`	no	Force GWAS pipeline (default if no pipeline flag provided).
`--ribosnitch/--ribosnitches`	`-`	no	Enable ribosnitch pipeline.
`--catalog_variants/--genomic_variants`	`-`	no	Enable catalog-variants pipeline (ClinVar/TCGA/1000G).

Catalog-variants-only arguments#

Argument	Type / choices	Required	Meaning
`--model_id_strategy`	`from_header` / `from_fasta_stem`	no	Catalog mode: choose checkpoint naming strategy.
`--k`	`int`	no	Catalog mode: forwarded to build_Transformer_embeddings.
`--pos_weight`	`float`	no	Catalog mode: pos_weight for BCEWithLogitsLoss.
`--strict_ref_match`	`-`	no	Catalog mode: skip if REF/ALT cannot be matched in the window.
`--disable_off_by_one`	`-`	no	Catalog mode: disable +/-1 fallback when locating SNV index.

Quickstart (the two most common commands)#

Below is an example using a 1000G FASTA file (typically you run the same FASTA twice: before and after):

python variant_aware.py \
  --catalog_variants \
  --variation_mode before \
  --fasta_sequence_path /path/to/1000genomes_diff.all.fa \
  --variant_out_file ./results/variant/mut_before_after_score/1000genomes.before.txt \
  --Transformer_path /path/to/RBPformer \
  --model_save_path ./results/model \
  --device cpu

python variant_aware.py \
  --catalog_variants \
  --variation_mode after \
  --fasta_sequence_path /path/to/1000genomes_diff.all.fa \
  --variant_out_file ./results/variant/mut_before_after_score/1000genomes.after.txt \
  --Transformer_path /path/to/RBPformer \
  --model_save_path ./results/model \
  --device cuda:0

What you should see in the terminal logs typically includes:

The number of FASTA records and the window length (if you enable the sanity-check cell below).
Per-record header parsing (errors/warnings will be printed if a header cannot be parsed).
A final message like “Results appended to …”.

Input FASTA requirements (Catalog variants)#

Header must include chr:start-end(strand) and POS:REF>ALT.
Default model id parsing expects ... <PROTEIN> in <CELL>; otherwise set --model_id_strategy from_fasta_stem.

Example input (recommended: copy and edit)#

A typical header should contain at least two required tokens:

Region token: chr:start-end(strand), e.g. chr1:100-200(+)
Variant token: POS:REF>ALT, e.g. 150:A>G

Optionally (but commonly), include protein and cell line information so the default checkpoint id can be inferred as model_id = <PROTEIN>_<CELL>:

>chr1:100-200(+) 150:A>G SRSF1 in K562
ACGTT... (a fixed-length window sequence, e.g. 101 nt)

What should the sequence length be?#

Many BRIDGE checkpoints (and the feature tensors created by this script) are designed around a 101 nt window (you will see 101 in the code).
If your trained/downloaded checkpoint uses a different window length, make sure your input sequence length matches the checkpoint; otherwise the embedding/feature dimensions may not align.

What happens when `strand` is `-`?#

When strand == '-', the script complements REF/ALT before locating/replacing the variant so that the sequence orientation matches what the model expects (consistent with how the model was trained).

Output format (Catalog variants)#

<header_without_>\tmodel_id=<...>\tmode=<before|after>\tPrediction_score:<float>

Example output (single line)#

chr1:100-200(+) 150:A>G SRSF1 in K562	model_id=SRSF1_K562	mode=after	Prediction_score:0.123456

mode=before: score the as-is window sequence (if the input already contains ALT at the variant position, the script will log it and still score the sequence as provided).
mode=after: locate the variant position in the window and replace it with ALT before scoring.
Prediction_score: a floating-point model score (in the current implementation, this is the logits/score returned by validate_without_sigmoid).

0. Imports + repo bootstrap#

from __future__ import annotations

import sys
from pathlib import Path

def find_repo_root(start: Path | None = None) -> Path:
    # Heuristically locate the BRIDGE repo root so we can import `variant_aware.py` and `utils/`.
    # Works when this notebook lives under docs/tutorials/notebooks/.
    start = (start or Path().resolve())
    for p in [start, *start.parents]:
        if (p / "variant_aware.py").exists() and (p / "utils").exists():
            return p
    raise FileNotFoundError(
        "Cannot locate repo root (expected to find variant_aware.py and utils/). "
        "Run this notebook from within the BRIDGE repository."
    )

REPO_ROOT = find_repo_root()
if str(REPO_ROOT) not in sys.path:
    sys.path.insert(0, str(REPO_ROOT))
print("Repo root:", REPO_ROOT)

import variant_aware as va
from utils import variant as var_utils

print("Imported:", va.__file__)

Repo root: /fs1/private/user/wangyubo/code/BRIDGE
Imported: /fs1/private/user/wangyubo/code/BRIDGE/variant_aware.py

1. Minimal runnable demo data (catalog-variant header)#

This branch is for ClinVar / TCGA / 1000G-style FASTA batches. Headers include a region token (chr:start-end(strand)) and an SNV token (POS:REF>ALT).

from pathlib import Path

def find_repo_root(start: Path) -> Path:
    for p in [start] + list(start.parents):
        if (p / ".git").exists():
            return p
    return start

repo_root = find_repo_root(Path.cwd())
fasta_path = repo_root / "dataset_variant" / "1000genomes_diff.all.fa"

n_records = 3
count = 0

with fasta_path.open("r") as f:
    while count < n_records:
        header = f.readline()
        if not header:
            break
        seq = f.readline()
        if not seq:
            break
        header = header.strip()
        seq = seq.strip()
        print(f"[{count+1}] {header}")
        print(f"    len={len(seq)} seq[:60]={seq[:60]}{'...' if len(seq)>60 else ''}")
        count += 1

[1] >variant_10 chrM:542-642(+)[]{Benign} 592:C>T gene MTPAP in K562
    len=101 seq[:60]=CCCGAACCAACCAAACCCCAAAGACACCCCCCACAGTTTATGTAGCTTACCTCCTCAAAG...
[2] >variant_11 chrM:542-642(+)[]{Benign} 592:C>T gene ILF3 in K562
    len=101 seq[:60]=CCCGAACCAACCAAACCCCAAAGACACCCCCCACAGTTTATGTAGCTTACCTCCTCAAAG...
[3] >variant_12 chrM:542-642(+)[]{Benign} 592:C>T gene RPS3 in K562
    len=101 seq[:60]=CCCGAACCAACCAAACCCCAAAGACACCCCCCACAGTTTATGTAGCTTACCTCCTCAAAG...

2. Load FASTA#

headers, seqs = va.read_fasta(fasta_path)
header0, seq0 = headers[0], seqs[0]
print("header0:", header0)
print("seq0 len:", len(seq0))

header0: >variant_10 chrM:542-642(+)[]{Benign} 592:C>T gene MTPAP in K562
seq0 len: 101

3. Parse header into structured fields (`parse_header_catalog`)#

ph = va.parse_header_catalog(header0)
print(ph)
print("protein/cell:", ph.protein, ph.cell_line)

ParsedHeader(header_raw='>variant_10 chrM:542-642(+)[]{Benign} 592:C>T gene MTPAP in K562', chrom='chrM', start=542, end=642, strand='+', var_pos=592, ref='C', alt='T', protein='MTPAP', cell_line='K562')
protein/cell: MTPAP K562

4. Find the variant index inside the window + demonstrate off-by-one fallback#

# Strand-aware handling (same as in the pipeline)
ref, alt = ph.ref, ph.alt
if ph.strand == "-":
    ref, alt = va.apply_complement(ref), va.apply_complement(alt)

idx0, state = va.find_variant_index(
    seq=seq0,
    seq_start=ph.start,
    var_pos=ph.var_pos,
    ref=ref,
    alt=alt,
    try_off_by_one=True,
)
print("idx0:", idx0, "| state:", state)

L = 8
ctx = seq0[max(0, idx0-L): idx0] + "[" + seq0[idx0] + "]" + seq0[idx0+1: idx0+1+L]
print("Context:", ctx)

idx0: 50 | state: ref
Context: TAGCTTAC[C]TCCTCAAA

5. Apply mutation according to `variation_mode`#

variation_mode = "after"  # try "before" vs "after"

if variation_mode == "before":
    modified_seq = seq0
else:
    modified_seq = va.substitute_base(seq0, idx0, alt) if state == "ref" else seq0

ctx_after = modified_seq[max(0, idx0-L): idx0] + "[" + modified_seq[idx0] + "]" + modified_seq[idx0+1: idx0+1+L]
print("mode:", variation_mode)
print("after:", ctx_after)

mode: after
after: TAGCTTAC[T]TCCTCAAA

6. Run full catalog-variants scoring function#

This calls process_sequences_catalog_variants(...) directly (same logic as CLI). Requires real checkpoints; model-id strategy is explained in the CLI help.

If you already have a working BRIDGE environment (dependencies installed and valid Transformer/checkpoint paths), you can run an end-to-end example as follows:

Make sure these three paths exist:
- FASTA_PATH: your FASTA file (ClinVar/TCGA/1000G catalog variants)
- TRANSFORMER_PATH: the Transformer directory (e.g. RBPformer)
- MODEL_SAVE_PATH: the BRIDGE checkpoint directory (containing <model_id>.pth or equivalent naming)
Start with a small file and run --device cpu as a smoke test, then switch to cuda:0 for large-scale scoring.

The output file is appended to. If a file with the same name already exists, new results will be appended at the end.

import argparse
from pathlib import Path
import torch
import os

Transformer_path = Path(os.environ.get("BRIDGE_TRANSFORMER_PATH", "../../../RBPformer"))
model_save_path = Path(os.environ.get("BRIDGE_MODEL_SAVE_PATH", "../../../results/model"))
variant_out_file = Path("catalog_scores.txt")

args = argparse.Namespace(
    variation_mode="after",
    fasta_sequence_path=fasta_path,
    variant_out_file=variant_out_file,
    Transformer_path=Transformer_path,
    model_save_path=model_save_path,
    catalog_variants=True,
    model_id_strategy="from_header",
    k=1,
    pos_weight=2.0,
    strict_ref_match=False,
    disable_off_by_one=False,
)

device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
hub = va.ModelHub(Transformer_path, device)

if Transformer_path.exists() and model_save_path.exists():
    va.process_sequences_catalog_variants(headers, seqs, args, hub)
    print("Wrote scores to:", variant_out_file)
    print(variant_out_file.read_text().splitlines()[:3])
else:
    print("Skip: set BRIDGE_TRANSFORMER_PATH and BRIDGE_MODEL_SAVE_PATH.")

Wrote scores to: catalog_scores.txt
['>variant_10 chrM:542-642(+)[]{Benign} 592:C>T gene MTPAP in K562\tmodel_id=MTPAP_K562\tmode=after\tPrediction_score:-13.840032', '>variant_11 chrM:542-642(+)[]{Benign} 592:C>T gene ILF3 in K562\tmodel_id=ILF3_K562\tmode=after\tPrediction_score:3.613627', '>variant_12 chrM:542-642(+)[]{Benign} 592:C>T gene RPS3 in K562\tmodel_id=RPS3_K562\tmode=after\tPrediction_score:-20.866917']

End-to-end run (optional): load the model and write scores#