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.

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

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())