Bio Tcr Bcr Analysis Immcantation Analysis

Version Compatibility

Reference examples tested with: alakazam 1.3+, shazam 1.2+, scoper 1.3+, dowser 2.x, tigger 1.1+ (Immcantation R suite), plus IgBLAST, Change-O, and PHYLIP/IgPhyML as external dependencies.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Note: createGermlines now lives in dowser (not shazam); BASELINe selection uses calcBaseline/groupBaseline (the old estimateBaseline name is gone); mutation R/S classification is set by regionDefinition, not a fake mutationDefinition=MUTATION_SCHEMES$S5F (that has no S5F member); the clonal threshold must come from findThreshold, never a literature constant.

Immcantation Analysis

"Find the B-cell clones and measure their affinity maturation" -> partition SHM-diverged sequences into clonal families, quantify somatic hypermutation and selection against a reconstructed germline, and build antibody lineage trees.

• R: shazam::distToNearest() + shazam::findThreshold() (threshold), scoper::hierarchicalClones()/scoper::spectralClones() (clones), dowser::createGermlines() + shazam::observedMutations() (SHM), shazam::calcBaseline() (selection), dowser::getTrees() (lineage trees)

The governing principle: the clonal threshold is derived, not assumed

Every downstream number in a BCR analysis -- clone counts, diversity, selection strength, tree topology -- inherits its error from one quantity: the nucleotide-distance cutoff used to group sequences into clonal families. That cutoff is NOT a literature constant. distToNearest computes each sequence's Hamming distance to its nearest neighbor within the same V gene, J gene, and junction length; because unrelated rearrangements almost never share V/J plus a near-identical junction by chance while clonally related sequences differ only by SHM, the resulting dist_nearest distribution is bimodal. findThreshold locates the VALLEY between the clonally-related mode (small distances) and the unrelated mode (large distances). That valley is the per-dataset threshold. A hardcoded threshold = 0.15 is the exact anti-pattern to avoid: the valley shifts with subject, locus, sequencing depth, and chemistry, and a wrong threshold silently merges independent lineages or shatters one clone into many (Gupta 2015 Bioinformatics 31:3356; Nouri 2018 Bioinformatics 34:i341).

If the dist_nearest histogram is UNIMODAL (no clear valley), a fixed threshold is undefined -- switch to spectralClones(method="novj"), whose adaptive local threshold does not require findThreshold.

Why BCR needs a different clonotype definition than TCR

TCR does not hypermutate, so all progeny of a founding T cell share the exact CDR3 nucleotide sequence and exact-CDR3 matching is correct. BCR hypermutates: members of one lineage are NOT identical, so exact-CDR3 shatters a single clone into hundreds of fragments. The field-standard BCR clone groups sequences sharing the same V gene, same J gene, and same junction LENGTH, then clusters within that partition by junction nucleotide distance at the derived threshold. Use nucleotide (not amino-acid) junction distance -- SHM is a nucleotide process and codon degeneracy would blur it.

Verify current best practice against the SCOPer vignette before committing to a method; the spectral vj model is the reason spectral clustering holds diverged clones together where a fixed threshold fragments them.

Pipeline order (load-bearing)

This order is not interchangeable; getting it wrong silently corrupts mutation and selection counts.

0. TIGGER genotype FIRST. An unrecorded personal germline polymorphism otherwise reads as recurrent SHM at a fixed position -- it inflates mutation and selection counts AND adds spurious junction distance that corrupts distToNearest.
1. createGermlines (per-sequence) to reconstruct the D-masked germline BEFORE any mutation counting (mutation = observed vs inferred germline).
2. distToNearest -> findThreshold to derive the threshold.
3. Clonal clustering (hierarchicalClones/spectralClones).
4. createGermlines again per-clone (clone consensus germline), then observedMutations with the CDR3/junction MASKED (the D-masked germline handles this; junctional N/P bases have no template).
5. BASELINe selection (calcBaseline -> groupBaseline) with a codon+motif-aware null -- raw R/S is biased by germline codon structure and SHM hotspot/transition bias, so naive R/S is not selection.
6. Dowser lineage trees.

