GPTomics/bio-tcr-bcr-analysis-scirpy-analysis

Version Compatibility

Reference examples tested with: scirpy 0.24+, scanpy 1.10+, anndata 0.10+, mudata 0.3+, awkward 2+

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

• Python: pip show <package> then help(module.function) to check signatures

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

Note: since scirpy 0.13 the AIRR receptor data lives as an awkward array in adata.obsm['airr'], NOT in per-chain adata.obs['IR_VJ_1_*'] columns. Fields are read with scirpy.get.airr(...) after scirpy.pp.index_chains(...). Legacy (pre-0.13) objects must be migrated with scirpy.io.upgrade_schema(). Paired GEX+AIRR is held in a MuData with modalities gex and airr; tool/plot functions take the MuData and namespace their output obs columns as airr:<column>.

scirpy Analysis

"Analyze my single-cell paired TCR/BCR alongside gene expression" -> Ingest AIRR receptor records, QC chain pairing, define clonotypes, and overlay clonality on the transcriptomic embedding, all in one AnnData/MuData object.

• Python: scirpy.io.read_10x_vdj() / read_airr() / from_dandelion(), scirpy.pp.index_chains(), scirpy.tl.chain_qc(), scirpy.tl.define_clonotypes() (TCR) or scirpy.tl.define_clonotype_clusters() (BCR)

The governing principles

Two choices silently bias every downstream expansion, diversity, and overlap number. State both explicitly whenever reporting a result.

Principle 1 - the clonotype definition is a choice, and TCR and BCR need different ones. TCR does not somatically hypermutate, so all progeny of a founding T cell share the exact CDR3 nucleotide sequence: tl.define_clonotypes (implicit metric='identity', sequence='nt') on CDR3-nt plus V and J is correct and defensible. B cells DO hypermutate during affinity maturation, so lineage members are NOT identical - tl.define_clonotypes shatters one true BCR lineage into dozens of fake singletons and destroys expansion and diversity estimates (Gupta 2015 Bioinformatics 31:3356). BCR requires tl.define_clonotype_clusters with a distance metric (normalized_hamming), sequence='nt' (SHM acts on nucleotides), and same_v_gene=True, same_j_gene=True to approximate clonal lineages. Even then scirpy returns clonal CLUSTERS, not germline-rooted phylogenies - hand off to Immcantation/dandelion/Dowser for true lineages, mutation calling, and selection.

Principle 2 - single-cell chain QC is the domain-specific hard part, and blanket filtering biases clonality upward. tl.chain_qc labels each cell (single pair, orphan VJ/VDJ, extra VJ/VDJ, two full chains, multichain, ambiguous). Multichain and TCR+BCR-ambiguous cells are likely doublets and are excluded from clonotype definition regardless. But dropping ALL orphan and extra-chain cells is not free: large clones capture both chains more often, so orphans are enriched for singletons, and deleting them preferentially removes small clones - inflating apparent clonal expansion and deflating diversity. Match the filter to the question, and report it.

Clonotype definition: which function

Goal: Pick the clonotyping approach that matches the receptor's biology.

| Question | Function | metric / sequence | Best when | Fails when | |----------|----------|-------------------|-----------|------------| | TCR clonal identity | define_clonotypes | identity / nt (implicit) | TCR (no SHM); exact founder-lineage identity | Applied to BCR - SHM fragments lineages | | BCR clonal lineage (approx.) | define_clonotype_clusters | normalized_hamming / nt + same_v_gene + same_j_gene | BCR; group SHM-diverged members of one lineage | Threshold not data-derived -> chains/merges clones | | Convergent/functional TCR clusters | define_clonotype_clusters | tcrdist or alignment / aa | Antigen-convergent TCRs (different nt, same specificity) | Interpreted as recombination-event counts | | Reconcile with bulk beta/heavy-only | define_clonotype_clusters | receptor_arms='VDJ' | Matching single-cell to bulk TRB/IGH repertoires | Paired-chain specificity is discarded |

