RNA Annotation¶

GTF-based transcript annotation for RNA-seq variant counting (v5.0.0).

With --gtf, GTF mode appends 17 columns total: exon_boundary_dist (1), transcript_read_counts + transcript_fragment_counts (2), and the 14 ASJD fields (7 metrics × 2 alleles).

GTF Requirements¶

gbcms supports Ensembl and GENCODE GTF annotation files.

Supported Formats¶

Download¶

# Ensembl (recommended for GRCh38)
wget https://ftp.ensembl.org/pub/release-112/gtf/homo_sapiens/Homo_sapiens.GRCh38.112.gtf.gz

# GENCODE
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_46/gencode.v46.annotation.gtf.gz

Parsing Details¶

The GTF loader (rust/src/annotation/) reads plain or gzip/BGZF files (detected from the content, not the extension) and keeps the exon rows of the variant chromosomes, extracting:

• seqname (column 1) → chromosome (with chr prefix normalization)
• start, end (columns 4-5) → 1-based exon coordinates, stored 0-based half-open
• strand (column 7) → +, - or unstranded .
• transcript_id attribute → transcript identifier (a GENCODE .N version is stripped)

Each exon row is checked column by column with noodles-gtf's grammar (the parser used before 6.6.0): positive integer coordinates, a numeric or . score, +/-/. strand, a ./0/1/2 frame, and key value; attributes (an unquoted value ends at its ;). Whitespace runs may separate a key from its value. A row with start after end, a coordinate past 2,147,483,647, a missing or empty transcript_id (after the version is stripped), or bytes that are not UTF-8 text is not loaded. Rejected rows are counted in one warning that names the first one's line number and reason.

A transcript is an ID on one chromosome and strand. An ID that the GTF uses on another chromosome or strand (the same RefSeq NM_ on chrX and chrY in the PARs, alternate loci) is another transcript: its exons, introns and per-transcript counts are its own, and are not sorted in with the first one's. The parser warns once, naming how many IDs are reused and the first, when it finds any among the loaded chromosomes. Exons of one ID at the same chromosome and strand are one transcript, including versions that differ only by suffix (NM_x.1 and NM_x.2).

Variant-Guided Filtering¶

To minimize memory usage, the GTF parser only loads chromosomes that contain variants. If your variant file has variants on 3 chromosomes out of 25 in the GTF, only those 3 are indexed.

Annotation Index Architecture¶

The annotation index is built once at pipeline startup and shared across all counting threads via Arc<AnnotationIndex>.

flowchart LR
    GTF["GTF File"] --> Parse["parse_gtf()"]
    Variants["Variant Chroms"] --> Parse
    Parse --> Index["AnnotationIndex"]
    Index --> COI["COITree per Chrom"]
    Index --> SpliceMask["Splice Sites per Transcript"]
    COI --> ExonDist["exon_boundary_dist()"]
    COI --> TxLookup["Overlapping Transcripts"]
    SpliceMask --> ASJD["Splice Junction Comparison"]

Feature 1: Exon Boundary Distance¶

