gptomics/bioskills/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())

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