# Bio Applied Vdj Biology

> Compute TCR/BCR clonotype diversity (Shannon, Simpson, clonality, Chao1, D50) from AIRR/10x VDJ tables; explains IMGT V/D/J nomenclature and CDR3 junctions. Use for repertoire diversity, clonal expansion, or scTCR-seq analysis.

- Skill: `pavel-kravchenko/bio-applied-vdj-biology` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-vdj-biology`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-vdj-biology/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/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

| Locus prefix | Chain | Receptor | Segments |
|---|---|---|---|
| `TRAV`, `TRAJ` | α | TCR αβ | V, J (no D) |
| `TRBV`, `TRBD`, `TRBJ` | β | TCR αβ | V, D, J |
| `TRGV`, `TRGJ` | γ | TCR γδ | V, J (no D) |
| `TRDV`, `TRDD`, `TRDJ` | δ | TCR γδ | V, D, J |
| `IGHV`, `IGHD`, `IGHJ` | heavy | BCR/Ab | V, D, J |
| `IGKV`, `IGKJ` / `IGLV`, `IGLJ` | κ/λ light | BCR/Ab | V, J (no D) |

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

| Strategy | Definition | Use case |
|---|---|---|
| Strict | Identical V + J + CDR3 nucleotide | Tracking exact clonal lineages |
| AA | Same V + J + CDR3 amino acid | Convergent recombination |
| CDR3-only | Identical CDR3 amino acid | Cross-individual "public" clones |
| Distance-based | CDR3 Hamming/Levenshtein distance ≤ 1 | Clonal family clustering (SHM) |

## 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.

```python
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.

```python
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.

```python
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