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.

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.

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

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

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.

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-A02) vs 4-digit (HLA-A02: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