Python Dependencies

immLynx uses basilisk to manage Python dependencies automatically. The first time you call a function that invokes Python, basilisk will create an isolated conda environment containing all required packages (tcrdist3, olga, soNNia, clusTCR, metaclonotypist, and PyTorch with the ESM-2 model). This initial setup may take several minutes but only needs to occur once per installation. Subsequent calls reuse the existing environment with minimal overhead.

Loading the Package and Example Data

We begin by loading immLynx along with scran and scater, which provide utilities for single-cell analysis and visualization. The immLynx_example dataset is a SingleCellExperiment object containing scRNA-seq data from multiple patients, with TCR repertoire information added via scRepertoire.

library(immLynx)
library(scran)
library(scater)

data("immLynx_example")

immLynx_example
#> class: SingleCellExperiment 
#> dim: 2000 500 
#> metadata(0):
#> assays(2): counts logcounts
#> rownames(2000): TRBV7-2 HSPA1B ... AC119396.1 BHLHE40
#> rowData names(0):
#> colnames(500): P17B_AAACCTGCAACACGCC-1 P17B_AAACGGGTCTCCAGGG-1 ...
#>   P20L_TCAACGATCCCTCAGT-1 P20L_TCTGGAAGTGGAAAGA-1
#> colData names(16): orig.ident nCount_RNA ... ident clusters
#> reducedDimNames(1): UMAP
#> mainExpName: RNA
#> altExpNames(0):

Extracting TCR Data

Before running downstream analyses, it is often useful to extract TCR data into a data.frame format for inspection and quality control. The extractTCRdata() function wraps immApex::getIR() and supports both long format (one row per chain) and wide format (one row per cell with columns for each chain). The validateTCRdata() function checks for common issues such as missing CDR3 sequences, non-standard amino acids, and inconsistent chain annotations.

tcr_data <- extractTCRdata(immLynx_example,
                           chains = "TRB")
head(tcr_data)
#>             cdr3_aa        v    d       j     c                 barcode chain
#> 1     CAISEQGKGELFF TRBV10-3 <NA> TRBJ2-2 TRBC2 P17B_AAACCTGCAACACGCC-1   TRB
#> 2 CASSVRRERANTGELFF    TRBV9 <NA> TRBJ2-2 TRBC2 P17B_AAACGGGTCTCCAGGG-1   TRB
#> 4 CASSVRRERANTGELFF    TRBV9 <NA> TRBJ2-2 TRBC2 P17B_AACTTTCTCAGCGACC-1   TRB
#> 6   CASSSRQDSTDTQYF   TRBV27 <NA> TRBJ2-3 TRBC2 P17B_AATCGGTCAAGCCCAC-1   TRB
#> 8 CASSVRRERANTGELFF    TRBV9 <NA> TRBJ2-2 TRBC2 P17B_ACGAGCCTCATTTGGG-1   TRB
#> 9  CASSPRAGTPNTEAFF  TRBV4-1 <NA> TRBJ1-1 TRBC1 P17B_ACGCCAGGTTCATGGT-1   TRB

tcr_wide <- extractTCRdata(immLynx_example,
                           chains = "both",
                           format = "wide")

validation <- validateTCRdata(tcr_data)
print(validation)
#> $valid
#> [1] TRUE
#> 
#> $errors
#> character(0)
#> 
#> $warnings
#> character(0)
#> 
#> $summary
#> $summary$n_rows
#> [1] 306
#> 
#> $summary$n_unique_barcodes
#> [1] 306
#> 
#> $summary$n_unique_cdr3
#> [1] 275
#> 
#> $summary$n_na_cdr3
#> [1] 0
#> 
#> $summary$chains
#> [1] "TRB"

Summarizing TCR Repertoire

The summarizeTCRrepertoire() function computes standard repertoire statistics including clonotype frequencies, diversity indices (Shannon entropy, Simpson index, clonality), CDR3 length distributions, and V/J gene usage. These metrics provide a high-level overview of repertoire structure before detailed analysis.

summary <- summarizeTCRrepertoire(immLynx_example,
                                  chains = "TRB")
