# Immunogenomics

> Analyze scTCR/BCR-seq with scirpy on 10x VDJ contigs, type HLA with OptiType, and score neoantigens with NetMHCpan/pVACseq. Use when doing clonotype/repertoire analysis, HLA typing, or building neoantigen pipelines.

- Skill: `pavel-kravchenko/immunogenomics` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/immunogenomics`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/immunogenomics/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/immunogenomics

---


# Immunogenomics

## When to Use

- Analyzing 10x Genomics scTCR-seq/scBCR-seq data (`filtered_contig_annotations.csv`) jointly with paired gene-expression data
- Defining clonotypes, computing clonal expansion, or plotting clonotype networks from single-cell V(D)J data
- Chaining a full somatic-mutation-to-neoantigen pipeline: VEP annotation to peptide extraction to NetMHCpan binding to pVACseq prioritization
- Typing patient HLA alleles from NGS reads before neoantigen or immunotherapy-response analysis
- Reconstructing bulk TCR/BCR repertoires from tumor RNA-seq with TRUST4 when no single-cell data exists

## Version Compatibility

scirpy ≥0.13, muon ≥0.1, anndata ≥0.10, scanpy ≥1.10, Python ≥3.10; TRUST4 ≥1.0; OptiType 1.3.5; NetMHCpan 4.1; pVACtools ≥4.0.

## Prerequisites

- `pip install scirpy muon scanpy pandas numpy scipy`
- 10x Cell Ranger `vdj` output (`filtered_contig_annotations.csv`) and, for joint analysis, matched `filtered_feature_bc_matrix.h5`
- Basic immunology: V(D)J recombination, CDR3 structure, MHC class I/II presentation (see `bio-applied-vdj-biology`, `bio-applied-hla-typing` for the deep dive on each)

## scirpy: Single-Cell V(D)J Analysis

**Goal:** define clonotypes, quantify clonal expansion, and visualize repertoire structure from paired scTCR/BCR-seq + scRNA-seq.
**Approach:** load VDJ contigs into an AnnData, pair it with the GEX AnnData in a MuData object, compute CDR3 sequence distances, then call scirpy's clonotype/diversity/plotting functions on the combined object.

```python
import scanpy as sc
import scirpy as ir
import muon as mu


def load_paired_vdj_gex(gex_h5, vdj_csv):
    """Build a MuData combining 10x gene expression and V(D)J contigs.

    gex_h5: path to Cell Ranger filtered_feature_bc_matrix.h5
    vdj_csv: path to Cell Ranger filtered_contig_annotations.csv
    Returns a MuData with 'gex' and 'airr' modalities aligned on cell barcode.
    """
    adata_gex = sc.read_10x_h5(gex_h5)
    adata_gex.var_names_make_unique()
    adata_vdj = ir.io.read_10x_vdj(vdj_csv)
    mdata = mu.MuData({"gex": adata_gex, "airr": adata_vdj})
    return mdata


mdata = load_paired_vdj_gex("filtered_feature_bc_matrix.h5", "filtered_contig_annotations.csv")

ir.pp.index_chains(mdata)
ir.tl.chain_qc(mdata)
# drop cells with only an orphan chain (no valid receptor pair)
mdata = mdata[mdata.obs["airr:chain_pairing"] != "orphan VDJ"].copy()

# clonotype definition: CDR3 amino-acid identity across both receptor arms
ir.pp.ir_dist(mdata, metric="identity", sequence="aa")
ir.tl.define_clonotypes(mdata, receptor_arms="all", dual_ir="primary_only")
ir.tl.clonal_expansion(mdata)

ir.tl.alpha_diversity(mdata, groupby="condition", target_col="clone_id")
ir.pl.vdj_usage(mdata, full_combination=False)
ir.pl.clonotype_network(mdata, color="condition")
```

## Bulk Repertoire from RNA-seq (TRUST4)

**Goal:** reconstruct TCR/BCR CDR3 sequences directly from tumor bulk RNA-seq when no dedicated V(D)J library exists.
**Approach:** TRUST4 scans reads against a V/J/C reference and IMGT sequences to assemble CDR3s and estimate clonal abundance.

```bash
# Extract TCR/BCR-derived reads from an aligned tumor RNA-seq BAM and assemble CDR3s
run-trust4 \
    -b tumor_rna.bam \
    -f hg38_bcrtcr.fa \
    --ref human_IMGT+C.fa \
    --thread 8 \
    -o trust4_output

# Output: trust4_output_report.tsv (CDR3, V, J, C, abundance)
# Downstream diversity/overlap metrics: see bio-applied-immune-repertoire
```

## HLA Typing (OptiType)

**Goal:** determine a patient's 4-digit HLA-A/B/C genotype from WES/WGS/RNA-seq reads, required before any binding prediction.
**Approach:** align reads against an HLA reference, keep only HLA-mapping read pairs, then let OptiType's ILP solver call the most likely allele pair per locus.

```bash
# 1. Extract reads mapping to HLA loci
bwa mem hla_reference.fa sample_R1.fastq.gz sample_R2.fastq.gz \
    | samtools view -b -F 4 > hla_reads.bam
