From 2e153946398be78f0df0b97c29f3741b1fb71c30 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 07:41:41 -0500 Subject: [PATCH 01/11] Add scXpand --- DESCRIPTION | 13 +- NAMESPACE | 7 + NEWS.md | 35 + R/basilisk.R | 31 + R/immLynx-package.R | 8 + R/runScXpand.R | 1232 ++++++++++++++++++++++++++++++ inst/WORDLIST | 21 + man/immLynx-package.Rd | 13 + man/listScXpandModels.Rd | 34 + man/runScXpand.Rd | 209 +++++ tests/testthat/helper-immLynx.R | 86 +++ tests/testthat/test-runScXpand.R | 624 +++++++++++++++ vignettes/advanced_analysis.Rmd | 78 ++ vignettes/immLynx_vignette.Rmd | 6 + 14 files changed, 2393 insertions(+), 4 deletions(-) create mode 100644 R/runScXpand.R create mode 100644 man/listScXpandModels.Rd create mode 100644 man/runScXpand.Rd create mode 100644 tests/testthat/test-runScXpand.R diff --git a/DESCRIPTION b/DESCRIPTION index 605456e..59a3252 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: immLynx Title: Linking Advanced TCR Python Pipelines and Hugging Face Models in R -Version: 1.1.2 +Version: 1.3.1 Authors@R: c( person(given = "Nick", family = "Borcherding", role = c("aut", "cre"), email = "ncborch@gmail.com", @@ -10,13 +10,14 @@ Description: A comprehensive toolkit that bridges popular Python-based immune R environment. Provides unified interfaces for TCR distance calculations (tcrdist3), sequence generation probability (OLGA), selection inference (soNNia), clustering (clusTCR), protein embeddings (ESM-2), metaclone - discovery (metaclonotypist). Fully compatible with the scRepertoire and immApex ecosystem for single-cell + discovery (metaclonotypist), and clonal expansion prediction from gene + expression (scXpand). Fully compatible with the scRepertoire and immApex ecosystem for single-cell immune repertoire analysis. License: MIT + file LICENSE Encoding: UTF-8 -RoxygenNote: 7.3.3 biocViews: Software, ImmunoOncology, SingleCell, Classification, Annotation, - Sequencing, MotifAnnotation, Clustering, DimensionReduction + Sequencing, MotifAnnotation, Clustering, DimensionReduction, + GeneExpression Depends: R (>= 4.5.0) Imports: @@ -28,15 +29,18 @@ Imports: SingleCellExperiment, stats, SummarizedExperiment, + tools, utils, zellkonverter Suggests: + AnnotationDbi, BiocStyle, data.table, ggplot2, knitr, markdown, Matrix, + org.Hs.eg.db, rmarkdown, scater, scran, @@ -51,3 +55,4 @@ URL: https://github.com/BorchLab/immLynx/ BugReports: https://github.com/BorchLab/immLynx/issues StagedInstall: no Config/testthat/edition: 3 +Config/roxygen2/version: 8.0.0 diff --git a/NAMESPACE b/NAMESPACE index 7f6bb0a..799e9eb 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,12 +5,14 @@ export(exportToScanpy) export(extractTCRdata) export(generateOLGA) export(huggingModel) +export(listScXpandModels) export(proteinEmbeddings) export(runClustTCR) export(runEmbeddings) export(runHLAassociation) export(runMetaclonotypist) export(runOLGA) +export(runScXpand) export(runSoNNia) export(runTCRdist) export(summarizeTCRrepertoire) @@ -19,18 +21,23 @@ export(validateTCRdata) exportClasses(TCR_summary) exportMethods(show) importFrom(SummarizedExperiment,"colData<-") +importFrom(SummarizedExperiment,assay) importFrom(SummarizedExperiment,assayNames) importFrom(SummarizedExperiment,colData) +importFrom(SummarizedExperiment,rowData) importFrom(immApex,getIR) +importFrom(methods,"slot<-") importFrom(methods,is) importFrom(methods,new) importFrom(methods,setClass) importFrom(methods,setMethod) importFrom(methods,show) +importFrom(methods,slot) importFrom(reticulate,array_reshape) importFrom(reticulate,import) importFrom(reticulate,py_to_r) importFrom(stats,median) importFrom(stats,sd) +importFrom(stats,setNames) importFrom(utils,head) importFrom(utils,write.csv) diff --git a/NEWS.md b/NEWS.md index d1aac7f..68f4172 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,38 @@ +# immLynx 1.3.1 + +* Added `runScXpand()`, which predicts T-cell clonal expansion from gene + expression alone using scXpand's pretrained pan-cancer models. Unlike the + other wrapped tools it needs no receptor sequences, so it works on + datasets with no paired TCR sequencing. Inference only; training and + hyperparameter optimization stay in Python. +* Added `listScXpandModels()` to enumerate the available pretrained models + without building the Python environment. +* Added a dedicated `scXpandEnv` basilisk environment (Python 3.11, CPU + PyTorch, scxpand 0.4.6). It cannot share the existing environments + because scxpand needs Python 3.11 and torch 2.5. The first call builds + several gigabytes and downloads the selected model from figshare. +* `runScXpand()` resolves gene identifiers to the Ensembl IDs scXpand's + models are indexed by, optionally mapping symbols through `org.Hs.eg.db`. + Ambiguous symbols (`HLA-DRA` alone maps to eight Ensembl IDs) and + collapsed duplicates are counted and reported rather than resolved + silently, because scXpand zero-fills genes it cannot find. +* When scRepertoire clone calls are present, `runScXpand()` derives + `clone_id_size`, `median_clone_size` and `expansion` per sample using + scXpand's 1.5x-median rule, so a gene-expression-only prediction can be + scored against the observed repertoire. Clone sizes are tabulated fresh + rather than read from `clonalFrequency`, which `combineExpression()` + computes under whatever grouping was in effect. +* Pretrained models are cached under `tools::R_user_dir("immLynx", "cache")`. + scXpand's own default would write a `.scxpand_cache` directory into the + current working directory, so `runScXpand()` downloads the model as an + explicit step with an explicit cache location. +* Worked around an upstream download bug: scXpand's registry points at + `figshare.com/ndownloader/articles/...`, which answers HTTP 202 with an + empty body, and the failed download is cached so retries keep failing. + The same archive on `ndownloader.figshare.com` serves correctly, so + `runScXpand()` rewrites the host. The rewrite becomes a no-op once + upstream fixes its URLs. + # immLynx 1.1.2 * Added `exportToScanpy()` to write a `SingleCellExperiment` or `Seurat` diff --git a/R/basilisk.R b/R/basilisk.R index 5cec0a0..d2a2e67 100644 --- a/R/basilisk.R +++ b/R/basilisk.R @@ -64,3 +64,34 @@ scanpyExportEnv <- basilisk::BasiliskEnvironment( "scirpy" ) ) + +# Separate environment for scXpand (clonal expansion prediction from gene +# expression). Cannot share scanpyExportEnv (python 3.10) or immLynxEnv +# (python 3.9, torch 2.1.2, numpy 1.23): scxpand requires python >= 3.11 +# and torch >= 2.5. +# +# pytorch-cpu is pulled from conda-forge rather than letting pip resolve +# `torch`, because the default PyPI torch wheel on linux-x86_64 bundles +# CUDA and runs to several gigabytes. basilisk's `pip` entries are bare +# specifiers, so there is no way to inject +# --index-url https://download.pytorch.org/whl/cpu. Installing the CPU +# build through conda first satisfies scxpand's torch>=2.5 requirement. +# +# Only scxpand itself is listed under `pip`; it resolves the rest of the +# stack (scanpy, anndata, scirpy, lightgbm, optuna, pooch, shap, ...). +# Enumerating those here would only create version-pin drift. +scXpandEnv <- basilisk::BasiliskEnvironment( + envname = "scXpandEnv", + pkgname = "immLynx", + packages = c( + "python=3.11", + "pytorch-cpu>=2.5", + "numpy", + "pandas", + "scipy", + "h5py" + ), + pip = c( + "scxpand==0.4.6" + ) +) diff --git a/R/immLynx-package.R b/R/immLynx-package.R index 5ce9643..dc0e5e7 100644 --- a/R/immLynx-package.R +++ b/R/immLynx-package.R @@ -34,6 +34,13 @@ #' \item \code{\link{runSoNNia}}: Infer selection with soNNia #' } #' +#' **Expansion Prediction:** +#' \itemize{ +#' \item \code{\link{runScXpand}}: Predict T-cell clonal expansion from gene +#' expression alone with scXpand, no paired TCR sequencing required +#' \item \code{\link{listScXpandModels}}: List the scXpand pretrained models +#' } +#' #' **Utility Functions:** #' \itemize{ #' \item \code{\link{extractTCRdata}}: Extract TCR data from SingleCellExperiment objects @@ -65,6 +72,7 @@ #' \item metaclonotypist - Metaclone discovery #' \item transformers - Hugging Face models #' \item torch - PyTorch for GPU support +#' \item scxpand - Clonal expansion prediction from gene expression #' } #' #' @seealso diff --git a/R/runScXpand.R b/R/runScXpand.R new file mode 100644 index 0000000..82af9fc --- /dev/null +++ b/R/runScXpand.R @@ -0,0 +1,1232 @@ +# Predict T-cell clonal expansion from gene expression using scXpand. +# Public entry points: runScXpand(), listScXpandModels(). +# +# Unlike every other tool wrapped by immLynx, scXpand consumes gene +# expression rather than receptor sequences: it infers whether a T cell +# belongs to an expanded clone without paired TCR sequencing. When the +# object does carry scRepertoire clone calls, those are used to derive +# ground-truth labels so the prediction can be benchmarked. + +#' @keywords internal +# Run a function inside the scXpandEnv basilisk env. Mirrors +# .run_in_scanpy_env in R/exportToScanpy.R but targets the torch/scxpand +# stack, which needs its own python (>= 3.11). +.run_in_scxpand_env <- function(FUN, ...) { + proc <- basilisk::basiliskStart(scXpandEnv) + on.exit(basilisk::basiliskStop(proc)) + basilisk::basiliskRun(proc, FUN, ...) +} + +# Static mirror of scxpand.pretrained.model_registry.PRETRAINED_MODELS. +# Kept in R so listScXpandModels() can answer without building a +# multi-gigabyte conda environment; listScXpandModels(refresh = TRUE) +# queries the installed package and is authoritative. +.SCXPAND_MODELS <- data.frame( + model_name = c("pan_cancer_autoencoder", "pan_cancer_mlp", + "pan_cancer_lightgbm", "pan_cancer_logistic", + "pan_cancer_svm"), + model_type = c("autoencoder", "mlp", "lightgbm", "logistic", "svm"), + version = rep("1.0.0", 5L), + description = c( + "Pan-cancer autoencoder (scXpand default)", + "Pan-cancer multi-layer perceptron", + "Pan-cancer LightGBM gradient boosting", + "Pan-cancer logistic regression", + "Pan-cancer support vector machine" + ), + stringsAsFactors = FALSE +) + +# =========================================================================== +# Gene identifier resolution +# =========================================================================== + +#' @keywords internal +.isEnsembl <- function(x) { + out <- grepl("^ENSG[0-9]{11}$", x) + out[is.na(x)] <- FALSE + out +} + +#' @keywords internal +# Strip trailing version suffixes, but only from strings that already look +# like versioned Ensembl gene IDs. A blanket sub("\\.[0-9]+$", "", x) would +# mangle legitimate symbols such as MARCH1.2 or 7SK.2. +.stripEnsemblVersion <- function(x) { + versioned <- grepl("^ENSG[0-9]{11}\\.[0-9]+$", x) + versioned[is.na(x)] <- FALSE + x[versioned] <- sub("\\.[0-9]+$", "", x[versioned]) + x +} + +#' @keywords internal +# Map gene symbols to Ensembl gene IDs via org.Hs.eg.db, trying SYMBOL +# first and falling back to ALIAS for unmatched keys. +# +# One symbol routinely maps to several ENSG IDs (HLA-DRA -> 8, HSPA1A -> 5, +# alt haplotypes and patch scaffolds). scXpand silently zero-fills any gene +# it cannot find in its panel, so an unlucky pick deletes that gene rather +# than erroring. Every ambiguous case is therefore counted and reported. +.mapSymbolsToEnsembl <- function(symbols, multi_map = c("first", "expand", + "drop")) { + multi_map <- match.arg(multi_map) + + empty <- list(map = data.frame(SYMBOL = character(0), + ENSEMBL = character(0), + stringsAsFactors = FALSE), + n_query = 0L, n_symbol = 0L, n_alias = 0L, + n_ambiguous = 0L, n_unmapped = 0L) + + uq <- unique(symbols[!is.na(symbols) & nzchar(symbols)]) + if (!length(uq)) return(empty) + + if (!requireNamespace("org.Hs.eg.db", quietly = TRUE) || + !requireNamespace("AnnotationDbi", quietly = TRUE)) { + stop("Symbol to Ensembl mapping requires org.Hs.eg.db and AnnotationDbi. ", + "Install with: BiocManager::install(c('org.Hs.eg.db', ", + "'AnnotationDbi')), or pass map_symbols = 'never' and supply ", + "Ensembl IDs via gene_ids.", call. = FALSE) + } + + db <- org.Hs.eg.db::org.Hs.eg.db + + .lookup <- function(keys, keytype) { + if (!length(keys)) return(NULL) + res <- tryCatch( + suppressWarnings(suppressMessages( + AnnotationDbi::select(db, keys = keys, keytype = keytype, + columns = "ENSEMBL"))), + error = function(e) NULL) + if (is.null(res) || !nrow(res)) return(NULL) + res <- res[!is.na(res$ENSEMBL), , drop = FALSE] + if (!nrow(res)) return(NULL) + data.frame(SYMBOL = as.character(res[[keytype]]), + ENSEMBL = as.character(res$ENSEMBL), + stringsAsFactors = FALSE) + } + + m <- .lookup(uq, "SYMBOL") + n_symbol <- if (is.null(m)) 0L else length(unique(m$SYMBOL)) + + unmatched <- setdiff(uq, if (is.null(m)) character(0) else m$SYMBOL) + a <- .lookup(unmatched, "ALIAS") + n_alias <- if (is.null(a)) 0L else length(unique(a$SYMBOL)) + + m <- rbind(m, a) + if (is.null(m) || !nrow(m)) { + empty$n_query <- length(uq) + empty$n_unmapped <- length(uq) + return(empty) + } + + m <- m[!duplicated(m), , drop = FALSE] + # Lexicographic order makes multi_map = "first" reproducible across + # machines and org.Hs.eg.db releases; AnnotationDbi's own + # multiVals = "first" follows database row order, which is not stable. + m <- m[order(m$SYMBOL, m$ENSEMBL), , drop = FALSE] + + per_symbol <- table(m$SYMBOL) + ambiguous <- names(per_symbol)[per_symbol > 1L] + + m <- switch(multi_map, + first = m[!duplicated(m$SYMBOL), , drop = FALSE], + drop = m[!(m$SYMBOL %in% ambiguous), , drop = FALSE], + expand = m) + + list(map = m, + n_query = length(uq), + n_symbol = n_symbol, + n_alias = n_alias, + n_ambiguous = length(ambiguous), + n_unmapped = length(setdiff(uq, m$SYMBOL))) +} + +#' @keywords internal +# Turn whatever identifiers the object carries into Ensembl gene IDs, and +# work out which rows of the object survive. +# +# Returns: +# ids final rownames for the staged matrix +# keep row indices into the input object; may repeat when +# multi_map = "expand", and is longer than `ids` when +# collapse_groups is non-NULL +# collapse_groups NULL, or a factor over `keep` whose levels are `ids`, +# meaning those rows must be summed together +# report named list summarising the mapping, surfaced to the +# user and stored in metadata()$scXpand$gene_mapping +.resolveGeneIDs <- function(sce, gene_ids = "rownames", + map_symbols = c("auto", "always", "never"), + ensembl_min_frac = 0.5, + multi_map = c("first", "expand", "drop"), + collapse = c("sum", "first", "drop"), + verbose = TRUE) { + map_symbols <- match.arg(map_symbols) + multi_map <- match.arg(multi_map) + collapse <- match.arg(collapse) + + n_in <- nrow(sce) + + # --- Extract the raw identifiers ------------------------------------- + if (is.character(gene_ids) && length(gene_ids) == 1L && + identical(gene_ids, "rownames")) { + raw <- rownames(sce) + if (is.null(raw)) { + stop("gene_ids = 'rownames' but the object has no rownames; pass a ", + "rowData column name or a character vector instead.", + call. = FALSE) + } + src <- "rownames" + } else if (is.character(gene_ids) && length(gene_ids) == 1L) { + rd <- SummarizedExperiment::rowData(sce) + if (!(gene_ids %in% colnames(rd))) { + stop("gene_ids must be \"rownames\", the name of a rowData column, ", + "or a character vector of length nrow(input). ", + "'", gene_ids, "' is not a rowData column", + if (ncol(rd)) paste0(" (available: ", + paste(colnames(rd), collapse = ", "), ")") + else " (the object has no rowData columns)", ".", + call. = FALSE) + } + raw <- as.character(rd[[gene_ids]]) + src <- paste0("rowData$", gene_ids) + } else if (is.character(gene_ids) && length(gene_ids) == n_in) { + raw <- as.character(gene_ids) + src <- "user-supplied vector" + } else { + stop("gene_ids must be \"rownames\", the name of a rowData column, or ", + "a character vector of length nrow(input) (", n_in, "); got ", + class(gene_ids)[1], " of length ", length(gene_ids), ".", + call. = FALSE) + } + + raw <- trimws(raw) + raw[!nzchar(raw)] <- NA_character_ + raw <- .stripEnsemblVersion(raw) + + is_ens <- .isEnsembl(raw) + frac <- mean(is_ens) + + # --- Decide whether to map ------------------------------------------- + do_map <- switch(map_symbols, + always = TRUE, + never = FALSE, + auto = frac < ensembl_min_frac) + + if (map_symbols == "never" && frac < ensembl_min_frac) { + offenders <- utils::head(unique(raw[!is_ens & !is.na(raw)]), 5L) + stop("Only ", sprintf("%.1f%%", 100 * frac), " of gene identifiers from ", + src, " look like Ensembl gene IDs (need at least ", + sprintf("%.1f%%", 100 * ensembl_min_frac), "). ", + "Examples: ", paste(offenders, collapse = ", "), ". ", + "scXpand's pretrained models are indexed by Ensembl ID. ", + "Use map_symbols = 'auto' to map symbols via org.Hs.eg.db, or ", + "supply Ensembl IDs through gene_ids.", call. = FALSE) + } + + # --- Ensembl rows pass straight through ------------------------------- + keep <- which(is_ens) + ids <- raw[keep] + + mapinfo <- NULL + if (do_map) { + mapinfo <- .mapSymbolsToEnsembl(raw[!is_ens], multi_map = multi_map) + m <- mapinfo$map + if (nrow(m)) { + by_sym <- split(m$ENSEMBL, m$SYMBOL) + n_hits <- lengths(by_sym)[raw] + n_hits[is.na(n_hits)] <- 0L + n_hits[is_ens] <- 0L # already handled above + hit_rows <- which(n_hits > 0L) + keep <- c(keep, rep(hit_rows, times = n_hits[hit_rows])) + ids <- c(ids, unlist(by_sym[raw[hit_rows]], use.names = FALSE)) + } + if (mapinfo$n_query > 0L && + (mapinfo$n_query - mapinfo$n_unmapped) / mapinfo$n_query < 0.25) { + warning("Fewer than 25% of gene symbols mapped to Ensembl IDs. ", + "Is this a human dataset, and are the identifiers really ", + "gene symbols?", call. = FALSE) + } + } + + ord <- order(keep) + keep <- keep[ord] + ids <- ids[ord] + + if (!length(ids)) { + stop("No gene identifiers could be resolved to Ensembl IDs from ", src, + ". scXpand cannot run without them.", call. = FALSE) + } + + # --- Post-mapping threshold re-check ---------------------------------- + post_frac <- length(unique(ids)) / max(1L, length(unique(raw[!is.na(raw)]))) + if (do_map && map_symbols == "auto" && post_frac < ensembl_min_frac) { + stop("After symbol mapping only ", sprintf("%.1f%%", 100 * post_frac), + " of input features resolved to Ensembl IDs (need at least ", + sprintf("%.1f%%", 100 * ensembl_min_frac), "). ", + "Lower ensembl_min_frac if this is expected (targeted panel), or ", + "supply Ensembl IDs through gene_ids.", call. = FALSE) + } + + # --- Collapse duplicate Ensembl IDs ----------------------------------- + n_precollapse <- length(ids) + dup_ids <- unique(ids[duplicated(ids)]) + collapse_groups <- NULL + + if (length(dup_ids)) { + if (collapse == "sum") { + collapse_groups <- factor(ids) + ids <- levels(collapse_groups) + } else if (collapse == "first") { + first <- !duplicated(ids) + keep <- keep[first] + ids <- ids[first] + } else { + okay <- !(ids %in% dup_ids) + keep <- keep[okay] + ids <- ids[okay] + } + } + + report <- list( + source = src, + n_input = n_in, + n_already_ensembl = sum(is_ens), + frac_ensembl_in = frac, + mapped = do_map, + n_mapped_symbol = if (is.null(mapinfo)) 0L else mapinfo$n_symbol, + n_mapped_alias = if (is.null(mapinfo)) 0L else mapinfo$n_alias, + n_unmapped = if (is.null(mapinfo)) 0L else mapinfo$n_unmapped, + n_ambiguous = if (is.null(mapinfo)) 0L else mapinfo$n_ambiguous, + multi_map = if (do_map) multi_map else NA_character_, + collapse = if (length(dup_ids)) collapse else NA_character_, + n_collapsed = length(dup_ids), + n_precollapse = n_precollapse, + n_output = length(ids) + ) + + if (verbose) message(.formatGeneMapping(report)) + + list(ids = ids, keep = keep, collapse_groups = collapse_groups, + report = report) +} + +#' @keywords internal +.formatGeneMapping <- function(r) { + lines <- c( + sprintf("scXpand gene ID resolution (%s):", r$source), + sprintf(" input features: %d", r$n_input), + sprintf(" already Ensembl: %d (%.1f%%)", r$n_already_ensembl, + 100 * r$frac_ensembl_in)) + if (isTRUE(r$mapped)) { + lines <- c(lines, + sprintf(" mapped via SYMBOL: %d", r$n_mapped_symbol), + sprintf(" mapped via ALIAS: %d", r$n_mapped_alias), + sprintf(" unmapped (dropped): %d", r$n_unmapped), + sprintf(" ambiguous 1:many: %d symbols (multi_map = '%s')", + r$n_ambiguous, r$multi_map)) + } + if (r$n_collapsed > 0L) { + lines <- c(lines, + sprintf(" collapsed many:1: %d Ensembl IDs received >1 row (collapse = '%s')", + r$n_collapsed, r$collapse)) + } + c(paste(c(lines, sprintf(" features written: %d", r$n_output)), + collapse = "\n")) +} + +# =========================================================================== +# Ground-truth expansion labels from scRepertoire clone calls +# =========================================================================== + +#' @keywords internal +# Derive scXpand's ground-truth obs fields from clone identity. +# +# scXpand defines a cell as expanded when its clone_id_size exceeds +# 1.5 x median_clone_size for that sample. Clone sizes are tabulated fresh +# here rather than read from scRepertoire's clonalFrequency column, because +# combineExpression() computes that under whatever group.by was in effect +# at combine time (or globally) and it is not guaranteed to be a per-sample +# count. +# +# Returns NULL when the object carries no usable clone data. +.deriveExpansionLabels <- function(sce, clone_col = NULL, sample_col = NULL, + median_basis = c("clone", "cell"), + verbose = TRUE) { + median_basis <- match.arg(median_basis) + cd <- SummarizedExperiment::colData(sce) + + if (is.null(clone_col)) { + candidates <- intersect(c("CTstrict", "CTaa"), colnames(cd)) + if (!length(candidates)) return(NULL) + clone_col <- candidates[1L] + } else if (!(clone_col %in% colnames(cd))) { + stop("clone_col '", clone_col, "' not found in colData.", call. = FALSE) + } + + cid <- as.character(cd[[clone_col]]) + # scRepertoire writes literal "NA" tokens into the CT* columns for cells + # with a missing chain; those are absent clone calls, not clone names. + cid[!nzchar(cid) | cid %in% c("NA", "None", "NA_NA")] <- NA_character_ + has <- !is.na(cid) + if (!any(has)) return(NULL) + + if (is.null(sample_col)) { + warning("sample_col is NULL; treating all ", ncol(sce), + " cells as a single sample. scXpand's expansion rule is defined ", + "per sample -- pass sample_col = '' (e.g. ", + "'orig.ident' or 'Patient') for correct labels.", call. = FALSE) + smp <- rep("__all__", ncol(sce)) + } else { + if (!(sample_col %in% colnames(cd))) { + stop("sample_col '", sample_col, "' not found in colData.", + call. = FALSE) + } + smp <- as.character(cd[[sample_col]]) + if (anyNA(smp)) { + warning(sum(is.na(smp)), " cells have a missing ", sample_col, + " value; they are pooled into their own stratum.", + call. = FALSE) + smp[is.na(smp)] <- "__NA__" + } + } + + size <- rep(NA_real_, ncol(sce)) + med <- rep(NA_real_, ncol(sce)) + samples <- unique(smp[has]) + median_by_sample <- stats::setNames(rep(NA_real_, length(samples)), samples) + + for (s in samples) { + in_s <- has & smp == s + tab <- table(cid[in_s]) + sz <- as.numeric(tab[cid[in_s]]) + size[in_s] <- sz + m <- if (median_basis == "clone") { + stats::median(as.numeric(tab)) # median over unique clones + } else { + stats::median(sz) # median over cells + } + med[in_s] <- m + median_by_sample[[s]] <- m + } + + expansion <- ifelse(is.na(size), NA_character_, + ifelse(size > 1.5 * med, "expanded", "non-expanded")) + + if (verbose) { + message("Derived expansion labels from '", clone_col, "'", + if (is.null(sample_col)) "" else paste0(" within '", sample_col, "'"), + " (median_basis = '", median_basis, "'): ", + sum(has), "/", ncol(sce), " cells labelled, ", + sum(expansion == "expanded", na.rm = TRUE), " expanded.") + } + + list(clone_id_size = size, + median_clone_size = med, + expansion = expansion, + clone_col = clone_col, + sample_col = sample_col, + median_basis = median_basis, + median_by_sample = median_by_sample, + n_labelled = sum(has), + n_clones = length(unique(cid[has])), + n_samples = length(samples)) +} + +#' @keywords internal +# Rank-based AUROC (Mann-Whitney U). Computed in R so the reported number +# never depends on how scXpand's evaluator handles NAs or ties. +.aurocR <- function(prob, label) { + ok <- !is.na(prob) & !is.na(label) + prob <- prob[ok] + label <- as.logical(label[ok]) + n1 <- sum(label) + n0 <- length(label) - n1 + if (n1 == 0L || n0 == 0L) return(NA_real_) + r <- rank(prob) + (sum(r[label]) - n1 * (n1 + 1) / 2) / (n1 * n0) +} + +# =========================================================================== +# H5AD staging +# =========================================================================== + +#' @keywords internal +.validateCountsAssay <- function(X, assay) { + normalized <- c("logcounts", "lognorm", "lognormcounts", "data", + "normcounts", "scaledata", "scale.data") + if (tolower(assay) %in% tolower(normalized)) { + stop("assay '", assay, "' looks like normalized or log-transformed data. ", + "scXpand requires raw UMI counts; pass assay = 'counts'.", + call. = FALSE) + } + + v <- if (methods::is(X, "dgCMatrix")) X@x else as.numeric(X) + if (!length(v) || all(v == 0)) { + stop("assay '", assay, "' contains no non-zero counts.", call. = FALSE) + } + if (min(v) < 0) { + stop("assay '", assay, "' contains negative values; scXpand requires ", + "raw UMI counts.", call. = FALSE) + } + # Integrality on a subsample: full checks on millions of cells are wasteful + # and a single non-integer is enough to reject. + if (length(v) > 1e5L) v <- v[seq_len(1e5L)] + if (any(abs(v - round(v)) > 1e-8)) { + stop("assay '", assay, "' contains non-integer values (it looks ", + "normalized or log-transformed). scXpand requires raw UMI counts; ", + "pass assay = 'counts'.", call. = FALSE) + } + invisible(TRUE) +} + +#' @keywords internal +# Sum rows sharing an Ensembl ID. Raw UMI counts are additive, so summing +# preserves the count semantics scXpand expects. +.collapseRowsSum <- function(X, groups) { + groups <- droplevels(as.factor(groups)) + if (requireNamespace("Matrix", quietly = TRUE)) { + g <- as.integer(groups) + agg <- Matrix::sparseMatrix(i = g, j = seq_along(g), x = 1, + dims = c(nlevels(groups), length(g))) + out <- agg %*% X + rownames(out) <- levels(groups) + colnames(out) <- colnames(X) + return(out) + } + out <- rowsum(as.matrix(X), group = groups, reorder = TRUE) + rownames(out) <- levels(groups) + out +} + +#' @keywords internal +# Build a minimal SCE holding only what scXpand reads (X and obs) and write +# it to H5AD. zellkonverter runs in its own basilisk env, so this happens on +# the R side and only the file path crosses into scXpandEnv. +.stageScXpandH5AD <- function(sce, res_ids, labels, obs_columns, assay, + h5ad_path, verbose = TRUE) { + + X <- SummarizedExperiment::assay(sce, assay) + .validateCountsAssay(X, assay) + + X <- X[res_ids$keep, , drop = FALSE] + if (!is.null(res_ids$collapse_groups)) { + X <- .collapseRowsSum(X, res_ids$collapse_groups) + } else { + rownames(X) <- res_ids$ids + } + colnames(X) <- colnames(sce) + + stage <- SingleCellExperiment::SingleCellExperiment( + assays = list(counts = X), + colData = SummarizedExperiment::colData(sce)) + + # drop_ct = TRUE: the CT* columns are long strings scXpand never reads, + # and their content is already distilled into clone_id_size. + stage <- .pruneColData(stage, obs_columns = obs_columns, drop_ct = TRUE) + + # Attached after pruning so obs_columns can never remove them. Names come + # from scxpand/data_util/constants.py. All three are written together or + # not at all: scXpand reads them as a set during evaluation, and a partial + # set is a plausible way to make its evaluator raise. + if (!is.null(labels) && !is.null(labels$expansion)) { + SummarizedExperiment::colData(stage)$clone_id_size <- + labels$clone_id_size + SummarizedExperiment::colData(stage)$median_clone_size <- + labels$median_clone_size + SummarizedExperiment::colData(stage)$expansion <- labels$expansion + } + + if (verbose) message("Writing staged H5AD: ", h5ad_path) + zellkonverter::writeH5AD(stage, file = h5ad_path, X_name = "counts") + + invisible(h5ad_path) +} + +# =========================================================================== +# Result handling +# =========================================================================== + +#' @keywords internal +# Match predictions back to R cells by barcode rather than trusting +# positional alignment through writeH5AD -> h5py -> AnnData -> DataLoader. +.alignPredictions <- function(preds, obs_names, expected_names) { + if (length(preds) != length(obs_names)) { + stop("scXpand returned ", length(preds), " predictions for ", + length(obs_names), " cells in the staged H5AD.", call. = FALSE) + } + if (length(preds) != length(expected_names) || + !setequal(obs_names, expected_names)) { + warning("H5AD obs_names do not match the staged cell names; falling ", + "back to positional alignment.", call. = FALSE) + names(preds) <- expected_names[seq_along(preds)] + return(preds) + } + names(preds) <- obs_names + preds[expected_names] +} + +#' @keywords internal +# Flatten scXpand's nested metrics dict to a named numeric vector with +# dot-joined keys. Non-scalar leaves (arrays, confusion matrices) are +# dropped here; the unflattened list is retained as metrics_raw. +.flattenMetrics <- function(x, prefix = "") { + if (is.null(x) || !length(x)) return(stats::setNames(numeric(0), character(0))) + nms <- names(x) + if (is.null(nms)) return(stats::setNames(numeric(0), character(0))) + out <- stats::setNames(numeric(0), character(0)) + for (i in seq_along(x)) { + nm <- nms[i] + if (is.null(nm) || !nzchar(nm)) next + v <- x[[i]] + key <- if (nzchar(prefix)) paste(prefix, nm, sep = ".") else nm + if (is.list(v)) { + out <- c(out, .flattenMetrics(v, key)) + } else if (length(v) == 1L && (is.numeric(v) || is.logical(v)) && + !is.na(v) && is.finite(as.numeric(v))) { + out <- c(out, stats::setNames(as.numeric(v), key)) + } + } + out +} + +#' @keywords internal +# Fill a per-cell column, leaving NA for cells that were not scored. +.writeCellColumn <- function(obj, col_name, values, cell_names) { + col_vec <- rep(NA, ncol(obj)) + names(col_vec) <- colnames(obj) + col_vec[cell_names] <- values + names(col_vec) <- NULL + if (methods::is(obj, "SingleCellExperiment")) { + SummarizedExperiment::colData(obj)[[col_name]] <- col_vec + } else { + obj[[col_name]] <- col_vec + } + obj +} + +#' @keywords internal +# Stash the run summary. SingleCellExperiment has metadata(); Seurat keeps +# the equivalent in the misc slot. +.writeObjMetadata <- function(obj, key, value) { + if (methods::is(obj, "SingleCellExperiment")) { + md <- S4Vectors::metadata(obj) + md[[key]] <- value + S4Vectors::metadata(obj) <- md + } else { + misc <- methods::slot(obj, "misc") + if (!is.list(misc)) misc <- list() + misc[[key]] <- value + methods::slot(obj, "misc") <- misc + } + obj +} + +# =========================================================================== +# Exported functions +# =========================================================================== + +#' List available scXpand pretrained models +#' +#' @description +#' Returns the pan-cancer models published with scXpand. By default this +#' reads a static table shipped with immLynx so it never triggers a +#' multi-gigabyte Python environment build. Set \code{refresh = TRUE} to +#' query the installed \code{scxpand} package, which is authoritative if +#' upstream has added models. +#' +#' @param refresh Logical. Query the installed \code{scxpand} package rather +#' than the static table. Requires the scXpand basilisk environment and +#' will build it on first use. Default \code{FALSE}. +#' +#' @return A \code{data.frame} with columns \code{model_name}, +#' \code{model_type}, \code{version} and \code{description}. When +#' \code{refresh = TRUE}, \code{model_type} and \code{description} may be +#' \code{NA} for models absent from the static table. +#' +#' @export +#' +#' @examples +#' listScXpandModels() +#' +#' \donttest{ +#' # Authoritative, but builds the scXpand environment on first call. +#' listScXpandModels(refresh = TRUE) +#' } +listScXpandModels <- function(refresh = FALSE) { + if (!isTRUE(refresh)) return(.SCXPAND_MODELS) + + live <- .run_in_scxpand_env(function() { + sx <- reticulate::import("scxpand", convert = FALSE) + bi <- reticulate::import_builtins(convert = FALSE) + reg <- sx$PRETRAINED_MODELS + nms <- as.character(reticulate::py_to_r(bi$list(reg$keys()))) + vers <- vapply(nms, function(n) tryCatch( + as.character(reticulate::py_to_r( + reticulate::py_get_attr(reg[[n]], "version"))), + error = function(e) NA_character_), character(1)) + data.frame(model_name = nms, version = unname(vers), + stringsAsFactors = FALSE) + }) + + idx <- match(live$model_name, .SCXPAND_MODELS$model_name) + data.frame(model_name = live$model_name, + model_type = .SCXPAND_MODELS$model_type[idx], + version = live$version, + description = .SCXPAND_MODELS$description[idx], + stringsAsFactors = FALSE) +} + +#' Predict T-cell clonal expansion from gene expression with scXpand +#' +#' @description +#' Runs inference with a scXpand pretrained pan-cancer model to estimate, +#' for each cell, the probability that it belongs to an expanded T-cell +#' clone. scXpand uses gene expression only, so this works on datasets with +#' no paired TCR sequencing. When the object does carry scRepertoire clone +#' calls, \code{runScXpand} also derives ground-truth expansion labels so +#' the prediction can be benchmarked against the observed repertoire. +#' +#' @details +#' \strong{Input requirements.} scXpand's pretrained models were trained on +#' raw UMI counts from T cells indexed by Ensembl gene ID. All three matter: +#' \itemize{ +#' \item \emph{Raw counts.} Normalized or log-transformed values are +#' rejected. Use \code{assay = "counts"}. +#' \item \emph{Ensembl IDs.} Most Seurat and scRepertoire objects carry +#' gene symbols; see \code{gene_ids} and \code{map_symbols}. +#' \item \emph{T cells only.} Filter to T cells before calling. There is +#' no reliable programmatic check, and predictions on non-T cells are +#' meaningless. +#' } +#' Genes the model expects but cannot find are zero-filled by scXpand, and +#' extra genes are ignored, so a partial overlap degrades quietly rather +#' than erroring. That is why the gene mapping summary is worth reading. +#' +#' \strong{Symbol mapping.} When identifiers are symbols, they are mapped +#' through \code{org.Hs.eg.db} (SYMBOL first, then ALIAS). One symbol often +#' maps to several Ensembl IDs because of alternate haplotypes: +#' \code{HLA-DRA} maps to eight. \code{multi_map = "first"} takes the +#' lexicographically first, which is reproducible but may pick a scaffold +#' the model does not know; \code{multi_map = "expand"} emits every +#' candidate, which maximizes overlap with the model's panel at the cost of +#' a larger staged matrix. Several input rows can also land on one Ensembl +#' ID, in which case \code{collapse} decides what happens (summing raw +#' counts is the default and preserves UMI semantics). +#' +#' \strong{Expansion labels.} scXpand calls a cell expanded when its clone +#' size exceeds 1.5 times the median clone size for its sample. Clone sizes +#' are tabulated fresh from \code{clone_col} within \code{sample_col}; +#' scRepertoire's \code{clonalFrequency} column is deliberately not used, +#' because \code{combineExpression()} computes it under whatever +#' \code{group.by} was in effect and it is not necessarily a per-sample +#' count. "Median clone size" is genuinely ambiguous, so +#' \code{median_basis} selects the reading: \code{"clone"} takes the median +#' over unique clones, which in typical 10x data is 1 (most clones are +#' singletons) and therefore reduces the rule to \emph{clone size at least +#' 2}; \code{"cell"} takes the median over cells, which is dominated by +#' large clones and is much stricter. +#' +#' \strong{Cost.} The first call builds a dedicated basilisk environment +#' (Python 3.11, PyTorch CPU, scxpand) that occupies several gigabytes, and +#' downloads the selected model (roughly 35 MB) from figshare into +#' \code{cache_dir}. Both are reused afterwards. If a download does fail, +#' clear \code{cache_dir} before retrying: the downloader caches the failed +#' response and will otherwise keep reusing it. +#' +#' @param input A \code{SingleCellExperiment} or \code{Seurat} object. +#' SingleCellExperiment is the native format and is used without +#' conversion; Seurat objects are converted in and returned as Seurat. +#' @param model_name Pretrained model to use. See +#' \code{\link{listScXpandModels}}. Default +#' \code{"pan_cancer_autoencoder"}. +#' @param assay Assay holding raw UMI counts. Default \code{"counts"}. +#' @param gene_ids Where to find gene identifiers: \code{"rownames"} +#' (default), the name of a \code{rowData} column, or a character vector +#' of length \code{nrow(input)}. +#' @param map_symbols One of \code{"auto"} (map symbols to Ensembl IDs only +#' if too few identifiers already look like Ensembl IDs), \code{"always"}, +#' or \code{"never"}. Mapping requires \code{org.Hs.eg.db} and +#' \code{AnnotationDbi}. +#' @param ensembl_min_frac Minimum fraction of identifiers that must resolve +#' to Ensembl IDs. Default \code{0.5}. Lower it for targeted panels. +#' @param multi_map How to resolve one symbol mapping to several Ensembl +#' IDs: \code{"first"} (default, lexicographically first), +#' \code{"expand"} (emit all), or \code{"drop"}. +#' @param collapse How to resolve several rows mapping to one Ensembl ID: +#' \code{"sum"} (default), \code{"first"}, or \code{"drop"}. +#' @param derive_labels Logical. Derive ground-truth expansion labels from +#' clone data when present. Default \code{TRUE}. +#' @param clone_col Column holding clone identity. Default \code{NULL}: +#' use \code{"CTstrict"} if present, else \code{"CTaa"}. +#' @param sample_col Column identifying the sample or patient that clone +#' sizes are tabulated within. Default \code{NULL}, which pools all cells +#' into one sample and warns. +#' @param median_basis \code{"clone"} (default) or \code{"cell"}. See +#' Details. +#' @param label_cells \code{"clonal"} (default) scores only cells with a +#' clone call, which guarantees the ground-truth column has no missing +#' values; \code{"all"} scores every cell. +#' @param threshold Probability cutoff for the binary call. Default +#' \code{0.5}. +#' @param obs_columns Character vector of \code{colData} columns to carry +#' into the staged H5AD. Default \code{NULL}: all of them, minus +#' list-columns and the scRepertoire \code{CT*} columns. +#' @param batch_size Inference batch size. Default \code{1024}. +#' @param num_workers DataLoader workers. Default \code{0}, because worker +#' processes forking out of embedded Python inside a basilisk child can +#' hang on macOS. Raise it on Linux for speed. +#' @param work_dir Directory for the staged H5AD and scXpand's own output. +#' Default \code{NULL}: a temporary directory, removed on exit. +#' @param cache_dir Directory for downloaded pretrained models. Default +#' \code{NULL}: a \code{scxpand} subdirectory of +#' \code{tools::R_user_dir("immLynx", "cache")}. Passing this explicitly +#' matters because scXpand's own default would write a +#' \code{.scxpand_cache} directory into the current working directory. +#' @param keep_files Logical. Keep the staged H5AD and scXpand outputs. +#' Default \code{FALSE}. +#' @param overwrite Logical. Replace an existing \code{adata.h5ad} in an +#' explicitly supplied \code{work_dir}. Default \code{FALSE}. +#' @param column_prefix Prefix for the columns written back. Default +#' \code{"scXpand"}. +#' @param return_object Logical. Return the input object with new columns +#' (default), or a \code{data.frame} of per-cell results. +#' @param verbose Logical. Emit progress messages. Default \code{TRUE}. +#' +#' @return If \code{return_object = TRUE}, the input object (same class) +#' with \code{_expansion_prob} and \code{_expansion_pred} +#' added to cell metadata, plus \code{_clone_id_size}, +#' \code{_median_clone_size} and \code{_expansion_truth} +#' when labels were derived. A run summary is stored under +#' \code{metadata(x)$scXpand} for SingleCellExperiment or +#' \code{x@misc$scXpand} for Seurat, holding the model info, the gene +#' mapping report, scXpand's metrics and an independently computed AUROC. +#' If \code{return_object = FALSE}, a \code{data.frame} with one row per +#' scored cell and the run summary attached as +#' \code{attr(x, "scXpand")}. +#' +#' @export +#' @importFrom methods is slot slot<- +#' @importFrom SummarizedExperiment colData colData<- rowData assay assayNames +#' @importFrom stats median setNames +#' @importFrom utils head +#' +#' @examples +#' data(immLynx_example) +#' +#' # The label derivation is pure R and runs without Python. Note that the +#' # derived clone sizes are tabulated per sample and so need not match the +#' # clonalFrequency column scRepertoire wrote. +#' labs <- immLynx:::.deriveExpansionLabels( +#' immLynx_example, +#' sample_col = "Patient" +#' ) +#' table(labs$expansion, useNA = "ifany") +#' +#' \donttest{ +#' # Full inference. The example object carries gene symbols and is not +#' # filtered to T cells, so this demonstrates the mechanics rather than a +#' # meaningful biological result. +#' sce <- runScXpand( +#' immLynx_example, +#' sample_col = "Patient", +#' assay = "counts" +#' ) +#' summary(SummarizedExperiment::colData(sce)$scXpand_expansion_prob) +#' S4Vectors::metadata(sce)$scXpand$auroc_immLynx +#' } +runScXpand <- function(input, + model_name = "pan_cancer_autoencoder", + assay = "counts", + gene_ids = "rownames", + map_symbols = c("auto", "always", "never"), + ensembl_min_frac = 0.5, + multi_map = c("first", "expand", "drop"), + collapse = c("sum", "first", "drop"), + derive_labels = TRUE, + clone_col = NULL, + sample_col = NULL, + median_basis = c("clone", "cell"), + label_cells = c("clonal", "all"), + threshold = 0.5, + obs_columns = NULL, + batch_size = 1024L, + num_workers = 0L, + work_dir = NULL, + cache_dir = NULL, + keep_files = FALSE, + overwrite = FALSE, + column_prefix = "scXpand", + return_object = TRUE, + verbose = TRUE) { + + map_symbols <- match.arg(map_symbols) + multi_map <- match.arg(multi_map) + collapse <- match.arg(collapse) + median_basis <- match.arg(median_basis) + label_cells <- match.arg(label_cells) + + # --- Input dispatch --------------------------------------------------- + # SingleCellExperiment passes through untouched; Seurat is converted in + # here and the results are written back onto the original object below. + sce <- .coerceToSCE(input) + + if (ncol(sce) == 0L) { + stop("input has no cells.", call. = FALSE) + } + if (is.null(colnames(sce))) { + stop("input must have cell names (colnames); predictions are matched ", + "back by barcode.", call. = FALSE) + } + if (anyDuplicated(colnames(sce))) { + stop("input has duplicated cell names; predictions cannot be matched ", + "back by barcode. Make colnames unique first.", call. = FALSE) + } + + # --- Scalar argument validation --------------------------------------- + if (!is.character(model_name) || length(model_name) != 1L || + is.na(model_name) || !nzchar(model_name)) { + stop("model_name must be a single non-empty string.", call. = FALSE) + } + if (!(model_name %in% .SCXPAND_MODELS$model_name)) { + warning("Unknown model_name '", model_name, "'; known models: ", + paste(.SCXPAND_MODELS$model_name, collapse = ", "), + ". Passing it through to scXpand anyway.", call. = FALSE) + } + if (!is.character(assay) || length(assay) != 1L || is.na(assay)) { + stop("assay must be a single string.", call. = FALSE) + } + if (!(assay %in% SummarizedExperiment::assayNames(sce))) { + stop("assay '", assay, "' not found. Available: ", + paste(SummarizedExperiment::assayNames(sce), collapse = ", "), ".", + call. = FALSE) + } + if (!is.numeric(threshold) || length(threshold) != 1L || is.na(threshold) || + threshold < 0 || threshold > 1) { + stop("threshold must be a single number in [0, 1].", call. = FALSE) + } + if (!is.numeric(ensembl_min_frac) || length(ensembl_min_frac) != 1L || + is.na(ensembl_min_frac) || ensembl_min_frac < 0 || + ensembl_min_frac > 1) { + stop("ensembl_min_frac must be a single number in [0, 1].", call. = FALSE) + } + .checkCount <- function(x, nm) { + if (!is.numeric(x) || length(x) != 1L || is.na(x) || x < 0 || + abs(x - round(x)) > 1e-8) { + stop(nm, " must be a single non-negative whole number.", call. = FALSE) + } + as.integer(round(x)) + } + batch_size <- .checkCount(batch_size, "batch_size") + num_workers <- .checkCount(num_workers, "num_workers") + if (batch_size < 1L) { + stop("batch_size must be at least 1.", call. = FALSE) + } + if (!is.character(column_prefix) || length(column_prefix) != 1L || + is.na(column_prefix)) { + stop("column_prefix must be a single string.", call. = FALSE) + } + + cd_names <- colnames(SummarizedExperiment::colData(sce)) + if (!is.null(clone_col) && + (!is.character(clone_col) || length(clone_col) != 1L || + !(clone_col %in% cd_names))) { + stop("clone_col '", paste(clone_col, collapse = ", "), + "' not found in colData.", call. = FALSE) + } + if (!is.null(sample_col) && + (!is.character(sample_col) || length(sample_col) != 1L || + !(sample_col %in% cd_names))) { + stop("sample_col '", paste(sample_col, collapse = ", "), + "' not found in colData.", call. = FALSE) + } + + # --- Work directory --------------------------------------------------- + transient <- is.null(work_dir) + if (transient) { + work_dir <- tempfile("scxpand_") + } else if (!is.character(work_dir) || length(work_dir) != 1L || + is.na(work_dir) || !nzchar(work_dir)) { + stop("work_dir must be a single non-empty path string.", call. = FALSE) + } + h5ad_path <- file.path(work_dir, "adata.h5ad") + out_dir <- file.path(work_dir, "scxpand_out") + + if (file.exists(h5ad_path)) { + if (!isTRUE(overwrite)) { + stop("Output exists (use overwrite = TRUE to replace): ", h5ad_path, + call. = FALSE) + } + unlink(h5ad_path) + } + if (!dir.exists(work_dir)) { + ok <- dir.create(work_dir, recursive = TRUE, showWarnings = FALSE) + if (!ok || !dir.exists(work_dir)) { + stop("Failed to create work_dir: ", work_dir, call. = FALSE) + } + } + if (transient && !isTRUE(keep_files)) { + on.exit(unlink(work_dir, recursive = TRUE), add = TRUE) + } + + # --- Model cache ------------------------------------------------------ + # scXpand would otherwise write ".scxpand_cache" into the current working + # directory. Keep downloaded models in the standard per-user R cache so + # they persist across sessions without touching the user's project. + if (is.null(cache_dir)) { + cache_dir <- file.path(tools::R_user_dir("immLynx", which = "cache"), + "scxpand") + } else if (!is.character(cache_dir) || length(cache_dir) != 1L || + is.na(cache_dir) || !nzchar(cache_dir)) { + stop("cache_dir must be a single non-empty path string.", call. = FALSE) + } + if (!dir.exists(cache_dir)) { + ok <- dir.create(cache_dir, recursive = TRUE, showWarnings = FALSE) + if (!ok || !dir.exists(cache_dir)) { + stop("Failed to create cache_dir: ", cache_dir, call. = FALSE) + } + } + cache_dir <- normalizePath(cache_dir, mustWork = TRUE) + + # --- Gene identifiers ------------------------------------------------- + res_ids <- .resolveGeneIDs(sce, gene_ids = gene_ids, + map_symbols = map_symbols, + ensembl_min_frac = ensembl_min_frac, + multi_map = multi_map, collapse = collapse, + verbose = verbose) + + # --- Ground-truth labels ---------------------------------------------- + labels <- NULL + if (isTRUE(derive_labels)) { + labels <- .deriveExpansionLabels(sce, clone_col = clone_col, + sample_col = sample_col, + median_basis = median_basis, + verbose = verbose) + if (is.null(labels)) { + warning("No clone data found in colData (looked for CTstrict, CTaa). ", + "Running inference without ground-truth labels; AUROC will be ", + "unavailable.", call. = FALSE) + } + } + + # --- Cell subset ------------------------------------------------------ + scored <- colnames(sce) + if (!is.null(labels) && label_cells == "clonal") { + has <- !is.na(labels$expansion) + if (!any(has)) { + stop("label_cells = 'clonal' but no cell has a clone call.", + call. = FALSE) + } + if (!all(has) && verbose) { + message("Scoring ", sum(has), "/", ncol(sce), + " cells with clone calls (label_cells = 'clonal').") + } + sce <- sce[, has, drop = FALSE] + labels$clone_id_size <- labels$clone_id_size[has] + labels$median_clone_size <- labels$median_clone_size[has] + labels$expansion <- labels$expansion[has] + scored <- colnames(sce) + } + + # A single-class ground-truth column makes scXpand's evaluator raise + # inside roc_auc_score, so drop it rather than lose the whole run. + if (!is.null(labels)) { + classes <- unique(labels$expansion[!is.na(labels$expansion)]) + if (length(classes) < 2L) { + warning("Derived expansion labels have only one class ('", + paste(classes, collapse = ", "), + "'); omitting the ground-truth column so scXpand's evaluation ", + "does not fail. Predictions are unaffected.", call. = FALSE) + labels$expansion_dropped <- labels$expansion + labels$expansion <- NULL + } + } + + # --- Stage the H5AD --------------------------------------------------- + .stageScXpandH5AD(sce, res_ids = res_ids, labels = labels, + obs_columns = obs_columns, assay = assay, + h5ad_path = h5ad_path, verbose = verbose) + + # --- Inference -------------------------------------------------------- + if (verbose) { + message("Running scXpand inference with '", model_name, "' on ", + length(scored), " cells x ", length(res_ids$ids), " features...") + } + + py_out <- tryCatch( + .run_in_scxpand_env(function(h5ad, model, save_dir, cache, bsz, nw) { + # convert = FALSE throughout: run_inference() returns an + # InferenceResults dataclass whose fields are a numpy array, a nested + # dict and a custom ModelInfo object. Auto-converting the container + # yields a half-converted structure, so we hold a proxy and convert + # each field explicitly, guarding them individually against upstream + # field renames. + sx <- reticulate::import("scxpand", convert = FALSE) + ad <- reticulate::import("anndata", convert = FALSE) + + # Download as an explicit step rather than letting run_inference do + # it, for two reasons. + # + # First, run_inference() calls download_pretrained_model() without a + # cache_dir, and that defaults to ".scxpand_cache" in the *current + # working directory* -- which would drop a few hundred megabytes into + # whatever project the user happens to be sitting in. + # + # Second, scXpand's registry points at + # https://figshare.com/ndownloader/articles/..., which answers HTTP + # 202 with an empty body; pooch writes that empty file out and the + # unzip then fails with "File is not a zip file". The same archive on + # the ndownloader.figshare.com host returns 200 and a valid zip, so + # we take the registry URL and rewrite the host. The substitution is + # a no-op once upstream fixes its URLs. + url <- as.character(reticulate::py_to_r(reticulate::py_get_attr( + sx$get_pretrained_model_info(model), "url"))) + url <- sub("^https://figshare\\.com/ndownloader/", + "https://ndownloader.figshare.com/", url) + + model_path <- as.character(reticulate::py_to_r( + sx$download_pretrained_model(model_url = url, cache_dir = cache))) + + res <- sx$run_inference(data_path = h5ad, + model_path = model_path, + save_path = save_dir, + batch_size = as.integer(bsz), + num_workers = as.integer(nw)) + + preds <- as.numeric(reticulate::py_to_r(res$predictions)) + + metrics <- tryCatch({ + if (isTRUE(reticulate::py_to_r(res$has_metrics))) { + reticulate::py_to_r(res$metrics) + } else NULL + }, error = function(e) NULL) + + version <- tryCatch( + as.character(reticulate::py_to_r( + reticulate::py_get_attr(sx, "__version__"))), + error = function(e) NA_character_) + + # Read obs_names back so R can match by barcode rather than trusting + # positional alignment through the file round-trip. + obs_names <- as.character(reticulate::py_to_r( + ad$read_h5ad(h5ad, backed = "r")$obs_names$to_list())) + + list(predictions = preds, metrics = metrics, model_path = model_path, + version = version, obs_names = obs_names) + }, + h5ad = h5ad_path, model = model_name, save_dir = out_dir, + cache = cache_dir, bsz = batch_size, nw = num_workers), + error = function(e) e + ) + + if (inherits(py_out, "error")) { + msg <- conditionMessage(py_out) + hint <- if (grepl("not a zip file|Failed to download", msg)) paste0( + "\nThis is a model download failure, not a problem with your data. ", + "The failed response is cached, so clear it before retrying:\n ", + "unlink(\"", cache_dir, "\", recursive = TRUE)") else "" + + stop("scXpand inference failed: ", msg, "\n", + "Staged H5AD: ", + if (isTRUE(keep_files) || !transient) h5ad_path else + "(deleted; re-run with keep_files = TRUE to inspect)", "\n", + "Model cache: ", cache_dir, hint, call. = FALSE) + } + + preds <- .alignPredictions(py_out$predictions, py_out$obs_names, scored) + + # --- Assemble the run summary ------------------------------------------ + truth <- if (!is.null(labels)) { + if (is.null(labels$expansion)) labels$expansion_dropped else labels$expansion + } else NULL + + auroc_immLynx <- if (!is.null(truth)) { + .aurocR(as.numeric(preds), truth == "expanded") + } else NA_real_ + + metrics_flat <- .flattenMetrics(py_out$metrics) + auroc_scxpand <- { + hit <- grep("(^|\\.)AUROC$", names(metrics_flat), ignore.case = TRUE) + if (length(hit)) unname(metrics_flat[[hit[1]]]) else NA_real_ + } + + # run_inference() only populates model_info for its registry branch, and + # we deliberately take the local-path branch to control the cache. Rebuild + # the same fields from the registry table plus the resolved path. + reg <- .SCXPAND_MODELS[.SCXPAND_MODELS$model_name == model_name, , + drop = FALSE] + model_info <- list( + model_name = model_name, + model_type = if (nrow(reg)) reg$model_type[1] else NA_character_, + version = if (nrow(reg)) reg$version[1] else NA_character_, + source = "registry", + path = py_out$model_path) + + summary_list <- list( + model = model_name, + model_info = model_info, + cache_dir = cache_dir, + scxpand_version = py_out$version, + n_cells = length(preds), + n_features = length(res_ids$ids), + threshold = threshold, + metrics = metrics_flat, + metrics_raw = py_out$metrics, + auroc = auroc_scxpand, + auroc_immLynx = auroc_immLynx, + labels = if (is.null(labels)) NULL else list( + clone_col = labels$clone_col, + sample_col = labels$sample_col, + median_basis = labels$median_basis, + median_by_sample = labels$median_by_sample, + n_labelled = labels$n_labelled, + n_clones = labels$n_clones, + n_samples = labels$n_samples), + gene_mapping = res_ids$report, + work_dir = if (isTRUE(keep_files) || !transient) work_dir else NULL, + timestamp = Sys.time()) + + if (verbose && !is.na(auroc_immLynx)) { + message("AUROC against derived expansion labels: ", + sprintf("%.3f", auroc_immLynx)) + } + + pred_call <- ifelse(is.na(preds), NA_character_, + ifelse(preds >= threshold, "expanded", "non-expanded")) + + # --- Return ------------------------------------------------------------ + if (isTRUE(return_object)) { + p <- function(x) if (nzchar(column_prefix)) + paste0(column_prefix, "_", x) else x + + input <- .writeCellColumn(input, p("expansion_prob"), + as.numeric(preds), scored) + input <- .writeCellColumn(input, p("expansion_pred"), pred_call, scored) + if (!is.null(labels)) { + input <- .writeCellColumn(input, p("clone_id_size"), + labels$clone_id_size, scored) + input <- .writeCellColumn(input, p("median_clone_size"), + labels$median_clone_size, scored) + input <- .writeCellColumn(input, p("expansion_truth"), truth, scored) + } + input <- .writeObjMetadata(input, "scXpand", summary_list) + + if (verbose) { + message("Added ", p("expansion_prob"), " and ", p("expansion_pred"), + " to cell metadata.") + } + return(input) + } + + df <- data.frame(barcode = scored, + expansion_prob = as.numeric(preds), + expansion_pred = pred_call, + stringsAsFactors = FALSE) + if (!is.null(labels)) { + df$clone_id_size <- labels$clone_id_size + df$median_clone_size <- labels$median_clone_size + df$expansion_truth <- truth + } + rownames(df) <- NULL + attr(df, "scXpand") <- summary_list + df +} diff --git a/inst/WORDLIST b/inst/WORDLIST index 9925c46..48aa6e1 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,4 +1,5 @@ AIRR +AUROC AutoModel BCR BioC @@ -8,8 +9,10 @@ CMD Clonotypes Codecov DBSCAN +DataLoader ESM Embeddings +Ensembl GEX GLIPH HLA @@ -19,6 +22,7 @@ IGK IGL IMGT Isacchini +LightGBM MCL Metaclone Metaclones @@ -27,6 +31,8 @@ MuData OLGA's ORCID Pgen +Pretrained +PyTorch SCE SCEPTR Scalability @@ -44,12 +50,15 @@ TSV Tokenize Tokenizer UMAP +UMI VDJ Valkiers aa al +autoencoder barcode barcodes +benchmarked benchmarking bfloat cdr @@ -68,12 +77,14 @@ cuda dbscan deduplicates disambiguated +downloader embeddings erroring esm et extractTCRdata facebook +figshare generateOLGA heterodimeric https @@ -84,6 +95,8 @@ immApex immGLIPH immunoGlobulin leiden +lexicographically +lightgbm louvain mcl metaclone @@ -100,7 +113,9 @@ nt numpy olga pgen +pooch pre +pretrained proteinEmbeddings pubmed pw @@ -112,12 +127,15 @@ runMetaclonotypist scRNA scRepertoire scRepertoire's +scXpand +scXpand's scanpy scater sceptr scirpy scirpy's scran +scxpand soNNia sonia sonnia @@ -127,4 +145,7 @@ tf thymic tokenized tokenizer +torch +transcriptome +unmapped unselected diff --git a/man/immLynx-package.Rd b/man/immLynx-package.Rd index 2d869c6..3696017 100644 --- a/man/immLynx-package.Rd +++ b/man/immLynx-package.Rd @@ -40,6 +40,13 @@ several popular Python-based tools and makes them accessible within R workflows. \item \code{\link{runSoNNia}}: Infer selection with soNNia } +**Expansion Prediction:** +\itemize{ + \item \code{\link{runScXpand}}: Predict T-cell clonal expansion from gene + expression alone with scXpand, no paired TCR sequencing required + \item \code{\link{listScXpandModels}}: List the scXpand pretrained models +} + **Utility Functions:** \itemize{ \item \code{\link{extractTCRdata}}: Extract TCR data from SingleCellExperiment objects @@ -75,6 +82,7 @@ Python packages are managed automatically by basilisk. The following are include \item metaclonotypist - Metaclone discovery \item transformers - Hugging Face models \item torch - PyTorch for GPU support + \item scxpand - Clonal expansion prediction from gene expression } } @@ -88,5 +96,10 @@ Useful links: \author{ \strong{Maintainer}: Nick Borcherding \email{ncborch@gmail.com} (\href{https://orcid.org/0000-0003-1427-6342}{ORCID}) +Authors: +\itemize{ + \item Nick Borcherding \email{ncborch@gmail.com} (\href{https://orcid.org/0000-0003-1427-6342}{ORCID}) +} + } \keyword{internal} diff --git a/man/listScXpandModels.Rd b/man/listScXpandModels.Rd new file mode 100644 index 0000000..3e20ac4 --- /dev/null +++ b/man/listScXpandModels.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/runScXpand.R +\name{listScXpandModels} +\alias{listScXpandModels} +\title{List available scXpand pretrained models} +\usage{ +listScXpandModels(refresh = FALSE) +} +\arguments{ +\item{refresh}{Logical. Query the installed \code{scxpand} package rather +than the static table. Requires the scXpand basilisk environment and +will build it on first use. Default \code{FALSE}.} +} +\value{ +A \code{data.frame} with columns \code{model_name}, + \code{model_type}, \code{version} and \code{description}. When + \code{refresh = TRUE}, \code{model_type} and \code{description} may be + \code{NA} for models absent from the static table. +} +\description{ +Returns the pan-cancer models published with scXpand. By default this +reads a static table shipped with immLynx so it never triggers a +multi-gigabyte Python environment build. Set \code{refresh = TRUE} to +query the installed \code{scxpand} package, which is authoritative if +upstream has added models. +} +\examples{ +listScXpandModels() + +\donttest{ +# Authoritative, but builds the scXpand environment on first call. +listScXpandModels(refresh = TRUE) +} +} diff --git a/man/runScXpand.Rd b/man/runScXpand.Rd new file mode 100644 index 0000000..4a117bf --- /dev/null +++ b/man/runScXpand.Rd @@ -0,0 +1,209 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/runScXpand.R +\name{runScXpand} +\alias{runScXpand} +\title{Predict T-cell clonal expansion from gene expression with scXpand} +\usage{ +runScXpand( + input, + model_name = "pan_cancer_autoencoder", + assay = "counts", + gene_ids = "rownames", + map_symbols = c("auto", "always", "never"), + ensembl_min_frac = 0.5, + multi_map = c("first", "expand", "drop"), + collapse = c("sum", "first", "drop"), + derive_labels = TRUE, + clone_col = NULL, + sample_col = NULL, + median_basis = c("clone", "cell"), + label_cells = c("clonal", "all"), + threshold = 0.5, + obs_columns = NULL, + batch_size = 1024L, + num_workers = 0L, + work_dir = NULL, + cache_dir = NULL, + keep_files = FALSE, + overwrite = FALSE, + column_prefix = "scXpand", + return_object = TRUE, + verbose = TRUE +) +} +\arguments{ +\item{input}{A \code{SingleCellExperiment} or \code{Seurat} object. +SingleCellExperiment is the native format and is used without +conversion; Seurat objects are converted in and returned as Seurat.} + +\item{model_name}{Pretrained model to use. See +\code{\link{listScXpandModels}}. Default +\code{"pan_cancer_autoencoder"}.} + +\item{assay}{Assay holding raw UMI counts. Default \code{"counts"}.} + +\item{gene_ids}{Where to find gene identifiers: \code{"rownames"} +(default), the name of a \code{rowData} column, or a character vector +of length \code{nrow(input)}.} + +\item{map_symbols}{One of \code{"auto"} (map symbols to Ensembl IDs only +if too few identifiers already look like Ensembl IDs), \code{"always"}, +or \code{"never"}. Mapping requires \code{org.Hs.eg.db} and +\code{AnnotationDbi}.} + +\item{ensembl_min_frac}{Minimum fraction of identifiers that must resolve +to Ensembl IDs. Default \code{0.5}. Lower it for targeted panels.} + +\item{multi_map}{How to resolve one symbol mapping to several Ensembl +IDs: \code{"first"} (default, lexicographically first), +\code{"expand"} (emit all), or \code{"drop"}.} + +\item{collapse}{How to resolve several rows mapping to one Ensembl ID: +\code{"sum"} (default), \code{"first"}, or \code{"drop"}.} + +\item{derive_labels}{Logical. Derive ground-truth expansion labels from +clone data when present. Default \code{TRUE}.} + +\item{clone_col}{Column holding clone identity. Default \code{NULL}: +use \code{"CTstrict"} if present, else \code{"CTaa"}.} + +\item{sample_col}{Column identifying the sample or patient that clone +sizes are tabulated within. Default \code{NULL}, which pools all cells +into one sample and warns.} + +\item{median_basis}{\code{"clone"} (default) or \code{"cell"}. See +Details.} + +\item{label_cells}{\code{"clonal"} (default) scores only cells with a +clone call, which guarantees the ground-truth column has no missing +values; \code{"all"} scores every cell.} + +\item{threshold}{Probability cutoff for the binary call. Default +\code{0.5}.} + +\item{obs_columns}{Character vector of \code{colData} columns to carry +into the staged H5AD. Default \code{NULL}: all of them, minus +list-columns and the scRepertoire \code{CT*} columns.} + +\item{batch_size}{Inference batch size. Default \code{1024}.} + +\item{num_workers}{DataLoader workers. Default \code{0}, because worker +processes forking out of embedded Python inside a basilisk child can +hang on macOS. Raise it on Linux for speed.} + +\item{work_dir}{Directory for the staged H5AD and scXpand's own output. +Default \code{NULL}: a temporary directory, removed on exit.} + +\item{cache_dir}{Directory for downloaded pretrained models. Default +\code{NULL}: a \code{scxpand} subdirectory of +\code{tools::R_user_dir("immLynx", "cache")}. Passing this explicitly +matters because scXpand's own default would write a +\code{.scxpand_cache} directory into the current working directory.} + +\item{keep_files}{Logical. Keep the staged H5AD and scXpand outputs. +Default \code{FALSE}.} + +\item{overwrite}{Logical. Replace an existing \code{adata.h5ad} in an +explicitly supplied \code{work_dir}. Default \code{FALSE}.} + +\item{column_prefix}{Prefix for the columns written back. Default +\code{"scXpand"}.} + +\item{return_object}{Logical. Return the input object with new columns +(default), or a \code{data.frame} of per-cell results.} + +\item{verbose}{Logical. Emit progress messages. Default \code{TRUE}.} +} +\value{ +If \code{return_object = TRUE}, the input object (same class) + with \code{_expansion_prob} and \code{_expansion_pred} + added to cell metadata, plus \code{_clone_id_size}, + \code{_median_clone_size} and \code{_expansion_truth} + when labels were derived. A run summary is stored under + \code{metadata(x)$scXpand} for SingleCellExperiment or + \code{x@misc$scXpand} for Seurat, holding the model info, the gene + mapping report, scXpand's metrics and an independently computed AUROC. + If \code{return_object = FALSE}, a \code{data.frame} with one row per + scored cell and the run summary attached as + \code{attr(x, "scXpand")}. +} +\description{ +Runs inference with a scXpand pretrained pan-cancer model to estimate, +for each cell, the probability that it belongs to an expanded T-cell +clone. scXpand uses gene expression only, so this works on datasets with +no paired TCR sequencing. When the object does carry scRepertoire clone +calls, \code{runScXpand} also derives ground-truth expansion labels so +the prediction can be benchmarked against the observed repertoire. +} +\details{ +\strong{Input requirements.} scXpand's pretrained models were trained on +raw UMI counts from T cells indexed by Ensembl gene ID. All three matter: +\itemize{ + \item \emph{Raw counts.} Normalized or log-transformed values are + rejected. Use \code{assay = "counts"}. + \item \emph{Ensembl IDs.} Most Seurat and scRepertoire objects carry + gene symbols; see \code{gene_ids} and \code{map_symbols}. + \item \emph{T cells only.} Filter to T cells before calling. There is + no reliable programmatic check, and predictions on non-T cells are + meaningless. +} +Genes the model expects but cannot find are zero-filled by scXpand, and +extra genes are ignored, so a partial overlap degrades quietly rather +than erroring. That is why the gene mapping summary is worth reading. + +\strong{Symbol mapping.} When identifiers are symbols, they are mapped +through \code{org.Hs.eg.db} (SYMBOL first, then ALIAS). One symbol often +maps to several Ensembl IDs because of alternate haplotypes: +\code{HLA-DRA} maps to eight. \code{multi_map = "first"} takes the +lexicographically first, which is reproducible but may pick a scaffold +the model does not know; \code{multi_map = "expand"} emits every +candidate, which maximizes overlap with the model's panel at the cost of +a larger staged matrix. Several input rows can also land on one Ensembl +ID, in which case \code{collapse} decides what happens (summing raw +counts is the default and preserves UMI semantics). + +\strong{Expansion labels.} scXpand calls a cell expanded when its clone +size exceeds 1.5 times the median clone size for its sample. Clone sizes +are tabulated fresh from \code{clone_col} within \code{sample_col}; +scRepertoire's \code{clonalFrequency} column is deliberately not used, +because \code{combineExpression()} computes it under whatever +\code{group.by} was in effect and it is not necessarily a per-sample +count. "Median clone size" is genuinely ambiguous, so +\code{median_basis} selects the reading: \code{"clone"} takes the median +over unique clones, which in typical 10x data is 1 (most clones are +singletons) and therefore reduces the rule to \emph{clone size at least +2}; \code{"cell"} takes the median over cells, which is dominated by +large clones and is much stricter. + +\strong{Cost.} The first call builds a dedicated basilisk environment +(Python 3.11, PyTorch CPU, scxpand) that occupies several gigabytes, and +downloads the selected model (roughly 35 MB) from figshare into +\code{cache_dir}. Both are reused afterwards. If a download does fail, +clear \code{cache_dir} before retrying: the downloader caches the failed +response and will otherwise keep reusing it. +} +\examples{ +data(immLynx_example) + +# The label derivation is pure R and runs without Python. Note that the +# derived clone sizes are tabulated per sample and so need not match the +# clonalFrequency column scRepertoire wrote. +labs <- immLynx:::.deriveExpansionLabels( + immLynx_example, + sample_col = "Patient" +) +table(labs$expansion, useNA = "ifany") + +\donttest{ +# Full inference. The example object carries gene symbols and is not +# filtered to T cells, so this demonstrates the mechanics rather than a +# meaningful biological result. +sce <- runScXpand( + immLynx_example, + sample_col = "Patient", + assay = "counts" +) +summary(SummarizedExperiment::colData(sce)$scXpand_expansion_prob) +S4Vectors::metadata(sce)$scXpand$auroc_immLynx +} +} diff --git a/tests/testthat/helper-immLynx.R b/tests/testthat/helper-immLynx.R index a5b30d3..c54a981 100644 --- a/tests/testthat/helper-immLynx.R +++ b/tests/testthat/helper-immLynx.R @@ -125,6 +125,92 @@ skip_if_no_scanpy_env <- function() { if (!ok) testthat::skip("scanpyExportEnv not available") } +# Skip test if the scXpand basilisk env is not available. This env carries +# python 3.11 + PyTorch and runs to several gigabytes, so its first-use +# install cost is much larger than the other two. +skip_if_no_scxpand_env <- function() { + skip_on_bioc_build() + ok <- tryCatch({ + proc <- basilisk::basiliskStart(immLynx:::scXpandEnv) + on.exit(basilisk::basiliskStop(proc)) + TRUE + }, error = function(e) FALSE) + if (!ok) testthat::skip("scXpandEnv not available") +} + +# Inference additionally downloads a pretrained model from figshare on +# first use. Gate that behind an explicit opt-in so it never fires +# unattended in CI. +skip_if_no_scxpand_model <- function() { + skip_if_no_scxpand_env() + if (!nzchar(Sys.getenv("IMMLYNX_TEST_SCXPAND"))) { + testthat::skip("set IMMLYNX_TEST_SCXPAND=1 to run scXpand model tests") + } +} + +# Mock SCE indexed by Ensembl IDs with raw integer counts — the shape +# runScXpand expects. Clone calls follow the scRepertoire "_" contract. +mock_ensembl_sce <- function(n_genes = 50, n_cells = 20, seed = 1) { + set.seed(seed) + counts <- Matrix::Matrix( + matrix(stats::rpois(n_genes * n_cells, lambda = 3), n_genes, n_cells), + sparse = TRUE) + rownames(counts) <- sprintf("ENSG%011d", seq_len(n_genes)) + colnames(counts) <- paste0("cell", seq_len(n_cells)) + + # Skewed clone sizes: one dominant clone, a couple of medium ones, the + # rest singletons. A flat distribution would make every cell fall on the + # same side of the 1.5x median cutoff and leave nothing to score. + n_big <- max(2L, n_cells %/% 4L) + n_mid <- max(1L, n_cells %/% 10L) + clone_idx <- c(rep(1L, n_big), rep(2L, n_mid), rep(3L, n_mid), + seq.int(4L, length.out = max(0L, n_cells - n_big - + 2L * n_mid))) + clone_idx <- clone_idx[seq_len(n_cells)] + clone <- paste0("CAS", clone_idx, "F_CAS", clone_idx, "F") + + SingleCellExperiment::SingleCellExperiment( + assays = list(counts = counts), + colData = S4Vectors::DataFrame( + sample = rep(c("S1", "S2"), length.out = n_cells), + CTstrict = clone, + row.names = colnames(counts) + ) + ) +} + +# Mock SCE whose clone layout makes the two median_basis readings diverge, +# so the tests can pin the arithmetic by hand. +# +# sample A: c1 x 4, c2 x 1, c3 x 1 +# median over unique clones = median(4, 1, 1) = 1 -> cutoff 1.5 +# -> only c1 expanded +# median over cells = median(4, 4, 4, 4, 1, 1) = 4 -> cutoff 6 +# -> nothing expanded +# sample B: d1 x 2, d2 x 2, plus one cell with no clone call +mock_clonal_sce <- function() { + clone <- c(rep("c1", 4), "c2", "c3", rep("d1", 2), rep("d2", 2), NA) + sample <- c(rep("A", 6), rep("B", 5)) + n <- length(clone) + + # Dimnames set up front, and values varied: Matrix() collapses a square + # constant matrix to a symmetric dsCMatrix, which ties rownames to + # colnames and would drop the gene IDs. + counts <- Matrix::Matrix( + matrix(seq_len(4L * n), nrow = 4L, + dimnames = list(sprintf("ENSG%011d", seq_len(4L)), + paste0("cell", seq_len(n)))), + sparse = TRUE) + + SingleCellExperiment::SingleCellExperiment( + assays = list(counts = counts), + colData = S4Vectors::DataFrame( + CTstrict = clone, sample = sample, + row.names = colnames(counts) + ) + ) +} + # Mock SCE with hand-crafted scRepertoire CT* fields. Used by the # exportToScanpy / .buildAIRR test suite. scRepertoire uses "_" as the # chain-slot separator (TRA before, TRB after) with explicit "NA" tokens diff --git a/tests/testthat/test-runScXpand.R b/tests/testthat/test-runScXpand.R new file mode 100644 index 0000000..e16332a --- /dev/null +++ b/tests/testthat/test-runScXpand.R @@ -0,0 +1,624 @@ +library(testthat) + +# =========================================================================== +# Signatures and the static model table +# =========================================================================== + +test_that("runScXpand signature has correct defaults", { + f <- formals(runScXpand) + expect_equal(f$model_name, "pan_cancer_autoencoder") + expect_equal(f$assay, "counts") + expect_equal(f$gene_ids, "rownames") + expect_equal(eval(f$map_symbols), c("auto", "always", "never")) + expect_equal(f$ensembl_min_frac, 0.5) + expect_equal(eval(f$multi_map), c("first", "expand", "drop")) + expect_equal(eval(f$collapse), c("sum", "first", "drop")) + expect_true(f$derive_labels) + expect_null(eval(f$clone_col)) + expect_null(eval(f$sample_col)) + expect_equal(eval(f$median_basis), c("clone", "cell")) + expect_equal(eval(f$label_cells), c("clonal", "all")) + expect_equal(f$threshold, 0.5) + expect_equal(f$batch_size, 1024L) + expect_equal(f$num_workers, 0L) + expect_null(eval(f$work_dir)) + expect_null(eval(f$cache_dir)) + expect_false(f$keep_files) + expect_false(f$overwrite) + expect_equal(f$column_prefix, "scXpand") + expect_true(f$return_object) + expect_true(f$verbose) +}) + +test_that("listScXpandModels returns the registry without touching Python", { + m <- listScXpandModels() + expect_s3_class(m, "data.frame") + expect_equal(nrow(m), 5L) + expect_equal(colnames(m), + c("model_name", "model_type", "version", "description")) + expect_true("pan_cancer_autoencoder" %in% m$model_name) + expect_setequal(m$model_type, + c("autoencoder", "mlp", "lightgbm", "logistic", "svm")) + expect_false(formals(listScXpandModels)$refresh) +}) + +# =========================================================================== +# Ensembl identifier helpers +# =========================================================================== + +test_that(".isEnsembl matches only well-formed Ensembl gene IDs", { + expect_true(immLynx:::.isEnsembl("ENSG00000198851")) + expect_false(immLynx:::.isEnsembl("ENSG0000019885")) # too short + expect_false(immLynx:::.isEnsembl("ENST00000198851")) # transcript + expect_false(immLynx:::.isEnsembl("CD3E")) + expect_false(immLynx:::.isEnsembl(NA_character_)) +}) + +test_that(".stripEnsemblVersion strips only Ensembl versions", { + expect_equal(immLynx:::.stripEnsemblVersion("ENSG00000198851.12"), + "ENSG00000198851") + # A blanket sub() would mangle these. + expect_equal(immLynx:::.stripEnsemblVersion("MARCH1.2"), "MARCH1.2") + expect_equal(immLynx:::.stripEnsemblVersion("7SK.2"), "7SK.2") + expect_equal(immLynx:::.stripEnsemblVersion(NA_character_), NA_character_) +}) + +# =========================================================================== +# .resolveGeneIDs +# =========================================================================== + +.gene_sce <- function(rn, id_col = NULL) { + n <- length(rn) + # Values must vary and dimnames must be set up front: a square matrix of + # constant values makes Matrix() return a symmetric dsCMatrix, which ties + # rownames to colnames and silently discards the gene identifiers. + counts <- Matrix::Matrix( + matrix(seq_len(n * 4L), nrow = n, + dimnames = list(rn, paste0("c", 1:4))), + sparse = TRUE) + sce <- SingleCellExperiment::SingleCellExperiment( + assays = list(counts = counts)) + if (!is.null(id_col)) SummarizedExperiment::rowData(sce)$ID <- id_col + sce +} + +test_that(".resolveGeneIDs accepts all three gene_ids forms", { + ens <- sprintf("ENSG%011d", 1:3) + + r1 <- immLynx:::.resolveGeneIDs(.gene_sce(ens), verbose = FALSE) + expect_equal(r1$ids, ens) + expect_equal(r1$keep, 1:3) + expect_null(r1$collapse_groups) + + r2 <- immLynx:::.resolveGeneIDs(.gene_sce(c("a", "b", "c"), ens), + gene_ids = "ID", verbose = FALSE) + expect_equal(r2$ids, ens) + + r3 <- immLynx:::.resolveGeneIDs(.gene_sce(c("a", "b", "c")), + gene_ids = ens, verbose = FALSE) + expect_equal(r3$ids, ens) +}) + +test_that(".resolveGeneIDs strips Ensembl versions from rownames", { + rn <- paste0(sprintf("ENSG%011d", 1:3), ".", 1:3) + r <- immLynx:::.resolveGeneIDs(.gene_sce(rn), verbose = FALSE) + expect_equal(r$ids, sprintf("ENSG%011d", 1:3)) +}) + +test_that(".resolveGeneIDs rejects malformed gene_ids", { + sce <- .gene_sce(c("a", "b", "c"), sprintf("ENSG%011d", 1:3)) + expect_error(immLynx:::.resolveGeneIDs(sce, gene_ids = "nope"), + "not a rowData column") + expect_error(immLynx:::.resolveGeneIDs(sce, gene_ids = c("x", "y")), + "character vector of length nrow") + expect_error(immLynx:::.resolveGeneIDs(sce, gene_ids = 1:3), + "character vector of length nrow") +}) + +test_that("map_symbols = 'never' errors on symbol identifiers", { + sce <- .gene_sce(c("CD3E", "CD8A", "GZMB")) + expect_error( + immLynx:::.resolveGeneIDs(sce, map_symbols = "never", verbose = FALSE), + "0\\.0% of gene identifiers") +}) + +test_that(".resolveGeneIDs collapses duplicate Ensembl IDs", { + ids <- c("ENSG00000000001", "ENSG00000000001", "ENSG00000000002") + sce <- .gene_sce(c("g1", "g2", "g3"), ids) + + r_sum <- immLynx:::.resolveGeneIDs(sce, gene_ids = "ID", collapse = "sum", + verbose = FALSE) + expect_equal(r_sum$ids, c("ENSG00000000001", "ENSG00000000002")) + expect_equal(r_sum$keep, 1:3) + expect_equal(as.character(r_sum$collapse_groups), ids) + + r_first <- immLynx:::.resolveGeneIDs(sce, gene_ids = "ID", + collapse = "first", verbose = FALSE) + expect_equal(r_first$keep, c(1L, 3L)) + expect_null(r_first$collapse_groups) + + r_drop <- immLynx:::.resolveGeneIDs(sce, gene_ids = "ID", collapse = "drop", + verbose = FALSE) + expect_equal(r_drop$ids, "ENSG00000000002") + expect_equal(r_drop$keep, 3L) +}) + +test_that(".resolveGeneIDs errors when nothing resolves", { + sce <- .gene_sce(c("aa", "bb", "cc")) + expect_error( + immLynx:::.resolveGeneIDs(sce, map_symbols = "never", + ensembl_min_frac = 0, verbose = FALSE), + "No gene identifiers could be resolved") +}) + +# =========================================================================== +# Symbol to Ensembl mapping +# =========================================================================== + +test_that(".mapSymbolsToEnsembl maps symbols and counts what it dropped", { + skip_if_not_installed("org.Hs.eg.db") + skip_if_not_installed("AnnotationDbi") + + m <- immLynx:::.mapSymbolsToEnsembl( + c("CD3E", "HLA-DRA", "NOTAGENEATALL"), multi_map = "first") + + expect_equal(m$map$ENSEMBL[m$map$SYMBOL == "CD3E"], "ENSG00000198851") + expect_equal(m$n_query, 3L) + expect_equal(m$n_unmapped, 1L) + # HLA-DRA sits on several alt haplotypes. + expect_gte(m$n_ambiguous, 1L) +}) + +test_that(".mapSymbolsToEnsembl multi_map modes behave differently", { + skip_if_not_installed("org.Hs.eg.db") + skip_if_not_installed("AnnotationDbi") + + first <- immLynx:::.mapSymbolsToEnsembl("HLA-DRA", multi_map = "first") + expand <- immLynx:::.mapSymbolsToEnsembl("HLA-DRA", multi_map = "expand") + dropped <- immLynx:::.mapSymbolsToEnsembl("HLA-DRA", multi_map = "drop") + + expect_equal(nrow(first$map), 1L) + expect_gt(nrow(expand$map), 1L) + expect_equal(nrow(dropped$map), 0L) + + # Lexicographic ordering makes "first" reproducible run to run. + again <- immLynx:::.mapSymbolsToEnsembl("HLA-DRA", multi_map = "first") + expect_identical(first$map, again$map) + expect_equal(first$map$ENSEMBL, min(expand$map$ENSEMBL)) +}) + +test_that(".resolveGeneIDs maps symbols end to end", { + skip_if_not_installed("org.Hs.eg.db") + skip_if_not_installed("AnnotationDbi") + + sce <- .gene_sce(c("CD3E", "CD8A", "GZMB", "NOTAGENEATALL")) + r <- immLynx:::.resolveGeneIDs(sce, map_symbols = "auto", verbose = FALSE) + + expect_true(all(immLynx:::.isEnsembl(r$ids))) + expect_equal(length(r$ids), 3L) # the bogus symbol is dropped + expect_true(r$report$mapped) + expect_equal(r$report$n_unmapped, 1L) + + # "expand" duplicates rows rather than dropping candidates. + sce2 <- .gene_sce(c("CD3E", "HLA-DRA")) + re <- immLynx:::.resolveGeneIDs(sce2, map_symbols = "auto", + multi_map = "expand", verbose = FALSE) + expect_gt(length(re$ids), 2L) + expect_true(any(duplicated(re$keep))) +}) + +# =========================================================================== +# Expansion label derivation +# =========================================================================== + +test_that(".deriveExpansionLabels reproduces scXpand's 1.5x median rule", { + sce <- mock_clonal_sce() + + l <- immLynx:::.deriveExpansionLabels(sce, sample_col = "sample", + median_basis = "clone", + verbose = FALSE) + + # Sample A: sizes 4,4,4,4,1,1; median over unique clones = 1; cutoff 1.5. + expect_equal(l$clone_id_size[1:6], c(4, 4, 4, 4, 1, 1)) + expect_equal(unname(l$median_clone_size[1:6]), rep(1, 6)) + expect_equal(l$expansion[1:6], + c(rep("expanded", 4), "non-expanded", "non-expanded")) + + # Sample B: two clones of size 2; median 2; cutoff 3; nothing expanded. + expect_equal(l$clone_id_size[7:10], rep(2, 4)) + expect_equal(l$expansion[7:10], rep("non-expanded", 4)) + + # The cell with no clone call stays NA in all three vectors. + expect_true(is.na(l$clone_id_size[11])) + expect_true(is.na(l$median_clone_size[11])) + expect_true(is.na(l$expansion[11])) + + expect_equal(l$n_labelled, 10L) + expect_equal(l$n_clones, 5L) + expect_equal(l$n_samples, 2L) + expect_equal(unname(l$median_by_sample[["A"]]), 1) +}) + +test_that("median_basis = 'cell' gives a stricter cutoff", { + sce <- mock_clonal_sce() + l <- immLynx:::.deriveExpansionLabels(sce, sample_col = "sample", + median_basis = "cell", + verbose = FALSE) + # Sample A median over cells = median(4,4,4,4,1,1) = 4; cutoff 6. + expect_equal(unname(l$median_by_sample[["A"]]), 4) + expect_equal(l$expansion[1:6], rep("non-expanded", 6)) + expect_false(any(l$expansion == "expanded", na.rm = TRUE)) +}) + +test_that(".deriveExpansionLabels warns when sample_col is NULL", { + sce <- mock_clonal_sce() + expect_warning( + l <- immLynx:::.deriveExpansionLabels(sce, verbose = FALSE), + "single sample") + # Pooled: sizes 4,1,1,2,2 -> median over unique clones = 2 -> cutoff 3. + expect_equal(unname(l$median_clone_size[1]), 2) + expect_equal(l$expansion[1], "expanded") + expect_equal(l$expansion[7], "non-expanded") +}) + +test_that(".deriveExpansionLabels falls back from CTstrict to CTaa", { + sce <- mock_clonal_sce() + SummarizedExperiment::colData(sce)$CTaa <- + SummarizedExperiment::colData(sce)$CTstrict + SummarizedExperiment::colData(sce)$CTstrict <- NULL + + l <- immLynx:::.deriveExpansionLabels(sce, sample_col = "sample", + verbose = FALSE) + expect_equal(l$clone_col, "CTaa") + expect_equal(l$n_labelled, 10L) +}) + +test_that(".deriveExpansionLabels returns NULL without usable clone data", { + sce <- mock_clonal_sce() + SummarizedExperiment::colData(sce)$CTstrict <- NULL + expect_null(immLynx:::.deriveExpansionLabels(sce, sample_col = "sample", + verbose = FALSE)) + + # scRepertoire writes literal "NA" tokens for cells with no clone call. + sce2 <- mock_clonal_sce() + SummarizedExperiment::colData(sce2)$CTstrict <- "NA" + expect_null(immLynx:::.deriveExpansionLabels(sce2, sample_col = "sample", + verbose = FALSE)) +}) + +test_that(".deriveExpansionLabels rejects unknown columns", { + sce <- mock_clonal_sce() + expect_error(immLynx:::.deriveExpansionLabels(sce, clone_col = "nope"), + "clone_col") + expect_error(immLynx:::.deriveExpansionLabels(sce, sample_col = "nope"), + "sample_col") +}) + +test_that("clone sizes are tabulated fresh from the clone column", { + data(immLynx_example, envir = environment()) + + l <- immLynx:::.deriveExpansionLabels(immLynx_example, + sample_col = "Patient", + verbose = FALSE) + cd <- SummarizedExperiment::colData(immLynx_example) + ok <- !is.na(l$clone_id_size) + expect_true(any(ok)) + + # Pin the arithmetic against a hand-rolled per-patient tabulation. + key <- paste(cd$Patient, cd$CTstrict, sep = "||") + expected <- as.numeric(table(key[ok])[key[ok]]) + expect_equal(l$clone_id_size[ok], expected) +}) + +test_that("per-sample tabulation diverges from a pooled count", { + # The shipped example has no clone shared across patients, so per-patient + # and pooled counts happen to coincide there. Construct a case where a + # clone does span samples: that is the situation in which reading + # scRepertoire's clonalFrequency (computed under whatever group.by was in + # effect) would give the wrong per-sample size. + clone <- c("cX", "cX", "cX", "cY", "cZ", "cX", "cX", "cW", "cV", "cU") + sample <- c(rep("A", 5), rep("B", 5)) + n <- length(clone) + counts <- Matrix::Matrix( + matrix(seq_len(3L * n), nrow = 3L, + dimnames = list(sprintf("ENSG%011d", 1:3), + paste0("cell", seq_len(n)))), + sparse = TRUE) + sce <- SingleCellExperiment::SingleCellExperiment( + assays = list(counts = counts), + colData = S4Vectors::DataFrame(CTstrict = clone, sample = sample, + row.names = colnames(counts))) + + per_sample <- immLynx:::.deriveExpansionLabels(sce, sample_col = "sample", + verbose = FALSE) + pooled <- suppressWarnings( + immLynx:::.deriveExpansionLabels(sce, verbose = FALSE)) + + # cX has 3 cells in A and 2 in B, but 5 overall. + expect_equal(per_sample$clone_id_size[1], 3) + expect_equal(per_sample$clone_id_size[6], 2) + expect_equal(pooled$clone_id_size[1], 5) + expect_false(isTRUE(all.equal(per_sample$clone_id_size, + pooled$clone_id_size))) + + # Within B the unique clone sizes are 2, 1, 1, 1, so the median is 1 and + # the cutoff 1.5: cX's two cells there clear it. Pooled, cX is counted as + # a single clone of 5 and B's cells inherit A's count. + expect_equal(per_sample$median_clone_size[6], 1) + expect_equal(per_sample$expansion[6], "expanded") + expect_equal(per_sample$median_clone_size[1], 1) +}) + +# =========================================================================== +# Counts validation and row collapsing +# =========================================================================== + +test_that(".validateCountsAssay accepts raw counts and rejects everything else", { + X <- Matrix::Matrix(matrix(as.numeric(rpois(40, 3)), 10, 4), sparse = TRUE) + expect_true(immLynx:::.validateCountsAssay(X, "counts")) + + expect_error(immLynx:::.validateCountsAssay(X, "logcounts"), + "normalized or log-transformed") + + Xn <- X; Xn@x <- Xn@x + 0.5 + expect_error(immLynx:::.validateCountsAssay(Xn, "counts"), + "non-integer values") + + Xneg <- X; Xneg@x[1] <- -1 + expect_error(immLynx:::.validateCountsAssay(Xneg, "counts"), + "negative values") + + expect_error( + immLynx:::.validateCountsAssay(Matrix::Matrix(0, 4, 4, sparse = TRUE), + "counts"), + "no non-zero counts") +}) + +test_that(".collapseRowsSum adds colliding rows together", { + Y <- Matrix::Matrix(matrix(c(1, 2, 3, 4, 5, 6), nrow = 3), sparse = TRUE) + colnames(Y) <- c("a", "b") + out <- immLynx:::.collapseRowsSum(Y, factor(c("E1", "E1", "E2"))) + + expect_equal(rownames(out), c("E1", "E2")) + expect_equal(colnames(out), c("a", "b")) + expect_equal(as.matrix(unname(as.matrix(out))), + matrix(c(3, 3, 9, 6), nrow = 2)) +}) + +# =========================================================================== +# Result handling +# =========================================================================== + +test_that(".aurocR matches the Mann-Whitney definition", { + expect_equal(immLynx:::.aurocR(c(0.1, 0.2, 0.8, 0.9), c(0, 0, 1, 1)), 1) + expect_equal(immLynx:::.aurocR(c(0.9, 0.8, 0.2, 0.1), c(0, 0, 1, 1)), 0) + expect_equal(immLynx:::.aurocR(rep(0.5, 4), c(0, 0, 1, 1)), 0.5) + expect_true(is.na(immLynx:::.aurocR(c(0.1, 0.2), c(0, 0)))) + # NAs are dropped pairwise rather than poisoning the result. + expect_equal(immLynx:::.aurocR(c(0.1, NA, 0.8, 0.9), c(0, 1, 1, 1)), 1) +}) + +test_that(".flattenMetrics flattens scalars and drops the rest", { + f <- immLynx:::.flattenMetrics( + list(overall = list(AUROC = 0.81, n = 10L), + arr = c(1, 2, 3), txt = "x", empty = list())) + expect_equal(names(f), c("overall.AUROC", "overall.n")) + expect_equal(unname(f), c(0.81, 10)) + expect_length(immLynx:::.flattenMetrics(NULL), 0L) + expect_length(immLynx:::.flattenMetrics(list(1, 2)), 0L) # unnamed +}) + +test_that(".alignPredictions matches by barcode", { + p <- c(0.1, 0.9, 0.5) + a <- immLynx:::.alignPredictions(p, c("c", "a", "b"), c("a", "b", "c")) + expect_equal(names(a), c("a", "b", "c")) + expect_equal(unname(a), c(0.9, 0.5, 0.1)) + + expect_error(immLynx:::.alignPredictions(p, c("a", "b"), c("a", "b")), + "3 predictions for 2 cells") + expect_warning( + immLynx:::.alignPredictions(p, c("x", "y", "z"), c("a", "b", "c")), + "positional alignment") +}) + +test_that(".writeCellColumn leaves NA for unscored cells", { + sce <- mock_clonal_sce() + out <- immLynx:::.writeCellColumn(sce, "prob", c(0.2, 0.8), + c("cell3", "cell1")) + v <- SummarizedExperiment::colData(out)$prob + expect_equal(v[1], 0.8) + expect_equal(v[3], 0.2) + expect_true(is.na(v[2])) + expect_length(v, ncol(sce)) +}) + +# =========================================================================== +# runScXpand argument validation (no Python) +# =========================================================================== + +test_that("runScXpand rejects malformed input objects", { + expect_error(runScXpand(list(a = 1)), + "SingleCellExperiment or Seurat") + + sce <- mock_ensembl_sce(n_genes = 5, n_cells = 4) + colnames(sce) <- c("a", "a", "b", "c") + expect_error(runScXpand(sce), "duplicated cell names") +}) + +test_that("runScXpand validates its scalar arguments", { + sce <- mock_ensembl_sce(n_genes = 5, n_cells = 6) + + expect_error(runScXpand(sce, assay = "logcounts"), "not found") + expect_error(runScXpand(sce, threshold = 1.5), "threshold") + expect_error(runScXpand(sce, ensembl_min_frac = 2), "ensembl_min_frac") + expect_error(runScXpand(sce, batch_size = -1), "batch_size") + expect_error(runScXpand(sce, num_workers = 1.5), "num_workers") + expect_error(runScXpand(sce, model_name = c("a", "b")), "model_name") + expect_error(runScXpand(sce, clone_col = "nope"), "clone_col") + expect_error(runScXpand(sce, sample_col = "nope"), "sample_col") +}) + +test_that("runScXpand warns rather than errors on an unknown model", { + sce <- mock_ensembl_sce(n_genes = 5, n_cells = 6) + # An unrecognised name must not hard-stop: upstream may add models. Pair + # it with a bad assay so validation aborts before any Python work. + expect_warning( + expect_error(runScXpand(sce, model_name = "future_model", + assay = "nope", verbose = FALSE), + "not found"), + "Unknown model_name") +}) + +test_that("runScXpand refuses to clobber an existing staged H5AD", { + sce <- mock_ensembl_sce(n_genes = 5, n_cells = 6) + wd <- tempfile("scxpand_test_") + dir.create(wd) + writeLines("placeholder", file.path(wd, "adata.h5ad")) + + expect_error(runScXpand(sce, work_dir = wd, verbose = FALSE), + "overwrite = TRUE") + unlink(wd, recursive = TRUE) +}) + +# =========================================================================== +# H5AD staging (zellkonverter only, no scXpand env) +# =========================================================================== + +test_that(".stageScXpandH5AD writes what scXpand expects", { + skip_if_not_installed("zellkonverter") + + sce <- mock_clonal_sce() + SummarizedExperiment::colData(sce)$junk <- + I(replicate(ncol(sce), 1:2, simplify = FALSE)) + + res <- immLynx:::.resolveGeneIDs(sce, verbose = FALSE) + labs <- immLynx:::.deriveExpansionLabels(sce, sample_col = "sample", + verbose = FALSE) + f <- tempfile(fileext = ".h5ad") + immLynx:::.stageScXpandH5AD(sce, res, labs, NULL, "counts", f, + verbose = FALSE) + expect_true(file.exists(f)) + + back <- zellkonverter::readH5AD(f) + expect_true(all(grepl("^ENSG", rownames(back)))) + expect_equal(colnames(back), colnames(sce)) + expect_true(all(c("clone_id_size", "median_clone_size", "expansion") %in% + colnames(SummarizedExperiment::colData(back)))) + # CT* and list-columns must not survive into the H5AD. + expect_false(any(c("CTstrict", "CTaa", "junk") %in% + colnames(SummarizedExperiment::colData(back)))) + unlink(f) +}) + +test_that(".stageScXpandH5AD sums collapsed rows", { + skip_if_not_installed("zellkonverter") + + sce <- mock_ensembl_sce(n_genes = 4, n_cells = 5) + SummarizedExperiment::rowData(sce)$ID <- c("ENSG00000000001", + "ENSG00000000001", + "ENSG00000000002", + "ENSG00000000003") + res <- immLynx:::.resolveGeneIDs(sce, gene_ids = "ID", collapse = "sum", + verbose = FALSE) + f <- tempfile(fileext = ".h5ad") + immLynx:::.stageScXpandH5AD(sce, res, NULL, NULL, "counts", f, + verbose = FALSE) + + back <- zellkonverter::readH5AD(f) + expect_equal(nrow(back), 3L) + raw <- SummarizedExperiment::assay(sce, "counts") + expect_equal( + as.numeric(SummarizedExperiment::assay(back)["ENSG00000000001", ]), + as.numeric(raw[1, ] + raw[2, ])) + unlink(f) +}) + +# =========================================================================== +# Python-dependent tests +# +# These build the scXpand basilisk env (several GB) and download a +# pretrained model from figshare. Opt in with IMMLYNX_TEST_SCXPAND=1. +# =========================================================================== + +test_that("runScXpand adds predictions to a SingleCellExperiment", { + skip_if_no_scxpand_model() + + sce <- mock_ensembl_sce(n_genes = 200, n_cells = 40) + out <- runScXpand(sce, sample_col = "sample", verbose = FALSE) + + expect_s4_class(out, "SingleCellExperiment") + prob <- SummarizedExperiment::colData(out)$scXpand_expansion_prob + expect_length(prob, ncol(sce)) + expect_true(all(prob >= 0 & prob <= 1, na.rm = TRUE)) + + pred <- SummarizedExperiment::colData(out)$scXpand_expansion_pred + expect_true(all(pred %in% c("expanded", "non-expanded", NA))) + + md <- S4Vectors::metadata(out)$scXpand + expect_equal(md$model, "pan_cancer_autoencoder") + expect_equal(md$threshold, 0.5) + expect_true(!is.null(md$gene_mapping)) +}) + +test_that("runScXpand returns a data.frame when asked", { + skip_if_no_scxpand_model() + + sce <- mock_ensembl_sce(n_genes = 200, n_cells = 40) + df <- runScXpand(sce, sample_col = "sample", return_object = FALSE, + verbose = FALSE) + + expect_s3_class(df, "data.frame") + expect_true(all(c("barcode", "expansion_prob", "expansion_pred") %in% + colnames(df))) + expect_true(all(df$barcode %in% colnames(sce))) + expect_false(is.null(attr(df, "scXpand"))) +}) + +test_that("runScXpand computes AUROC against derived labels", { + skip_if_no_scxpand_model() + + sce <- mock_ensembl_sce(n_genes = 200, n_cells = 60) + out <- runScXpand(sce, sample_col = "sample", verbose = FALSE) + md <- S4Vectors::metadata(out)$scXpand + + expect_true(!is.null(md$labels)) + expect_true(is.na(md$auroc_immLynx) || + (md$auroc_immLynx >= 0 && md$auroc_immLynx <= 1)) +}) + +test_that("runScXpand returns a Seurat object for Seurat input", { + skip_if_no_scxpand_model() + skip_if_not_installed("Seurat") + + sce <- mock_ensembl_sce(n_genes = 200, n_cells = 40) + # Seurat chatters about coercing the matrix and about empty layers while + # building the fixture; none of it comes from runScXpand. + seu <- suppressWarnings(Seurat::CreateSeuratObject( + counts = as.matrix(SummarizedExperiment::assay(sce, "counts")), + meta.data = as.data.frame(SummarizedExperiment::colData(sce)))) + + out <- suppressWarnings( + runScXpand(seu, sample_col = "sample", verbose = FALSE)) + expect_s4_class(out, "Seurat") + expect_true("scXpand_expansion_prob" %in% colnames(out[[]])) + expect_false(is.null(out@misc$scXpand)) +}) + +test_that("listScXpandModels(refresh = TRUE) matches the static table", { + skip_if_no_scxpand_env() + + live <- listScXpandModels(refresh = TRUE) + expect_true(all(immLynx:::.SCXPAND_MODELS$model_name %in% live$model_name)) +}) + +test_that("keep_files preserves the staged H5AD", { + skip_if_no_scxpand_model() + + sce <- mock_ensembl_sce(n_genes = 200, n_cells = 30) + wd <- tempfile("scxpand_keep_") + out <- runScXpand(sce, sample_col = "sample", work_dir = wd, + keep_files = TRUE, verbose = FALSE) + expect_true(file.exists(file.path(wd, "adata.h5ad"))) + expect_equal(S4Vectors::metadata(out)$scXpand$work_dir, wd) + unlink(wd, recursive = TRUE) +}) diff --git a/vignettes/advanced_analysis.Rmd b/vignettes/advanced_analysis.Rmd index 5bce2a0..c61ce9b 100644 --- a/vignettes/advanced_analysis.Rmd +++ b/vignettes/advanced_analysis.Rmd @@ -199,6 +199,84 @@ cat("Number of clustcr clusters:", nrow(comparison), "\n") cat("Number of unique clonotypes:", ncol(comparison), "\n") ``` +# Predicting Clonal Expansion from Gene Expression + +Every other tool in immLynx starts from receptor sequences. `runScXpand()` +goes the other way: it uses [scXpand](https://github.com/yizhak-lab-ccg/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: + +- **Raw UMI counts.** Normalized or log-transformed values are rejected. +- **Ensembl gene IDs.** Most Seurat and scRepertoire objects carry symbols. +- **T cells only.** Filter before calling. There is no reliable + programmatic check, and predictions on other cell types are meaningless. + +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. + +```{r scxpand-models, eval=FALSE} +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. + +```{r scxpand-run, eval=FALSE} +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") +``` + +## Benchmarking against the observed repertoire + +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. + +```{r scxpand-benchmark, eval=FALSE} +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_mapping +``` + # Analyzing Selection Pressure Generation probability (Pgen) from OLGA provides a window into diff --git a/vignettes/immLynx_vignette.Rmd b/vignettes/immLynx_vignette.Rmd index 515ce3c..cf4bf8a 100644 --- a/vignettes/immLynx_vignette.Rmd +++ b/vignettes/immLynx_vignette.Rmd @@ -100,6 +100,12 @@ 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. +Two functions use their own environments rather than that one. +`exportToScanpy()` needs the scanpy/scirpy stack, and `runScXpand()` needs +Python 3.11 with a current PyTorch, which conflicts with the pinned +versions above. Each builds on first use. The scXpand environment is the +largest, at several gigabytes. + # Quick Start ## Loading the Package and Example Data From e1f354e7ecd2bb08608eab826f6e6c5ce0af147c Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 07:45:22 -0500 Subject: [PATCH 02/11] deeptcr addition #8 --- R/basilisk.R | 59 ++++++++++++++++- R/runDeepTCR.R | 176 +++++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 234 insertions(+), 1 deletion(-) create mode 100644 R/runDeepTCR.R diff --git a/R/basilisk.R b/R/basilisk.R index d2a2e67..f2598c4 100644 --- a/R/basilisk.R +++ b/R/basilisk.R @@ -58,7 +58,12 @@ scanpyExportEnv <- basilisk::BasiliskEnvironment( "h5py" ), pip = c( - "anndata>=0.8", + # Pinned exactly rather than ">=0.8": basilisk hands the pip vector to + # system2(), which runs it through an unquoted shell, so ">=0.8" is + # parsed as a redirect to a file named "=0.8". That both dropped the + # version floor and littered the working directory. 0.11.4 is the + # version pip already resolved to under python=3.10. + "anndata==0.11.4", "scanpy", "muon", "scirpy" @@ -95,3 +100,55 @@ scXpandEnv <- basilisk::BasiliskEnvironment( "scxpand==0.4.6" ) ) + +# Separate environment for DeepTCR. DeepTCR pins its whole scientific stack +# (numpy 1.23.5, pandas 1.5.3, scipy 1.10.1, TensorFlow 2.12) with "==", so it +# cannot share immLynxEnv. +# +# Three constraints drove this layout, each verified by building the +# environment and training a VAE: +# +# 1. python=3.10, not 3.11. DeepTCR pins biopython==1.76, which ships no +# wheel past cp38 and none for macOS arm64, so it compiles from source. +# Its C extension assigns to Py_TYPE(), which CPython 3.11 made a hard +# error. It compiles cleanly on 3.10. +# +# 2. The stack comes from conda rather than pip so that TensorFlow arrives +# without Apple's metal plugin. On macOS, DeepTCR's requirements pull +# tensorflow-metal==0.8.0, which is built against TF 2.11 and aborts the +# process at import under TF 2.12 with "platform is already registered +# with name: METAL". conda-forge's tensorflow has no such plugin. +# +# 3. biopython stays on pip. conda-forge's oldest osx-arm64 build is 1.78, +# which removed Bio.Alphabet, and DeepTCR still imports it. +# +# "--no-deps" therefore applies to the whole pip step: every dependency is +# already satisfied by conda, and it also keeps DeepTCR's jupyterlab and +# notebook requirements out of the environment. +deepTCREnv <- basilisk::BasiliskEnvironment( + envname = "deepTCREnv", + pkgname = "immLynx", + packages = c( + "python=3.10", + "tensorflow=2.12", + "numpy=1.23.5", + "pandas=1.5.3", + "scipy=1.10.1", + "h5py=3.8.0", + "scikit-learn=1.2.2", + "matplotlib-base=3.7.2", + "seaborn=0.12.2", + "umap-learn", + "networkx", + "tqdm", + "psutil" + ), + pip = c( + "--no-deps", + "DeepTCR==2.1.29", + "biopython==1.76", + "logomaker==0.8", + "distinctipy==1.2.1", + "python-louvain==0.16" + ) +) diff --git a/R/runDeepTCR.R b/R/runDeepTCR.R new file mode 100644 index 0000000..fae9964 --- /dev/null +++ b/R/runDeepTCR.R @@ -0,0 +1,176 @@ +#' Featurize TCR Sequences with the DeepTCR Variational Autoencoder +#' +#' @description Extracts CDR3 amino acid sequences from a +#' SingleCellExperiment object and learns a low-dimensional representation +#' using DeepTCR's unsupervised variational autoencoder (DeepTCR_U). +#' +#' @param input A SingleCellExperiment object containing scRepertoire TCR data. +#' @param chains Which chain(s) to featurize: "TRB", "TRA", or "both". +#' Default is "TRB". +#' @param latent_dim Width of the VAE latent space. Default is 256. +#' @param use_genes Logical. If TRUE, includes V and J gene usage alongside the +#' CDR3 sequence. Default is FALSE, sequence only. +#' @param stop_criterion Training convergence threshold passed to Train_VAE. +#' Default is 0.01. +#' @param combine_chains Logical. If TRUE and chains="both", concatenates alpha +#' and beta sequences with "_". Default is FALSE. +#' @param reduction_name Name for the dimensional reduction. +#' Default is "tcr_deeptcr". +#' @param reduction_key Key prefix for the reduction columns. +#' Default is "DeepTCR_". +#' @param return_object Logical. If TRUE, adds features as a dimensional +#' reduction. If FALSE, returns the raw features. Default is TRUE. +#' @param seed Optional integer passed to Train_VAE as `graph_seed` for +#' reproducible training. Default is NULL. +#' @param verbose Logical. If FALSE, suppresses DeepTCR's training output. +#' Default is TRUE. +#' @param ... Additional arguments passed to calculate.deepTCR(). +#' +#' @return If return_object=TRUE, the input object with features added as the +#' `reduction_name` reduction. Cells with no sequence for the requested chain +#' get NA. If return_object=FALSE, a list with `features` (one row per unique +#' sequence), `sequences`, `barcodes`, and `explained_variance_ratio`. +#' +#' @details The VAE is trained on unique sequences rather than on cells. +#' DeepTCR's own documentation notes that `Load_Data` does not merge +#' identical amino acid sequences, so passing a redundant repertoire would +#' both waste training time and weight the loss toward expanded clones. +#' Features are expanded back onto cells afterward, so two cells sharing a +#' CDR3 share a feature vector. +#' +#' DeepTCR runs in its own basilisk environment (`deepTCREnv`) because it +#' pins TensorFlow 2.12 and an older numpy, pandas, and scipy stack that +#' cannot coexist with the main immLynx environment. The first call builds +#' that environment, which is a large download. +#' +#' `latent_dim` sets the requested width. DeepTCR prunes uninformative latent +#' features during training, so the returned matrix may be narrower. +#' +#' @export +#' @importFrom methods is +#' @importFrom SummarizedExperiment colData colData<- +#' +#' @examples +#' data(immLynx_example) +#' \donttest{ +#' # Learn a 64-dimensional representation of the TRB repertoire +#' sce <- runDeepTCR(immLynx_example, +#' chains = "TRB", +#' latent_dim = 64) +#' +#' # Include V and J gene usage +#' sce <- runDeepTCR(immLynx_example, +#' chains = "TRB", +#' use_genes = TRUE) +#' +#' # Get the raw feature matrix +#' res <- runDeepTCR(immLynx_example, +#' chains = "TRB", +#' return_object = FALSE) +#' dim(res$features) +#' } +runDeepTCR <- function(input, + chains = c("TRB", "TRA", "both"), + latent_dim = 256, + use_genes = FALSE, + stop_criterion = 0.01, + combine_chains = FALSE, + reduction_name = "tcr_deeptcr", + reduction_key = "DeepTCR_", + return_object = TRUE, + seed = NULL, + verbose = TRUE, + ...) { + + chains <- match.arg(chains) + + if (!methods::is(input, "SingleCellExperiment")) { + stop("Input must be a SingleCellExperiment object") + } + + if (!is.numeric(latent_dim) || length(latent_dim) != 1L || + is.na(latent_dim) || latent_dim < 1) { + stop("latent_dim must be a single positive integer") + } + latent_dim <- as.integer(latent_dim) + + message("Extracting TCR sequences from object...") + seq_map <- .extractChainSeqs(input, chains = chains, + combine_chains = combine_chains, + with_genes = use_genes) + + keep <- !duplicated(seq_map$sequences) + unique_seqs <- seq_map$sequences[keep] + + if (length(unique_seqs) == 0L) { + stop("No CDR3 sequences found for the requested chain(s).") + } + + message("Training DeepTCR VAE on ", length(unique_seqs), + " unique sequences (latent_dim = ", latent_dim, ")...") + + res <- calculate.deepTCR( + sequences = unique_seqs, + v_genes = if (use_genes) seq_map$v[keep] else NULL, + j_genes = if (use_genes) seq_map$j[keep] else NULL, + latent_dim = latent_dim, + stop_criterion = stop_criterion, + seed = seed, + verbose = verbose, + ... + ) + + features <- res$features + if (is.null(dim(features))) { + features <- matrix(features, nrow = length(unique_seqs)) + } + message("Learned ", ncol(features), " latent feature(s).") + + if (!return_object) { + return(list( + features = features, + sequences = unique_seqs, + barcodes = seq_map$barcodes, + explained_variance_ratio = res$explained_variance_ratio + )) + } + + # Expand unique-sequence features onto the cells that carry them. + cell_features <- .deeptcrFeaturesToCells(features, unique_seqs, + seq_map$sequences) + + full <- matrix(NA_real_, + nrow = ncol(input), + ncol = ncol(cell_features), + dimnames = list(colnames(input), + paste0(reduction_key, + seq_len(ncol(cell_features))))) + full[seq_map$barcodes, ] <- cell_features + + SingleCellExperiment::reducedDim(input, reduction_name) <- full + + message("Features added as '", reduction_name, "' reduction") + input +} + + +#' Expand per-sequence features onto cells +#' +#' @param features Matrix with one row per unique sequence. +#' @param unique_seqs Character vector matching the rows of `features`. +#' @param cell_seqs Character vector of the sequence carried by each cell. +#' @return A matrix with one row per entry of `cell_seqs`. +#' @keywords internal +.deeptcrFeaturesToCells <- function(features, unique_seqs, cell_seqs) { + + if (nrow(features) != length(unique_seqs)) { + stop("features must have as many rows as unique sequences") + } + + idx <- match(cell_seqs, unique_seqs) + if (anyNA(idx)) { + stop("every cell sequence must appear in the unique sequence set") + } + + features[idx, , drop = FALSE] +} From a613bc5bbcb3872f0e7490fc7bf74caefd0f43a2 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 07:46:00 -0500 Subject: [PATCH 03/11] XT support (symdel) #6 --- R/calculate_helpers.R | 141 ++++++++++++++++++++++++++++++ R/runSymdelNeighbors.R | 191 +++++++++++++++++++++++++++++++++++++++++ R/utils.R | 68 +++++++++++++++ 3 files changed, 400 insertions(+) create mode 100644 R/runSymdelNeighbors.R diff --git a/R/calculate_helpers.R b/R/calculate_helpers.R index a2a5f6d..2bb9639 100644 --- a/R/calculate_helpers.R +++ b/R/calculate_helpers.R @@ -334,3 +334,144 @@ calculate.sonia <- function(data_folder, return(results) } + +#' Find Near-Neighbor Sequences using pyrepseq symdel +#' +#' @description Internal function that calls pyrepseq's symmetric deletion +#' lookup via basilisk. Returns raw (i, j, distance) triplets using 0-based +#' indices into `sequences`; both directions of each pair are present. +#' @param sequences Character vector of unique CDR3 amino acid sequences +#' @param max_edits Integer maximum edit distance defining a neighbor +#' @param max_returns Integer maximum neighbors per sequence, or NULL +#' @param n_cpu Integer number of processes +#' @param ... Additional arguments passed to pyrepseq.nn.symdel +#' @return An n x 3 integer matrix of (i, j, distance), or NULL if no pairs +#' @keywords internal +calculate.symdel <- function(sequences, + max_edits = 1, + max_returns = NULL, + n_cpu = 1, + ...) { + + .run_in_basilisk(function(sequences, max_edits, max_returns, n_cpu, extra) { + nn <- reticulate::import("pyrepseq.nn", convert = FALSE) + np <- reticulate::import("numpy", convert = FALSE) + + args <- c( + list( + reticulate::r_to_py(as.character(sequences)), + max_edits = as.integer(max_edits), + max_returns = if (is.null(max_returns)) { + NULL + } else { + as.integer(max_returns) + }, + n_cpu = as.integer(n_cpu), + output_type = "triplets" + ), + extra + ) + + triplets <- do.call(nn$symdel, args) + + # Marshal in one array rather than element-by-element: a large repertoire + # can return millions of triplets and per-tuple conversion is very slow. + arr <- np$asarray(triplets, dtype = np$int64) + if (reticulate::py_to_r(arr$size) == 0) { + return(NULL) + } + reticulate::py_to_r(arr$reshape(reticulate::tuple(-1L, 3L))) + + }, sequences = sequences, max_edits = max_edits, + max_returns = max_returns, n_cpu = n_cpu, extra = list(...)) +} + +# Shared basilisk helper for the DeepTCR environment. +.run_in_deeptcr <- function(FUN, ...) { + proc <- basilisk::basiliskStart(deepTCREnv) + on.exit(basilisk::basiliskStop(proc)) + basilisk::basiliskRun(proc, FUN, ...) +} + +#' Train a DeepTCR Variational Autoencoder +#' +#' @description Internal function that trains DeepTCR_U via basilisk and +#' returns the learned latent features. +#' @param sequences Character vector of unique CDR3 amino acid sequences +#' @param v_genes Optional character vector of V genes, aligned to sequences +#' @param j_genes Optional character vector of J genes, aligned to sequences +#' @param latent_dim Integer width of the latent space +#' @param stop_criterion Numeric training convergence threshold +#' @param seed Optional integer graph seed for reproducible training +#' @param verbose Logical, whether to show DeepTCR training output +#' @param ... Additional arguments passed to Train_VAE +#' @return A list with `features` and `explained_variance_ratio` +#' @keywords internal +calculate.deepTCR <- function(sequences, + v_genes = NULL, + j_genes = NULL, + latent_dim = 256, + stop_criterion = 0.01, + seed = NULL, + verbose = TRUE, + ...) { + + .run_in_deeptcr(function(sequences, v_genes, j_genes, latent_dim, + stop_criterion, seed, verbose, extra) { + + np <- reticulate::import("numpy", convert = FALSE) + dt <- reticulate::import("DeepTCR.DeepTCR", convert = FALSE) + + # DeepTCR_U writes model checkpoints and result folders relative to the + # working directory, keyed on the object Name. Contain that in a temp dir + # and restore the caller's working directory afterwards. + workdir <- tempfile("deeptcr_") + dir.create(workdir, recursive = TRUE) + oldwd <- setwd(workdir) + on.exit({ + setwd(oldwd) + unlink(workdir, recursive = TRUE, force = TRUE) + }, add = TRUE) + + obj <- dt$DeepTCR_U("immLynx_vae") + + load_args <- list( + beta_sequences = np$array(reticulate::r_to_py(as.character(sequences))) + ) + if (!is.null(v_genes)) { + load_args$v_beta <- np$array(reticulate::r_to_py(as.character(v_genes))) + } + if (!is.null(j_genes)) { + load_args$j_beta <- np$array(reticulate::r_to_py(as.character(j_genes))) + } + do.call(obj$Load_Data, load_args) + + train_args <- c( + list( + latent_dim = as.integer(latent_dim), + stop_criterion = as.numeric(stop_criterion), + suppress_output = !isTRUE(verbose) + ), + extra + ) + if (!is.null(seed)) { + train_args$graph_seed <- as.integer(seed) + } + do.call(obj$Train_VAE, train_args) + + features <- reticulate::py_to_r( + np$asarray(obj$features, dtype = np$float64) + ) + + evr <- tryCatch( + reticulate::py_to_r(np$asarray(obj$explained_variance_ratio_, + dtype = np$float64)), + error = function(e) NULL + ) + + list(features = features, explained_variance_ratio = evr) + + }, sequences = sequences, v_genes = v_genes, j_genes = j_genes, + latent_dim = latent_dim, stop_criterion = stop_criterion, + seed = seed, verbose = verbose, extra = list(...)) +} diff --git a/R/runSymdelNeighbors.R b/R/runSymdelNeighbors.R new file mode 100644 index 0000000..707089f --- /dev/null +++ b/R/runSymdelNeighbors.R @@ -0,0 +1,191 @@ +#' Find Near-Neighbor CDR3 Sequences by Symmetric Deletion Lookup +#' +#' @description Extracts CDR3 amino acid sequences from a +#' SingleCellExperiment object and finds all pairs within a given edit +#' distance using the symmetric deletion lookup ("symdel") algorithm from +#' pyrepseq. This is the same algorithm described in the XT-neighbor +#' preprint, running on CPU rather than GPU. +#' +#' @param input A SingleCellExperiment object containing scRepertoire TCR data. +#' @param chains Which chain(s) to search: "TRB", "TRA", or "both". +#' Default is "TRB". +#' @param max_edits Maximum edit distance defining a neighbor. Default is 1. +#' @param max_returns Maximum number of neighbors to return per sequence. +#' Default is NULL, meaning no limit. +#' @param n_cpu Number of CPU processes for the search. Default is 1. +#' @param combine_chains Logical. If TRUE and chains="both", concatenates alpha +#' and beta sequences with "_" before searching. Default is FALSE. +#' @param return_object Logical. If TRUE, adds a per-cell neighbor count to the +#' input object. If FALSE, returns the neighbor edge list. Default is TRUE. +#' @param column_prefix Prefix for the new metadata column. Default is "symdel", +#' producing a "symdel_degree" column. +#' @param ... Additional arguments passed to calculate.symdel(). +#' +#' @return If return_object=TRUE, the input object with a `_degree` +#' column added to colData giving the number of neighbors found for each +#' cell's CDR3 sequence. Cells with no sequence for the requested chain get +#' NA. If return_object=FALSE, a data.frame with columns `from_seq`, +#' `to_seq`, and `distance`, one row per undirected neighbor pair. +#' +#' @details The search runs over unique sequences rather than over cells. +#' Repertoires carry heavy clonal redundancy, so this is both substantially +#' faster and what the symdel algorithm expects. Results are mapped back to +#' cells afterward. +#' +#' Degree is a neighborhood-density measure. High-degree sequences sit inside +#' dense clusters of similar receptors, which is the signal used to identify +#' convergent recombination and antigen-driven expansion. For explicit +#' cluster assignments, see \code{\link{runClustTCR}}. +#' +#' @export +#' @importFrom immApex getIR +#' @importFrom methods is +#' @importFrom SummarizedExperiment colData colData<- +#' +#' @examples +#' data(immLynx_example) +#' \donttest{ +#' # Neighbors within one edit on the TRB chain +#' sce <- runSymdelNeighbors(immLynx_example, chains = "TRB") +#' +#' # Widen the search to two edits +#' sce <- runSymdelNeighbors(immLynx_example, +#' chains = "TRB", +#' max_edits = 2) +#' +#' # Get the neighbor edge list instead +#' edges <- runSymdelNeighbors(immLynx_example, +#' chains = "TRB", +#' return_object = FALSE) +#' head(edges) +#' } +runSymdelNeighbors <- function(input, + chains = c("TRB", "TRA", "both"), + max_edits = 1, + max_returns = NULL, + n_cpu = 1, + combine_chains = FALSE, + return_object = TRUE, + column_prefix = "symdel", + ...) { + + chains <- match.arg(chains) + + if (!methods::is(input, "SingleCellExperiment")) { + stop("Input must be a SingleCellExperiment object") + } + + if (!is.numeric(max_edits) || length(max_edits) != 1L || + is.na(max_edits) || max_edits < 1) { + stop("max_edits must be a single positive integer") + } + max_edits <- as.integer(max_edits) + + message("Extracting TCR sequences from object...") + seq_map <- .extractChainSeqs(input, chains = chains, + combine_chains = combine_chains) + + unique_seqs <- unique(seq_map$sequences[!is.na(seq_map$sequences)]) + if (length(unique_seqs) == 0L) { + stop("No CDR3 sequences found for the requested chain(s).") + } + + message("Searching ", length(unique_seqs), " unique sequences ", + "(max_edits = ", max_edits, ")...") + + triplets <- calculate.symdel(unique_seqs, + max_edits = max_edits, + max_returns = max_returns, + n_cpu = n_cpu, + ...) + + edges <- .symdelTripletsToEdges(triplets, unique_seqs) + message("Found ", nrow(edges), " neighbor pair(s).") + + if (!return_object) { + return(edges) + } + + deg <- .symdelDegree(edges, unique_seqs) + + # Map degrees from unique sequences back onto cells. Cells with no sequence + # for the requested chain stay NA. + cell_deg <- unname(deg[match(seq_map$sequences, unique_seqs)]) + + col_vec <- rep(NA_integer_, ncol(input)) + names(col_vec) <- colnames(input) + col_vec[seq_map$barcodes] <- cell_deg + + colData(input)[[paste0(column_prefix, "_degree")]] <- unname(col_vec) + input +} + + +#' Convert symdel triplets into an undirected edge data.frame +#' +#' @description symdel returns (i, j, distance) using 0-based indices and emits +#' both directions of every pair. This collapses them to one row per +#' undirected pair and drops any self-pair. +#' +#' @param triplets A list of length-3 numeric vectors, or an n x 3 matrix. +#' @param unique_seqs Character vector the indices refer to. +#' @return A data.frame with `from_seq`, `to_seq`, and `distance`. +#' @keywords internal +.symdelTripletsToEdges <- function(triplets, unique_seqs) { + + empty <- data.frame(from_seq = character(0), to_seq = character(0), + distance = integer(0), stringsAsFactors = FALSE) + + if (is.null(triplets) || length(triplets) == 0L) { + return(empty) + } + + m <- if (is.matrix(triplets)) triplets else do.call(rbind, triplets) + if (is.null(m) || nrow(m) == 0L) { + return(empty) + } + + # 0-based Python indices to 1-based R indices. + i <- as.integer(m[, 1L]) + 1L + j <- as.integer(m[, 2L]) + 1L + d <- as.integer(m[, 3L]) + + # Drop self-pairs, then canonicalize orientation so the two directions of a + # pair collapse onto the same row. + keep <- i != j + i <- i[keep]; j <- j[keep]; d <- d[keep] + + lo <- pmin(i, j) + hi <- pmax(i, j) + dedup <- !duplicated(paste(lo, hi)) + + out <- data.frame( + from_seq = unique_seqs[lo[dedup]], + to_seq = unique_seqs[hi[dedup]], + distance = d[dedup], + stringsAsFactors = FALSE + ) + rownames(out) <- NULL + out +} + + +#' Neighbor count per unique sequence +#' +#' @param edges A data.frame from .symdelTripletsToEdges(). +#' @param unique_seqs Character vector of every sequence searched. +#' @return A named integer vector of degrees, in the order of `unique_seqs`. +#' @keywords internal +.symdelDegree <- function(edges, unique_seqs) { + + deg <- integer(length(unique_seqs)) + names(deg) <- unique_seqs + + if (nrow(edges) > 0L) { + # Each undirected edge contributes one to both endpoints. + counts <- table(c(edges$from_seq, edges$to_seq)) + deg[names(counts)] <- as.integer(counts) + } + + deg +} diff --git a/R/utils.R b/R/utils.R index 142a98f..a5c9bc1 100644 --- a/R/utils.R +++ b/R/utils.R @@ -558,3 +558,71 @@ summarizeTCRrepertoire <- function(input, } + +#' Extract CDR3 sequences, barcodes, and optionally V/J genes +#' +#' @description Shared chain-extraction logic for the neighbor-search and +#' featurization wrappers. Returns one entry per cell that carries a usable +#' sequence for the requested chain(s). +#' +#' @param input A SingleCellExperiment object. +#' @param chains "TRB", "TRA", or "both". +#' @param combine_chains Logical. When chains="both", concatenate TRA and TRB +#' with "_" and keep only cells carrying both. +#' @param with_genes Logical. If TRUE, also return `v` and `j` gene vectors +#' aligned to `sequences`. +#' @return A list with `sequences` and `barcodes`, plus `v` and `j` when +#' `with_genes = TRUE`. +#' @keywords internal +.extractChainSeqs <- function(input, chains, combine_chains = FALSE, + with_genes = FALSE) { + + if (chains != "both") { + dat <- immApex::getIR(input, chains = chains) + keep <- !is.na(dat$cdr3_aa) + out <- list(sequences = dat$cdr3_aa[keep], barcodes = dat$barcode[keep]) + if (with_genes) { + out$v <- dat$v[keep] + out$j <- dat$j[keep] + } + return(out) + } + + tra <- immApex::getIR(input, chains = "TRA") + trb <- immApex::getIR(input, chains = "TRB") + merged <- merge(tra, trb, by = "barcode", suffixes = c("_TRA", "_TRB"), + all = TRUE) + + if (combine_chains) { + # Only cells carrying both chains can form a paired sequence. + keep <- !is.na(merged$cdr3_aa_TRA) & !is.na(merged$cdr3_aa_TRB) + if (!any(keep)) { + stop("No cells found with both TRA and TRB chains.") + } + out <- list( + sequences = paste0(merged$cdr3_aa_TRA[keep], "_", + merged$cdr3_aa_TRB[keep]), + barcodes = merged$barcode[keep] + ) + if (with_genes) { + # Gene usage is only meaningful per chain; report the beta genes, which + # carry most of the specificity signal. + out$v <- merged$v_TRB[keep] + out$j <- merged$j_TRB[keep] + } + return(out) + } + + # Not concatenating: take beta where present, otherwise alpha, keeping one + # row per cell so results map back cleanly. + use_trb <- !is.na(merged$cdr3_aa_TRB) + seqs <- ifelse(use_trb, merged$cdr3_aa_TRB, merged$cdr3_aa_TRA) + keep <- !is.na(seqs) + + out <- list(sequences = seqs[keep], barcodes = merged$barcode[keep]) + if (with_genes) { + out$v <- ifelse(use_trb, merged$v_TRB, merged$v_TRA)[keep] + out$j <- ifelse(use_trb, merged$j_TRB, merged$j_TRA)[keep] + } + out +} From e2e81a80fd899cd47ee4d8d8c91e3dcc56945323 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 07:46:14 -0500 Subject: [PATCH 04/11] Update exportToScanpy.R --- R/exportToScanpy.R | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/R/exportToScanpy.R b/R/exportToScanpy.R index ec0f3d4..308c8e7 100644 --- a/R/exportToScanpy.R +++ b/R/exportToScanpy.R @@ -408,7 +408,18 @@ exportToScanpy <- function(input, stop("Seurat input requires the Seurat package. ", "Install with: install.packages('Seurat')", call. = FALSE) } - return(Seurat::as.SingleCellExperiment(input)) + # Seurat emits one warning per empty layer ("Layer 'data' is empty") when + # converting an object that only carries counts. That is the normal state + # for a freshly created object and is not actionable here, so muffle just + # those and let every other warning through. + return(withCallingHandlers( + Seurat::as.SingleCellExperiment(input), + warning = function(w) { + if (grepl("^Layer '.*' is empty", conditionMessage(w))) { + invokeRestart("muffleWarning") + } + } + )) } stop("input must be a SingleCellExperiment or Seurat object", call. = FALSE) From 1169043397f2d4bd164f5932bf01d1029aa02993 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 07:46:41 -0500 Subject: [PATCH 05/11] update documentation --- NAMESPACE | 2 + NEWS.md | 18 ++++ inst/WORDLIST | 147 ++++++++++++++++++------------- man/calculate.deepTCR.Rd | 42 +++++++++ man/calculate.symdel.Rd | 28 ++++++ man/dot-extractChainSeqs.Rd | 29 ++++++ man/dot-symdelDegree.Rd | 20 +++++ man/dot-symdelTripletsToEdges.Rd | 22 +++++ man/runDeepTCR.Rd | 102 +++++++++++++++++++++ man/runSymdelNeighbors.Rd | 85 ++++++++++++++++++ 10 files changed, 433 insertions(+), 62 deletions(-) create mode 100644 man/calculate.deepTCR.Rd create mode 100644 man/calculate.symdel.Rd create mode 100644 man/dot-extractChainSeqs.Rd create mode 100644 man/dot-symdelDegree.Rd create mode 100644 man/dot-symdelTripletsToEdges.Rd create mode 100644 man/runDeepTCR.Rd create mode 100644 man/runSymdelNeighbors.Rd diff --git a/NAMESPACE b/NAMESPACE index 799e9eb..80c5d41 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -8,12 +8,14 @@ export(huggingModel) export(listScXpandModels) export(proteinEmbeddings) export(runClustTCR) +export(runDeepTCR) export(runEmbeddings) export(runHLAassociation) export(runMetaclonotypist) export(runOLGA) export(runScXpand) export(runSoNNia) +export(runSymdelNeighbors) export(runTCRdist) export(summarizeTCRrepertoire) export(tokenizeSequences) diff --git a/NEWS.md b/NEWS.md index 68f4172..fa98fb2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -33,6 +33,24 @@ `runScXpand()` rewrites the host. The rewrite becomes a no-op once upstream fixes its URLs. +* Added `runSymdelNeighbors()` for near-neighbor CDR3 search by symmetric + deletion lookup, backed by `pyrepseq.nn.symdel` in the existing + `immLynxEnv`. Returns either a neighbor edge list or a per-cell neighbor + count. This closes the XT-neighbor request (#6) without a GPU dependency: + XT-neighbor is CUDA-only and its own documentation redirects users to the + CPU implementation of the same algorithm. + +* Added `runDeepTCR()` for unsupervised VAE featurization of CDR3 sequences + via DeepTCR (#8), writing features to a dimensional reduction. Runs in a new + `deepTCREnv` basilisk environment. The dependency conflict that originally + blocked this was resolved upstream in DeepTCR 2.1.29. + +* Fixed `scanpyExportEnv` declaring `anndata>=0.8`. basilisk passes the `pip` + vector through an unquoted shell, so `>=0.8` was parsed as a redirect: the + version floor was silently dropped and a stray file named `=0.8` was written + to the working directory. Now pinned to `anndata==0.11.4`, the version pip + already resolved to. + # immLynx 1.1.2 * Added `exportToScanpy()` to write a `SingleCellExperiment` or `Seurat` diff --git a/inst/WORDLIST b/inst/WORDLIST index 48aa6e1..7653e57 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,151 +1,174 @@ -AIRR -AUROC -AutoModel -BCR -BioC -Bioc -CDR -CMD -Clonotypes -Codecov -DBSCAN -DataLoader -ESM -Embeddings -Ensembl -GEX -GLIPH -HLA -Hyperexpanded -IGH -IGK -IGL -IMGT -Isacchini -LightGBM -MCL -Metaclone -Metaclones -Metaclonotypist -MuData -OLGA's -ORCID -Pgen -Pretrained -PyTorch -SCE -SCEPTR -Scalability -Sethna -SingleCellExperiment -Symdel -TCR -TCRdist -TCRs -TRA -TRB -TRD -TRG -TSV -Tokenize -Tokenizer -UMAP -UMI -VDJ -Valkiers aa +AIRR al +AUROC +Autoencoder autoencoder +AutoModel barcode barcodes +BCR benchmarked benchmarking bfloat +Bioc +BioC cdr +CDR clonality clonotype +Clonotypes clonotypes clonotyping cls +clustcr clusTCR clusTCR's -clustcr +CMD +Codecov +colData complementarity conda cpu cuda +DataLoader dbscan +DBSCAN deduplicates +DeepTCR +deeptcr +deepTCR +DeepTCR's +deepTCREnv disambiguated downloader +Embeddings embeddings +Ensembl erroring esm +ESM et extractTCRdata facebook +featurization +Featurize +featurize figshare generateOLGA +GEX +GLIPH heterodimeric +HLA https humanIGH humanTRA humanTRB +Hyperexpanded +IGH +IGK +IGL +IMGT immApex immGLIPH immunoGlobulin +Isacchini leiden lexicographically lightgbm +LightGBM louvain mcl +MCL +Metaclone metaclone +Metaclones metaclones metaclonotypist +Metaclonotypist minibatched mouseTRB mps +MuData muon ncbi nih nlm +nn nt numpy olga +OLGA's +ORCID pgen +Pgen pooch pre +preprint +Pretrained pretrained proteinEmbeddings pubmed pw +pyrepseq +pyrepseq's +PyTorch qimmuno quadratically reformats roxygen runMetaclonotypist -scRNA -scRepertoire -scRepertoire's -scXpand -scXpand's +Scalability scanpy scater +SCE sceptr +SCEPTR +scipy scirpy scirpy's scran +scRepertoire +scRepertoire's +scRNA scxpand -soNNia +scXpand +scXpand's +seqs +Sethna +SingleCellExperiment sonia +soNNia sonnia +Symdel +symdel +symdelTripletsToEdges +TCR tcr +TCRdist tcrdist +TCRs tf thymic +Tokenize tokenized tokenizer +Tokenizer torch +TRA transcriptome +TRB +TRD +TRG +TSV +UMAP +UMI +UniProt unmapped unselected +VAE +Valkiers +Variational +variational +VDJ +XT diff --git a/man/calculate.deepTCR.Rd b/man/calculate.deepTCR.Rd new file mode 100644 index 0000000..09d2007 --- /dev/null +++ b/man/calculate.deepTCR.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calculate_helpers.R +\name{calculate.deepTCR} +\alias{calculate.deepTCR} +\title{Train a DeepTCR Variational Autoencoder} +\usage{ +calculate.deepTCR( + sequences, + v_genes = NULL, + j_genes = NULL, + latent_dim = 256, + stop_criterion = 0.01, + seed = NULL, + verbose = TRUE, + ... +) +} +\arguments{ +\item{sequences}{Character vector of unique CDR3 amino acid sequences} + +\item{v_genes}{Optional character vector of V genes, aligned to sequences} + +\item{j_genes}{Optional character vector of J genes, aligned to sequences} + +\item{latent_dim}{Integer width of the latent space} + +\item{stop_criterion}{Numeric training convergence threshold} + +\item{seed}{Optional integer graph seed for reproducible training} + +\item{verbose}{Logical, whether to show DeepTCR training output} + +\item{...}{Additional arguments passed to Train_VAE} +} +\value{ +A list with `features` and `explained_variance_ratio` +} +\description{ +Internal function that trains DeepTCR_U via basilisk and + returns the learned latent features. +} +\keyword{internal} diff --git a/man/calculate.symdel.Rd b/man/calculate.symdel.Rd new file mode 100644 index 0000000..8a98bbf --- /dev/null +++ b/man/calculate.symdel.Rd @@ -0,0 +1,28 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/calculate_helpers.R +\name{calculate.symdel} +\alias{calculate.symdel} +\title{Find Near-Neighbor Sequences using pyrepseq symdel} +\usage{ +calculate.symdel(sequences, max_edits = 1, max_returns = NULL, n_cpu = 1, ...) +} +\arguments{ +\item{sequences}{Character vector of unique CDR3 amino acid sequences} + +\item{max_edits}{Integer maximum edit distance defining a neighbor} + +\item{max_returns}{Integer maximum neighbors per sequence, or NULL} + +\item{n_cpu}{Integer number of processes} + +\item{...}{Additional arguments passed to pyrepseq.nn.symdel} +} +\value{ +An n x 3 integer matrix of (i, j, distance), or NULL if no pairs +} +\description{ +Internal function that calls pyrepseq's symmetric deletion + lookup via basilisk. Returns raw (i, j, distance) triplets using 0-based + indices into `sequences`; both directions of each pair are present. +} +\keyword{internal} diff --git a/man/dot-extractChainSeqs.Rd b/man/dot-extractChainSeqs.Rd new file mode 100644 index 0000000..f388770 --- /dev/null +++ b/man/dot-extractChainSeqs.Rd @@ -0,0 +1,29 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{.extractChainSeqs} +\alias{.extractChainSeqs} +\title{Extract CDR3 sequences, barcodes, and optionally V/J genes} +\usage{ +.extractChainSeqs(input, chains, combine_chains = FALSE, with_genes = FALSE) +} +\arguments{ +\item{input}{A SingleCellExperiment object.} + +\item{chains}{"TRB", "TRA", or "both".} + +\item{combine_chains}{Logical. When chains="both", concatenate TRA and TRB +with "_" and keep only cells carrying both.} + +\item{with_genes}{Logical. If TRUE, also return `v` and `j` gene vectors +aligned to `sequences`.} +} +\value{ +A list with `sequences` and `barcodes`, plus `v` and `j` when + `with_genes = TRUE`. +} +\description{ +Shared chain-extraction logic for the neighbor-search and + featurization wrappers. Returns one entry per cell that carries a usable + sequence for the requested chain(s). +} +\keyword{internal} diff --git a/man/dot-symdelDegree.Rd b/man/dot-symdelDegree.Rd new file mode 100644 index 0000000..1e07e13 --- /dev/null +++ b/man/dot-symdelDegree.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/runSymdelNeighbors.R +\name{.symdelDegree} +\alias{.symdelDegree} +\title{Neighbor count per unique sequence} +\usage{ +.symdelDegree(edges, unique_seqs) +} +\arguments{ +\item{edges}{A data.frame from .symdelTripletsToEdges().} + +\item{unique_seqs}{Character vector of every sequence searched.} +} +\value{ +A named integer vector of degrees, in the order of `unique_seqs`. +} +\description{ +Neighbor count per unique sequence +} +\keyword{internal} diff --git a/man/dot-symdelTripletsToEdges.Rd b/man/dot-symdelTripletsToEdges.Rd new file mode 100644 index 0000000..72ae50c --- /dev/null +++ b/man/dot-symdelTripletsToEdges.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/runSymdelNeighbors.R +\name{.symdelTripletsToEdges} +\alias{.symdelTripletsToEdges} +\title{Convert symdel triplets into an undirected edge data.frame} +\usage{ +.symdelTripletsToEdges(triplets, unique_seqs) +} +\arguments{ +\item{triplets}{A list of length-3 numeric vectors, or an n x 3 matrix.} + +\item{unique_seqs}{Character vector the indices refer to.} +} +\value{ +A data.frame with `from_seq`, `to_seq`, and `distance`. +} +\description{ +symdel returns (i, j, distance) using 0-based indices and emits + both directions of every pair. This collapses them to one row per + undirected pair and drops any self-pair. +} +\keyword{internal} diff --git a/man/runDeepTCR.Rd b/man/runDeepTCR.Rd new file mode 100644 index 0000000..614f2bc --- /dev/null +++ b/man/runDeepTCR.Rd @@ -0,0 +1,102 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/runDeepTCR.R +\name{runDeepTCR} +\alias{runDeepTCR} +\title{Featurize TCR Sequences with the DeepTCR Variational Autoencoder} +\usage{ +runDeepTCR( + input, + chains = c("TRB", "TRA", "both"), + latent_dim = 256, + use_genes = FALSE, + stop_criterion = 0.01, + combine_chains = FALSE, + reduction_name = "tcr_deeptcr", + reduction_key = "DeepTCR_", + return_object = TRUE, + seed = NULL, + verbose = TRUE, + ... +) +} +\arguments{ +\item{input}{A SingleCellExperiment object containing scRepertoire TCR data.} + +\item{chains}{Which chain(s) to featurize: "TRB", "TRA", or "both". +Default is "TRB".} + +\item{latent_dim}{Width of the VAE latent space. Default is 256.} + +\item{use_genes}{Logical. If TRUE, includes V and J gene usage alongside the +CDR3 sequence. Default is FALSE, sequence only.} + +\item{stop_criterion}{Training convergence threshold passed to Train_VAE. +Default is 0.01.} + +\item{combine_chains}{Logical. If TRUE and chains="both", concatenates alpha +and beta sequences with "_". Default is FALSE.} + +\item{reduction_name}{Name for the dimensional reduction. +Default is "tcr_deeptcr".} + +\item{reduction_key}{Key prefix for the reduction columns. +Default is "DeepTCR_".} + +\item{return_object}{Logical. If TRUE, adds features as a dimensional +reduction. If FALSE, returns the raw features. Default is TRUE.} + +\item{seed}{Optional integer passed to Train_VAE as `graph_seed` for +reproducible training. Default is NULL.} + +\item{verbose}{Logical. If FALSE, suppresses DeepTCR's training output. +Default is TRUE.} + +\item{...}{Additional arguments passed to calculate.deepTCR().} +} +\value{ +If return_object=TRUE, the input object with features added as the + `reduction_name` reduction. Cells with no sequence for the requested chain + get NA. If return_object=FALSE, a list with `features` (one row per unique + sequence), `sequences`, `barcodes`, and `explained_variance_ratio`. +} +\description{ +Extracts CDR3 amino acid sequences from a + SingleCellExperiment object and learns a low-dimensional representation + using DeepTCR's unsupervised variational autoencoder (DeepTCR_U). +} +\details{ +The VAE is trained on unique sequences rather than on cells. + DeepTCR's own documentation notes that `Load_Data` does not merge + identical amino acid sequences, so passing a redundant repertoire would + both waste training time and weight the loss toward expanded clones. + Features are expanded back onto cells afterward, so two cells sharing a + CDR3 share a feature vector. + + DeepTCR runs in its own basilisk environment (`deepTCREnv`) because it + pins TensorFlow 2.12 and an older numpy, pandas, and scipy stack that + cannot coexist with the main immLynx environment. The first call builds + that environment, which is a large download. + + `latent_dim` sets the requested width. DeepTCR prunes uninformative latent + features during training, so the returned matrix may be narrower. +} +\examples{ +data(immLynx_example) +\donttest{ + # Learn a 64-dimensional representation of the TRB repertoire + sce <- runDeepTCR(immLynx_example, + chains = "TRB", + latent_dim = 64) + + # Include V and J gene usage + sce <- runDeepTCR(immLynx_example, + chains = "TRB", + use_genes = TRUE) + + # Get the raw feature matrix + res <- runDeepTCR(immLynx_example, + chains = "TRB", + return_object = FALSE) + dim(res$features) +} +} diff --git a/man/runSymdelNeighbors.Rd b/man/runSymdelNeighbors.Rd new file mode 100644 index 0000000..f0b034d --- /dev/null +++ b/man/runSymdelNeighbors.Rd @@ -0,0 +1,85 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/runSymdelNeighbors.R +\name{runSymdelNeighbors} +\alias{runSymdelNeighbors} +\title{Find Near-Neighbor CDR3 Sequences by Symmetric Deletion Lookup} +\usage{ +runSymdelNeighbors( + input, + chains = c("TRB", "TRA", "both"), + max_edits = 1, + max_returns = NULL, + n_cpu = 1, + combine_chains = FALSE, + return_object = TRUE, + column_prefix = "symdel", + ... +) +} +\arguments{ +\item{input}{A SingleCellExperiment object containing scRepertoire TCR data.} + +\item{chains}{Which chain(s) to search: "TRB", "TRA", or "both". +Default is "TRB".} + +\item{max_edits}{Maximum edit distance defining a neighbor. Default is 1.} + +\item{max_returns}{Maximum number of neighbors to return per sequence. +Default is NULL, meaning no limit.} + +\item{n_cpu}{Number of CPU processes for the search. Default is 1.} + +\item{combine_chains}{Logical. If TRUE and chains="both", concatenates alpha +and beta sequences with "_" before searching. Default is FALSE.} + +\item{return_object}{Logical. If TRUE, adds a per-cell neighbor count to the +input object. If FALSE, returns the neighbor edge list. Default is TRUE.} + +\item{column_prefix}{Prefix for the new metadata column. Default is "symdel", +producing a "symdel_degree" column.} + +\item{...}{Additional arguments passed to calculate.symdel().} +} +\value{ +If return_object=TRUE, the input object with a `_degree` + column added to colData giving the number of neighbors found for each + cell's CDR3 sequence. Cells with no sequence for the requested chain get + NA. If return_object=FALSE, a data.frame with columns `from_seq`, + `to_seq`, and `distance`, one row per undirected neighbor pair. +} +\description{ +Extracts CDR3 amino acid sequences from a + SingleCellExperiment object and finds all pairs within a given edit + distance using the symmetric deletion lookup ("symdel") algorithm from + pyrepseq. This is the same algorithm described in the XT-neighbor + preprint, running on CPU rather than GPU. +} +\details{ +The search runs over unique sequences rather than over cells. + Repertoires carry heavy clonal redundancy, so this is both substantially + faster and what the symdel algorithm expects. Results are mapped back to + cells afterward. + + Degree is a neighborhood-density measure. High-degree sequences sit inside + dense clusters of similar receptors, which is the signal used to identify + convergent recombination and antigen-driven expansion. For explicit + cluster assignments, see \code{\link{runClustTCR}}. +} +\examples{ +data(immLynx_example) +\donttest{ + # Neighbors within one edit on the TRB chain + sce <- runSymdelNeighbors(immLynx_example, chains = "TRB") + + # Widen the search to two edits + sce <- runSymdelNeighbors(immLynx_example, + chains = "TRB", + max_edits = 2) + + # Get the neighbor edge list instead + edges <- runSymdelNeighbors(immLynx_example, + chains = "TRB", + return_object = FALSE) + head(edges) +} +} From 4419d906850f52a99c43ab7892aaec0aae881c56 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 07:46:50 -0500 Subject: [PATCH 06/11] Update DESCRIPTION --- DESCRIPTION | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 59a3252..9ec4481 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -10,8 +10,10 @@ Description: A comprehensive toolkit that bridges popular Python-based immune R environment. Provides unified interfaces for TCR distance calculations (tcrdist3), sequence generation probability (OLGA), selection inference (soNNia), clustering (clusTCR), protein embeddings (ESM-2), metaclone - discovery (metaclonotypist), and clonal expansion prediction from gene - expression (scXpand). Fully compatible with the scRepertoire and immApex ecosystem for single-cell + discovery (metaclonotypist), clonal expansion prediction from gene + expression (scXpand), near-neighbor search by symmetric deletion lookup + (pyrepseq), and deep generative featurization (DeepTCR). Fully compatible + with the scRepertoire and immApex ecosystem for single-cell immune repertoire analysis. License: MIT + file LICENSE Encoding: UTF-8 From 40204bbc42b40c11b7c703f835d7ae5b1cd01b2e Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 07:46:53 -0500 Subject: [PATCH 07/11] Create dot-deeptcrFeaturesToCells.Rd --- man/dot-deeptcrFeaturesToCells.Rd | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 man/dot-deeptcrFeaturesToCells.Rd diff --git a/man/dot-deeptcrFeaturesToCells.Rd b/man/dot-deeptcrFeaturesToCells.Rd new file mode 100644 index 0000000..95be1ab --- /dev/null +++ b/man/dot-deeptcrFeaturesToCells.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/runDeepTCR.R +\name{.deeptcrFeaturesToCells} +\alias{.deeptcrFeaturesToCells} +\title{Expand per-sequence features onto cells} +\usage{ +.deeptcrFeaturesToCells(features, unique_seqs, cell_seqs) +} +\arguments{ +\item{features}{Matrix with one row per unique sequence.} + +\item{unique_seqs}{Character vector matching the rows of `features`.} + +\item{cell_seqs}{Character vector of the sequence carried by each cell.} +} +\value{ +A matrix with one row per entry of `cell_seqs`. +} +\description{ +Expand per-sequence features onto cells +} +\keyword{internal} From 3c751161e721320db4970a80c72968b54eeb9dcd Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 07:47:05 -0500 Subject: [PATCH 08/11] update unit tests --- tests/testthat/helper-immLynx.R | 20 +++ tests/testthat/test-exportToScanpy.R | 11 +- tests/testthat/test-runDeepTCR.R | 115 +++++++++++++++ tests/testthat/test-runSymdelNeighbors.R | 169 +++++++++++++++++++++++ tests/testthat/test-zz-basilisk-specs.R | 126 +++++++++++++++++ 5 files changed, 438 insertions(+), 3 deletions(-) create mode 100644 tests/testthat/test-runDeepTCR.R create mode 100644 tests/testthat/test-runSymdelNeighbors.R create mode 100644 tests/testthat/test-zz-basilisk-specs.R diff --git a/tests/testthat/helper-immLynx.R b/tests/testthat/helper-immLynx.R index c54a981..bdb3719 100644 --- a/tests/testthat/helper-immLynx.R +++ b/tests/testthat/helper-immLynx.R @@ -256,3 +256,23 @@ mock_bcr_sce <- function() { ) ) } + +# Skip if the DeepTCR basilisk env is unavailable. DeepTCR pins TensorFlow +# 2.12, which publishes no linux-aarch64 wheel, so that platform is skipped +# outright rather than attempting a doomed multi-hundred-MB build. +skip_if_no_deeptcr <- function() { + skip_on_bioc_build() + + si <- Sys.info() + if (identical(unname(si[["sysname"]]), "Linux") && + grepl("aarch64|arm64", unname(si[["machine"]]))) { + testthat::skip("DeepTCR unavailable: no TensorFlow 2.12 linux-aarch64 wheel") + } + + ok <- tryCatch({ + proc <- basilisk::basiliskStart(immLynx:::deepTCREnv) + on.exit(basilisk::basiliskStop(proc)) + TRUE + }, error = function(e) FALSE) + if (!ok) testthat::skip("deepTCREnv not available") +} diff --git a/tests/testthat/test-exportToScanpy.R b/tests/testthat/test-exportToScanpy.R index 1f06713..bfe800c 100644 --- a/tests/testthat/test-exportToScanpy.R +++ b/tests/testthat/test-exportToScanpy.R @@ -437,9 +437,14 @@ test_that("exportToScanpy accepts a Seurat input via .coerceToSCE", { skip_if_not_installed("Matrix") skip_on_bioc_build() - m <- matrix(rpois(40 * 25, lambda = 2), 40, 25, - dimnames = list(paste0("g", seq_len(40)), - paste0("c", seq_len(25)))) + # Build the counts as a sparse matrix up front; CreateSeuratObject warns + # when it has to coerce a dense matrix itself. + m <- Matrix::Matrix( + matrix(rpois(40 * 25, lambda = 2), 40, 25, + dimnames = list(paste0("g", seq_len(40)), + paste0("c", seq_len(25)))), + sparse = TRUE + ) sobj <- Seurat::CreateSeuratObject(counts = m) sobj$sample <- rep(c("A", "B"), length.out = 25) diff --git a/tests/testthat/test-runDeepTCR.R b/tests/testthat/test-runDeepTCR.R new file mode 100644 index 0000000..654ae54 --- /dev/null +++ b/tests/testthat/test-runDeepTCR.R @@ -0,0 +1,115 @@ +# Tests for DeepTCR VAE featurization (issue #8). +# +# DeepTCR's Load_Data explicitly does not merge identical amino acid +# sequences, so the wrapper trains on unique sequences and expands the +# resulting feature rows back onto cells. That expansion is the pure-R logic +# worth testing without Python. + +# =========================================================================== +# .deeptcrFeaturesToCells +# =========================================================================== + +test_that(".deeptcrFeaturesToCells expands unique-sequence rows onto cells", { + # Two unique sequences, three cells, the first sequence used twice. + features <- matrix(c(1, 2, + 3, 4), nrow = 2, byrow = TRUE) + unique_seqs <- c("A", "B") + cell_seqs <- c("A", "B", "A") + + out <- immLynx:::.deeptcrFeaturesToCells(features, unique_seqs, cell_seqs) + + expect_identical(dim(out), c(3L, 2L)) + expect_identical(out[1, ], c(1, 2)) + expect_identical(out[2, ], c(3, 4)) + # Cells sharing a sequence must share a feature vector. + expect_identical(out[3, ], out[1, ]) +}) + +test_that(".deeptcrFeaturesToCells preserves cell order", { + features <- matrix(c(10, 20, 30), nrow = 3, ncol = 1) + unique_seqs <- c("X", "Y", "Z") + cell_seqs <- c("Z", "X", "Y") + + out <- immLynx:::.deeptcrFeaturesToCells(features, unique_seqs, cell_seqs) + + expect_identical(as.vector(out), c(30, 10, 20)) +}) + +test_that(".deeptcrFeaturesToCells errors when a cell sequence is unmatched", { + features <- matrix(1:2, nrow = 2, ncol = 1) + + expect_error( + immLynx:::.deeptcrFeaturesToCells(features, c("A", "B"), c("A", "C")), + "unique" + ) +}) + +test_that(".deeptcrFeaturesToCells errors on a length mismatch", { + features <- matrix(1:2, nrow = 2, ncol = 1) + + expect_error( + immLynx:::.deeptcrFeaturesToCells(features, c("A", "B", "C"), c("A")), + "rows" + ) +}) + +# =========================================================================== +# Argument validation, no Python required +# =========================================================================== + +test_that("runDeepTCR rejects non-SingleCellExperiment input", { + expect_error(runDeepTCR(data.frame(x = 1)), "SingleCellExperiment") +}) + +test_that("runDeepTCR rejects an invalid latent_dim", { + data("immLynx_example", package = "immLynx") + + expect_error(runDeepTCR(immLynx_example, latent_dim = 0), "latent_dim") + expect_error(runDeepTCR(immLynx_example, latent_dim = -5), "latent_dim") +}) + +# =========================================================================== +# End-to-end, requires the DeepTCR environment +# =========================================================================== + +test_that("runDeepTCR adds a reduction of the requested width", { + skip_if_no_deeptcr() + data("immLynx_example", package = "immLynx") + + sce <- runDeepTCR(immLynx_example, chains = "TRB", + latent_dim = 16, return_object = TRUE) + + expect_true("tcr_deeptcr" %in% SingleCellExperiment::reducedDimNames(sce)) + rd <- SingleCellExperiment::reducedDim(sce, "tcr_deeptcr") + expect_identical(nrow(rd), ncol(immLynx_example)) + expect_lte(ncol(rd), 16L) + # Cells carrying a TRB sequence must have finite features. + expect_true(any(stats::complete.cases(rd))) +}) + +test_that("runDeepTCR returns raw features when return_object is FALSE", { + skip_if_no_deeptcr() + data("immLynx_example", package = "immLynx") + + res <- runDeepTCR(immLynx_example, chains = "TRB", + latent_dim = 16, return_object = FALSE) + + expect_type(res, "list") + expect_true(all(c("features", "sequences", "explained_variance_ratio") %in% + names(res))) + expect_true(is.matrix(res$features)) + # One feature row per unique sequence trained on. + expect_identical(nrow(res$features), length(res$sequences)) +}) + +test_that("runDeepTCR leaves no model directory behind", { + skip_if_no_deeptcr() + data("immLynx_example", package = "immLynx") + + before <- list.files(getwd(), all.files = TRUE) + invisible(runDeepTCR(immLynx_example, chains = "TRB", + latent_dim = 16, return_object = FALSE)) + after <- list.files(getwd(), all.files = TRUE) + + expect_identical(setdiff(after, before), character(0)) +}) diff --git a/tests/testthat/test-runSymdelNeighbors.R b/tests/testthat/test-runSymdelNeighbors.R new file mode 100644 index 0000000..ad8b132 --- /dev/null +++ b/tests/testthat/test-runSymdelNeighbors.R @@ -0,0 +1,169 @@ +# Tests for symdel neighbor search (issue #6). +# +# The pure-R helpers are tested without Python. symdel returns a list of +# (i, j, distance) tuples using 0-based indices, and it emits both directions +# of every pair, so the conversion helpers carry the real logic. + +# =========================================================================== +# .symdelTripletsToEdges +# =========================================================================== + +test_that(".symdelTripletsToEdges converts 0-based indices to sequences", { + seqs <- c("CASSA", "CASSB", "CASSC") + # Both directions of a single pair (0, 1), as symdel emits them. + triplets <- list(c(0, 1, 1), c(1, 0, 1)) + + edges <- immLynx:::.symdelTripletsToEdges(triplets, seqs) + + expect_s3_class(edges, "data.frame") + expect_identical(names(edges), c("from_seq", "to_seq", "distance")) + expect_identical(nrow(edges), 1L) + expect_identical(edges$from_seq, "CASSA") + expect_identical(edges$to_seq, "CASSB") + expect_identical(edges$distance, 1L) +}) + +test_that(".symdelTripletsToEdges collapses both directions into one edge", { + seqs <- c("A", "B", "C") + triplets <- list(c(0, 1, 1), c(1, 0, 1), c(1, 2, 1), c(2, 1, 1)) + + edges <- immLynx:::.symdelTripletsToEdges(triplets, seqs) + + expect_identical(nrow(edges), 2L) + # Every edge is stored in a canonical orientation, so no pair repeats. + keys <- paste(edges$from_seq, edges$to_seq) + expect_identical(anyDuplicated(keys), 0L) +}) + +test_that(".symdelTripletsToEdges drops self-pairs", { + seqs <- c("A", "B") + triplets <- list(c(0, 0, 0), c(0, 1, 1), c(1, 0, 1)) + + edges <- immLynx:::.symdelTripletsToEdges(triplets, seqs) + + expect_identical(nrow(edges), 1L) + expect_false(any(edges$from_seq == edges$to_seq)) +}) + +test_that(".symdelTripletsToEdges returns an empty frame for no neighbors", { + edges <- immLynx:::.symdelTripletsToEdges(list(), c("A", "B")) + + expect_s3_class(edges, "data.frame") + expect_identical(nrow(edges), 0L) + expect_identical(names(edges), c("from_seq", "to_seq", "distance")) +}) + +# =========================================================================== +# .symdelDegree +# =========================================================================== + +test_that(".symdelDegree counts neighbors for each unique sequence", { + seqs <- c("A", "B", "C", "D") + edges <- data.frame( + from_seq = c("A", "B"), + to_seq = c("B", "C"), + distance = c(1L, 1L), + stringsAsFactors = FALSE + ) + + deg <- immLynx:::.symdelDegree(edges, seqs) + + # B touches both edges, A and C touch one each, D is isolated. + expect_identical(deg[["A"]], 1L) + expect_identical(deg[["B"]], 2L) + expect_identical(deg[["C"]], 1L) + expect_identical(deg[["D"]], 0L) +}) + +test_that(".symdelDegree returns all zeros when there are no edges", { + seqs <- c("A", "B") + edges <- data.frame(from_seq = character(0), to_seq = character(0), + distance = integer(0), stringsAsFactors = FALSE) + + deg <- immLynx:::.symdelDegree(edges, seqs) + + expect_identical(unname(deg), c(0L, 0L)) + expect_identical(names(deg), seqs) +}) + +test_that(".symdelDegree names every unique sequence in input order", { + seqs <- c("Z", "Y", "X") + edges <- data.frame(from_seq = "Z", to_seq = "X", distance = 1L, + stringsAsFactors = FALSE) + + deg <- immLynx:::.symdelDegree(edges, seqs) + + expect_identical(names(deg), seqs) +}) + +# =========================================================================== +# Argument validation +# =========================================================================== + +test_that("runSymdelNeighbors rejects non-SingleCellExperiment input", { + expect_error(runSymdelNeighbors(data.frame(x = 1)), + "SingleCellExperiment") +}) + +test_that("runSymdelNeighbors rejects an invalid max_edits", { + data("immLynx_example", package = "immLynx") + + expect_error(runSymdelNeighbors(immLynx_example, max_edits = 0), + "max_edits") + expect_error(runSymdelNeighbors(immLynx_example, max_edits = -1), + "max_edits") +}) + +# =========================================================================== +# End-to-end, requires the Python environment +# =========================================================================== + +test_that("runSymdelNeighbors returns an edge list", { + skip_if_no_python() + data("immLynx_example", package = "immLynx") + + edges <- runSymdelNeighbors(immLynx_example, chains = "TRB", + max_edits = 1, return_object = FALSE) + + expect_s3_class(edges, "data.frame") + expect_identical(names(edges), c("from_seq", "to_seq", "distance")) + expect_true(all(edges$distance <= 1)) + expect_false(any(edges$from_seq == edges$to_seq)) +}) + +test_that("runSymdelNeighbors adds a degree column to colData", { + skip_if_no_python() + data("immLynx_example", package = "immLynx") + + sce <- runSymdelNeighbors(immLynx_example, chains = "TRB", + max_edits = 1, return_object = TRUE) + + expect_true("symdel_degree" %in% colnames(SummarizedExperiment::colData(sce))) + deg <- SummarizedExperiment::colData(sce)$symdel_degree + expect_identical(length(deg), ncol(immLynx_example)) + # Cells without a TRB chain get NA, everything else is a non-negative count. + expect_true(all(deg[!is.na(deg)] >= 0)) +}) + +test_that("runSymdelNeighbors honors column_prefix", { + skip_if_no_python() + data("immLynx_example", package = "immLynx") + + sce <- runSymdelNeighbors(immLynx_example, chains = "TRB", + column_prefix = "nn", return_object = TRUE) + + expect_true("nn_degree" %in% colnames(SummarizedExperiment::colData(sce))) +}) + +test_that("runSymdelNeighbors finds more neighbors at a larger max_edits", { + skip_if_no_python() + data("immLynx_example", package = "immLynx") + + e1 <- runSymdelNeighbors(immLynx_example, chains = "TRB", + max_edits = 1, return_object = FALSE) + e2 <- runSymdelNeighbors(immLynx_example, chains = "TRB", + max_edits = 2, return_object = FALSE) + + expect_gte(nrow(e2), nrow(e1)) + expect_true(all(e2$distance <= 2)) +}) diff --git a/tests/testthat/test-zz-basilisk-specs.R b/tests/testthat/test-zz-basilisk-specs.R new file mode 100644 index 0000000..80c51e4 --- /dev/null +++ b/tests/testthat/test-zz-basilisk-specs.R @@ -0,0 +1,126 @@ +# Regression guards for basilisk environment specifications. +# +# Named "zz" so testthat runs these last, after any test that may have +# triggered a first-use environment build. +# +# Background: basilisk::setupBasiliskEnv() hands the `pip` vector straight to +# system2(), which pastes the arguments into a single string and runs it via +# /bin/sh without quoting. A spec such as "anndata>=0.8" is therefore split by +# the shell into the word "anndata", a ">" redirect, and the target "=0.8". +# Two things go wrong at once: the version floor is silently dropped, and pip's +# stdout is written to a file literally named "=0.8" in the working directory. +# +# The conda `packages` vector takes a different route: it reaches +# reticulate::conda_install(), which applies maybe_shQuote() to every argument +# before calling system2(). Comparison operators are therefore safe there, and +# a floor such as "pytorch-cpu>=2.5" is deliberate — pinning conda specs +# exactly makes the solve brittle across platforms for no safety gain. So the +# two vectors get two different rules, matching the two code paths. + +# pip: characters that survive an unquoted /bin/sh word split. Everything in +# the current specs is covered — version pins ("==", "="), dotted versions, +# dashes, underscores, local-version "+", and the git+https://...@tag URL used +# for clusTCR (":", "/", "@", "+"). +SHELL_SAFE <- "^[A-Za-z0-9._+/@=:-]+$" + +# conda: quoting makes redirects harmless, but a spec should still never carry +# whitespace, quotes, or command-substitution characters. Those signal a typo +# ("numpy >= 1.2") or an injection rather than a version bound. +# Note "<" and ">" are intentionally absent: they are valid, quoted, and +# meaningful in a conda version bound. +CONDA_JUNK <- "[[:space:]\"'`$;&|()]" + +# Collect every BasiliskEnvironment object defined in the package namespace, +# so environments added later are covered without editing this test. +basilisk_envs <- function() { + ns <- asNamespace("immLynx") + found <- list() + for (nm in ls(ns, all.names = TRUE)) { + obj <- tryCatch(get(nm, envir = ns), error = function(e) NULL) + if (methods::is(obj, "BasiliskEnvironment")) { + found[[nm]] <- obj + } + } + found +} + +test_that("every basilisk environment is discoverable for inspection", { + envs <- basilisk_envs() + + # Guard the guard: if the collector silently returns nothing, the spec + # checks below would pass vacuously. + expect_gt(length(envs), 0) + expect_true("immLynxEnv" %in% names(envs)) + expect_true("scanpyExportEnv" %in% names(envs)) +}) + +test_that("pip specs contain no shell metacharacters", { + envs <- basilisk_envs() + + for (nm in names(envs)) { + specs <- envs[[nm]]@pip + if (!length(specs)) next + + bad <- specs[!grepl(SHELL_SAFE, specs)] + expect_identical( + bad, character(0), + info = paste0( + "Unsafe pip spec(s) in ", nm, ": ", + paste(sprintf("'%s'", bad), collapse = ", "), + ". basilisk passes the pip vector through an unquoted shell, so '>' ", + "or '<' becomes a redirect: the version bound is dropped and a stray ", + "file is created. Use an exact '==' pin instead." + ) + ) + } +}) + +test_that("conda package specs are well formed", { + envs <- basilisk_envs() + + for (nm in names(envs)) { + specs <- envs[[nm]]@packages + if (!length(specs)) next + + bad <- specs[grepl(CONDA_JUNK, specs)] + expect_identical( + bad, character(0), + info = paste0( + "Malformed conda spec(s) in ", nm, ": ", + paste(sprintf("'%s'", bad), collapse = ", "), + ". Version bounds such as 'pytorch-cpu>=2.5' are fine here, but a ", + "spec must not contain whitespace, quotes, or shell command ", + "characters." + ) + ) + } +}) + +test_that("no shell-redirect artefact files were left in the package tree", { + # An unquoted ">" redirect writes to the working directory of the R process + # that built the environment. Under devtools::test() that is tests/testthat; + # under R CMD check it is the corresponding directory inside .Rcheck. Walking + # up two levels covers both, plus the package root for interactive use. + here <- normalizePath(getwd(), mustWork = FALSE) + roots <- unique(c( + here, + normalizePath(file.path(here, ".."), mustWork = FALSE), + normalizePath(file.path(here, "..", ".."), mustWork = FALSE) + )) + roots <- roots[dir.exists(roots)] + + stray <- unlist(lapply(roots, function(r) { + list.files(r, pattern = "^=", all.files = TRUE, full.names = TRUE) + }), use.names = FALSE) + stray <- if (is.null(stray)) character(0) else stray + + expect_identical( + stray, character(0), + info = paste0( + "Found file(s) whose name begins with '=': ", + paste(stray, collapse = ", "), + ". These are shell redirect artefacts from an unquoted version spec. ", + "Delete them and fix the offending spec in R/basilisk.R." + ) + ) +}) From 22bbf36d6b7110303f82d320ae00e10a6a3b42fa Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Sun, 9 Aug 2026 07:47:15 -0500 Subject: [PATCH 09/11] update vignettes --- vignettes/advanced_analysis.Rmd | 53 +++++++++++++++++ vignettes/immLynx_vignette.Rmd | 101 +++++++++++++++++++++++++++++++- 2 files changed, 151 insertions(+), 3 deletions(-) diff --git a/vignettes/advanced_analysis.Rmd b/vignettes/advanced_analysis.Rmd index c61ce9b..4f93af5 100644 --- a/vignettes/advanced_analysis.Rmd +++ b/vignettes/advanced_analysis.Rmd @@ -176,6 +176,43 @@ sce_paired <- scater::runUMAP(sce_paired, scater::plotReducedDim(sce_paired, dimred = "UMAP") ``` +## Dataset-Specific Features with DeepTCR + +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. + +```{r deeptcr-compare, eval=FALSE} +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. + # Integration with scRepertoire Clonotypes immLynx is designed to complement `r Biocpkg("scRepertoire")`, which @@ -342,6 +379,13 @@ help manage memory and runtime: distances, which scales quadratically. For very large datasets, computing distances on a representative subset is a practical alternative. +- **Near-neighbor search instead of a distance matrix**: when the + question is "which sequences are similar" rather than "how far apart + is every pair", `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. ```{r large-data, eval=FALSE} sce_large <- runEmbeddings( @@ -358,6 +402,15 @@ sce_large <- runClustTCR( 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 +) ``` # HLA Association Analysis diff --git a/vignettes/immLynx_vignette.Rmd b/vignettes/immLynx_vignette.Rmd index cf4bf8a..348dca3 100644 --- a/vignettes/immLynx_vignette.Rmd +++ b/vignettes/immLynx_vignette.Rmd @@ -75,6 +75,11 @@ following analytical capabilities: structural properties - **Metaclone discovery** via metaclonotypist, identifying groups of related TCR clonotypes that may share antigen specificity +- **Near-neighbor search** via symmetric deletion lookup, finding all + CDR3 pairs within a given edit distance without computing a full + pairwise distance matrix +- **Deep generative featurization** via DeepTCR, learning a compact + latent representation of CDR3 sequences with a variational autoencoder This vignette demonstrates the core functionality of immLynx, walking through data preparation, individual analysis functions, and a combined @@ -96,9 +101,15 @@ immLynx uses `r Biocpkg("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. +pyrepseq, 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. + +A few tools pin dependency versions that cannot coexist with that main +environment, so they get their own. `exportToScanpy()` uses a +scanpy/scirpy environment, and `runDeepTCR()` uses a TensorFlow 2.12 +environment. Each is built on first use of the function that needs it, +so you only pay for the ones you actually call. Two functions use their own environments rather than that one. `exportToScanpy()` needs the scanpy/scirpy stack, and `runScXpand()` needs @@ -186,6 +197,50 @@ sce <- runClustTCR( table(sce$clustcr_TRB) ``` +## Near-Neighbor Search with symdel + +Symmetric deletion lookup finds every pair of CDR3 sequences within a +given edit distance without ever building a full pairwise distance +matrix. It indexes each sequence by its set of single-character +deletions, so two sequences within one edit share at least one deletion +variant and collide in the index. Cost scales with the number of +sequences rather than with their square, which makes it practical on +repertoires where a full distance matrix would not fit in memory. + +`runSymdelNeighbors()` searches over unique sequences and maps results +back onto cells. By default it adds a `symdel_degree` column to +`colData` giving the number of neighbors each cell's CDR3 has. High +degree marks sequences sitting in dense similarity neighborhoods, the +signature of convergent recombination and antigen-driven expansion. + +```{r symdel, eval=python_available} +sce <- runSymdelNeighbors( + immLynx_example, + chains = "TRB", + max_edits = 1 +) + +table(sce$symdel_degree, useNA = "ifany") +``` + +Setting `return_object = FALSE` returns the neighbor pairs themselves, +which is the more useful form when building a similarity graph. + +```{r symdel-edges, eval=python_available} +edges <- runSymdelNeighbors( + immLynx_example, + chains = "TRB", + max_edits = 2, + return_object = FALSE +) + +head(edges) +``` + +Note that this is the same algorithm benchmarked in the XT-neighbor +work. immLynx uses the CPU implementation from `pyrepseq`, so no GPU is +required. + ## TCR Distance Calculations with tcrdist3 tcrdist3 computes pairwise distances between TCR sequences using a @@ -262,6 +317,39 @@ sce <- scater::runUMAP(sce, scater::plotReducedDim(sce, dimred = "UMAP") ``` +## Deep Generative Featurization with DeepTCR + +DeepTCR learns a latent representation of CDR3 sequences with a +variational autoencoder. Where ESM-2 brings general protein knowledge +learned from UniProt, the DeepTCR VAE is trained on the repertoire in +front of it, so its features describe the structure of that specific +dataset. The two are complementary rather than interchangeable. + +`runDeepTCR()` trains on unique sequences and writes the learned +features to a dimensional reduction. DeepTCR prunes uninformative latent +dimensions during training, so the result may be narrower than +`latent_dim` requests. + +```{r deeptcr, eval=FALSE} +sce <- runDeepTCR( + immLynx_example, + chains = "TRB", + latent_dim = 64, + verbose = FALSE +) + +dim(SingleCellExperiment::reducedDim(sce, "tcr_deeptcr")) +``` + +This chunk is not evaluated when the vignette is built. DeepTCR pins +TensorFlow 2.12 and its own numpy, pandas, and scipy versions, so it +runs in a separate `deepTCREnv` environment that is a large first-use +download. Set `eval = TRUE` to run it locally. + +Passing `use_genes = TRUE` includes V and J gene usage alongside the +CDR3 sequence, which helps when gene usage carries specificity signal +that the CDR3 alone does not. + ## Metaclone Discovery with Metaclonotypist Metaclonotypist identifies metaclones---groups of TCR clonotypes that @@ -331,6 +419,13 @@ immLynx_example <- runClustTCR( method = "mcl" ) +# Count near neighbors within one edit +immLynx_example <- runSymdelNeighbors( + immLynx_example, + chains = "TRB", + max_edits = 1 +) + # Calculate generation probability immLynx_example <- runOLGA( immLynx_example, From 483de9ae0b16c256b95b4b9a4669bb7ceba3a6b9 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Mon, 10 Aug 2026 07:05:53 -0500 Subject: [PATCH 10/11] Add full Seurat Support --- .Rbuildignore | 1 + .gitignore | 1 + R/runClustTCR.R | 35 +----- R/runDeepTCR.R | 8 +- R/runEmbeddings.R | 20 ++- R/runMetaclonotypist.R | 5 +- R/runOLGA.R | 26 +--- R/runScXpand.R | 33 +---- R/runSoNNia.R | 8 +- R/runSymdelNeighbors.R | 14 +-- R/runTCRdist.R | 17 +-- R/utils.R | 115 +++++++++++++++++ man/dot-assertSCObject.Rd | 20 +++ man/dot-isSCObject.Rd | 18 +++ man/dot-writeCellColumn.Rd | 25 ++++ man/dot-writeObjMetadata.Rd | 23 ++++ man/dot-writeReduction.Rd | 24 ++++ man/runClustTCR.Rd | 2 +- man/runDeepTCR.Rd | 2 +- man/runEmbeddings.Rd | 2 +- man/runOLGA.Rd | 2 +- man/runSoNNia.Rd | 2 +- man/runSymdelNeighbors.Rd | 2 +- man/runTCRdist.Rd | 2 +- tests/testthat/test-object-accessors.R | 168 +++++++++++++++++++++++++ tests/testthat/test-runClustTCR.R | 6 +- tests/testthat/test-runEmbeddings.R | 8 +- tests/testthat/test-runOLGA.R | 2 +- tests/testthat/test-runSoNNia.R | 2 +- tests/testthat/test-runTCRdist.R | 4 +- tests/testthat/test-seurat-parity.R | 144 +++++++++++++++++++++ 31 files changed, 599 insertions(+), 142 deletions(-) create mode 100644 man/dot-assertSCObject.Rd create mode 100644 man/dot-isSCObject.Rd create mode 100644 man/dot-writeCellColumn.Rd create mode 100644 man/dot-writeObjMetadata.Rd create mode 100644 man/dot-writeReduction.Rd create mode 100644 tests/testthat/test-object-accessors.R create mode 100644 tests/testthat/test-seurat-parity.R diff --git a/.Rbuildignore b/.Rbuildignore index cc0c9b9..23c8784 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -11,3 +11,4 @@ ^cran-comments\.md$ ^doc$ ^Meta$ +^\.claude$ diff --git a/.gitignore b/.gitignore index cb7105d..25ef0f6 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,4 @@ .Rhistory /doc/ /Meta/ +.claude/ diff --git a/R/runClustTCR.R b/R/runClustTCR.R index dc706c9..4d19d43 100644 --- a/R/runClustTCR.R +++ b/R/runClustTCR.R @@ -4,7 +4,7 @@ #' SingleCellExperiment object with scRepertoire data and performs clustering #' using the clusTCR algorithm. #' -#' @param input A SingleCellExperiment object containing scRepertoire +#' @param input A SingleCellExperiment or Seurat object containing scRepertoire #' TCR data in the metadata. #' @param chains Character string specifying which chains to use: "TRA", "TRB", or "both". #' Default is "TRB". @@ -63,36 +63,13 @@ runClustTCR <- function(input, chains <- match.arg(chains) # Determine input type - .is_sce <- methods::is(input, "SingleCellExperiment") + .assertSCObject(input) - if (!.is_sce) { - stop("Input must be a SingleCellExperiment object") - } - - # Get cell names based on object type - .get_cells <- function(obj) { - if (methods::is(obj, "SingleCellExperiment")) { - colnames(obj) - } else { - colnames(obj) - } - } + .get_cells <- function(obj) colnames(obj) - # Add metadata based on object type - .add_metadata <- function(obj, col_name, values, cell_names) { - if (methods::is(obj, "SingleCellExperiment")) { - col_vec <- rep(NA, ncol(obj)) - names(col_vec) <- colnames(obj) - col_vec[cell_names] <- values - colData(obj)[[col_name]] <- col_vec - } else { - col_vec <- rep(NA, ncol(obj)) - names(col_vec) <- colnames(obj) - col_vec[cell_names] <- values - obj[[col_name]] <- col_vec - } - obj - } + # Delegates to the shared writer so SingleCellExperiment and Seurat produce + # identical columns, including dropping names from the stored vector. + .add_metadata <- .writeCellColumn # Extract TCR data using immApex message("Extracting TCR sequences from object...") diff --git a/R/runDeepTCR.R b/R/runDeepTCR.R index fae9964..0bf3983 100644 --- a/R/runDeepTCR.R +++ b/R/runDeepTCR.R @@ -4,7 +4,7 @@ #' SingleCellExperiment object and learns a low-dimensional representation #' using DeepTCR's unsupervised variational autoencoder (DeepTCR_U). #' -#' @param input A SingleCellExperiment object containing scRepertoire TCR data. +#' @param input A SingleCellExperiment or Seurat object containing scRepertoire TCR data. #' @param chains Which chain(s) to featurize: "TRB", "TRA", or "both". #' Default is "TRB". #' @param latent_dim Width of the VAE latent space. Default is 256. @@ -84,9 +84,7 @@ runDeepTCR <- function(input, chains <- match.arg(chains) - if (!methods::is(input, "SingleCellExperiment")) { - stop("Input must be a SingleCellExperiment object") - } + .assertSCObject(input) if (!is.numeric(latent_dim) || length(latent_dim) != 1L || is.na(latent_dim) || latent_dim < 1) { @@ -147,7 +145,7 @@ runDeepTCR <- function(input, seq_len(ncol(cell_features))))) full[seq_map$barcodes, ] <- cell_features - SingleCellExperiment::reducedDim(input, reduction_name) <- full + input <- .writeReduction(input, reduction_name, full, reduction_key) message("Features added as '", reduction_name, "' reduction") input diff --git a/R/runEmbeddings.R b/R/runEmbeddings.R index a9e7806..a877899 100644 --- a/R/runEmbeddings.R +++ b/R/runEmbeddings.R @@ -3,7 +3,7 @@ #' @description Extracts TCR CDR3 sequences from a SingleCellExperiment object #' and generates embeddings using a protein language model (e.g., ESM-2). #' -#' @param input A SingleCellExperiment object containing scRepertoire TCR data. +#' @param input A SingleCellExperiment or Seurat object containing scRepertoire TCR data. #' @param chains Which chain(s) to embed: "TRB", "TRA", or "both". Default is "TRB". #' @param model_name Hugging Face model name. Default is "facebook/esm2_t12_35M_UR50D". #' Other options: "facebook/esm2_t33_650M_UR50D", "facebook/esm2_t36_3B_UR50D" @@ -69,11 +69,7 @@ runEmbeddings <- function(input, chains <- match.arg(chains) # Determine input type - .is_sce <- methods::is(input, "SingleCellExperiment") - - if (!.is_sce) { - stop("Input must be a SingleCellExperiment object") - } + .assertSCObject(input) message("Loading Hugging Face model: ", model_name) hf_components <- huggingModel(model_name = model_name) @@ -165,13 +161,13 @@ runEmbeddings <- function(input, cell_embeddings[barcodes, ] <- embeddings # Add as reduced dimension - SingleCellExperiment::reducedDim(input, reduction_name) <- cell_embeddings + input <- .writeReduction(input, reduction_name, cell_embeddings, + reduction_key) - # Add chain info to colData - chain_meta <- rep(NA_character_, ncol(input)) - names(chain_meta) <- colnames(input) - chain_meta[barcodes] <- chain_info - colData(input)[[paste0(reduction_name, "_chain")]] <- chain_meta + # Add chain info to cell metadata + input <- .writeCellColumn(input, + paste0(reduction_name, "_chain"), + chain_info, barcodes) message("Embeddings added as '", reduction_name, "' reduction") message("Use RunUMAP(obj, reduction='", reduction_name, "') to visualize") diff --git a/R/runMetaclonotypist.R b/R/runMetaclonotypist.R index 1d65951..9d2a962 100644 --- a/R/runMetaclonotypist.R +++ b/R/runMetaclonotypist.R @@ -76,8 +76,9 @@ if (is.null(max_dist)) { max_dist <- if (method == "tcrdist") 20 else 1.5 } - # Extract TCR data - is_sc_object <- methods::is(input, "SingleCellExperiment") + # Extract TCR data. Anything that is not a single-cell container is + # treated as a pre-extracted data.frame. + is_sc_object <- .isSCObject(input) if (is_sc_object) { chain_code <- if (chains == "alpha") "TRA" else "TRB" diff --git a/R/runOLGA.R b/R/runOLGA.R index b2b3002..801a40a 100644 --- a/R/runOLGA.R +++ b/R/runOLGA.R @@ -3,7 +3,7 @@ #' @description Extracts TCR sequences from a SingleCellExperiment object #' and calculates their generation probability using OLGA. #' -#' @param input A SingleCellExperiment object containing scRepertoire TCR data. +#' @param input A SingleCellExperiment or Seurat object containing scRepertoire TCR data. #' @param chains Which chain to analyze: "TRA" or "TRB". Default is "TRB". #' @param model OLGA model to use. Options: "humanTRB", "humanTRA", "humanIGH", "mouseTRB". #' If NULL, will be inferred from organism and chains parameters. @@ -55,27 +55,11 @@ runOLGA <- function(input, chains <- match.arg(chains) # Determine input type - .is_sce <- methods::is(input, "SingleCellExperiment") + .assertSCObject(input) - if (!.is_sce) { - stop("Input must be a SingleCellExperiment object") - } - - # Helper to add metadata - .add_metadata <- function(obj, col_name, values, cell_names) { - if (methods::is(obj, "SingleCellExperiment")) { - col_vec <- rep(NA_real_, ncol(obj)) - names(col_vec) <- colnames(obj) - col_vec[cell_names] <- values - colData(obj)[[col_name]] <- col_vec - } else { - col_vec <- rep(NA_real_, ncol(obj)) - names(col_vec) <- colnames(obj) - col_vec[cell_names] <- values - obj[[col_name]] <- col_vec - } - obj - } + # Delegates to the shared writer so SingleCellExperiment and Seurat produce + # identical columns, including dropping names from the stored vector. + .add_metadata <- .writeCellColumn # Infer model if not specified if (is.null(model)) { diff --git a/R/runScXpand.R b/R/runScXpand.R index 82af9fc..ffc3b13 100644 --- a/R/runScXpand.R +++ b/R/runScXpand.R @@ -589,37 +589,8 @@ out } -#' @keywords internal -# Fill a per-cell column, leaving NA for cells that were not scored. -.writeCellColumn <- function(obj, col_name, values, cell_names) { - col_vec <- rep(NA, ncol(obj)) - names(col_vec) <- colnames(obj) - col_vec[cell_names] <- values - names(col_vec) <- NULL - if (methods::is(obj, "SingleCellExperiment")) { - SummarizedExperiment::colData(obj)[[col_name]] <- col_vec - } else { - obj[[col_name]] <- col_vec - } - obj -} - -#' @keywords internal -# Stash the run summary. SingleCellExperiment has metadata(); Seurat keeps -# the equivalent in the misc slot. -.writeObjMetadata <- function(obj, key, value) { - if (methods::is(obj, "SingleCellExperiment")) { - md <- S4Vectors::metadata(obj) - md[[key]] <- value - S4Vectors::metadata(obj) <- md - } else { - misc <- methods::slot(obj, "misc") - if (!is.list(misc)) misc <- list() - misc[[key]] <- value - methods::slot(obj, "misc") <- misc - } - obj -} +# .writeCellColumn() and .writeObjMetadata() now live in R/utils.R, shared +# with every other wrapper that accepts Seurat or SingleCellExperiment. # =========================================================================== # Exported functions diff --git a/R/runSoNNia.R b/R/runSoNNia.R index b595ffe..0e84d46 100644 --- a/R/runSoNNia.R +++ b/R/runSoNNia.R @@ -3,7 +3,7 @@ #' @description Infer selection pressures on TCRs using soNNia. Requires a #' background dataset of unselected sequences (generated by OLGA). #' -#' @param input A SingleCellExperiment object containing scRepertoire TCR data. +#' @param input A SingleCellExperiment or Seurat object containing scRepertoire TCR data. #' @param chains Which chain to analyze: "TRB" or "TRA". Default is "TRB". #' @param background_file Path to CSV file with background sequences (from generateOLGA). #' @param organism Organism: "human" or "mouse". Default is "human". @@ -54,11 +54,7 @@ runSoNNia <- function(input, chains <- match.arg(chains) # Determine input type - .is_sce <- methods::is(input, "SingleCellExperiment") - - if (!.is_sce) { - stop("Input must be a SingleCellExperiment object") - } + .assertSCObject(input) if (!file.exists(background_file)) { stop("Background file not found: ", background_file) diff --git a/R/runSymdelNeighbors.R b/R/runSymdelNeighbors.R index 707089f..fe81fe5 100644 --- a/R/runSymdelNeighbors.R +++ b/R/runSymdelNeighbors.R @@ -6,7 +6,7 @@ #' pyrepseq. This is the same algorithm described in the XT-neighbor #' preprint, running on CPU rather than GPU. #' -#' @param input A SingleCellExperiment object containing scRepertoire TCR data. +#' @param input A SingleCellExperiment or Seurat object containing scRepertoire TCR data. #' @param chains Which chain(s) to search: "TRB", "TRA", or "both". #' Default is "TRB". #' @param max_edits Maximum edit distance defining a neighbor. Default is 1. @@ -71,9 +71,7 @@ runSymdelNeighbors <- function(input, chains <- match.arg(chains) - if (!methods::is(input, "SingleCellExperiment")) { - stop("Input must be a SingleCellExperiment object") - } + .assertSCObject(input) if (!is.numeric(max_edits) || length(max_edits) != 1L || is.na(max_edits) || max_edits < 1) { @@ -112,12 +110,8 @@ runSymdelNeighbors <- function(input, # for the requested chain stay NA. cell_deg <- unname(deg[match(seq_map$sequences, unique_seqs)]) - col_vec <- rep(NA_integer_, ncol(input)) - names(col_vec) <- colnames(input) - col_vec[seq_map$barcodes] <- cell_deg - - colData(input)[[paste0(column_prefix, "_degree")]] <- unname(col_vec) - input + .writeCellColumn(input, paste0(column_prefix, "_degree"), + cell_deg, seq_map$barcodes) } diff --git a/R/runTCRdist.R b/R/runTCRdist.R index 0f93013..ca820fc 100644 --- a/R/runTCRdist.R +++ b/R/runTCRdist.R @@ -4,7 +4,7 @@ #' SingleCellExperiment object with scRepertoire data and calculates pairwise #' TCR distances using tcrdist3. #' -#' @param input A SingleCellExperiment object containing scRepertoire TCR data. +#' @param input A SingleCellExperiment or Seurat object containing scRepertoire TCR data. #' @param chains Character vector specifying chains: "alpha", "beta", or c("alpha", "beta"). #' Default is "beta". #' @param organism Organism: "human" or "mouse". Default is "human". @@ -50,11 +50,7 @@ runTCRdist <- function(input, add_to_object = FALSE) { # Determine input type - .is_sce <- methods::is(input, "SingleCellExperiment") - - if (!.is_sce) { - stop("Input must be a SingleCellExperiment object") - } + .assertSCObject(input) message("Extracting TCR sequences from object...") @@ -167,8 +163,13 @@ runTCRdist <- function(input, ) if (add_to_object) { - S4Vectors::metadata(input)$tcrdist <- output - message("TCR distances added to object at metadata(obj)$tcrdist") + input <- .writeObjMetadata(input, "tcrdist", output) + message("TCR distances added to object at ", + if (methods::is(input, "SingleCellExperiment")) { + "metadata(obj)$tcrdist" + } else { + "obj@misc$tcrdist" + }) return(input) } else { return(output) diff --git a/R/utils.R b/R/utils.R index a5c9bc1..766aeb8 100644 --- a/R/utils.R +++ b/R/utils.R @@ -626,3 +626,118 @@ summarizeTCRrepertoire <- function(input, } out } + +# =========================================================================== +# Single-cell object accessors +# +# immLynx accepts either a SingleCellExperiment or a Seurat object and returns +# whichever class it was given. These helpers are the only place that knows +# the difference. Extraction is already class-agnostic because +# immApex::getIR() handles both, so only the guards and the write-back paths +# need to branch. +# =========================================================================== + +#' Is this a supported single-cell container? +#' +#' @param x Any object. +#' @return TRUE for SingleCellExperiment or Seurat, FALSE otherwise. +#' @keywords internal +.isSCObject <- function(x) { + methods::is(x, "SingleCellExperiment") || methods::is(x, "Seurat") +} + +#' Stop unless the input is a supported single-cell container +#' +#' @param x Object to check. +#' @param arg Name to use in the error message. +#' @return Invisibly TRUE; called for the error. +#' @keywords internal +.assertSCObject <- function(x, arg = "input") { + if (!.isSCObject(x)) { + stop(arg, " must be a SingleCellExperiment or Seurat object.", + call. = FALSE) + } + invisible(TRUE) +} + +#' Write a per-cell column back to the object +#' +#' @description Values are placed by barcode, so `cell_names` may cover only a +#' subset of the object. Cells not named are left as NA. +#' +#' @param obj A SingleCellExperiment or Seurat object. +#' @param col_name Column name to write. +#' @param values Vector of values, aligned to `cell_names`. +#' @param cell_names Barcodes the values belong to. +#' @return `obj` with the column added, same class as supplied. +#' @keywords internal +.writeCellColumn <- function(obj, col_name, values, cell_names) { + col_vec <- rep(NA, ncol(obj)) + names(col_vec) <- colnames(obj) + col_vec[cell_names] <- values + names(col_vec) <- NULL + if (methods::is(obj, "SingleCellExperiment")) { + SummarizedExperiment::colData(obj)[[col_name]] <- col_vec + } else { + obj[[col_name]] <- col_vec + } + obj +} + +#' Write a dimensional reduction back to the object +#' +#' @param obj A SingleCellExperiment or Seurat object. +#' @param name Name of the reduction. +#' @param mat Matrix with one row per cell, rownames matching `colnames(obj)`. +#' @param key Column key prefix, used by Seurat for axis labels. +#' @return `obj` with the reduction added, same class as supplied. +#' @keywords internal +.writeReduction <- function(obj, name, mat, key) { + if (is.null(rownames(mat))) { + rownames(mat) <- colnames(obj) + } + + if (methods::is(obj, "SingleCellExperiment")) { + SingleCellExperiment::reducedDim(obj, name) <- mat + return(obj) + } + + # Seurat validates the key: it must start with a letter, contain only + # alphanumerics, and end in an underscore. + key <- gsub("[^A-Za-z0-9]", "", key) + if (!nzchar(key)) key <- "DIM" + if (!grepl("^[A-Za-z]", key)) key <- paste0("X", key) + key <- paste0(key, "_") + + colnames(mat) <- paste0(key, seq_len(ncol(mat))) + obj[[name]] <- Seurat::CreateDimReducObject( + embeddings = mat, + key = key, + assay = Seurat::DefaultAssay(obj) + ) + obj +} + +#' Stash a run summary on the object +#' +#' @description SingleCellExperiment has metadata(); Seurat keeps the +#' equivalent in the misc slot. +#' +#' @param obj A SingleCellExperiment or Seurat object. +#' @param key Name to store under. +#' @param value Value to store. +#' @return `obj` with the entry added, same class as supplied. +#' @keywords internal +.writeObjMetadata <- function(obj, key, value) { + if (methods::is(obj, "SingleCellExperiment")) { + md <- S4Vectors::metadata(obj) + md[[key]] <- value + S4Vectors::metadata(obj) <- md + } else { + misc <- methods::slot(obj, "misc") + if (!is.list(misc)) misc <- list() + misc[[key]] <- value + methods::slot(obj, "misc") <- misc + } + obj +} diff --git a/man/dot-assertSCObject.Rd b/man/dot-assertSCObject.Rd new file mode 100644 index 0000000..aef1ea5 --- /dev/null +++ b/man/dot-assertSCObject.Rd @@ -0,0 +1,20 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{.assertSCObject} +\alias{.assertSCObject} +\title{Stop unless the input is a supported single-cell container} +\usage{ +.assertSCObject(x, arg = "input") +} +\arguments{ +\item{x}{Object to check.} + +\item{arg}{Name to use in the error message.} +} +\value{ +Invisibly TRUE; called for the error. +} +\description{ +Stop unless the input is a supported single-cell container +} +\keyword{internal} diff --git a/man/dot-isSCObject.Rd b/man/dot-isSCObject.Rd new file mode 100644 index 0000000..40cd944 --- /dev/null +++ b/man/dot-isSCObject.Rd @@ -0,0 +1,18 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{.isSCObject} +\alias{.isSCObject} +\title{Is this a supported single-cell container?} +\usage{ +.isSCObject(x) +} +\arguments{ +\item{x}{Any object.} +} +\value{ +TRUE for SingleCellExperiment or Seurat, FALSE otherwise. +} +\description{ +Is this a supported single-cell container? +} +\keyword{internal} diff --git a/man/dot-writeCellColumn.Rd b/man/dot-writeCellColumn.Rd new file mode 100644 index 0000000..aada449 --- /dev/null +++ b/man/dot-writeCellColumn.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{.writeCellColumn} +\alias{.writeCellColumn} +\title{Write a per-cell column back to the object} +\usage{ +.writeCellColumn(obj, col_name, values, cell_names) +} +\arguments{ +\item{obj}{A SingleCellExperiment or Seurat object.} + +\item{col_name}{Column name to write.} + +\item{values}{Vector of values, aligned to `cell_names`.} + +\item{cell_names}{Barcodes the values belong to.} +} +\value{ +`obj` with the column added, same class as supplied. +} +\description{ +Values are placed by barcode, so `cell_names` may cover only a + subset of the object. Cells not named are left as NA. +} +\keyword{internal} diff --git a/man/dot-writeObjMetadata.Rd b/man/dot-writeObjMetadata.Rd new file mode 100644 index 0000000..ae8951f --- /dev/null +++ b/man/dot-writeObjMetadata.Rd @@ -0,0 +1,23 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{.writeObjMetadata} +\alias{.writeObjMetadata} +\title{Stash a run summary on the object} +\usage{ +.writeObjMetadata(obj, key, value) +} +\arguments{ +\item{obj}{A SingleCellExperiment or Seurat object.} + +\item{key}{Name to store under.} + +\item{value}{Value to store.} +} +\value{ +`obj` with the entry added, same class as supplied. +} +\description{ +SingleCellExperiment has metadata(); Seurat keeps the + equivalent in the misc slot. +} +\keyword{internal} diff --git a/man/dot-writeReduction.Rd b/man/dot-writeReduction.Rd new file mode 100644 index 0000000..e717b9d --- /dev/null +++ b/man/dot-writeReduction.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils.R +\name{.writeReduction} +\alias{.writeReduction} +\title{Write a dimensional reduction back to the object} +\usage{ +.writeReduction(obj, name, mat, key) +} +\arguments{ +\item{obj}{A SingleCellExperiment or Seurat object.} + +\item{name}{Name of the reduction.} + +\item{mat}{Matrix with one row per cell, rownames matching `colnames(obj)`.} + +\item{key}{Column key prefix, used by Seurat for axis labels.} +} +\value{ +`obj` with the reduction added, same class as supplied. +} +\description{ +Write a dimensional reduction back to the object +} +\keyword{internal} diff --git a/man/runClustTCR.Rd b/man/runClustTCR.Rd index 74acfc5..83509b0 100644 --- a/man/runClustTCR.Rd +++ b/man/runClustTCR.Rd @@ -15,7 +15,7 @@ runClustTCR( ) } \arguments{ -\item{input}{A SingleCellExperiment object containing scRepertoire +\item{input}{A SingleCellExperiment or Seurat object containing scRepertoire TCR data in the metadata.} \item{chains}{Character string specifying which chains to use: "TRA", "TRB", or "both". diff --git a/man/runDeepTCR.Rd b/man/runDeepTCR.Rd index 614f2bc..21b5fad 100644 --- a/man/runDeepTCR.Rd +++ b/man/runDeepTCR.Rd @@ -20,7 +20,7 @@ runDeepTCR( ) } \arguments{ -\item{input}{A SingleCellExperiment object containing scRepertoire TCR data.} +\item{input}{A SingleCellExperiment or Seurat object containing scRepertoire TCR data.} \item{chains}{Which chain(s) to featurize: "TRB", "TRA", or "both". Default is "TRB".} diff --git a/man/runEmbeddings.Rd b/man/runEmbeddings.Rd index 61344c0..481618e 100644 --- a/man/runEmbeddings.Rd +++ b/man/runEmbeddings.Rd @@ -17,7 +17,7 @@ runEmbeddings( ) } \arguments{ -\item{input}{A SingleCellExperiment object containing scRepertoire TCR data.} +\item{input}{A SingleCellExperiment or Seurat object containing scRepertoire TCR data.} \item{chains}{Which chain(s) to embed: "TRB", "TRA", or "both". Default is "TRB".} diff --git a/man/runOLGA.Rd b/man/runOLGA.Rd index a0e0004..5d720aa 100644 --- a/man/runOLGA.Rd +++ b/man/runOLGA.Rd @@ -15,7 +15,7 @@ runOLGA( ) } \arguments{ -\item{input}{A SingleCellExperiment object containing scRepertoire TCR data.} +\item{input}{A SingleCellExperiment or Seurat object containing scRepertoire TCR data.} \item{chains}{Which chain to analyze: "TRA" or "TRB". Default is "TRB".} diff --git a/man/runSoNNia.Rd b/man/runSoNNia.Rd index 439d810..5033058 100644 --- a/man/runSoNNia.Rd +++ b/man/runSoNNia.Rd @@ -15,7 +15,7 @@ runSoNNia( ) } \arguments{ -\item{input}{A SingleCellExperiment object containing scRepertoire TCR data.} +\item{input}{A SingleCellExperiment or Seurat object containing scRepertoire TCR data.} \item{chains}{Which chain to analyze: "TRB" or "TRA". Default is "TRB".} diff --git a/man/runSymdelNeighbors.Rd b/man/runSymdelNeighbors.Rd index f0b034d..9a31b18 100644 --- a/man/runSymdelNeighbors.Rd +++ b/man/runSymdelNeighbors.Rd @@ -17,7 +17,7 @@ runSymdelNeighbors( ) } \arguments{ -\item{input}{A SingleCellExperiment object containing scRepertoire TCR data.} +\item{input}{A SingleCellExperiment or Seurat object containing scRepertoire TCR data.} \item{chains}{Which chain(s) to search: "TRB", "TRA", or "both". Default is "TRB".} diff --git a/man/runTCRdist.Rd b/man/runTCRdist.Rd index 5c9c1e9..4a2954f 100644 --- a/man/runTCRdist.Rd +++ b/man/runTCRdist.Rd @@ -13,7 +13,7 @@ runTCRdist( ) } \arguments{ -\item{input}{A SingleCellExperiment object containing scRepertoire TCR data.} +\item{input}{A SingleCellExperiment or Seurat object containing scRepertoire TCR data.} \item{chains}{Character vector specifying chains: "alpha", "beta", or c("alpha", "beta"). Default is "beta".} diff --git a/tests/testthat/test-object-accessors.R b/tests/testthat/test-object-accessors.R new file mode 100644 index 0000000..d38a58e --- /dev/null +++ b/tests/testthat/test-object-accessors.R @@ -0,0 +1,168 @@ +# Tests for the shared single-cell object accessors. +# +# immLynx functions accept either a SingleCellExperiment or a Seurat object +# and return the same class they were given. These helpers are the only place +# that knows the difference, so they carry the burden of that contract. + +make_sce <- function() { + data("immLynx_example", package = "immLynx") + immLynx_example +} + +make_seurat <- function() { + skip_if_not_installed("Seurat") + skip_if_not_installed("Matrix") + sce <- make_sce() + cts <- SummarizedExperiment::assay(sce, "counts") + seu <- Seurat::CreateSeuratObject( + counts = Matrix::Matrix(as.matrix(cts), sparse = TRUE) + ) + md <- as.data.frame(SummarizedExperiment::colData(sce)) + for (nm in colnames(md)) seu[[nm]] <- md[[nm]] + seu +} + +# =========================================================================== +# .isSCObject / .assertSCObject +# =========================================================================== + +test_that(".isSCObject recognizes both supported classes", { + expect_true(immLynx:::.isSCObject(make_sce())) + expect_true(immLynx:::.isSCObject(make_seurat())) +}) + +test_that(".isSCObject rejects everything else", { + expect_false(immLynx:::.isSCObject(data.frame(x = 1))) + expect_false(immLynx:::.isSCObject(matrix(1))) + expect_false(immLynx:::.isSCObject(NULL)) +}) + +test_that(".assertSCObject names both accepted classes in its error", { + err <- tryCatch(immLynx:::.assertSCObject(data.frame(x = 1)), + error = function(e) conditionMessage(e)) + expect_match(err, "SingleCellExperiment") + expect_match(err, "Seurat") +}) + +test_that(".assertSCObject passes valid objects through silently", { + expect_silent(immLynx:::.assertSCObject(make_sce())) + expect_silent(immLynx:::.assertSCObject(make_seurat())) +}) + +# =========================================================================== +# .writeCellColumn +# =========================================================================== + +test_that(".writeCellColumn writes to colData for SingleCellExperiment", { + sce <- make_sce() + cells <- colnames(sce)[1:5] + + out <- immLynx:::.writeCellColumn(sce, "myscore", seq_len(5), cells) + + expect_true("myscore" %in% colnames(SummarizedExperiment::colData(out))) + vals <- SummarizedExperiment::colData(out)$myscore + expect_identical(length(vals), ncol(sce)) + expect_identical(vals[1:5], seq_len(5)) + # Cells outside the supplied set stay NA. + expect_true(all(is.na(vals[6:length(vals)]))) +}) + +test_that(".writeCellColumn writes to meta.data for Seurat", { + seu <- make_seurat() + cells <- colnames(seu)[1:5] + + out <- immLynx:::.writeCellColumn(seu, "myscore", seq_len(5), cells) + + expect_true("myscore" %in% colnames(out[[]])) + vals <- out[["myscore"]][, 1] + # Seurat's ncol() returns a double where SingleCellExperiment returns an + # integer, so compare by value rather than by type. + expect_equal(length(vals), ncol(seu)) + expect_identical(vals[1:5], seq_len(5)) + expect_true(all(is.na(vals[6:length(vals)]))) +}) + +test_that(".writeCellColumn returns the class it was given", { + sce <- make_sce() + seu <- make_seurat() + + expect_s4_class( + immLynx:::.writeCellColumn(sce, "x", 1, colnames(sce)[1]), + "SingleCellExperiment" + ) + expect_s4_class( + immLynx:::.writeCellColumn(seu, "x", 1, colnames(seu)[1]), + "Seurat" + ) +}) + +# =========================================================================== +# .writeReduction +# =========================================================================== + +test_that(".writeReduction stores a reduction on SingleCellExperiment", { + sce <- make_sce() + m <- matrix(rnorm(ncol(sce) * 4), nrow = ncol(sce), ncol = 4, + dimnames = list(colnames(sce), NULL)) + + out <- immLynx:::.writeReduction(sce, "myred", m, "MY_") + + expect_true("myred" %in% SingleCellExperiment::reducedDimNames(out)) + rd <- SingleCellExperiment::reducedDim(out, "myred") + expect_identical(dim(rd), c(ncol(sce), 4L)) +}) + +test_that(".writeReduction stores a reduction on Seurat", { + seu <- make_seurat() + m <- matrix(rnorm(ncol(seu) * 4), nrow = ncol(seu), ncol = 4, + dimnames = list(colnames(seu), NULL)) + + out <- immLynx:::.writeReduction(seu, "myred", m, "MY_") + + expect_true("myred" %in% Seurat::Reductions(out)) + emb <- Seurat::Embeddings(out, "myred") + expect_equal(dim(emb), c(ncol(seu), 4L)) + expect_identical(rownames(emb), colnames(seu)) +}) + +test_that(".writeReduction round-trips values unchanged", { + sce <- make_sce() + seu <- make_seurat() + m <- matrix(seq_len(ncol(sce) * 3) * 1.0, nrow = ncol(sce), ncol = 3, + dimnames = list(colnames(sce), NULL)) + + rd <- SingleCellExperiment::reducedDim( + immLynx:::.writeReduction(sce, "r", m, "R_"), "r") + emb <- Seurat::Embeddings( + immLynx:::.writeReduction(seu, "r", m, "R_"), "r") + + expect_equal(unname(rd), unname(m)) + expect_equal(unname(emb), unname(m)) +}) + +test_that(".writeReduction tolerates rows of NA for cells without sequences", { + sce <- make_sce() + m <- matrix(NA_real_, nrow = ncol(sce), ncol = 2, + dimnames = list(colnames(sce), NULL)) + m[1:10, ] <- 1 + + expect_silent(immLynx:::.writeReduction(sce, "r", m, "R_")) +}) + +# =========================================================================== +# .writeObjMetadata +# =========================================================================== + +test_that(".writeObjMetadata stores under metadata() for SingleCellExperiment", { + sce <- make_sce() + out <- immLynx:::.writeObjMetadata(sce, "mykey", list(a = 1)) + + expect_identical(S4Vectors::metadata(out)$mykey, list(a = 1)) +}) + +test_that(".writeObjMetadata stores under misc for Seurat", { + seu <- make_seurat() + out <- immLynx:::.writeObjMetadata(seu, "mykey", list(a = 1)) + + expect_identical(methods::slot(out, "misc")$mykey, list(a = 1)) +}) diff --git a/tests/testthat/test-runClustTCR.R b/tests/testthat/test-runClustTCR.R index 18e21db..dbea710 100644 --- a/tests/testthat/test-runClustTCR.R +++ b/tests/testthat/test-runClustTCR.R @@ -17,7 +17,7 @@ test_that("runClustTCR rejects non-SCE input", { expect_error( runClustTCR(tcr_data, chains = "TRB"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) @@ -26,7 +26,7 @@ test_that("runClustTCR rejects matrix input", { expect_error( runClustTCR(mat, chains = "TRB"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) @@ -35,7 +35,7 @@ test_that("runClustTCR rejects list input", { expect_error( runClustTCR(lst, chains = "TRB"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) diff --git a/tests/testthat/test-runEmbeddings.R b/tests/testthat/test-runEmbeddings.R index a9f5b8d..cc6ac65 100644 --- a/tests/testthat/test-runEmbeddings.R +++ b/tests/testthat/test-runEmbeddings.R @@ -21,7 +21,7 @@ test_that("runEmbeddings rejects non-SCE input", { expect_error( runEmbeddings(tcr_data, chains = "TRB"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) @@ -30,7 +30,7 @@ test_that("runEmbeddings rejects data.frame input", { expect_error( runEmbeddings(df, chains = "TRB"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) @@ -39,7 +39,7 @@ test_that("runEmbeddings rejects matrix input", { expect_error( runEmbeddings(mat, chains = "TRB"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) @@ -48,7 +48,7 @@ test_that("runEmbeddings rejects list input", { expect_error( runEmbeddings(lst, chains = "TRB"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) diff --git a/tests/testthat/test-runOLGA.R b/tests/testthat/test-runOLGA.R index 5d06f21..e60a92c 100644 --- a/tests/testthat/test-runOLGA.R +++ b/tests/testthat/test-runOLGA.R @@ -17,7 +17,7 @@ test_that("runOLGA rejects non-SCE input", { expect_error( runOLGA(tcr_data, chains = "TRB"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) diff --git a/tests/testthat/test-runSoNNia.R b/tests/testthat/test-runSoNNia.R index f0ae134..9f6087b 100644 --- a/tests/testthat/test-runSoNNia.R +++ b/tests/testthat/test-runSoNNia.R @@ -31,7 +31,7 @@ test_that("runSoNNia rejects non-SCE input", { expect_error( runSoNNia(tcr_data, chains = "TRB", background_file = temp_bg), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) diff --git a/tests/testthat/test-runTCRdist.R b/tests/testthat/test-runTCRdist.R index b19ed6e..f57093e 100644 --- a/tests/testthat/test-runTCRdist.R +++ b/tests/testthat/test-runTCRdist.R @@ -18,7 +18,7 @@ test_that("runTCRdist rejects non-SCE input", { expect_error( runTCRdist(tcr_data, chains = "beta"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) @@ -27,7 +27,7 @@ test_that("runTCRdist rejects data.frame input", { expect_error( runTCRdist(df, chains = "beta"), - "Input must be a SingleCellExperiment object" + "SingleCellExperiment or Seurat" ) }) diff --git a/tests/testthat/test-seurat-parity.R b/tests/testthat/test-seurat-parity.R new file mode 100644 index 0000000..9a29f3a --- /dev/null +++ b/tests/testthat/test-seurat-parity.R @@ -0,0 +1,144 @@ +# Every immLynx wrapper must accept a SingleCellExperiment or a Seurat object +# and return the class it was given. +# +# The extraction layer is already class-agnostic because immApex::getIR() +# handles both, so these tests target the guards and the write-back paths, +# which are where the two diverge. + +sce_fixture <- function() { + data("immLynx_example", package = "immLynx") + immLynx_example +} + +seurat_fixture <- function() { + skip_if_not_installed("Seurat") + skip_if_not_installed("Matrix") + sce <- sce_fixture() + cts <- SummarizedExperiment::assay(sce, "counts") + seu <- Seurat::CreateSeuratObject( + counts = Matrix::Matrix(as.matrix(cts), sparse = TRUE) + ) + md <- as.data.frame(SummarizedExperiment::colData(sce)) + for (nm in colnames(md)) seu[[nm]] <- md[[nm]] + seu +} + +# Read a per-cell column back out regardless of container. +cell_col <- function(obj, nm) { + if (methods::is(obj, "SingleCellExperiment")) { + SummarizedExperiment::colData(obj)[[nm]] + } else { + obj[[nm]][, 1] + } +} + +# =========================================================================== +# Guards: no wrapper may reject a Seurat object outright +# =========================================================================== + +test_that("no wrapper reports the old SingleCellExperiment-only error", { + # The historical guard read "Input must be a SingleCellExperiment object" + # and turned Seurat users away. Assert on behaviour rather than on source + # text, so this holds under R CMD check where the sources are not present. + wrappers <- list(runClustTCR, runOLGA, runEmbeddings, runTCRdist, + runSymdelNeighbors, runDeepTCR, runSoNNia) + + for (fn in wrappers) { + msg <- tryCatch(fn(data.frame(x = 1)), + error = function(e) conditionMessage(e)) + expect_false(grepl("must be a SingleCellExperiment object", msg, + fixed = TRUE)) + # It must still name Seurat as an accepted class. + expect_match(msg, "Seurat") + } +}) + +test_that("wrappers still reject objects that are neither class", { + for (fn in list(runClustTCR, runOLGA, runEmbeddings, runTCRdist, + runSymdelNeighbors, runDeepTCR, runSoNNia)) { + expect_error(fn(data.frame(x = 1)), + "SingleCellExperiment or Seurat") + } +}) + +test_that("the shared guard accepts a Seurat object", { + expect_silent(immLynx:::.assertSCObject(seurat_fixture())) +}) + +# =========================================================================== +# Round trip: class in, same class out +# =========================================================================== + +test_that("runSymdelNeighbors round-trips both classes", { + skip_if_no_python() + + sce <- runSymdelNeighbors(sce_fixture(), chains = "TRB", max_edits = 1) + seu <- runSymdelNeighbors(seurat_fixture(), chains = "TRB", max_edits = 1) + + expect_s4_class(sce, "SingleCellExperiment") + expect_s4_class(seu, "Seurat") + + a <- cell_col(sce, "symdel_degree") + b <- cell_col(seu, "symdel_degree") + expect_equal(a, b) +}) + +test_that("runClustTCR round-trips both classes", { + skip_if_no_python() + + sce <- runClustTCR(sce_fixture(), chains = "TRB") + seu <- runClustTCR(seurat_fixture(), chains = "TRB") + + expect_s4_class(sce, "SingleCellExperiment") + expect_s4_class(seu, "Seurat") + expect_true("clustcr_TRB" %in% colnames(seu[[]])) +}) + +test_that("runOLGA round-trips both classes", { + skip_if_no_python() + + sce <- runOLGA(sce_fixture(), chains = "TRB") + seu <- runOLGA(seurat_fixture(), chains = "TRB") + + expect_s4_class(sce, "SingleCellExperiment") + expect_s4_class(seu, "Seurat") + expect_equal(cell_col(sce, "olga_pgen_TRB"), + cell_col(seu, "olga_pgen_TRB")) +}) + +test_that("runEmbeddings writes a usable reduction on both classes", { + skip_if_no_python() + + sce <- runEmbeddings(sce_fixture(), chains = "TRB") + seu <- runEmbeddings(seurat_fixture(), chains = "TRB") + + expect_s4_class(sce, "SingleCellExperiment") + expect_s4_class(seu, "Seurat") + + expect_true("tcr_esm" %in% SingleCellExperiment::reducedDimNames(sce)) + expect_true("tcr_esm" %in% Seurat::Reductions(seu)) + + rd <- SingleCellExperiment::reducedDim(sce, "tcr_esm") + emb <- Seurat::Embeddings(seu, "tcr_esm") + expect_equal(dim(rd), dim(emb)) + expect_identical(rownames(emb), colnames(seu)) +}) + +test_that("runTCRdist stores results on both classes", { + skip_if_no_tcrdist() + + seu <- runTCRdist(seurat_fixture(), chains = "beta", add_to_object = TRUE) + + expect_s4_class(seu, "Seurat") + expect_true("tcrdist" %in% names(methods::slot(seu, "misc"))) +}) + +test_that("runDeepTCR writes a reduction on Seurat", { + skip_if_no_deeptcr() + + seu <- runDeepTCR(seurat_fixture(), chains = "TRB", + latent_dim = 16, verbose = FALSE) + + expect_s4_class(seu, "Seurat") + expect_true("tcr_deeptcr" %in% Seurat::Reductions(seu)) +}) From 5d0b29f11e5a15a2d47192a59a0f195d0aa31fa2 Mon Sep 17 00:00:00 2001 From: theHumanBorch Date: Mon, 10 Aug 2026 07:24:38 -0500 Subject: [PATCH 11/11] update unit tests --- tests/testthat/helper-immLynx.R | 14 ++++++++++++++ tests/testthat/test-runSymdelNeighbors.R | 8 ++++---- tests/testthat/test-seurat-parity.R | 2 +- tests/testthat/test-zz-basilisk-specs.R | 23 ++++++++++++++++++++--- 4 files changed, 39 insertions(+), 8 deletions(-) diff --git a/tests/testthat/helper-immLynx.R b/tests/testthat/helper-immLynx.R index bdb3719..484daa7 100644 --- a/tests/testthat/helper-immLynx.R +++ b/tests/testthat/helper-immLynx.R @@ -276,3 +276,17 @@ skip_if_no_deeptcr <- function() { }, error = function(e) FALSE) if (!ok) testthat::skip("deepTCREnv not available") } + +# Skip if pyrepseq's neighbor-search module cannot be imported. +# +# pyrepseq.nn imports pwseqdist at module level, which pulls numba and then +# llvmlite. llvmlite loads a compiled shared object, and on some CI images that +# load fails with "OSError: Could not find/load shared object file" even though +# the package itself installed cleanly. That is an environment problem rather +# than a defect in the wrapper, so the symdel tests skip instead of failing. +skip_if_no_symdel <- function() { + skip_if_no_python() + if (!can_import_module("pyrepseq.nn")) { + testthat::skip("pyrepseq.nn not importable (llvmlite/numba shared library)") + } +} diff --git a/tests/testthat/test-runSymdelNeighbors.R b/tests/testthat/test-runSymdelNeighbors.R index ad8b132..9c6f784 100644 --- a/tests/testthat/test-runSymdelNeighbors.R +++ b/tests/testthat/test-runSymdelNeighbors.R @@ -119,7 +119,7 @@ test_that("runSymdelNeighbors rejects an invalid max_edits", { # =========================================================================== test_that("runSymdelNeighbors returns an edge list", { - skip_if_no_python() + skip_if_no_symdel() data("immLynx_example", package = "immLynx") edges <- runSymdelNeighbors(immLynx_example, chains = "TRB", @@ -132,7 +132,7 @@ test_that("runSymdelNeighbors returns an edge list", { }) test_that("runSymdelNeighbors adds a degree column to colData", { - skip_if_no_python() + skip_if_no_symdel() data("immLynx_example", package = "immLynx") sce <- runSymdelNeighbors(immLynx_example, chains = "TRB", @@ -146,7 +146,7 @@ test_that("runSymdelNeighbors adds a degree column to colData", { }) test_that("runSymdelNeighbors honors column_prefix", { - skip_if_no_python() + skip_if_no_symdel() data("immLynx_example", package = "immLynx") sce <- runSymdelNeighbors(immLynx_example, chains = "TRB", @@ -156,7 +156,7 @@ test_that("runSymdelNeighbors honors column_prefix", { }) test_that("runSymdelNeighbors finds more neighbors at a larger max_edits", { - skip_if_no_python() + skip_if_no_symdel() data("immLynx_example", package = "immLynx") e1 <- runSymdelNeighbors(immLynx_example, chains = "TRB", diff --git a/tests/testthat/test-seurat-parity.R b/tests/testthat/test-seurat-parity.R index 9a29f3a..5aa9faf 100644 --- a/tests/testthat/test-seurat-parity.R +++ b/tests/testthat/test-seurat-parity.R @@ -70,7 +70,7 @@ test_that("the shared guard accepts a Seurat object", { # =========================================================================== test_that("runSymdelNeighbors round-trips both classes", { - skip_if_no_python() + skip_if_no_symdel() sce <- runSymdelNeighbors(sce_fixture(), chains = "TRB", max_edits = 1) seu <- runSymdelNeighbors(seurat_fixture(), chains = "TRB", max_edits = 1) diff --git a/tests/testthat/test-zz-basilisk-specs.R b/tests/testthat/test-zz-basilisk-specs.R index 80c51e4..fbe537f 100644 --- a/tests/testthat/test-zz-basilisk-specs.R +++ b/tests/testthat/test-zz-basilisk-specs.R @@ -32,11 +32,16 @@ CONDA_JUNK <- "[[:space:]\"'`$;&|()]" # Collect every BasiliskEnvironment object defined in the package namespace, # so environments added later are covered without editing this test. +# +# inherits = FALSE matters: a namespace's parent chain reaches the imports +# environment and eventually the search path, so the default lookup can return +# objects that do not belong to immLynx at all. basilisk_envs <- function() { ns <- asNamespace("immLynx") found <- list() for (nm in ls(ns, all.names = TRUE)) { - obj <- tryCatch(get(nm, envir = ns), error = function(e) NULL) + obj <- tryCatch(get(nm, envir = ns, inherits = FALSE), + error = function(e) NULL) if (methods::is(obj, "BasiliskEnvironment")) { found[[nm]] <- obj } @@ -44,6 +49,18 @@ basilisk_envs <- function() { found } +# Read a spec vector off an environment without assuming the slot exists. +# BasiliskEnvironment has gained and renamed slots across basilisk versions, +# and the package is built and checked against whatever the platform ships. +env_specs <- function(env, slot_name) { + if (!slot_name %in% methods::slotNames(env)) { + return(character(0)) + } + out <- tryCatch(methods::slot(env, slot_name), + error = function(e) character(0)) + if (is.null(out)) character(0) else as.character(out) +} + test_that("every basilisk environment is discoverable for inspection", { envs <- basilisk_envs() @@ -58,7 +75,7 @@ test_that("pip specs contain no shell metacharacters", { envs <- basilisk_envs() for (nm in names(envs)) { - specs <- envs[[nm]]@pip + specs <- env_specs(envs[[nm]], "pip") if (!length(specs)) next bad <- specs[!grepl(SHELL_SAFE, specs)] @@ -79,7 +96,7 @@ test_that("conda package specs are well formed", { envs <- basilisk_envs() for (nm in names(envs)) { - specs <- envs[[nm]]@packages + specs <- env_specs(envs[[nm]], "packages") if (!length(specs)) next bad <- specs[grepl(CONDA_JUNK, specs)]