pp.ir_dist computes and caches the VJ/VDJ distance matrices; the subsequent define_* call MUST use the SAME metric and sequence or the cached distances silently mismatch the grouping. For BCR the cutoff for normalized_hamming is a PERCENT distance, not a nucleotide count: cutoff=15 means 15% mismatch (~85% identity). Set it from the bimodal distance-to-nearest-neighbor histogram (within-clone mode near 0 vs between-clone mode).

Clustering parameters and their biology

| Parameter | Options | Meaning / when to change | |-----------|---------|--------------------------| | receptor_arms | all / any / VJ / VDJ | all (default): BOTH VJ (alpha/light) and VDJ (beta/heavy) must match - stringent, high specificity. any rescues single-arm dropout but can merge distinct clones sharing only a beta (beta convergence is real). VDJ mimics bulk beta/heavy-only clonotyping. | | dual_ir | any / primary_only / all | Handles two chains of one arm. ~30% of T cells carry two productive TRA (allelic inclusion, Padovan 1993 Science 262:422) - so extra-VJ is real dual-TCR, not junk. primary_only uses the highest-UMI chain; any links cells sharing any chain; all requires both to correspond. | | same_v_gene / same_j_gene | False / True | Require identical V (and J) gene, not just CDR3. Two cells can convergently share a CDR3 from different V genes; requiring same V/J enforces common ancestry. Turn ON for BCR lineage stringency. | | within_group | 'receptor_type' (default) / obs col | Never merge clonotypes across this grouping. Default stops a B cell and a T cell joining one clonotype; set to sample/patient to forbid cross-sample clonotypes. |

Load VDJ and build the joint object

Goal: Ingest receptor contigs and pair them with gene expression in one MuData.

Approach: read_10x_vdj (or read_airr / from_dandelion) returns an AIRR AnnData; wrap it with the GEX AnnData in a MuData keyed gex/airr, then index chains before any QC or clonotyping.

import scirpy as ir
import scanpy as sc
import mudata as mu

adata_gex = sc.read_10x_h5('filtered_feature_bc_matrix.h5')
adata_airr = ir.io.read_10x_vdj('filtered_contig_annotations.csv')  # returns an AnnData, does NOT modify in place
# ir.io.read_airr(['tra.tsv', 'trb.tsv'])  # AIRR TSV from dandelion/Immcantation/airrflow
# ir.io.from_dandelion(dandelion_obj)      # round-trip a dandelion Dandelion object
# ir.io.upgrade_schema(legacy_adata)       # migrate a pre-0.13 obs-column object first

mdata = mu.MuData({'gex': adata_gex, 'airr': adata_airr})
ir.pp.index_chains(mdata)  # REQUIRED before QC/clonotyping; builds obsm['chain_indices']

Chain QC and question-aware filtering

Goal: Categorize chain pairing and remove doublets without silently biasing clonality.

Approach: Run chain_qc, always drop multichain and TCR+BCR-ambiguous doublets, and decide orphan/extra retention by the downstream question - keep orphans for pure GEX overlay, drop them only for paired-clonotype/specificity work.

ir.tl.chain_qc(mdata)  # writes obs: airr:receptor_type, airr:receptor_subtype, airr:chain_pairing
print(mdata.obs['airr:chain_pairing'].value_counts())

# Always exclude likely doublets from clonotype definition.
drop = ['multichain']
keep_types = mdata.obs['airr:receptor_type'].isin(['TCR', 'BCR'])  # exclude 'ambiguous' (TCR+BCR doublet)
paired = mdata[keep_types & ~mdata.obs['airr:chain_pairing'].isin(drop)].copy()

# For paired-clonotype/specificity analysis also require a complete receptor (drop orphans),
# but note this preferentially deletes small clones -> inflates expansion, deflates diversity.
complete = paired[paired.obs['airr:chain_pairing'].isin(['single pair', 'extra VJ', 'extra VDJ'])].copy()

Define clonotypes - TCR (identity)

Goal: Group T cells sharing an exact CDR3-nucleotide founder rearrangement.

Approach: Cache identity distances, then partition; identity on CDR3-nt plus matching arms is the correct TCR clonotype.

ir.pp.ir_dist(complete, metric='identity', sequence='nt', cutoff=0)
ir.tl.define_clonotypes(complete, receptor_arms='all', dual_ir='primary_only')  # writes airr:clone_id
print('TCR clonotypes:', complete.obs['airr:clone_id'].nunique())

