GPTomics/bio-tcr-bcr-analysis-vdjtools-analysis

Version Compatibility

Reference examples tested with: VDJtools 1.2.1+, Java (JRE 8+), R 4.x with ggplot2/reshape2/gridExtra (for Plot* modules), immunarch 0.9+/1.0+

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

• CLI: java -jar vdjtools.jar prints the current routine list; <routine> with no args prints its flags
• R: packageVersion('immunarch') then ?repDiversity / ?repOverlap to confirm .method strings

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

Note: routine names are CamelCase and case-sensitive, and the depth-resampling routine is DownSample (capital S). Run RInstall once so the Plot* modules can call R. VDJtools is post-analysis only: it consumes clonotype tables (from MiXCR etc.), not FASTQ.

VDJtools Analysis

"Compute diversity and compare my TCR/BCR repertoires" -> summarize each repertoire's clonal structure, compare samples at equal depth, and quantify overlap.

• CLI: java -jar vdjtools.jar CalcDiversityStats | CalcPairwiseDistances | TrackClonotypes
• R alternative: immunarch repDiversity(), repOverlap(), repClonality()

The governing principle: diversity is sampling-depth-dependent

A repertoire is a sample of an enormous, unevenly expanded clonal population with a long tail of rare clonotypes, so observed richness never saturates: deeper sequencing keeps discovering new clonotypes. Observed richness, Shannon entropy, clonality, Jaccard, and shared-clonotype counts are all functions of read depth. Comparing raw values across libraries of unequal depth measures depth, not biology -- this is the field's single most common and most invalidating error.

The fix is mandatory before any cross-sample claim: bring all samples to a common depth. Two routes:

• DownSample -x <reads> every sample to a shared depth, then analyze; or
• read CalcDiversityStats at a common depth from its resampled table. CalcDiversityStats emits two tables, diversity.<i>.txt (original) and diversity.<i>.resampled.txt (downsampled to the smallest sample or -x). Use the resampled/normalized values for between-sample comparison; the original table is for within-sample description only.

Choosing the normalization depth is itself a decision: downsampling every sample to the cohort minimum discards data and can leave everyone underpowered if one library is tiny. Set the common depth near the cohort's lower quartile, and EXCLUDE (do not drag everyone down to) any sample far below it -- a sample whose rarefaction curve is still steeply climbing well below the chosen depth is under-sampled and cannot support a diversity claim at all. Report the chosen depth and any excluded samples. PlotQuantileStats and the rarefaction curves show which samples are safe to include.

Rarefaction makes the problem visible: RarefactionPlot draws interpolated + extrapolated diversity-vs-depth curves. Compare samples at a common x, never at curve endpoints of different depth. Extrapolation is reliable only to ~2-3x observed depth and degrades for q=0 (Chao 2014).

Report a Hill profile, not one number

A single index misleads because indices weight the abundance distribution differently. Report the Hill profile -- effective number of clonotypes at orders q=0, 1, 2 -- whose shape (steep drop from q=0 to q=2 = a few dominant clones over a large rare tail) is the informative object (Greiff 2015 Genome Med 7:49; Chao 2014 Ecol Monogr 84:45-67). Two repertoires can share richness yet have opposite clonality.

CalcDiversityStats emits these columns (each with _mean/_std); Gini is NOT among them (it is an immunarch option, not a VDJtools output):

