TCR/BCR Repertoire Analysis

The worked walkthrough is TCR-focused: it follows filtered single-cell V(D)J contigs from Cell Ranger into scRepertoire, summarizes clone-size structure, and attaches receptor labels to a Seurat object. A separate, explicitly synthetic BCR example below demonstrates productive-chain filtering, a heavy/light-chain presence check, and a heavy-chain-only summary. Both are descriptive; neither establishes antigen specificity, malignancy, or cell function.

What you need and what you will produce

Use one filtered contig-annotation CSV per sample, produced by a compatible V(D)J workflow, plus a Seurat object containing the corresponding cells and a sample label in its metadata. This example assumes matched TCR and RNA libraries, compatible cell barcodes, and one sample per input directory. It uses R and scRepertoire; package APIs change, so check the linked reference for your installed version.

Work from your project directory (set it explicitly with setwd() or use an RStudio project). Check the runtime with R --version in a shell, then in R record R.version.string and packageVersion("scRepertoire"), packageVersion("Seurat"), and packageVersion("ggplot2"). Install packages using current official instructions for your R version. These snippets are templates; no specific package-version combination is claimed as tested.

Before loading, verify that each sample directory contains filtered_contig_annotations.csv with expected barcode, chain, productive, and CDR3 columns. Cell Ranger describes this as contig-level rows for high-confidence cells, so multiple chain rows per barcode are expected. Record the Cell Ranger run/version and reference configuration.

At the end, record a clone-size summary and a UMAP colored by clone-size category, alongside the exact clone definition, chain handling, sample grouping, and number of receptor-bearing cells included. These are summaries to interpret, not biological findings supplied in advance.

1. Load contigs and make sample-aware clone calls

Cell Ranger's filtered_contig_annotations.csv contains contigs associated with high-confidence cell barcodes; it is not a table with one row per cell because a cell can have multiple receptor-chain contigs. Keep sample identity explicit: barcodes can recur between libraries, and a sample prefix prevents accidental merging of identical barcode strings.

library(scRepertoire)
library(Seurat)
library(ggplot2)

# Replace these directories with one directory per V(D)J sample.
sample_dirs <- c("data/sample_A", "data/sample_B")
sample_ids <- c("sample_A", "sample_B")
stopifnot(length(sample_dirs) == length(sample_ids))

contig_list <- lapply(sample_dirs, function(path) {
  loadContigs(input = path, format = "10X")[[1]]
})
combined_tcr <- combineTCR(
  contig_list,
  samples = sample_ids,
  ID = sample_ids
)

The directories must contain Cell Ranger's named filtered-contig CSV. Check that each list element loaded the intended sample, that barcode prefixes distinguish samples, and that productive chains and missing/extra chains are handled as intended. The defaults for removing or retaining incomplete and multi-chain cells are analysis choices, not neutral cleanup; inspect the package documentation before changing them.

Barcode check before attachment: samples and ID add sample strings as prefixes in combineTCR(); with both set to sample_A, do not assume the result remains the raw 10x barcode. Inspect head(combined_tcr[[1]]$barcode) against head(colnames(seurat_obj)). If RNA barcodes are not sample-qualified, apply one documented convention consistently to both assays before combining; do not strip suffixes blindly, which can merge cells across samples. Require exact overlap and inspect unmatched counts; stop if overlap is unexpectedly low.

2. Compare clone-size structure

Choose a clonotype definition before interpreting an expansion plot. The example uses aa, amino-acid CDR3 calls. Amino-acid matching can group distinct nucleotide rearrangements that encode the same sequence; nucleotide-based definitions are more stringent for sequence identity. The package's available definitions and chain rules depend on current version and receptor type.

# Summarize the fraction of each repertoire in clone-size bins.
clonalHomeostasis(
  combined_tcr,
  clone.call = "aa",
  group.by = "sample"
)

# Show clone contribution by rank; this is a different view of abundance.
clonalProportion(
  combined_tcr,
  clone.call = "aa",
  group.by = "sample"
)

Homeostasis groups clones by their proportion of the selected repertoire; proportion plots rank clones into abundance groups. Neither is a diversity estimate by itself. Read the plot with its denominator and grouping in mind: the same clone may be assigned differently when the compared group changes. Record the clone-call definition, chain selection, grouping variable, and whether cells or reads underpin the input data.

3. Attach clone labels to Seurat metadata

Make the RNA object's sample metadata agree with the sample names used above. The attachment relies on matching cell barcodes. The example groups frequencies within each sample and requests amino-acid clonotypes; inspect the resulting metadata columns rather than assuming a particular label or bin exists.

# Load your matched Seurat object, e.g.:
# seurat_obj <- readRDS("data/combined_seurat.rds")
# Metadata column "sample" must contain values in sample_ids.
stopifnot(exists("seurat_obj"), "sample" %in% colnames(seurat_obj[[]]))
stopifnot(all(unique(as.character(seurat_obj$sample)) %in% sample_ids))
contig_barcodes <- unique(unlist(lapply(combined_tcr, `[[`, "barcode")))
rna_barcodes <- colnames(seurat_obj)
barcode_overlap <- intersect(contig_barcodes, rna_barcodes)
message(length(barcode_overlap), " exact barcode matches; ",
        length(setdiff(contig_barcodes, rna_barcodes)), " contig barcodes unmatched")
if (length(barcode_overlap) == 0L) stop("No exact barcode overlap; resolve prefixes before attachment")
seurat_obj <- combineExpression(
  combined_tcr,
  seurat_obj,
  clone.call = "aa",
  group.by = "sample"
)