Define clonotypes - BCR (distance clusters)

Goal: Group SHM-diverged B cells of one lineage that identity clonotyping would shatter.

Approach: Use normalized Hamming distance on nucleotides within same-V/same-J partitions - this approximates a clonal lineage; identity clonotyping is WRONG for BCR.

# cutoff=15 is a PERCENT distance for normalized_hamming (15% mismatch ~= 85% identity), NOT 15 nt;
# confirm from the distance-to-nearest-neighbor histogram (bimodal trough).
ir.pp.ir_dist(complete, metric='normalized_hamming', sequence='nt', cutoff=15)
ir.tl.define_clonotype_clusters(
    complete,
    sequence='nt', metric='normalized_hamming',
    receptor_arms='all', dual_ir='any',
    same_v_gene=True, same_j_gene=True,   # enforce common ancestry for lineage-grade clones
)  # writes airr:cc_nt_normalized_hamming (a clonotype-cluster id column)
# For true germline-rooted lineages, SHM, and selection: hand off to immcantation-analysis / dandelion / Dowser.

Clonal expansion and diversity

Goal: Quantify how expanded each clone is and how diverse each group's repertoire is.

Approach: Bin cells by clone size, then compute per-group diversity - but remember both numbers depend entirely on the QC filter and clonotype definition above, so report them alongside.

# target_col is resolved WITHIN the airr modality, so pass the bare name 'clone_id', not 'airr:clone_id'.
ir.tl.clonal_expansion(mdata, target_col='clone_id')  # bins per cell: singleton / 2 / >= 3 (breakpoints=(1, 2))
ir.pl.clonal_expansion(mdata, target_col='clone_id', groupby='airr:receptor_subtype')

# Alpha diversity per group; groupby names a full mdata.obs column, so it keeps its modality prefix.
ir.tl.alpha_diversity(mdata, groupby='gex:sample', target_col='clone_id', metric='normalized_shannon_entropy')

# Pairwise repertoire sharing (public/expanded clones, trafficking) - depth-sensitive; compare at equal depth.
ir.tl.repertoire_overlap(mdata, groupby='gex:sample', target_col='clone_id')
ir.pl.repertoire_overlap(mdata, groupby='gex:sample')

Overlay clonality on the transcriptome

Goal: See which cell states the expanded clones occupy.

Approach: Cluster on GEX independently (never on receptor sequence), then color the transcriptomic UMAP by a clonality column pushed into the GEX modality.

# GEX pipeline lives on mdata['gex']: normalize -> HVG -> PCA -> neighbors -> leiden -> umap (see single-cell/clustering).
mdata['gex'].obs['clonal_expansion'] = mdata.obs['airr:clonal_expansion']
sc.pl.umap(mdata['gex'], color='clonal_expansion')

# clonotype_modularity tests whether a clone's cells are more transcriptionally connected than random
# (needs sc.pp.neighbors on the GEX modality first) - distinguishes a coherent functional clone from scatter.
ir.tl.clonotype_modularity(mdata, target_col='clone_id')

Gene usage and specificity

Goal: Summarize V(D)J segment usage and annotate antigen specificity.

Approach: Plot usage/spectratype directly; for specificity, match receptors against a reference database by sequence distance (not ML prediction).

ir.pl.vdj_usage(mdata, full_combination=False)      # V-D-J segment flow (Sankey/ribbon)
ir.pl.spectratype(mdata, chain='VDJ_1', color='airr:receptor_subtype')  # CDR3-length distribution

# Antigen specificity by sequence match to a reference DB (reuses the ir_dist machinery).
vdjdb = ir.datasets.vdjdb()
ir.tl.ir_query(mdata, vdjdb, metric='identity', sequence='aa')
ir.tl.ir_query_annotate(mdata, vdjdb, include_ref_cols=['antigen.species', 'antigen.epitope'])
# For deeper TCR specificity modelling leave scirpy for tcrdist3 / CoNGA (see specificity-annotation).

Export AIRR

Goal: Hand the receptor table to a bulk/interchange tool.

Approach: Write the AIRR modality as a standard rearrangement TSV; do NOT reconstruct it from stale per-chain obs columns (they no longer exist).