For each variant, gbcms computes the distance (bp) from its REF span to the nearest annotated exon boundary on its contig, across all transcripts: the least distance from any REF base, and 0 when a boundary lies inside the span (as in Ensembl VEP's overlap test). The span is the REF as gbcms normalizes it, left-aligned and trimmed as in a VCF, so it can differ from the row's coordinates. For an SNV or an insertion (a one-base REF) it is the distance from that base. A pure deletion's span starts at its VCF anchor base, so its distance to a boundary on its left is one less than VEP's, which drops the anchor. The distance is unsigned: exonic and intronic positions both count up from the edge. A boundary is an exon's first base or the first intron base after it, so an exon's last base is at distance 1. The contig is matched in any naming (chr1 ~ 1, chrM ~ MT).

Before 6.6.0 the distance was measured from the variant's first base only (#106).

Output column: exon_boundary_dist (MAF)

Feature 2: Per-Transcript Counting¶

When --gtf is provided, the engine counts ALT reads and fragments per overlapping transcript. This resolves ambiguity when a variant overlaps multiple transcripts with different exon structures.

Algorithm¶

1. For each variant, query the COITree to find all transcripts whose exons hold the variant position (an exon's first through last base; the intron bases beside it are not the exon's).
2. For each overlapping transcript, extract the splice site mask.
3. During counting, reads are attributed to transcripts based on splice junction compatibility:
4. A read's CIGAR N operations (splice junctions) are compared against the transcript's annotated splice sites.
5. Reads with matching splice junctions are counted toward that transcript.
6. A read with no junction fits every transcript it overlaps (no junction rules any out): it counts in the DP and RD (or AD) of each overlapping spliced transcript, and once in the variant's counts (register RJ-32 in Read Judgment).
7. Only spliced transcripts are counted. A single-exon transcript has no splice structure to check reads against, so it gets no per-transcript entry; a missing entry says nothing about the reads.
8. Alleles are classified under the same base-quality rules as the main counts, including the exon-boundary BAQ exception (RNA Splice-Junction Handling); a splice lowers no base's quality in either. For the transcript every read is compatible with, the per-transcript counts match the variant's counts. The exception is when the main counts took a decomposed form of the variant or were adjusted for sibling alleles, since per-transcript counting classifies the variant as written.

Output Columns¶

Feature 3: Aberrant Splice Junction Detection (ASJD)¶

ASJD compares the splice junctions observed in reads against the annotated splice sites from the GTF to detect potential aberrant splicing.

How It Works¶

flowchart TD
    Read(["📖 Splice-Spanning Read"]) --> Extract["Extract CIGAR N junctions"]
    Extract --> Compare{"Junction in annotation?"}
    Compare -->|"Yes"| Annotated["Annotated junction count++"]
    Compare -->|"No"| Novel["Novel junction count++"]
    Annotated --> Summary["Per-variant ASJD summary"]
    Novel --> Summary

For each variant:

1. All reads with CIGAR N operations (splice junctions) are collected.
2. Each observed junction (chrom, start, end) is looked up in the GTF splice mask.
3. Junctions present in the annotation are annotated; absent ones are novel.
4. Counts are stratified by allele (REF vs ALT) for differential analysis. Only reads judged REF or ALT vote: a read carrying another allele is in neither partition, as its junction is not the given ALT's, so LOW_ALT_JUNC can fire where partial_alt is high.
5. Counts are deduped per fragment (by QNAME): a molecule whose R1 and R2 both span the same junction votes once, so the junction totals and the strand-discordance test reflect independent fragments, not mates.
6. Each allele's dominant junction is the one with the most fragments. A tie is not a divergence. When a REF top junction and an ALT top junction are the same splice event, each allele reports its own of that pair and no test is run. "Same" means both ends within 5bp, the tolerance ASJD uses throughout; an exact match is preferred. Otherwise each allele reports, among its tied junctions, the one the other allele's fragments use most (the leftmost only on a further tie), and Fisher's exact test compares the two. The verdict therefore follows the reads, not the coordinates. The choice is deterministic: the same input always gives the same junctions and p-value.
7. Alleles are classified under the main counts' base-quality rules, including the exon-boundary BAQ exception, so the spliced reads at an exon-edge variant are ASJD evidence, as they are in the main counts.

Output Columns (14 ASJD)¶

The 14 ASJD fields are part of the 17 columns GTF mode adds (alongside exon_boundary_dist and the two per-transcript count columns). See Output Formats → ASJD for the complete column reference.

Diagnostic Flags¶

The asjd_diagnostic column provides semicolon-separated QC flags. All junction counts are per fragment (a molecule's R1 and R2 are deduped to one vote):

The floors are ASJD's own junction-evidence minimums (LOW_REF_JUNC / LOW_ALT_JUNC). When RETENTION_DOMINANT is present, any junction flags alongside it (NOVEL_ALT_JUNC, MULTI_JUNCTION, …) describe the allele-classified reads — by the marker's own condition the junction-bearing minority of an intron-retaining population — not the spliced majority, which only the ASJD-2 markers report.

Performance and Memory¶

Memory¶

Time¶

Related¶

• RNA Command — --gtf CLI option
• Output Formats — GTF-aware column reference
• RNA Splice Handling — the splice-aware evidence rule and BAQ
• Architecture — System overview with annotation layer