Bio Applied Vdj Biology

V(D)J Recombination and Adaptive Immune Receptor Biology

When to Use

• Explaining V(D)J recombination mechanics, IMGT gene nomenclature (TRBV20-1*01), or CDR3 junction structure to interpret TCR/BCR sequencing results
• Deciding how to define a clonotype (nucleotide-identical, amino-acid-identical, CDR3-only, or distance-based) before comparing repertoires
• Parsing an AIRR-format rearrangement TSV or a 10x Genomics filtered_contig_annotations.csv VDJ output into a clonotype table
• Computing repertoire diversity/clonality (Shannon entropy, Simpson index, Chao1, D50) to compare conditions (tumor vs blood, pre- vs post-treatment)
• Detecting clonal expansion, spectratyping CDR3 length, or flagging public clonotypes against a reference (e.g. VDJdb)

Version Compatibility

• Python ≥3.10, pandas ≥2.0, numpy ≥1.24, scipy ≥1.11, matplotlib ≥3.7
• AIRR Community Rearrangement schema ≥1.4 (https://docs.airr-community.org/)
• 10x Genomics Cell Ranger vdj ≥7.0 output format (filtered_contig_annotations.csv)
• For scRNA-seq-paired VDJ objects use scirpy ≥0.13 or Dandelion ≥0.3 (see bio-tcr-bcr-analysis-scirpy-analysis)

Prerequisites

• pip install pandas numpy scipy matplotlib
• Basic immunology: T-cell/B-cell receptor structure, MHC/peptide presentation
• A clonotype table (AIRR TSV or 10x CSV) with v_call/j_call/junction_aa/duplicate_count columns, or raw contigs to build one

Core biology reference

CDR3 (IMGT positions 105–117) spans V-end Cys → junction (N-nucleotides + D + N-nucleotides) → J-start Phe/Trp — it is the most variable part of the receptor and the basis of clonotype identity. B cells additionally undergo AID-mediated somatic hypermutation after antigen activation (~10⁻³ mutations/bp/generation), further diversifying CDRs — this does not apply to T cells.

Clonotype definition strategies (pick one and state it explicitly before comparing samples):

Building an AIRR Clonotype Table

Goal: turn raw V(D)J calls (as would come from IgBlast/MiXCR/10x Cell Ranger) into an AIRR-style clonotype table with per-clone frequency. Approach: one row per unique clone with v_call/d_call/j_call/junction_aa/duplicate_count; frequency is duplicate_count normalized by total reads/UMIs in the sample.

import numpy as np
import pandas as pd


def build_airr_clonotype_table(n_clones=200, seed=0):
    """Construct an AIRR-format clonotype table (stand-in for parsed
    10x `filtered_contig_annotations.csv` or MiXCR/IgBlast output).

    Returns a DataFrame with clone_id/v_call/d_call/j_call/junction_aa/
    duplicate_count/cdr3_length/frequency, sorted by clone size.
    """
    rng = np.random.default_rng(seed)
    trbv_genes = ["TRBV2", "TRBV5-1", "TRBV6-5", "TRBV9", "TRBV12-3",
                  "TRBV20-1", "TRBV28", "TRBV29-1", "TRBV7-2", "TRBV30"]
    trbj_genes = ["TRBJ1-1", "TRBJ1-2", "TRBJ2-1", "TRBJ2-3", "TRBJ2-7"]
    amino_acids = list("ACDEFGHIKLMNPQRSTVWY")

    def random_cdr3_aa(min_len=8, max_len=18):
        """CDR3 amino acid sequence: starts with conserved Cys, ends with Phe."""
        length = rng.integers(min_len, max_len + 1)
        middle = "".join(rng.choice(amino_acids, length - 2))
        return "C" + middle + "F"

    clone_sizes = np.sort(rng.zipf(1.8, n_clones))[::-1]  # power-law clone sizes
    clone_sizes = np.minimum(clone_sizes, 5000)

    table = pd.DataFrame({
        "clone_id": [f"clonotype{i + 1}" for i in range(n_clones)],
        "v_call": rng.choice(trbv_genes, n_clones),
        "d_call": rng.choice(["TRBD1", "TRBD2"], n_clones),
        "j_call": rng.choice(trbj_genes, n_clones),
        "junction_aa": [random_cdr3_aa() for _ in range(n_clones)],
        "duplicate_count": clone_sizes,
    })
    table["cdr3_length"] = table["junction_aa"].str.len()
    table["frequency"] = table["duplicate_count"] / table["duplicate_count"].sum()
    return table.sort_values("duplicate_count", ascending=False).reset_index(drop=True)


clonotype_table = build_airr_clonotype_table()
print(f"Total clones: {len(clonotype_table)}")
print(f"Top clone frequency: {clonotype_table['frequency'].max():.1%}")
print(clonotype_table.head(10).to_string(index=False))

Repertoire Diversity Metrics

Goal: quantify how dominated a repertoire is by expanded clones (tumor infiltration, acute infection, leukemia) vs how diverse/polyclonal it is. Approach: compute Shannon entropy, normalized Shannon, Simpson index, clonality (1 − normalized Shannon), D50 (clones needed to reach 50% of reads), and Chao1 (estimated true richness from singletons/doubletons) directly from clone frequencies/counts.

def shannon_entropy(freqs):
    """Shannon entropy H = -sum(p_i * ln(p_i)) over nonzero clone frequencies."""
    freqs = np.asarray(freqs)
    freqs = freqs[freqs > 0]
    return -np.sum(freqs * np.log(freqs))


def normalized_shannon(freqs):
    """H' = H / ln(N); 0-1 scale, 1.0 = perfectly uniform repertoire."""
    n = int(np.sum(np.asarray(freqs) > 0))
    h = shannon_entropy(freqs)
    return h / np.log(n) if n > 1 else 0.0


def simpson_d(freqs):
    """Simpson index D = sum(p_i^2); D -> 1 as a single clone dominates."""
    return float(np.sum(np.asarray(freqs) ** 2))


def clonality(freqs):
    """Clonality = 1 - H'; approaches 1 for monoclonal expansions (e.g. leukemia)."""
    return 1 - normalized_shannon(freqs)


def d50_index(freqs):
    """Number of top clones needed to reach 50% of total reads."""
    sorted_freqs = np.sort(np.asarray(freqs))[::-1]
    cum = np.cumsum(sorted_freqs)
    return int(np.searchsorted(cum, 0.5) + 1)


def chao1_estimator(duplicate_counts):
    """Chao1 richness estimator: true clone count inferred from singletons (f1)
    and doubletons (f2) among observed clones. duplicate_counts: per-clone
    read/UMI counts (S_obs = number of observed clones).
    """
    counts = np.asarray(duplicate_counts)
    s_obs = len(counts)
    f1 = np.sum(counts == 1)
    f2 = np.sum(counts == 2)
    if f2 == 0:
        return s_obs + f1 * (f1 - 1) / 2  # bias-corrected form when f2 = 0
    return s_obs + (f1 ** 2) / (2 * f2)


freqs = clonotype_table["frequency"].values
print(f"Shannon H:      {shannon_entropy(freqs):.3f}")
print(f"Normalized H':  {normalized_shannon(freqs):.3f}")
print(f"Simpson D:      {simpson_d(freqs):.4f}")
print(f"Clonality:      {clonality(freqs):.3f}")
print(f"D50:            {d50_index(freqs)}")
print(f"Chao1 richness: {chao1_estimator(clonotype_table['duplicate_count']):.1f}")

Comparing V-Gene Usage and CDR3 Length Between Conditions

Goal: compare TRBV gene usage, CDR3 length distribution (spectratype), and clonal expansion between two samples (e.g. healthy blood vs tumor-infiltrating lymphocytes). Approach: build a per-sample V-gene usage histogram, overlay CDR3-length distributions, and plot rank-abundance curves on a log scale to visualize clonal skewing.

import matplotlib.pyplot as plt


def compare_vgene_usage(sample_a_counts, sample_b_counts, gene_names,
                         label_a="Sample A", label_b="Sample B"):
    """Bar-plot V-gene usage (%) side by side for two repertoires.

    sample_a_counts/sample_b_counts: arrays of per-gene clone/read counts,
    aligned to gene_names.
    """
    x = np.arange(len(gene_names))
    w = 0.35
    fig, ax = plt.subplots(figsize=(9, 4))
    ax.bar(x - w / 2, np.asarray(sample_a_counts) / np.sum(sample_a_counts) * 100,
           w, label=label_a, color="steelblue", alpha=0.8)
    ax.bar(x + w / 2, np.asarray(sample_b_counts) / np.sum(sample_b_counts) * 100,
           w, label=label_b, color="firebrick", alpha=0.8)
    ax.set_xticks(x)
    ax.set_xticklabels(gene_names, rotation=45, ha="right", fontsize=8)
    ax.set_ylabel("Usage (%)")
    ax.set_title("V-Gene Usage Comparison")
    ax.legend()
    plt.tight_layout()
    return fig


rng = np.random.default_rng(42)
trbv_all = clonotype_table["v_call"].unique().tolist()
healthy_counts = rng.multinomial(2000, np.ones(len(trbv_all)) / len(trbv_all))
til_probs = np.ones(len(trbv_all))
til_probs[trbv_all.index("TRBV20-1")] = 8.0  # simulate tumor-reactive expansion
til_probs /= til_probs.sum()
til_counts = rng.multinomial(1200, til_probs)

compare_vgene_usage(healthy_counts, til_counts, trbv_all,
                     label_a="Healthy Blood", label_b="Tumor TIL")
plt.savefig("vgene_usage_comparison.png", dpi=120, bbox_inches="tight")

Pitfalls

• Clonotype definition matters: nucleotide-identical, CDR3-amino-acid-identical, and distance-based (Hamming ≤ 1) definitions give different clone counts — pick one and apply it consistently across all samples being compared.
• Sequencing depth confounds diversity metrics: Shannon/Simpson/Chao1 are not comparable across samples with different total read/UMI counts without rarefaction (subsample to equal depth first).
• Use UMI-collapsed counts, not raw reads: bulk TCR/BCR-seq without UMIs inflates apparent clone size via PCR duplication; duplicate_count should reflect unique molecules.
• Not all chains have a D segment: TRA/TRG/IGK/IGL rearrange V+J only — don't apply TRB/IGH-style junctional-diversity or CDR3-length assumptions to light chains.
• 10x filtered_contig_annotations.csv is one row per contig/chain, not per cell: aggregate by barcode and resolve multi-chain/multiplet cells before building clonotypes.
• Multiple testing: apply FDR correction (Benjamini-Hochberg) when testing diversity or clonotype frequency across many samples/timepoints simultaneously.

See Also

• bio-tcr-bcr-analysis-scirpy-analysis — scRNA-seq-paired VDJ analysis with scirpy/AnnData
• bio-tcr-bcr-analysis-vdjtools-analysis — command-line repertoire diversity and overlap statistics
• bio-tcr-bcr-analysis-immcantation-analysis — BCR lineage reconstruction and somatic hypermutation
• bio-tcr-bcr-analysis-repertoire-visualization — clonotype network and repertoire plotting