utils.variant

utils.variant#

Variant-aware inference utilities for BRIDGE (GWAS / ribosnitches-style workflows).

This module implements I/O and parsing helpers used by BRIDGE variant scoring pipelines. It focuses on FASTA inputs whose headers encode variant coordinates and alleles, and provides utilities to (1) reconstruct the alternate-allele sequence window, and (2) cache heavy Transformer/BRIDGE models for high-throughput scoring.

Key ideas#

  1. FASTA parsing with wrapped sequences read_fasta() supports standard multi-line/wrapped FASTA sequences. Each record begins with a header line starting with ‘>’ and is followed by one or more sequence lines. Sequence lines are concatenated and returned in upper-case.

  2. Variant metadata encoded in headers

    Two parsers are provided:

    • parse_variant_block() (legacy / fixed token positions) Assumes the region token is at fields[1] and the variant token is at fields[-2]. This matches the original GWAS implementation and is kept for backward compatibility.

    • parse_variant_block_flexible() (robust / token search) Searches the header tokens for:

      • a region token that contains :, -, (, and )

      • a variant token matching ^\d+:[ACGT]>[ACGT]$ (case-insensitive)

      This is intended for headers where trailing tokens vary (for example, extra annotations or cell-line suffixes), and is often used by ribosnitches-style inputs.

    Both parsers return the same 5-tuple:

    (variant_pos, ref_base, alt_base, strand, seq_start)

    where:

    • variant_pos: genomic coordinate of the SNV

    • ref_base: reference allele base (A/C/G/T)

    • alt_base: alternate allele base (A/C/G/T)

    • strand: '+' or '-' parsed from the region token

    • seq_start: genomic start of the provided sequence window

  3. Coordinate conversion: genomic -> window index

    Given:

    idx0 = variant_pos - seq_start

    the index is 0-based into the sequence window returned by read_fasta().

    Most pipelines then validate that:

    seq[idx0] == ref_base

    (or its complement for the '-' strand) before substituting the alternate base.

  4. Strand handling and complements

    This module provides:

    • COMPLEMENT mapping for DNA letters {A, T, C, G, N}

    • apply_complement() to map one base to its Watson-Crick complement

    • substitute_base() to write an alternate allele at a 0-based window index

    Important:

    • If your sequence window is given on the '-' strand, typical pipelines usually do one of the following:

      • store the window already reverse-complemented (then ref_base / alt_base can be used directly), or

      • store the window in genomic '+' orientation (then allele complementation may be required)

    • This module only provides the complement primitive. The exact policy should be enforced by the caller (that is, the caller decides whether to complement ref_base / alt_base when strand == '-').

    • The complement mapping uses 'T' (DNA). If your windows are RNA ('U'), consider extending COMPLEMENT with {"U": "A"} and adjusting parsing/validation accordingly.

  5. High-throughput model reuse via ModelHub

    BRIDGE variant scoring typically requires:

    • a tokenizer + Transformer encoder (BERT-like) for k-mer embeddings

    • a BRIDGE checkpoint per experiment/model name

    ModelHub caches these heavy components:

    • loads tokenizer and Transformer once from transformer_path

    • caches BRIDGE checkpoints by filename_stem to avoid repeated disk I/O when a FASTA file contains many records spanning multiple models

Constants#

COMPLEMENT

DNA Watson–Crick complement mapping used by apply_complement().

RIBOSNITCHES_MAX_LEN

Default fixed window length (101) used by ribosnitches-derived pipelines. Many downstream feature builders assume length 101; if you deviate, ensure you also update shape-dependent modules.

I/O helpers#

read_fasta(fasta_path)

Read headers and sequences from a FASTA file. Supports wrapped sequences.

open_output(out_path)

Create parent directories and return a Path suitable for writing/appending.

Variant utilities#

parse_variant_block(fasta_header)

Fixed-position parser (legacy GWAS rule).

parse_variant_block_flexible(fasta_header)

Search-based parser (robust to header token drift).

apply_complement(base)

Complement a single base (A/T/C/G), returning unchanged for unknown letters.

substitute_base(seq, pos0, alt)

Replace the base at 0-based index pos0 with alt and return the new string.

ModelHub#

ModelHub(transformer_path, device)

Loads tokenizer/Transformer once, and caches BRIDGE checkpoints.

ModelHub.load_bridge(model_dir, filename_stem)

Load <model_dir>/<filename_stem>.pth into a BRIDGE model on the hub device. Returns None if the checkpoint is missing.

Common failure modes and recommendations#

  • Header format drift:

    If parse_variant_block() raises ValueError or yields wrong tokens, switch to parse_variant_block_flexible() or call it as a fallback.

  • Variant token alphabet:

    The default regex accepts only A/C/G/T. If your headers can contain ‘U’ (e.g., A>U), extend _VARIANT_TOKEN_RE.

  • Bounds checking:

    Always check 0 <= (variant_pos - seq_start) < len(window_seq) before indexing.

  • Ref allele validation:

    Before writing alt allele, validate that the observed base matches the expected ref allele (possibly after complementing for ‘-’ strand depending on your policy).

Logging#

This module uses the standard logging module. The caller should configure logging handlers/levels (e.g., via logging.basicConfig) if runtime diagnostics are desired.

Functions

apply_complement(base)

Return Watson-Crick complement for A/T/C/G; otherwise return base unchanged.

open_output(out_path)

Create parent directories and return a Path for appending outputs.

parse_variant_block(fasta_header)

Parse a FASTA header and extract variant coordinates.

parse_variant_block_flexible(fasta_header)

Parse a FASTA header and extract variant coordinates (flexible token search rule).