# Inspect the values actually added before plotting.
colnames(seurat_obj[[]])
DimPlot(seurat_obj, group.by = "cloneSize")

If cells are unexpectedly unannotated, compare the barcode strings on both sides, including any sample prefixes, and verify that the sample metadata column has the intended values. Also check whether filtering removed receptor-negative cells. The plot describes where labeled cells lie in the embedding; it does not demonstrate that a clone caused a transcriptional state or that the cells are malignant.

Interpreting clone-associated expression

A clone-size category is an abundance grouping, not a malignancy label. In a lymphoma study, determining whether a clone is malignant requires disease-specific evidence and study design beyond this repertoire workflow. A cell-level contrast such as FindMarkers() can be exploratory, but cells from the same donor are not independent biological replicates. For condition-level claims, account for between-sample variation, for example with an appropriate pseudobulk design, and validate findings independently.

TCRs and BCRs share V(D)J-based diversity but differ in chain organization and biology. The main code above is TCR-focused. For a real BCR analysis, use BCR-specific handling such as combineBCR() on filtered contig annotations and verify productive IGH plus IGK and/or IGL chains and CDR3 fields. For combineExpression(), choose a BCR-appropriate clone.call and chain (documented choices include IGH, IGL, and both); the scRepertoire documentation describes its heavy/light-chain handling and related-clone options. A heavy-chain-only call does not represent paired receptors. Check barcode overlap and resulting metadata before plotting. Interpret pairing and somatic hypermutation with BCR-specific methods. The small synthetic example below is a teaching summary, not a validated BCR pipeline.

BCR: reproducible synthetic example

This base R exercise starts from invented, normalized contig rows with one row per chain call. It needs no external data or packages. It checks whether each barcode has a productive heavy chain and at least one productive light chain, then counts distinct cells by an explicitly simplified exact heavy-chain key: sample, V gene, J gene, and CDR3 amino-acid sequence. This key is not a scRepertoire clone call and may split related sequences that differ through somatic hypermutation.

bcr <- data.frame(
  sample_id = rep("toy_A", 8),
  cell_barcode = c("cell_01", "cell_01", "cell_02", "cell_02",
                   "cell_03", "cell_03", "cell_04", "cell_04"),
  chain = c("IGH", "IGK", "IGH", "IGK", "IGH", "IGL", "IGH", "IGK"),
  productive = c(TRUE, TRUE, TRUE, TRUE, TRUE, TRUE, FALSE, TRUE),
  v_gene = c("IGHV1-2", NA_character_, "IGHV1-2", NA_character_,
             "IGHV3-7", NA_character_, "IGHV1-2", NA_character_),
  j_gene = c("IGHJ4", NA_character_, "IGHJ4", NA_character_,
             "IGHJ6", NA_character_, "IGHJ4", NA_character_),
  junction_aa = c("CARGGYYGMDV", NA_character_, "CARGGYYGMDV", NA_character_,
                  "CARRDYGMDV", NA_character_, "CARGGYYGMDV", NA_character_),
  stringsAsFactors = FALSE
)

productive_bcr <- subset(bcr, productive)
chains_by_cell <- split(
  productive_bcr$chain,
  paste(productive_bcr$sample_id, productive_bcr$cell_barcode, sep = ":")
)
has_paired_chains <- vapply(chains_by_cell, function(chains) {
  "IGH" %in% chains && any(chains %in% c("IGK", "IGL"))
}, logical(1))
stopifnot(sum(has_paired_chains) == 3L, sum(!has_paired_chains) == 1L)
cat("cells_with_IGH_and_light=", sum(has_paired_chains),
    "; without_productive_pair=", sum(!has_paired_chains), "\n", sep = "")

productive_igh <- subset(productive_bcr, chain == "IGH")
bcr_clones <- aggregate(
  cell_barcode ~ sample_id + v_gene + j_gene + junction_aa,
  data = productive_igh,
  FUN = function(x) length(unique(x))
)
names(bcr_clones)[names(bcr_clones) == "cell_barcode"] <- "cells"
bcr_clones <- bcr_clones[order(-bcr_clones$cells, bcr_clones$junction_aa), ]
stopifnot(identical(as.integer(bcr_clones$cells), c(2L, 1L)))
write.table(bcr_clones, row.names = FALSE, quote = FALSE, sep = "\t")

Expected output (tab-separated after the first line):

cells_with_IGH_and_light=3; without_productive_pair=1
sample_id	v_gene	j_gene	junction_aa	cells
toy_A	IGHV1-2	IGHJ4	CARGGYYGMDV	2
toy_A	IGHV3-7	IGHJ6	CARRDYGMDV	1

In this toy table, three cells have a productive IGH and at least one productive IGK or IGL; one has only productive IGK after its nonproductive IGH row is excluded. The summary counts invented cell barcodes only. Real data can contain multiple productive contigs, doublets, or barcode mismatches, so inspect chain multiplicity and sample-aware barcode handling before calling paired receptors or interpreting counts. This exact amino-acid grouping does not infer lineage, affinity, antigen specificity, or biological function. For BCR clone calling and somatic-hypermutation analysis, use a version-appropriate BCR method and report its parameters.

Decisions to report with your analysis

• Input and denominator: assay/software output, included samples, and whether abundance counts cells or another unit.
• Clonotype definition: nucleotide or amino-acid sequence, gene requirements, selected chains, and treatment of dual or incomplete chains.
• Comparison: sample/group variable, normalization or subsampling, and receptor-negative cell handling.
• Interpretation: describe expansion as observed abundance. State what independent evidence would be needed for antigen specificity, function, or malignancy.