This vignette builds on the foundations established in
vignette("immLynx_vignette", package = "immLynx") to explore advanced
analytical workflows for T-cell receptor (TCR) repertoire data. While
the introductory vignette demonstrates each analysis function
individually, real-world studies often require combining multiple methods
to gain complementary views of repertoire biology. For instance,
sequence clustering identifies groups of structurally related TCRs,
while generation probability highlights sequences that are unlikely to
arise from random recombination and may therefore reflect
antigen-driven selection. Protein language model embeddings offer yet
another perspective by capturing learned biochemical features that do
not depend on predefined distance metrics.
This vignette covers the following advanced topics:
If you use immLynx in your research, please cite Belean et al. (2026),
Frontiers in Immunology
doi:10.3389/fimmu.2026.1916535,
along with the original publications for each Python tool you use. See
citation("immLynx") for the full reference.
We begin by loading the required packages and the example dataset. The
immLynx_example object is a SingleCellExperiment containing
scRNA-seq data with TCR annotations from scRepertoire.
library(immLynx)
library(scran)
library(scater)
library(ggplot2)
data("immLynx_example")The MCL algorithm used by clusTCR has a key parameter, inflation, that
controls the granularity of cluster assignments. Lower inflation values
(e.g., 2.0) favor larger, more inclusive clusters, while higher values
(e.g., 3.0) produce tighter, more specific groupings. Choosing the
appropriate inflation parameter depends on the biological question: broad
clusters may capture groups of TCRs with shared structural motifs, while
tight clusters more closely approximate antigen-specificity groups.
Here we run clusTCR with two different inflation parameters on the same
dataset and compare the resulting cluster counts. By using distinct
column_prefix values, both sets of assignments are stored in the same
SingleCellExperiment object for side-by-side comparison.
sce_mcl_2 <- runClustTCR(immLynx_example,
inflation = 2.0,
column_prefix = "mcl_2")
sce_mcl_3 <- runClustTCR(sce_mcl_2,
inflation = 3.0,
column_prefix = "mcl_3")
cat("MCL (inflation=2):",
length(unique(na.omit(sce_mcl_3$mcl_2_TRB))),
"clusters\n")
cat("MCL (inflation=3):",
length(unique(na.omit(sce_mcl_3$mcl_3_TRB))),
"clusters\n")The ESM-2 family of protein language models ranges from 8 million to 15 billion parameters. Larger models generally produce embeddings that better capture structural and functional properties of protein sequences, but they require substantially more memory and computation time. For TCR CDR3 sequences, which are short (typically 10-20 amino acids), smaller models often perform well and may be preferable for large-scale analyses.
The following example compares embeddings from the 35M and 650M
parameter variants. Both sets of embeddings are stored as separate
dimensionality reductions in the SingleCellExperiment and can be
independently visualized with UMAP.
sce_small <- runEmbeddings(
immLynx_example,
model_name = "facebook/esm2_t12_35M_UR50D",
reduction_name = "esm_small"
)
sce_med <- runEmbeddings(
sce_small,
model_name = "facebook/esm2_t33_650M_UR50D",
reduction_name = "esm_medium"
)
sce_med <- scater::runUMAP(sce_med,
dimred = "esm_small",
name = "umap_small")
sce_med <- scater::runUMAP(sce_med,
dimred = "esm_medium",
name = "umap_medium")
p1 <- scater::plotReducedDim(sce_med,
dimred = "umap_small") +
ggtitle("ESM-2 35M")
p2 <- scater::plotReducedDim(sce_med,
dimred = "umap_medium") +
ggtitle("ESM-2 650M")
p1 + p2For paired alpha-beta TCR data, embedding both chains together captures
the joint receptor structure. The chains = "both" option concatenates
the alpha and beta CDR3 sequences before computing embeddings, producing
a single representation per cell that reflects the full receptor
identity. This is particularly useful when alpha chain pairing
contributes to antigen specificity.
sce_paired <- runEmbeddings(
immLynx_example,
chains = "both",
reduction_name = "tcr_paired"
)
sce_paired <- scater::runUMAP(sce_paired,
dimred = "tcr_paired")
scater::plotReducedDim(sce_paired, dimred = "UMAP")ESM-2 and DeepTCR answer different questions. ESM-2 is pretrained on UniProt, so its embeddings encode general protein biochemistry and transfer to any repertoire without fitting anything. DeepTCR’s variational autoencoder is trained on the repertoire you give it, so its latent space describes the structure of that dataset specifically. That makes DeepTCR features sharper for within-dataset structure and unsuitable for comparing across datasets trained separately.
Because both write to a dimensional reduction, they can live in the same object and be compared directly.
sce <- runEmbeddings(sce, chains = "TRB",
reduction_name = "tcr_esm")
sce <- runDeepTCR(sce, chains = "TRB",
latent_dim = 64,
reduction_name = "tcr_deeptcr",
verbose = FALSE)
# Project each into UMAP and compare how they organize the repertoire
sce <- scater::runUMAP(sce, dimred = "tcr_esm", name = "umap_esm")
sce <- scater::runUMAP(sce, dimred = "tcr_deeptcr", name = "umap_deeptcr")
scater::plotReducedDim(sce, dimred = "umap_esm",
colour_by = "clustcr_TRB")
scater::plotReducedDim(sce, dimred = "umap_deeptcr",
colour_by = "clustcr_TRB")A caveat on training size: the VAE needs enough unique sequences to learn anything meaningful. On a few hundred sequences the latent space is dominated by noise, and the reconstruction accuracy reported during training is the quickest way to see it.
immLynx is designed to complement scRepertoire, which defines clonotypes based on exact CDR3 sequence matching. By comparing clusTCR cluster assignments with scRepertoire’s clonotype definitions, we can assess how sequence-similarity-based clustering relates to strict clonotype identity. Ideally, each cluster should contain one or a few closely related clonotypes, and cells sharing the same clonotype should be assigned to the same cluster.
sce_clust <- runClustTCR(immLynx_example,
chains = "TRB")
comparison <- table(
clustcr = sce_clust$clustcr_TRB,
clonotype = sce_clust$CTstrict
)
cat("Number of clustcr clusters:", nrow(comparison), "\n")
cat("Number of unique clonotypes:", ncol(comparison), "\n")Every other tool in immLynx starts from receptor sequences. runScXpand()
goes the other way: it uses scXpand
to predict, from gene expression alone, whether a T cell belongs to an
expanded clone. That makes it usable on the many scRNA-seq datasets with no
paired TCR arm, and on datasets that do have TCR data it gives you a way to
ask how much of clonal expansion is written into the transcriptome.
scXpand’s pretrained models impose three hard requirements on the input:
Genes the model expects but cannot find are silently zero-filled, so a poor
identifier match degrades the prediction quietly instead of erroring. Read
the gene mapping summary that runScXpand() prints.
The first call builds a dedicated basilisk environment (Python 3.11, CPU PyTorch, scxpand) that occupies several gigabytes, and downloads the chosen model from figshare. Both are cached.
listScXpandModels()If the object is indexed by gene symbols, map_symbols = "auto" routes them
through org.Hs.eg.db. One symbol often maps to several Ensembl IDs
(HLA-DRA maps to eight, across alternate haplotypes). multi_map = "first"
takes the lexicographically first, which is reproducible;
multi_map = "expand" emits every candidate, which maximizes overlap with
the model’s panel at the cost of a larger staged matrix.
sce_exp <- runScXpand(
immLynx_example,
model_name = "pan_cancer_autoencoder",
assay = "counts",
sample_col = "Patient",
map_symbols = "auto"
)
summary(sce_exp$scXpand_expansion_prob)
table(sce_exp$scXpand_expansion_pred, useNA = "ifany")When the object carries scRepertoire clone calls, runScXpand() also writes
scXpand’s ground-truth fields. It calls a cell expanded when its clone size
exceeds 1.5 times the median clone size for its sample, tabulating clone
sizes fresh from CTstrict within sample_col. Note that
clonalFrequency is deliberately not used: combineExpression() computes
it under whatever group.by was in effect, so it is not necessarily a
per-sample count.
“Median clone size” is ambiguous, and the two readings give different
labels. median_basis = "clone" (the default) takes the median over unique
clones, which in typical 10x data is 1 because most clones are singletons,
so the rule reduces to clone size at least 2. median_basis = "cell"
takes the median over cells, which is dominated by large clones and is much
stricter.
run <- S4Vectors::metadata(sce_exp)$scXpand
# AUROC of the gene-expression-only prediction against the observed clones.
run$auroc_immLynx
# Per-sample median clone sizes behind the label cutoff.
run$labels$median_by_sample
# How many features actually reached the model.
run$gene_mappingGeneration probability (Pgen) from OLGA provides a window into selection pressures acting on the TCR repertoire. Sequences with low Pgen are rarely produced by V(D)J recombination; their presence in the observed repertoire suggests positive selection, potentially driven by antigen encounter. By comparing Pgen distributions across experimental conditions or sample types, we can identify conditions associated with enrichment of rare, potentially antigen-selected sequences.
sce_pgen <- runOLGA(immLynx_example,
chains = "TRB")
ggplot(as.data.frame(colData(sce_pgen)), aes(x = Type,
y = olga_pgen_log10_TRB)) +
geom_boxplot() +
labs(title = "TCR Generation Probability by Sample Type",
x = "Sample Type",
y = "log10(Pgen)")
low_pgen_cells <- which(colData(sce_pgen)$olga_pgen_log10_TRB < -15)
cat("Cells with Pgen < 10^-15:", length(low_pgen_cells), "\n")tcrdist3 distance matrices provide a continuous measure of TCR
similarity that can serve as input to a variety of unsupervised learning
methods. Here we demonstrate hierarchical clustering of the distance
matrix, which complements the MCL-based approach used by clusTCR. The
complete linkage method ensures that all members of a cluster are
within a maximum distance of each other, producing compact groups. Note
that because tcrdist3 deduplicates input sequences, the distance matrix
dimensions correspond to unique TCR sequences rather than individual
cells.
dist_results <- runTCRdist(immLynx_example, chains = "beta")
d <- as.dist(dist_results$distances$pw_beta)
hc <- hclust(d, method = "complete")
n_unique <- nrow(dist_results$distances$pw_beta)
clusters <- cutree(hc, k = min(50, n_unique))
cat("Distance matrix:", n_unique, "x", n_unique, "\n")
cat("Hierarchical clusters:", length(unique(clusters)), "\n")For datasets with thousands or tens of thousands of cells, computational resources may become a limiting factor. The following strategies can help manage memory and runtime:
chunk_size parameter in runEmbeddings()
controls how many sequences are processed in each batch. Smaller
chunks reduce peak memory usage at the cost of slightly longer runtime.runSymdelNeighbors() avoids the quadratic step
entirely. It indexes sequences by their single-character deletions, so
cost scales with sequence count rather than its square, and it returns
a sparse edge list rather than a dense matrix. This is the practical
option past a few tens of thousands of unique sequences.sce_large <- runEmbeddings(
large_sce,
chunk_size = 16,
pool = "mean"
)
sce_large <- runClustTCR(
large_sce,
inflation = 2.0
)
sample_cells <- sample(colnames(large_sce), 5000)
subset_obj <- large_sce[, sample_cells]
dist_results <- runTCRdist(subset_obj)
# Sparse neighbor search over the full dataset, no subsampling needed
edges <- runSymdelNeighbors(
large_sce,
chains = "TRB",
max_edits = 1,
n_cpu = 4,
return_object = FALSE
)When HLA genotyping data is available, immLynx can test for associations
between metaclone assignments and specific HLA alleles using
runHLAassociation(). This analysis identifies metaclones whose
frequency differs significantly between HLA-positive and HLA-negative
individuals, suggesting potential HLA-restricted antigen recognition.
The function performs Fisher’s exact test for each metaclone-HLA
combination and returns adjusted p-values.
metaclones <- runMetaclonotypist(immLynx_example,
return_input = FALSE)
hla_data <- data.frame(
barcode = metaclones$barcode,
HLA_A_01_01 = sample(c(TRUE, FALSE),
nrow(metaclones),
replace = TRUE),
HLA_A_02_01 = sample(c(TRUE, FALSE),
nrow(metaclones),
replace = TRUE),
HLA_B_07_02 = sample(c(TRUE, FALSE),
nrow(metaclones),
replace = TRUE)
)
hla_results <- runHLAassociation(metaclones,
hla_data)
head(hla_results)All immLynx results can be exported to standard file formats for use
with external tools or for archiving. The convertToTcrdist() function
reformats TCR data into the tabular format expected by standalone
tcrdist3, facilitating interoperability with Python-based analysis
pipelines.
tcr_data <- extractTCRdata(immLynx_example,
chains = "both",
format = "wide")
tcrdist_format <- convertToTcrdist(tcr_data,
chains = "both")
write.csv(tcrdist_format, "tcr_for_tcrdist.csv", row.names = FALSE)
clusters <- data.frame(
barcode = colnames(sce_clust),
clustcr = sce_clust$clustcr_TRB
)
write.csv(clusters, "cluster_assignments.csv", row.names = FALSE)
embeddings <- reducedDim(sce_small, "esm_small")
write.csv(embeddings, "tcr_embeddings.csv")Data Quality: Always validate TCR data before analysis using
validateTCRdata(). This catches common issues such as missing CDR3
sequences, non-standard amino acids, and inconsistent chain
annotations that can cause downstream errors or misleading results.
validation <- validateTCRdata(tcr_data, check_sequences = TRUE)
if (!validation$valid) {
warning(validation$errors)
}Memory Management: For large datasets, use smaller chunk sizes in
runEmbeddings() and consider subsampling strategies for
distance-based methods.
Reproducibility: Set random seeds before running stochastic algorithms such as MCL clustering to ensure reproducible results.
set.seed(42)
sce <- runClustTCR(sce)Multiple Methods: Use multiple clustering or grouping methods and look for consensus. Sequences that are consistently co-clustered across methods are more likely to represent biologically meaningful groups.
Biological Validation: Computational clusters should be validated against biological evidence such as antigen specificity, phenotypic markers, or clinical outcomes whenever possible.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] scater_1.41.3 ggplot2_4.0.3
#> [3] scran_1.41.1 scuttle_1.23.2
#> [5] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
#> [7] Biobase_2.73.2 GenomicRanges_1.65.4
#> [9] Seqinfo_1.3.2 IRanges_2.47.5
#> [11] S4Vectors_0.51.10 BiocGenerics_0.59.12
#> [13] generics_0.1.4 MatrixGenerics_1.25.0
#> [15] matrixStats_1.5.0 immLynx_1.1.3
#> [17] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 viridisLite_0.4.3 vipor_0.4.7
#> [4] dplyr_1.2.1 farver_2.1.2 viridis_0.6.5
#> [7] filelock_1.0.3 S7_0.2.2 fastmap_1.2.0
#> [10] bluster_1.23.1 digest_0.6.39 rsvd_1.0.6
#> [13] lifecycle_1.0.5 cluster_2.1.8.3 statmod_1.5.2
#> [16] magrittr_2.0.5 compiler_4.6.1 rlang_1.3.0
#> [19] sass_0.4.10 tools_4.6.1 igraph_2.3.4
#> [22] yaml_2.3.12 knitr_1.52 S4Arrays_1.13.2
#> [25] dqrng_0.4.1 reticulate_1.47.0 DelayedArray_0.39.8
#> [28] RColorBrewer_1.1-3 abind_1.4-8 BiocParallel_1.47.0
#> [31] withr_3.0.3 grid_4.6.1 beachmat_2.29.3
#> [34] edgeR_4.99.10 scales_1.4.0 dichromat_2.0-1
#> [37] cli_3.6.6 rmarkdown_2.32 otel_0.2.0
#> [40] metapod_1.21.0 immApex_1.7.2 ggbeeswarm_0.7.3
#> [43] cachem_1.1.0 stringr_1.6.0 parallel_4.6.1
#> [46] BiocManager_1.30.27 XVector_0.53.0 basilisk_1.25.0
#> [49] vctrs_0.7.3 Matrix_1.7-6 jsonlite_2.0.0
#> [52] dir.expiry_1.21.0 bookdown_0.48 BiocSingular_1.29.2
#> [55] BiocNeighbors_2.7.3 ggrepel_0.9.8 beeswarm_0.4.0
#> [58] irlba_2.4.1 locfit_1.5-9.12 limma_3.99.3
#> [61] jquerylib_0.1.4 glue_1.8.1 codetools_0.2-20
#> [64] stringi_1.8.9 gtable_0.3.6 ScaledMatrix_1.21.0
#> [67] tibble_3.3.1 pillar_1.11.1 htmltools_0.5.9
#> [70] R6_2.6.1 evaluate_1.0.5 lattice_0.23-1
#> [73] png_0.1-9 bslib_0.12.0 Rcpp_1.1.2
#> [76] gridExtra_2.3.1 SparseArray_1.13.4 xfun_0.62
#> [79] pkgconfig_2.0.3