| Column | Hill order | Question it answers | Depth-robustness | |--------|-----------|---------------------|------------------| | observedDiversity | q=0 | How many distinct clonotypes were seen | Poor -- must downsample | | chao1 | q=0 | Nonparametric richness lower bound (uses singletons f1, doubletons f2) | Poor; breaks without count data (f2=0) or if rare clones were pre-filtered; PCR error inflates it | | chaoE | q=0 | Chao richness extrapolated, normalized for cross-sample use | Moderate (VDJtools' preferred richness proxy) | | efronThisted | q=0 | Efron-Thisted lower-bound total diversity | Poor; a lower bound, not the truth | | shannonWienerIndex | q=1 | exp(Shannon), effective number weighting by frequency | Moderate | | normalizedShannonWienerIndex | -- | Pielou evenness H'/ln(S), range 0-1 | Depth-dependent through ln(S) | | inverseSimpson | q=2 | 1/sum(p^2), dominated by abundant clones | Best -- most depth-robust | | d50 | -- | Fewest top clones covering 50% of reads | Poor; coarse descriptor |

Default: report q=0 (chaoE or downsampled observedDiversity), q=1 (shannonWienerIndex), and q=2 (inverseSimpson) together. Feed chao1/efronThisted only genuine count data with singletons and doubletons; on non-UMI, non-error-corrected data, PCR/sequencing errors manufacture singletons and inflate them arbitrarily.

Clonality: the field default and its three flaws

Clonality = 1 - normalizedShannonWienerIndex = 1 - H'/ln(S). It runs 0 (even/polyclonal) to 1 (one clone dominates) and is the near-universal one-number summary because it is bounded and intuitive. State its flaws in any report:

1. Depth-dependent through the ln(S) denominator -- compare clonality only on depth-normalized samples.
2. It discards richness (it is a rescaled evenness): two repertoires with identical clonality can differ 100x in richness.
3. It is dominated by the middle of the abundance distribution, not the top clones a clinician cares about.

Fix: report clonality alongside a q=2 Hill number (inverseSimpson) and a rarefaction curve, never alone.

Overlap: pick a depth-robust metric and hold the match key fixed

CalcPairwiseDistances builds an N x N matrix; OverlapPair compares two samples. Any count-of-shared-clonotypes or set index is dominated by the shallower sample's depth: a clone can only be shared if sampled in both, so the shallow sample caps the intersection. Downsample both samples to a common depth first, and prefer abundance-weighted metrics.

| Metric | Basis | Best when | Fails when | |--------|-------|-----------|------------| | MorisitaHorn | Abundance, size-normalized | Unequal depth; the default choice | -- (near-invariant to depth; dominated by abundant shared clones) | | F2 | Sum of per-clonotype geometric-mean frequencies | Frequency-weighted overlap robust to a single dominant shared clone | -- (preferred VDJtools frequency metric) | | F | Geometric mean of summed shared frequencies | Quick frequency overlap | One large shared clone dominates it | | R | Pearson of log-frequencies over shared clones only | Concordance of abundances among shared clones | Ignores private clones entirely | | D | Shared count / geometric-mean diversities | Descriptive | Numerator (shared count) still depth-biased | | Jaccard | Presence/absence | Equal-depth, denoised samples only | Dominated by the shallower sample's depth |

Overlap magnitude also swings by orders of magnitude with the clonotype match key (-i): nt (strict, few coincidental shares) vs aa (convergent recombination inflates sharing), and whether V/J must match (ntV, ntVJ, aaVJ, ...). Fix one key and hold it constant across every comparison in a study; state it in every figure. TrackClonotypes does ordered all-vs-all intersection for time courses -- a clone scoring "absent" at a timepoint is often a sampling zero, so downsample timepoints to common depth before declaring contraction.

An overlap number is only interpretable against a null: some sharing is expected by chance from convergent recombination of high-Pgen clonotypes. To claim overlap EXCEEDS chance, compare the observed statistic to a background of unrelated-donor pairs, or to shuffled/label-permuted repertoires at the same depth, and for public-clonotype claims condition on generation probability (specificity-annotation). Two related individuals or two timepoints from one host will always overlap more than two random donors regardless of biology.

Public is not antigen-driven

"Public" (a clonotype shared across individuals) is largely an artifact of generation probability, not shared antigen selection. High-Pgen CDR3s -- short, few insertions, near-germline (and fetal-generated) -- are independently produced by many donors, and convergent recombination compounds this at the aa level (Venturi 2006 PNAS 103:18691). So a public/shared count is enriched for stochastic high-Pgen sequences, not evidence of a shared response. To argue antigen association, condition on Pgen (OLGA/IGoR) or intersect with an antigen database (ScanDatabase against VDJdb; immunarch dbAnnotate against VDJdb/McPAS-TCR) -- and even a database hit is a sequence match, not proof of binding. Hand Pgen-aware interpretation off to specificity-annotation.

Segment usage and spectratype

CalcSegmentUsage yields per-sample V/J frequency vectors; CalcSpectratype yields the CDR3-length histogram (PlotFancySpectratype overlays the top-N clones; PlotSpectratypeV stacks by V family; PlotFancyVJUsage is the V-J chord plot).

• Spectratype shape: a Gaussian/bell length distribution indicates a diverse polyclonal (naive-like) repertoire; skew or spikes at particular lengths indicate clonal expansion(s). Weighting by reads shows expansions; weighting by unique clonotypes shows underlying diversity.
• Confound: multiplex-PCR primer sets have V-gene-specific amplification bias, so apparent V/J usage differences between platforms or batches are frequently primer artifacts, not biology (Barennes 2021 Nat Biotechnol 39:236). Compare usage only within one protocol, or use 5'-RACE/UMI data. Usage vectors are compositional (sum to 1): CLR-transform before PCA and check that PC1 is not just depth/batch.

VDJtools vs immunarch

Both consume the same clonotype tables; pick by pipeline, not by metric.

| | VDJtools | immunarch | |--|----------|-----------| | Language | Java CLI (calls R for plots) | R / tidyverse | | Maintenance | Stable, low activity (~1.2.1) | Actively maintained (v1.0 adds airr_*) | | Ingestion | Convert -S <fmt> | repLoad() auto-detects MiXCR/Adaptive/10x/AIRR/VDJtools | | Plotting | Fixed Plot* PDFs | vis() returns editable ggplot objects | | Strengths | Reference F/F2/chaoE + resampled tables; legacy reproducibility; pairs with MiXCR/VDJdb | 10x single-cell, k-mer/motif, publication plots, ML feature matrices |

Prefer VDJtools for CLI/legacy MiXCR pipelines and its exact resampled diversity tables; prefer immunarch for R, single-cell, k-mer/motif, or editable figures. Many groups convert with VDJtools and analyze/plot with immunarch.

Core operations in immunarch (verify .method strings on the installed version):

library(immunarch)
data <- repLoad('samples_dir/')                       # metadata + tidy clonotype tables

repDiversity(data$data, .method = 'raref')            # rarefaction/extrapolation curves (the depth control)
repDiversity(data$data, .method = 'hill')             # Hill profile across q
repDiversity(data$data, .method = 'inv.simp')         # q=2, depth-robust
repClonality(data$data, .method = 'homeo')            # clonal-space homeostasis (Rare..Hyperexpanded bins)
repOverlap(data$data, .method = 'morisita')           # depth-robust overlap; 'jaccard'/'public' are depth-biased
geneUsage(data$data[[1]])                             # V/J usage vector
trackClonotypes(data$data, list('Sample1', 1:10))     # longitudinal tracking
dbAnnotate(data$data, vdjdb, 'CDR3.aa', 'cdr3')       # antigen-database annotation

Prepare and normalize (CLI)

Convert upstream output, drop nonfunctional clones, and downsample to a shared depth before any comparison.

# Import MiXCR clonotypes to VDJtools format (also: -S migec/immunoseq/imgt/vidjil ...)
java -jar vdjtools.jar Convert -S mixcr mixcr_clones.txt converted/

# Keep only functional (in-frame, no-stop) clonotypes for functional-repertoire analysis
java -jar vdjtools.jar FilterNonFunctional -m metadata.txt filtered/

# Optional: remove cross-sample contamination (barcode switching / chimeras) before cross-sample work
java -jar vdjtools.jar Decontaminate -m metadata.txt decontaminated/

# Downsample every sample to a common read depth (capital S; -x/--size = target reads)
java -jar vdjtools.jar DownSample -x 100000 -m metadata.txt downsampled/

The metadata file is tab-delimited: a #file.name sample.id <covariate...> header row, then one row per sample. Most multi-sample routines consume it via -m.

Diversity and overlap (CLI)

# Diversity: emits diversity.<i>.txt (original) AND diversity.<i>.resampled.txt (depth-normalized)
# Compare across samples using the RESAMPLED table only.
java -jar vdjtools.jar CalcDiversityStats -m metadata.txt diversity/

# Rarefaction curves -- the correct visual for comparing diversity across depths (read at common x)
java -jar vdjtools.jar RarefactionPlot -m metadata.txt rarefaction/

# Pairwise overlap; -i sets the clonotype match key (hold constant study-wide).
# Report MorisitaHorn / F2 columns; treat Jaccard as depth-biased.
java -jar vdjtools.jar CalcPairwiseDistances -i aa -m metadata.txt overlap/
java -jar vdjtools.jar ClusterSamples -e MorisitaHorn overlap/ clustered/

# Longitudinal tracking across an ordered sample set (downsample timepoints first)
java -jar vdjtools.jar TrackClonotypes -m metadata_timecourse.txt tracking/

Parse VDJtools output in Python

import pandas as pd

def load_resampled_diversity(prefix):
    '''Load the depth-normalized diversity table for cross-sample comparison.'''
    return pd.read_csv(f'{prefix}.strict.resampled.txt', sep='\t')

def load_overlap_matrix(path, metric='MorisitaHorn'):
    '''Load one depth-robust overlap metric from the pairwise-distance output.'''
    df = pd.read_csv(path, sep='\t')
    return df.pivot(index='1_sample_id', columns='2_sample_id', values=metric)

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | Deeper libraries look 'more diverse' every time | Comparing raw richness/Shannon/clonality across unequal depth -- measuring depth, not biology | DownSample to a common depth or read the .resampled.txt table; compare rarefaction curves at a common x | | Overlap flips when samples are swapped or re-sequenced | Jaccard / public counts dominated by the shallower sample's depth | Downsample both, and report MorisitaHorn or F2 | | Overlap magnitude differs wildly between studies | Different clonotype match key (-i aa vs nt, +/-V/J) | Fix one -i value and state it in every figure | | chao1/efronThisted are NaN or absurdly large | Fed frequency-only or rare-clone-filtered data, or PCR errors created singletons | Provide genuine count data (singletons/doubletons); UMI/error-correct first; or use inverseSimpson | | One clonality number reported as 'the diversity' | Clonality is a rescaled evenness that discards richness and is depth-dependent | Report Hill q=0/1/2 (add inverseSimpson) plus a rarefaction curve | | 'Public' clones claimed as antigen-specific | Publicity is mostly high-Pgen convergent recombination | Condition on Pgen (OLGA) or intersect VDJdb via ScanDatabase; hand off to specificity-annotation | | V/J usage differs between cohorts by platform | Multiplex-PCR primer bias, not biology | Compare usage only within a protocol; CLR-transform before PCA and check PC1 is not depth/batch | | Plot* routine errors on start | R plotting dependencies missing | Run java -jar vdjtools.jar RInstall once |

Related Skills

• mixcr-analysis - Generate input clonotype tables
• repertoire-visualization - Rarefaction, spectratype and overlap figures
• immcantation-analysis - BCR-aware diversity and clonal analysis
• specificity-annotation - Pgen-aware interpretation of public clonotypes
• experimental-design/sample-size - Sequencing depth and power planning
• workflows/tcr-pipeline - End-to-end orchestration

References

• Shugay M, et al. VDJtools: unifying post-analysis of T cell receptor repertoires. PLoS Comput Biol 2015; 11(11):e1004503.
• Chao A, et al. Rarefaction and extrapolation with Hill numbers: a framework for sampling and estimation in species diversity studies. Ecol Monogr 2014; 84(1):45-67.
• Greiff V, et al. A bioinformatic framework for immune repertoire diversity profiling enables detection of immunological status. Genome Med 2015; 7:49.
• Chao A. Nonparametric estimation of the number of classes in a population. Scand J Stat 1984; 11:265-270.
• Venturi V, et al. Sharing of T cell receptors in antigen-specific responses is driven by convergent recombination. PNAS 2006; 103(49):18691-18696.
• Barennes P, et al. Benchmarking of T cell receptor repertoire profiling methods reveals large systematic biases. Nat Biotechnol 2021; 39:236-245.
• ImmunoMind Team. immunarch: an R package for painless analysis of T-cell and B-cell immune repertoires. CRAN / immunarch.com (v1.x).