Immcantation reads and writes one AIRR TSV. Expected columns: sequence_id, v_call, j_call, junction, junction_length, sequence_alignment, germline_alignment_d_mask, clone_id (plus locus and cell_id for single-cell). These are lowercase snake_case; legacy UPPERCASE Change-O names (V_CALL, JUNCTION, CLONE) are deprecated and mixing schemas is a silent failure.

Personalize the germline with TIGGER

Goal: Build the subject's own V-gene genotype so germline polymorphisms are not miscounted as somatic mutations.

Approach: Detect novel alleles from the mutation-frequency-vs-position signature, infer the personal genotype, and re-call V alleles against it before anything downstream.

library(tigger)

ighv <- readIgFasta('IMGT_Human_IGHV.fasta')             # named vector of germline V alleles
novel <- findNovelAlleles(db, germline_db = ighv, v_call = 'v_call', nproc = 1)
genotype <- inferGenotypeBayesian(db, germline_db = ighv, novel = novel, find_unmutated = TRUE)
gt_seqs <- genotypeFasta(genotype, germline_db = ighv, novel = novel)
db <- reassignAlleles(db, genotype_db = gt_seqs)         # collapse ambiguous calls to alleles the subject carries

Derive the clonal threshold

Goal: Obtain the per-dataset nucleotide-distance cutoff that separates clonally related from unrelated sequences.

Approach: Compute each sequence's distance to its nearest same-V/J/length neighbor, then find the valley of the bimodal distribution. Inspect the histogram before trusting the value.

library(shazam)

db <- distToNearest(db, sequenceColumn = 'junction', vCallColumn = 'v_call',
                    jCallColumn = 'j_call', model = 'ham', normalize = 'len', nproc = 1)
# Single-cell: add cellIdColumn='cell_id', locusColumn='locus', onlyHeavy=TRUE
#   (light chains lack the junction diversity to define clones alone)

thr_obj <- findThreshold(db$dist_nearest, method = 'density')   # 'gmm' makes the FP/FN tradeoff explicit
threshold <- thr_obj@threshold                                   # S4 slot; NA/unimodal -> use spectralClones('novj')
plot(thr_obj)                                                    # confirm bimodality before proceeding

Cluster sequences into clonal families

Goal: Group SHM-diverged sequences descended from one naive B cell into clones.

Approach: Cluster within V/J/junction-length partitions at the derived threshold; for single-cell paired data, cluster on heavy chains, then resolve light chains as a separate step.

library(scoper)

results <- hierarchicalClones(db, threshold = threshold, method = 'nt', linkage = 'single')
db <- as.data.frame(results)                       # adds clone_id

# Single-cell paired BCR: cluster on heavy only, then split clones by light-chain V/J.
# The scoper only_heavy/split_light args are DEPRECATED; use dowser::resolveLightChains:
# db <- dowser::resolveLightChains(db)

# Unimodal repertoire (no clear threshold): adaptive, SHM-aware alternative
# db <- as.data.frame(spectralClones(db, method = 'vj',
#     germline = 'germline_alignment', sequence = 'sequence_alignment'))

Reconstruct germline and quantify SHM

Goal: Measure somatic hypermutation as replacement (R) and silent (S) frequency by region, the signal of affinity maturation.

Approach: Rebuild the D-masked clonal germline, then compare each observed V-region to it. Use frequency (not raw counts) when coverage varies, and restrict to the V segment so the untemplated junction is excluded.

library(dowser)

references <- readIMGT('imgt/human/vdj')           # IMGT-gapped V/D/J reference dir
db <- createGermlines(db, references)              # per-clone germline; adds germline_alignment_d_mask

db <- observedMutations(db, sequenceColumn = 'sequence_alignment',
                        germlineColumn = 'germline_alignment_d_mask',
                        regionDefinition = IMGT_V,             # V only; stops before CDR3/junction
                        frequency = TRUE, nproc = 1)
# Adds mu_freq_cdr_r, mu_freq_cdr_s, mu_freq_fwr_r, mu_freq_fwr_s
# For property-based R/S use mutationDefinition = CHARGE_MUTATIONS (or HYDROPATHY/POLARITY/VOLUME).
# S5F is a TARGETING model (HH_S5F) for selection, NOT a mutationDefinition.