print(summary)
#> === TCR Repertoire Summary ===
#> Chain(s): TRB 
#> 
#> --- Basic Statistics ---
#> Total cells with TCR: 306 
#> Unique clonotypes: 275 
#> Clonotype ratio: 0.8987 
#> 
#> --- Diversity Metrics ---
#> Shannon entropy: 5.511 
#> Clonality: 0.0188 
#> Simpson index: 0.0057 
#> Inverse Simpson: 174.69 
#> 
#> --- Top Clonotypes ---
#>            cdr3_aa count  proportion
#>  CASSVRRERANTGELFF    14 0.045751634
#>      CASSPTVAGEQFF     3 0.009803922
#>   CASSQAPFSTSGELFF     3 0.009803922
#>   CASSQDRTGLDYEQYF     3 0.009803922
#>   CASSRLRTGALYEQYF     3 0.009803922
#> 
#> --- CDR3 Length Distribution ---
#> Mean: 14.8 
#> Median: 15 
#> Range: 10 - 21

TCR Clustering with clusTCR

clusTCR groups TCR sequences by CDR3 similarity using the Markov Cluster (MCL) algorithm. The inflation parameter controls cluster granularity: higher values produce more, smaller clusters, while lower values yield fewer, larger clusters. Cluster assignments are stored directly in the colData of the SingleCellExperiment object, making them immediately available for downstream visualization and statistical testing.

sce <- runClustTCR(
  immLynx_example,
  chains = "TRB",
  method = "mcl",
  inflation = 2.0
)

table(sce$clustcr_TRB)

TCR Distance Calculations with tcrdist3

tcrdist3 computes pairwise distances between TCR sequences using a substitution-matrix-based metric that accounts for CDR3 amino acid sequence, V gene, and J gene identity. The resulting distance matrix can be used for hierarchical clustering, nearest-neighbor analysis, or as input to dimensionality reduction methods. Note that tcrdist3 deduplicates sequences internally, so the distance matrix dimensions correspond to the number of unique TCR sequences rather than the total number of cells.

dist_results <- runTCRdist(
  immLynx_example,
  chains = "beta",
  organism = "human"
)

dim(dist_results$distances$pw_beta)

Generation Probability with OLGA

OLGA (Optimized Likelihood estimate of immunoGlobulin Amino-acid sequences) computes the generation probability (Pgen) of a TCR sequence—the probability that V(D)J recombination produces that specific CDR3 amino acid sequence. Low Pgen values indicate sequences that are unlikely to arise by chance, which may suggest antigen-driven selection. The runOLGA() function adds log10-transformed Pgen values to the colData of the input object.

sce <- runOLGA(
  immLynx_example,
  chains = "TRB",
  model = "humanTRB"
)

hist(log10(sce$olga_pgen_TRB), breaks = 50)

The generateOLGA() function can also produce synthetic TCR sequences sampled from the recombination model, which is useful for constructing background distributions or benchmarking analyses.

random_seqs <- generateOLGA(n = 100,
                            model = "humanTRB")
head(random_seqs)

Protein Language Model Embeddings

Protein language models such as ESM-2 learn dense vector representations of amino acid sequences from large-scale protein databases. When applied to CDR3 sequences, these embeddings capture biochemical and structural properties that complement sequence-based distance metrics. The runEmbeddings() function generates per-sequence embeddings and stores them as a dimensionality reduction in the SingleCellExperiment object, enabling visualization with standard tools such as UMAP or PCA. Multiple pooling strategies (mean, cls, max) are available to aggregate per-residue representations into a fixed-length vector.

sce <- runEmbeddings(
  immLynx_example,
  chains = "TRB",
  model_name = "facebook/esm2_t12_35M_UR50D",
  pool = "mean"
)

sce <- scater::runUMAP(sce,
                       dimred = "tcr_esm")
scater::plotReducedDim(sce, dimred = "UMAP")

Metaclone Discovery with Metaclonotypist

Metaclonotypist identifies metaclones—groups of TCR clonotypes that share enough sequence similarity to potentially recognize the same antigen. This approach extends traditional exact-match clonotyping by allowing controlled levels of CDR3 edit distance and tcrdist-based similarity. The max_edits and max_dist parameters control the stringency of metaclone grouping.

sce <- runMetaclonotypist(
  immLynx_example,
  chains = "beta",
  method = "tcrdist",
  max_edits = 2,
  max_dist = 20
)

table(sce$metaclone)

Selection Inference with soNNia

soNNia uses a neural network to model post-thymic selection by comparing observed TCR sequences against a background distribution generated by V(D)J recombination. The model estimates a selection factor for each sequence, quantifying how much more (or less) likely a sequence is to appear in the observed repertoire relative to the recombination baseline. This requires a background sequence file, which can be generated using generateOLGA().

background <- generateOLGA(n = 10000, model = "humanTRB")
write.csv(background, "background.csv", row.names = FALSE)

sce <- runSoNNia(
  immLynx_example,
  chains = "TRB",
  background_file = "background.csv"
)