Create publication-quality visualizations of immune repertoire data including circos plots, clone tracking, diversity plots, and network graphs. Use when generating figures for repertoire comparisons, clonal dynamics, or V(D)J gene usage.
Reference examples tested with: matplotlib 3.8+, seaborn 0.13+, pandas 2.2+, numpy 1.26+; R circlize 0.4+, iNEXT 3.0+; VDJtools 1.2.1+
Before using code patterns, verify installed versions match. If versions differ:
pip show <package> then help(module.function) to check signaturespackageVersion('<pkg>') then ?function_name to verify parameters<tool> --version then <tool> --help to confirm flagsIf code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.
Note: matplotlib 3.9 removed plt.cm.get_cmap(name, N); use plt.get_cmap(name) and sample it, or matplotlib.colormaps[name].resampled(N). seaborn 0.13 deprecated bare palette= without hue=; assign hue= and legend=False.
"Plot my TCR/BCR repertoire" -> Render V-J usage, CDR3-length spectratype, clonal-space, clonal tracking, diversity/rarefaction, overlap, and similarity-network figures, and read each one correctly.
vdjtools PlotFancyVJUsage, vdjtools RarefactionPlot (delegates plotting to R)circlize::chordDiagram (V-J), iNEXT::iNEXT + ggiNEXT (rarefaction/extrapolation)matplotlib/seaborn for bespoke figures; numpy multinomial resampling for rarefaction; networkx + rapidfuzz for similarity networksA repertoire figure is only as valid as the numbers behind it, and those numbers depend on two decisions made before any plot is drawn:
| Figure | What it reveals | How to compare correctly | |--------|-----------------|--------------------------| | V-J chord/circos | Combinatorial V-J pairing bias within one sample | Descriptive per sample; never compare raw usage across primer sets/platforms (primer bias masquerades as biology) | | CDR3 spectratype | Clonal structure: Gaussian length distribution = polyclonal/naive, spikes = expansion | Weight by frequency to see expansions; by unique clonotypes to see underlying diversity | | Clonal-space / proportion | Fraction of repertoire held by rare vs expanded clones | Bin by clone frequency (strata), not raw count; robust to depth if frequency-based | | Clonal tracking | Expansion/contraction/persistence of clones over time | Downsample timepoints to common depth first; an "absent" clone is often a sampling zero | | Rarefaction/extrapolation | Diversity comparison done right (interpolate to common depth) | Read all curves at a shared x; extrapolate at most 2-3x observed depth | | Overlap heatmap | Pairwise repertoire similarity | Use Morisita-Horn (depth-robust) on depth-normalized samples; Jaccard is depth-biased | | Similarity network | Clusters of related CDR3s (candidate specificity groups) | Structure depends entirely on the distance threshold; state it and test sensitivity |
The chord/circos shows which V and J segments pair, weighted by clonotype count, within one sample.
VDJtools route (delegates plotting to R; requires RInstall once):
vdjtools PlotFancyVJUsage -m metadata.txt output_dir/
R with circlize - a V-by-J count matrix becomes a chord diagram:
library(circlize)
plot_vj_chord <- function(clone_df) {
vj_matrix <- table(clone_df$v_gene, clone_df$j_gene)
chordDiagram(vj_matrix, transparency = 0.5, annotationTrack = c('grid', 'name'))
}
Python heatmap alternative (samples x V gene, or V x J for one sample) avoids a chord dependency and reads more quantitatively for many segments:
import seaborn as sns
import matplotlib.pyplot as plt
def plot_vj_heatmap(clone_df):
vj = clone_df.pivot_table(index='v_gene', columns='j_gene', values='frequency', aggfunc='sum', fill_value=0)
fig, ax = plt.subplots(figsize=(8, 6))
sns.heatmap(vj, cmap='viridis', ax=ax)
ax.set_title('V-J pairing frequency')
return fig
The spectratype is a histogram of CDR3 length. A roughly Gaussian, bell-shaped distribution indicates a diverse, polyclonal (naive-like) repertoire; skew or discrete spikes at particular lengths indicate clonal expansion(s). The classic immunoscope spectratype (and VDJtools CalcSpectratype) bins nucleotide length, where in-frame clones sit 3 nt apart and the periodicity is part of the readout; amino-acid length is a common simplification - state which is plotted. Weighting by read/UMI frequency shows expansions; weighting by unique clonotypes shows the underlying diversity - the two views can look opposite, so label which is plotted.
def plot_spectratype(clone_df, length_col='cdr3_length'):
fig, ax = plt.subplots(figsize=(9, 5))
bins = range(clone_df[length_col].min(), clone_df[length_col].max() + 2)
ax.hist(clone_df[length_col], bins=bins, weights=clone_df['frequency'], color='steelblue')
ax.set_xlabel('CDR3 length (aa)')
ax.set_ylabel('Frequency') # Gaussian = polyclonal; spikes = clonal expansion
return fig
Bin clones into frequency strata (Rare / Small / Medium / Large / Hyperexpanded) and stack their summed frequency. Because strata are defined on frequency, this view is comparatively depth-robust and shows clonal-space homeostasis at a glance. A treemap of top clones is an alternative when individual dominant clones matter.
import numpy as np
def plot_clonal_space(clone_df, sample_col='sample'):
# Frequency-based strata (immunarch homeo convention); frequency makes this depth-robust
edges = [0, 1e-4, 1e-3, 1e-2, 1e-1, 1.0]
labels = ['Rare', 'Small', 'Medium', 'Large', 'Hyperexpanded']
clone_df = clone_df.copy()
clone_df['stratum'] = pd.cut(clone_df['frequency'], bins=edges, labels=labels)
space = clone_df.groupby([sample_col, 'stratum'], observed=True)['frequency'].sum().unstack(fill_value=0)
fig, ax = plt.subplots(figsize=(8, 5))
space[labels].plot(kind='bar', stacked=True, ax=ax, colormap='viridis')
ax.set_ylabel('Fraction of repertoire')
return fig
Track individual clone frequencies over an ordered sample set (vaccination, infection, therapy). A line plot of the top clones, or an alluvial/stream for the same data, shows expansion, contraction, and persistence. Downsample timepoints to a common depth before declaring contraction - a clone scored absent at one timepoint is frequently below the sampling floor, not truly gone.
def plot_clone_tracking(clone_df, top_n=10, clone_col='cdr3_aa', time_col='timepoint'):
top = clone_df.groupby(clone_col)['frequency'].sum().nlargest(top_n).index
fig, ax = plt.subplots(figsize=(9, 5))
for clone in top:
d = clone_df[clone_df[clone_col] == clone].sort_values(time_col)
ax.plot(d[time_col], d['frequency'], marker='o', label=clone[:12])
ax.set_xlabel('Timepoint'); ax.set_ylabel('Clone frequency')
ax.legend(bbox_to_anchor=(1.02, 1), loc='upper left', fontsize=7)
return fig
A bare bar of Shannon (or observed richness) across samples of unequal depth is misleading - it plots depth as much as biology. The defensible comparison is a rarefaction/extrapolation curve: interpolate each sample down to (and extrapolate modestly above) a shared depth, then read diversity at a common x. iNEXT computes this for Hill numbers q=0/1/2 with confidence intervals (Hsieh et al. 2016 Methods Ecol Evol 7:1451; Chao et al. 2014 Ecol Monogr 84:45).
R route (gold standard, gives CIs):
library(iNEXT)
# count_list: named list of per-sample integer clonotype-count vectors
out <- iNEXT(count_list, q = c(0, 1, 2), datatype = 'abundance')
ggiNEXT(out, type = 1) # diversity vs sample size, read at a common x
VDJtools route: vdjtools RarefactionPlot -m metadata.txt output_dir/.
Python resampling route (interpolation by multinomial subsampling) when iNEXT is unavailable:
def rarefaction_curve(counts, depths, reps=20, rng=None):
# Interpolate observed richness by drawing 'm' reads without-replacement-like via multinomial
rng = rng or np.random.default_rng(0)
counts = np.asarray(counts, dtype=float)
p = counts / counts.sum()
total = int(counts.sum())
richness = []
for m in depths:
if m > total: # interpolation only; do not extrapolate past observed depth here
richness.append(np.nan); continue
obs = [np.count_nonzero(rng.multinomial(m, p)) for _ in range(reps)]
richness.append(np.mean(obs))
return richness
An N-by-N similarity heatmap summarizes pairwise repertoire overlap. The metric choice is the decision: Morisita-Horn is abundance-weighted and near-invariant to sample size, so it is the default across unequal depths; Jaccard (presence/absence) is dominated by the shallower sample's depth and should not be compared across unequal-depth pairs. Compute the matrix in vdjtools-analysis on depth-normalized samples, then render it here. Label the metric in the title.
def plot_overlap_heatmap(overlap_matrix, metric='Morisita-Horn'):
fig, ax = plt.subplots(figsize=(7, 6))
sns.heatmap(overlap_matrix, annot=True, fmt='.2f', cmap='YlOrRd', vmin=0, vmax=1, square=True, ax=ax)
ax.set_title(f'Repertoire overlap ({metric})') # state the metric; Jaccard is depth-biased
return fig
Nodes are CDR3s, edges connect sequences within a chosen distance, and clusters are candidate specificity groups. The network structure depends entirely on the threshold: too loose chains distinct clones together, too tight fragments real groups. State the threshold and metric, and test sensitivity before interpreting clusters. Same-length CDR3s with Hamming distance is the conservative default; Levenshtein allows indels.
import networkx as nx
from rapidfuzz.distance import Levenshtein
def build_similarity_network(clone_df, max_norm_dist=0.15, clone_col='cdr3_aa'):
# normalized_similarity in [0,1]; edge when 1 - similarity <= threshold. Structure is threshold-dependent.
clones = clone_df[clone_col].unique()
g = nx.Graph()
g.add_nodes_from(clones)
for i, a in enumerate(clones):
for b in clones[i + 1:]:
if 1 - Levenshtein.normalized_similarity(a, b) <= max_norm_dist:
g.add_edge(a, b)
return g
| Symptom | Cause | Fix |
|---------|-------|-----|
| Diversity "differs" between groups but tracks library size | Bare Shannon/richness bar across unequal depth | Downsample to common depth or plot rarefaction curves read at a shared x (iNEXT/RarefactionPlot) |
| Overlap heatmap shows huge differences driven by one shallow sample | Jaccard/shared-count on raw, unequal-depth counts | Use Morisita-Horn on depth-normalized samples; label the metric |
| Network clusters change completely on re-run with a new cutoff | Arbitrary distance threshold; single-linkage chaining | Fix and report the metric + threshold; run a sensitivity sweep; prefer Hamming on same-length CDR3s |
| Two figures disagree on sharing/diversity | Built on different clonotype definitions (nt vs aa, +/- V/J) | Hold one clonotype definition constant across all figures in a comparison and state it |
| Spectratype "expansion" vanishes when re-plotted | Switched between frequency-weighted and clonotype-weighted histogram | Choose one weighting per figure and label it (frequency shows expansions) |
| V-usage differences between batches look biological | Multiplex-PCR primer bias | Compare usage only within one protocol, or use UMI/5'RACE data |
| plt.cm.get_cmap(name, N) raises AttributeError | Removed in matplotlib 3.9+ | Use plt.get_cmap(name) or matplotlib.colormaps[name].resampled(N) |
npx skills add GPTomics/bio-tcr-bcr-analysis-repertoire-visualization下载完整 Skill 目录，包含 SKILL.md 及所有相关文件
Category:science-education