How It Works¶

ResistanceProfiler compares a new sample against a curated project database and reports the resistance mutations it finds. This page explains what happens behind the scenes, from input file to final report.

The core idea: ResPro stores curated rules against internal project references, then maps new FASTA- or VCF-based inputs back into that internal reference space before comparing amino-acid mutations against the rules.

Pipeline overview¶

flowchart TD
    A[Input: FASTA or VCF + reference FASTA] --> B[Reference matching]
    B --> C[Mutation detection]
    C --> D[Rule evaluation]
    D --> E[Report generation]

Build and use one consistent project database¶

The project database stores internal references, CDS annotations, and curated resistance rules. All later profiling steps are anchored to this database. References and rules are matched during database creation to ensure internal consistency.

Input parsing¶

• VCF mode uses a variant file plus a reference FASTA.
• FASTA mode uses a consensus sequence directly.

In both modes, input is converted into a common internal representation before rule matching.

FASTA-mode specifics¶

• IUPAC ambiguity codes in the consensus sequence are expanded into all possible alternative bases. Each mutation ALT receives a fractional allele frequency of 1 / len(options), where options is the full IUPAC base set (the reference base is one of the equally-likely possibilities = no mutation, so it is included in the denominator but not emitted as an ALT). For example, ref=A, query=R (R={A,G}) → G=0.5; ref=A, query=Y (Y={C,T}) → C=0.5, T=0.5; ref=A, query=N (N={A,C,G,T}) → C=0.25, G=0.25, T=0.25.
• N-stretch coverage gaps — full-NNN codons are treated as non-covered positions and reported as coverage gaps rather than variants. Partial-N codons (1–2 N bases) remain assessable and emit expanded IUPAC variants for the non-N positions.
• FASTA-mode AF bins use adjusted thresholds for the discrete IUPAC-derived frequencies: high (0.75–1.0), intermediate (0.35–0.74), low (0.01–0.34).
• VCF-mode AF bins use: high (0.75–1.0), intermediate (0.25–0.7499), low (0.01–0.2499).
• VCF-mode AF source — for VCF input, per-allele allele frequency is resolved from a fixed INFO-then-FORMAT precedence (INFO/AF → INFO/VAF → INFO/FREQ → FORMAT/AF → FORMAT/AD-derived), reading only the first sample for FORMAT-level values. Missing entries and short allele-specific arrays use a residual fallback (max(0, 1 - sum(known)) split equally among missing alleles) rather than assuming a missing allele is fully present. See CLI Reference — Allele-frequency source for the full contract.
• Insertions and deletions in FASTA mode are detected from alignment gaps (CIGAR-based) between the aligned query and the reference CDS.

Reference matching and coordinate mapping¶

Query sequence context is aligned to the internal project references and features. Only rule-backed features (has_rules = 1, i.e. features carrying at least one resistance rule) participate in alignment; a FASTA record that aligns to no rule-backed feature is dropped with a warning rather than profiled. The reference is determined automatically from minimap2-based (mappy) CDS matching, and the sequence with the highest identity is selected.

CIGAR-based coordinate remapping¶

For VCF input, variants are defined in user-provided reference coordinates. ResPro uses the CIGAR string from the minimap2 alignment to build a bidirectional coordinate map between query positions and internal CDS positions:

• Each VCF variant position is projected from the user reference into the internal CDS coordinate space.
• REF alleles are verified against the user query FASTA; mismatches produce a warning and the variant is skipped.
• For reverse-strand features, alleles are reverse-complemented to the internal forward strand.
• Anchor-changed indels (where the VCF anchor base differs from the internal reference) are automatically split into a substitution plus a canonical indel before annotation.

Spliced genes and unspliced queries¶

When a project feature is a spliced CDS (more than one segment in feature_segment) and the user supplies an unspliced whole-genome query, minimap2 reports the inter-exon intron as a single large I (FASTA) or D (raw mappy) operation in the CIGAR. Without special handling this intron would be misinterpreted as a giant coding insertion (e.g. a multi-kilobase frameshift) and would crash per-exon identity.

ResPro classifies such intron operations:

• The exon-junction CDS offsets are derived from feature.segments in genomic 5'→3' order (matching the normalized CIGAR's walking order for both strands).
• A CIGAR I op is classified as an intron when both hold: (1) its CDS position coincides with a known junction offset within alignment.intron_junction_tolerance (default 5 nt; configurable in defaults.toml), and (2) its length is strictly greater than that same tolerance.
• Classified intron I ops are removed from the CIGAR stored on the FeatureMatch (producing an exon-only CIGAR) and recorded as IntronIntervals (CDS junction position, query span, length).
• Identity and CDS coverage are recomputed over exons only, so a perfect exons match reports ~99% identity rather than the genomic span including the intron.
• In FASTA mode the intron query span is stripped from the region before codon walking, so no intron insertion variant is emitted. Real coding insertions (any I op not at a junction within tolerance, or of length ≤ tolerance) are still emitted.
• In VCF mode the query-to-CDS coordinate map skips intron query positions and offsets exon-2 (and later) CDS positions past the intron span, so variants inside the intron are excluded and exon-2 variants remap to the correct CDS offset.

Amino-acid consequence interpretation¶

Per-feature nucleotide changes are translated into amino-acid consequences. Supported classes include:

• Synonymous
• Missense
• Stop changes
• Frameshift
• Insertion
• Deletion
• Unknown

Multiple SNPs in one codon¶

When two or more SNPs fall within the same codon, they are evaluated together as a combined codon event. Each member SNP is emitted as its own per-SNP annotation carrying the amino-acid outcomes it can participate in — both the single-exchange effect (this SNP alone, co-codon SNPs absent) and the combined-state effects (this SNP plus co-occurring SNPs), each with its amino-acid frequency.

There are two paths, selected by whether an alignment file is supplied:

• No BAM (--bam omitted or the mode is fasta) — lower-bound frequencies. ResPro uses Fréchet (probability) bounds to determine which combined codon states are guaranteed to exist, and reports the conservative lower-bound frequency (assuming no linkage information). These frequencies are tagged lower bound in the report.
• With BAM (--bam provided) — observed frequencies. ResPro slices the BAM for each affected codon, counts the spanning reads that carry each candidate codon state (after MAPQ filtering and paired-end dedup/agreement), and reports the exact observed co-occurrence frequency. The Fréchet path is not taken at all. These frequencies are tagged observed in the report.

A single-exchange that is impossible (guaranteed absent, lower bound 0 on the no-BAM path; observed frequency 0 on the BAM path) is omitted from the report and cannot match a single-amino-acid rule — only its combined states can fire. On the no-BAM (Fréchet) path, a combined state is accepted only when its forced fraction reaches 2/3 (min_cooccurrence_codon_fraction in the config). On the BAM path there is no per-state frequency threshold: every candidate state observed in at least one spanning molecule is accepted, and the only gate is the spanning-molecule count (min_depth) — below it the codon falls back to single events.

The report formats each amino-acid effect as <mutation> | <frequency> | <tag>, where the tag is observed (the value is the frequency itself — variant allele frequency or BAM read count) or lower bound (a guaranteed minimum that the true frequency may exceed, for combined codons without read-level data).

For the full mathematical derivation, worked examples, the multiallelic-site generalisation, and the BAM co-occurrence algorithm, see Technical Reference §5.3 — Multiple SNPs in one codon and §5.4 — BAM-based exact co-occurrence.

Rule matching¶

• Single-mutation rules are matched directly against the observed amino-acid events.
• Combination rules evaluate boolean expressions (AND, OR, NOT, XOR) over their atomic member IDs. A member is considered present only when its matched effect's amino-acid frequency lower bound exceeds eps. An AND clause fires only when the joint Fréchet lower bound's forced fraction reaches min_cooccurrence_combination_fraction (default 2/3) — a partial combination whose members do not co-occur often enough does not fire. The reported hit frequency is the joint Fréchet lower bound: a guaranteed minimum on how often the members co-occur, which can be lower than any individual member's frequency. When multiple OR branches match, the branch with the highest member lower bound is selected deterministically (lexical tiebreak on member IDs).
• Interpretation algorithms extend rule evaluation with additional logic such as phenotype counting, score-based thresholds, and IC50-based drug interpretation. See Interpretation Algorithms for details.

Reporting and exports¶

• The HTML report is always generated.
• Optional JSON, PDF, and TSV exports are available via --export (repeatable). See Output Interpretation for details.

Coverage assessment¶

Non-covered codon positions are identified differently depending on the input mode:

• FASTA mode — any stretch of N characters spanning a full codon (NNN) is treated as a coverage gap. Each contiguous N-run is reported as a CoverageGap entry with the affected feature, codon start, and codon end. Rule positions that fall within a gap are not evaluated.
• VCF mode with BAM — an optional BAM file (aligned against the same query reference as the VCF) can be provided via --bam. Per-base depth is extracted from the BAM and projected through the CIGAR map onto internal CDS coordinates. Codons where any nucleotide falls below --min-depth (default 10) are reported as coverage gaps. This allows the report to flag positions that technically lack sufficient read support even if no variant was called. In the web app's batch VCF mode, an optional BAM can be attached to each sample: BAMs are auto-paired to VCFs by filename stem (e.g. sample1.vcf ↔ sample1.bam) on multi-select upload, and any pairing can be overridden per row. Pairing is positional per uploaded VCF row, not by on-disk filename (uploaded files are stored under anonymous names). A sample without a BAM skips coverage-gap analysis, exactly like a single-VCF run without --bam.
• VCF mode without BAM — coverage gaps are not assessed; the report assumes all positions targeted by variants are covered.