ir.io.write_airr(mdata['airr'], 'scirpy_airr.tsv')
# Pull specific fields for a custom table with the get accessor, not obs indexing:
junction_vj = ir.get.airr(mdata, 'junction_aa', 'VJ_1')   # pandas Series
with ir.get.airr_context(mdata, 'junction_aa', ['VJ_1', 'VDJ_1']):
    pass  # AIRR fields temporarily materialized into obs for grouping/plotting

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | BCR lineages appear as hundreds of singletons; no expansion | Identity define_clonotypes used on B cells; SHM makes members non-identical | Use define_clonotype_clusters with metric='normalized_hamming', sequence='nt', same_v_gene=True, same_j_gene=True | | KeyError: 'IR_VJ_1_junction_aa' / obs receptor columns missing | Pre-0.13 schema assumed; AIRR now lives in obsm['airr'] | Access via ir.get.airr(...) after pp.index_chains; migrate legacy objects with io.upgrade_schema() | | Expansion looks high, diversity looks low vs a collaborator | Blanket-filtered all orphan/extra-chain cells, deleting small clones | Keep orphans for GEX overlay; only drop them for paired-clonotype work, and report the filter | | define_* gives grouping that ignores the chosen metric | pp.ir_dist metric/sequence differ from the define_* call | Match metric and sequence between ir_dist and define_clonotype_clusters | | Clonotype/QC functions error or return nothing | pp.index_chains not run before QC/clonotyping | Run ir.pp.index_chains(mdata) immediately after building the MuData | | Real T/B cells look receptor-negative | GEX-only cells (contig dropout) treated as VDJ-negative | Keep GEX-only cells with NaN clonotype for cell-state analysis; only drop VDJ-only cells failing GEX QC | | Spurious shared/secondary chains in a hyperexpanded sample | Ambient VDJ mRNA from a dominant clone mis-assigned to droplets | Start from CellRanger filtered_contig_annotations (is_cell/high_confidence/productive), then drop secondary chains with very low UMI support (duplicate_count/consensus_count, e.g. < 2-3) before clonotyping | | BCR distance clusters look degraded even with correct settings | CellRanger BCR contigs are not IMGT-numbered and include partial/nonproductive contigs | Reannotate with IgBLAST (dandelion/airrflow) before define_clonotype_clusters, or hand off to Immcantation |

Related Skills

• mixcr-analysis - Process raw single-cell VDJ FASTQ
• immcantation-analysis - Proper BCR clonal lineages and SHM downstream
• specificity-annotation - Antigen-specificity clustering on single-cell clonotypes
• single-cell/data-io - Load and manage the GEX AnnData/MuData
• single-cell/clustering - Cell-state clustering to overlay clonality
• single-cell/doublet-detection - Corroborate multichain doublet calls

References

• Sturm G, Szabo T, Fotakis G, Haider M, Rieder D, Trajanoski Z, Finotello F. Scirpy: a Scanpy extension for analyzing single-cell T-cell receptor-sequencing data. Bioinformatics 2020;36(18):4817-4818. doi:10.1093/bioinformatics/btaa611.
• Suo C, Polanski K, Dann E, et al. Dandelion uses the single-cell adaptive immune receptor repertoire to explore lymphocyte developmental origins. Nature Biotechnology 2024;42:40-51. doi:10.1038/s41587-023-01734-7.
• Gupta NT, Vander Heiden JA, Uduman M, Gadala-Maria D, Yaari G, Kleinstein SH. Change-O: a toolkit for analyzing large-scale B cell immunoglobulin repertoire sequencing data. Bioinformatics 2015;31(20):3356-3358. doi:10.1093/bioinformatics/btv359.
• Padovan E, Casorati G, Dellabona P, Meyer S, Brockhaus M, Lanzavecchia A. Expression of two T cell receptor alpha chains: dual receptor T cells. Science 1993;262:422-424. doi:10.1126/science.8211163.
• Vander Heiden JA, Marquez S, Marthandan N, et al. AIRR Community standardized representations for annotated immune repertoires. Frontiers in Immunology 2018;9:2206. doi:10.3389/fimmu.2018.02206.