samtools sort -n hla_reads.bam | samtools fastq -1 hla_R1.fq -2 hla_R2.fq

# 2. OptiType HLA class I typing
OptiTypePipeline.py \
    -i hla_R1.fq hla_R2.fq \
    --dna --verbose \
    --outdir hla_typing/ \
    --prefix sample
# Output: sample_result.tsv  ->  HLA-A, HLA-B, HLA-C 4-digit alleles
```

## Neoantigen Prediction Pipeline

**Goal:** rank tumor somatic mutations as candidate neoantigens for vaccine/TCR-therapy design.
**Approach:** annotate somatic variants (VEP), extract 9-11 mer mutant peptide windows around each mutation, score binding with NetMHCpan against the patient's HLA alleles, then classify by %Rank_EL / IC50 thresholds (a simplified stand-in for the full pVACseq pipeline, which also weighs expression and clonal fraction).

```python
def extract_peptides(mut_aa_seq, position, lengths=(9, 10, 11)):
    """Extract all peptide windows of each length that cover a mutated residue.

    mut_aa_seq: full mutant protein sequence
    position: 0-based index of the mutated residue
    """
    peptides = []
    for length in lengths:
        for start in range(max(0, position - length + 1), position + 1):
            end = start + length
            if end <= len(mut_aa_seq):
                peptides.append(mut_aa_seq[start:end])
    return peptides


def classify_binders(df, rank_col="Rank_EL", ic50_col="IC50_nM"):
    """Classify NetMHCpan-4.1 output rows into binding tiers.

    Thresholds (NetMHCpan 4.1 convention): %Rank_EL < 0.5 or IC50 < 50 nM
    is a strong binder; %Rank_EL < 2.0 or IC50 < 500 nM is a weak binder;
    everything else is a non-binder.
    """
    def _level(rank, ic50):
        if rank < 0.5 or ic50 < 50:
            return "Strong Binder"
        if rank < 2.0 or ic50 < 500:
            return "Weak Binder"
        return "Non-binder"

    out = df.copy()
    out["Binding_Level"] = [_level(r, i) for r, i in zip(out[rank_col], out[ic50_col])]
    return out


# 1. Annotate somatic variants:
#    vep -i somatic.vcf -o annotated.vcf --cache --everything
# 2. Extract mutant peptide windows with extract_peptides() per variant
# 3. Score binding with NetMHCpan against the patient's OptiType alleles:
#    netMHCpan -p peptides.txt -a HLA-A02:01,HLA-B07:02 -l 9,10,11 -BA > binding.txt
# 4. Parse binding.txt into a DataFrame with Rank_EL/IC50_nM columns, then:
#    ranked = classify_binders(binding_df).sort_values(rank_col_name)
# 5. For full prioritization (expression, clonal fraction, agretopicity), use pVACseq.
```

## B-Cell Lineage Trees (R / dowser)

**Goal:** reconstruct a BCR somatic-hypermutation lineage tree from a clone's sequence variants.
**Approach:** partition sequences into clones with SCOPer, build germline-rooted lineages with `dowser`, and plot with `ggtree`.

```r
library(dowser)
library(alakazam)

# db: an Immcantation-format AIRR table with clone_id, sequence_alignment,
# germline_alignment_d_mask columns (from IgBlast + SCOPer clonal assignment)
clones <- formatClones(db, traits = "c_call", num_fields = "duplicate_count")

# Build maximum-parsimony trees per clone (requires IgPhyML or a phangorn backend)
trees <- getTrees(clones, build = "pratchet")
plotTrees(trees)[[1]]
```

## Key Databases

- **IMGT** — V/J/D/C gene segment reference sequences and nomenclature
- **VDJdb** — antigen-specific TCR/BCR sequences with HLA restrictions
- **McPAS-TCR** — manually curated pathology-associated TCRs
- **IEDB** — immune epitope database for T-cell/B-cell epitopes

## Pitfalls

- **Chain pairing**: 10x gives paired alpha/beta by default, but some cells carry two TCR-alpha chains or only an orphan chain — run `chain_qc` and filter before defining clonotypes.
- **HLA resolution**: 2-digit (HLA-A*02) vs 4-digit (HLA-A*02:01) typing changes which binding predictions are even valid — always type to 4 digits before NetMHCpan.
- **Neoantigen filtering**: binding affinity alone over-predicts immunogenicity — also require tumor expression, clonal fraction, and antigen processing (TAP/proteasome cleavage) before prioritizing.
- **TRUST4 sensitivity**: needs >50M reads for reliable reconstruction; low tumor purity further reduces sensitivity.
- **Diversity metrics need equal depth**: Shannon/Simpson/D50 are not comparable across samples with different sequencing depth without rarefaction.

## See Also

- `bio-applied-vdj-biology` — clonotype definitions and diversity/richness metrics in depth
- `bio-applied-immune-repertoire` — bulk TRUST4/MiXCR repertoire overlap and clonal tracking
- `bio-applied-hla-typing` — HLA typing and NetMHCpan binding prediction in depth
- `bio-applied-single-cell-scanpy` — scRNA-seq preprocessing for the paired GEX modality