read_fasta(fasta_path)

Read a FASTA file (supports wrapped / multi-line sequences).

substitute_base(seq, pos0, alt)

Return a new sequence where seq[pos0] is replaced by alt.

Classes

ModelHub(transformer_path, device)

Caches heavy models & tokenizers for the GWAS workflow.
utils.variant.read_fasta(fasta_path)[source]#

Read a FASTA file (supports wrapped / multi-line sequences).

This reader supports multi-line (wrapped) FASTA sequences. Each record begins with a header line starting with ‘>’ and is followed by one or more sequence lines. Sequence lines are concatenated and returned in upper-case.

Parameters:

fasta_path (Path) – Path to a FASTA file on disk.

Return type:

Tuple[List[str], List[str]]

Returns:

Tuple[List[str], List[str]] – A tuple (headers, seqs) where:

  • headers: List of header lines (including the leading ‘>’), one per record.

  • seqs: List of concatenated, upper-cased sequences, one per record.

Raises:

Notes

  • Empty/blank lines are ignored.

  • This function does not validate alphabet (A/C/G/T/U/N). If you need strict validation, do it downstream.

utils.variant.open_output(out_path)[source]#

Create parent directories and return a Path for appending outputs.

Return type:

Path

utils.variant.parse_variant_block(fasta_header)[source]#

Parse a FASTA header and extract variant coordinates.

This is the original parsing rule used by the GWAS branch, kept intact for backward compatibility.

Expected header (example):
>variant_1 chr1:27891903-27892003(-)[...]{NA} 27891953:T>A ...
Token usage in the original implementation:
fields = fasta_header.lstrip('>').split()

  • fields[1] is the region token like: chr1:27891903-27892003(-)[...] We parse:

    • strand: text between ‘(’ and ‘)’, e.g. ‘+’ or ‘-’

    • seq_start: window start coordinate, the first number after ‘:’

  • fields[-2] is the variant token like: 27891953:T>A We parse:

    • variant_pos: genomic position (int)

    • ref_base: reference base (str)

    • alt_base: alternate base (str)

Parameters:

fasta_header (str) – FASTA header line including the leading ‘>’.

Return type:

Tuple[int, str, str, str, int]

Returns:

Tuple[int, str, str, str, int] – (variant_pos, ref_base, alt_base, strand, seq_start) where:

  • variant_pos (int): Genomic coordinate of the variant.

  • ref_base (str): Reference allele base (A/C/G/T).

  • alt_base (str): Alternate allele base (A/C/G/T).

  • strand (str): ‘+’ or ‘-’ parsed from the region token.

  • seq_start (int): Genomic coordinate of the window start (used to compute 0-based index into the sequence).

Raises:

ValueError – If the header does not contain enough tokens to parse with this rule (e.g., fewer than 3 whitespace-separated fields).

Notes

  • This parser assumes fixed token positions. If your headers contain extra trailing tokens (e.g., cell line names), consider using parse_variant_block_flexible as a fallback.

utils.variant.parse_variant_block_flexible(fasta_header)[source]#

Parse a FASTA header and extract variant coordinates (flexible token search rule).

This parser is designed for headers where the variant token is not necessarily at a fixed index (e.g. when the last two tokens are cell-line names). It is used by the ribosnitches-after branch, but can also serve as a fallback when parse_variant_block() fails.

Parsing strategy:

  1. Split header into tokens: fields = fasta_header.lstrip('>').split()

  2. Locate: - region token: a token containing ‘:’, ‘-’, ‘(’ and ‘)’ - variant token: matches ^\d+:[ACGT]>[ACGT]$ (case-insensitive)

  3. Extract: - strand and seq_start from region token - variant_pos/ref/alt from variant token

Parameters:

fasta_header (str) – FASTA header line including the leading ‘>’.

Return type:

Tuple[int, str, str, str, int]

Returns:

Tuple[int, str, str, str, int] – (variant_pos, ref_base, alt_base, strand, seq_start) where each field has the same meaning as in parse_variant_block.

utils.variant.apply_complement(base)[source]#

Return Watson-Crick complement for A/T/C/G; otherwise return base unchanged.

Return type:

str

utils.variant.substitute_base(seq, pos0, alt)[source]#

Return a new sequence where seq[pos0] is replaced by alt.

Parameters:

  • seq (str) – Input sequence (window).

  • pos0 (int) – 0-based index into the window.

  • alt (str) – Alternate allele to write at pos0.

Return type:

str

Notes

  • If seq[pos0] already equals alt, we return the original string.

class utils.variant.ModelHub(transformer_path, device)[source]#

Bases: object

Caches heavy models & tokenizers for the GWAS workflow.

The tokenizer/transformer are loaded once and held for reuse. BRIDGE checkpoints are cached by filename stem to avoid repeated disk loads in long FASTA batches.

device#

Inference device (CPU/CUDA).

Type:

torch.device

tokenizer#

Tokenizer loaded from transformer_path.

Type:

BertTokenizer

transformer#

Transformer encoder loaded from transformer_path and set to eval().

Type:

BertModel

bridge_cache#

Cache mapping filename_stem to loaded BRIDGE models.

Type:

Dict[str, BRIDGE]

load_bridge(model_dir, filename_stem)[source]#

Load (or reuse cached) BRIDGE checkpoint: <model_dir>/<filename_stem>.pth.

Parameters:

  • model_dir (Path) – Directory containing .pth checkpoints.

  • filename_stem (str) – Stem used to construct checkpoint name.

Return type:

Optional[BRIDGE]

Returns:

Optional[BRIDGE] Loaded BRIDGE model in .eval() mode, or None if the file does not exist.

Contents