From 89d36ac167707376dad110136d50ff17e4c6d00c Mon Sep 17 00:00:00 2001 From: Robrecht Cannoodt Date: Wed, 20 Aug 2025 16:52:51 +0200 Subject: [PATCH 1/4] add suerat components --- src/methods/seurat_cca/config.vsh.yaml | 58 ++++++++++ src/methods/seurat_cca/script.R | 112 +++++++++++++++++++ src/methods/seurat_rpca/config.vsh.yaml | 59 ++++++++++ src/methods/seurat_rpca/script.R | 136 ++++++++++++++++++++++++ 4 files changed, 365 insertions(+) create mode 100644 src/methods/seurat_cca/config.vsh.yaml create mode 100644 src/methods/seurat_cca/script.R create mode 100644 src/methods/seurat_rpca/config.vsh.yaml create mode 100644 src/methods/seurat_rpca/script.R diff --git a/src/methods/seurat_cca/config.vsh.yaml b/src/methods/seurat_cca/config.vsh.yaml new file mode 100644 index 000000000..f689d6faa --- /dev/null +++ b/src/methods/seurat_cca/config.vsh.yaml @@ -0,0 +1,58 @@ +__merge__: /src/api/comp_method.yaml +name: seurat_cca +label: Seurat CCA +summary: Integration using Seurat's anchor-based CCA integration +description: | + Seurat's Canonical Correlation Analysis (CCA) integration method identifies shared + correlation structures across datasets to perform batch correction. This method is + effective for datasets with shared cell types across conditions/batches. + + The method works by: + 1. Finding highly variable features for each dataset + 2. Identifying integration anchors using CCA + 3. Using anchors to harmonize datasets + 4. Generating integrated low-dimensional embedding +references: + # Stuart, T., Butler, A., Hoffman, P. et al. + # Comprehensive Integration of Single-Cell Data. + # Cell 177, 1888-1902.e21 (2019). https://doi.org/10.1016/j.cell.2019.05.031 + doi: 10.1016/j.cell.2019.05.031 +links: + repository: https://github.com/satijalab/seurat + documentation: https://satijalab.org/seurat/articles/seurat5_integration.html +info: + method_types: [embedding] + preferred_normalization: log_cp10k +arguments: + - name: --dims + type: integer + description: Number of dimensions to use for integration. + default: 30 + - name: --k_anchor + type: integer + description: Number of neighbors to use when picking anchors. + default: 5 + - name: --k_filter + type: integer + description: Number of neighbors to use when filtering anchors. + default: 200 + - name: --k_score + type: integer + description: Number of neighbors to use when scoring anchors. + default: 30 +resources: + - type: r_script + path: script.R +engines: + - type: docker + image: openproblems/base_r:1 + setup: + - type: r + cran: + - Seurat + - SeuratObject +runners: + - type: executable + - type: nextflow + directives: + label: [lowcpu, highmem, hightime] diff --git a/src/methods/seurat_cca/script.R b/src/methods/seurat_cca/script.R new file mode 100644 index 000000000..d57df0ee1 --- /dev/null +++ b/src/methods/seurat_cca/script.R @@ -0,0 +1,112 @@ +cat("Loading dependencies\n") +suppressPackageStartupMessages({ + requireNamespace("anndata", quietly = TRUE) + library(Matrix, warn.conflicts = FALSE) + library(Seurat, warn.conflicts = FALSE) + library(SeuratObject, warn.conflicts = FALSE) +}) + +## VIASH START +par <- list( + input = 'resources_test/task_batch_integration/cxg_immune_cell_atlas/dataset.h5ad', + output = 'output.h5ad', + dims = 30L, + k_anchor = 5L, + k_filter = 200L, + k_score = 30L +) +meta <- list( + name = "seurat_cca" +) +## VIASH END + +cat("Read input\n") +adata <- anndata::read_h5ad(par$input) + +cat("Create Seurat object using precomputed data\n") +# Extract preprocessed data +norm_data <- t(adata$layers[["normalized"]]) +obs <- adata$obs +var <- adata$var + +# Convert to dgCMatrix if needed (Seurat v5 compatibility) +if (inherits(norm_data, "dgRMatrix")) { + dense_temp <- as.matrix(norm_data) + norm_data <- as(dense_temp, "dgCMatrix") +} + +# Ensure proper dimnames for other matrix types +rownames(norm_data) <- rownames(var) +colnames(norm_data) <- rownames(obs) + +# Create Seurat object +seurat_obj <- CreateSeuratObject( + counts = norm_data, + meta.data = obs, + assay = "RNA" +) + +# In Seurat v5, we need to set the data layer for normalized data +seurat_obj[["RNA"]]$data <- norm_data + +cat("Set highly variable genes from input\n") +hvg_genes <- rownames(adata$var)[adata$var$hvg] +cat("Using", length(hvg_genes), "HVGs from input dataset\n") +VariableFeatures(seurat_obj) <- hvg_genes + +cat("Split by batch and perform CCA integration\n") +# Split the object by batch +seurat_list <- SplitObject(seurat_obj, split.by = "batch") + +# Find integration anchors using CCA +anchors <- FindIntegrationAnchors( + object.list = seurat_list, + anchor.features = hvg_genes, + dims = seq_len(par$dims), + k.anchor = par$k_anchor, + k.filter = par$k_filter, + k.score = par$k_score, + verbose = FALSE +) + +# Integrate the data +integrated <- IntegrateData( + anchorset = anchors, + dims = seq_len(par$dims), + verbose = FALSE +) + +cat("Scale integrated data and run PCA\n") +DefaultAssay(integrated) <- "integrated" +integrated <- ScaleData(integrated, verbose = FALSE) +integrated <- RunPCA(integrated, npcs = par$dims, verbose = FALSE) + +cat("Generate UMAP embedding\n") +integrated <- RunUMAP( + integrated, + reduction = "pca", + dims = seq_len(par$dims), + verbose = FALSE +) + +cat("Extract embedding\n") +embedding <- Embeddings(integrated, reduction = "umap") + +cat("Store outputs\n") +output <- anndata::AnnData( + obs = adata$obs, + var = adata$var, + obsm = list( + X_emb = embedding + ), + uns = list( + dataset_id = adata$uns[["dataset_id"]], + normalization_id = adata$uns[["normalization_id"]], + method_id = meta$name + ) +) + +cat("Write output to file\n") +zzz <- output$write_h5ad(par$output, compression = "gzip") + +cat("Finished\n") diff --git a/src/methods/seurat_rpca/config.vsh.yaml b/src/methods/seurat_rpca/config.vsh.yaml new file mode 100644 index 000000000..fa9df05c6 --- /dev/null +++ b/src/methods/seurat_rpca/config.vsh.yaml @@ -0,0 +1,59 @@ +__merge__: /src/api/comp_method.yaml +name: seurat_rpca +label: Seurat RPCA +summary: Integration using Seurat's anchor-based RPCA integration +description: | + Seurat's Reciprocal PCA (RPCA) integration method represents a faster and more + conservative (less correction) approach compared to CCA integration. This method + is especially useful for large datasets or datasets with similar biological + composition across batches. + + The method works by: + 1. Finding highly variable features for each dataset + 2. Running PCA on each dataset + 3. Identifying integration anchors using RPCA + 4. Using anchors to harmonize datasets with more conservative correction +references: + # Stuart, T., Butler, A., Hoffman, P. et al. + # Comprehensive Integration of Single-Cell Data. + # Cell 177, 1888-1902.e21 (2019). https://doi.org/10.1016/j.cell.2019.05.031 + doi: 10.1016/j.cell.2019.05.031 +links: + repository: https://github.com/satijalab/seurat + documentation: https://satijalab.org/seurat/articles/seurat5_integration.html +info: + method_types: [embedding] + preferred_normalization: log_cp10k +arguments: + - name: --dims + type: integer + description: Number of dimensions to use for integration. + default: 30 + - name: --k_anchor + type: integer + description: Number of neighbors to use when picking anchors. + default: 5 + - name: --k_filter + type: integer + description: Number of neighbors to use when filtering anchors. + default: 200 + - name: --k_score + type: integer + description: Number of neighbors to use when scoring anchors. + default: 30 +resources: + - type: r_script + path: script.R +engines: + - type: docker + image: openproblems/base_r:1 + setup: + - type: r + cran: + - Seurat + - SeuratObject +runners: + - type: executable + - type: nextflow + directives: + label: [lowcpu, highmem, hightime] diff --git a/src/methods/seurat_rpca/script.R b/src/methods/seurat_rpca/script.R new file mode 100644 index 000000000..302db7924 --- /dev/null +++ b/src/methods/seurat_rpca/script.R @@ -0,0 +1,136 @@ +cat("Loading dependencies\n") +suppressPackageStartupMessages({ + requireNamespace("anndata", quietly = TRUE) + library(Matrix, warn.conflicts = FALSE) + library(Seurat, warn.conflicts = FALSE) + library(SeuratObject, warn.conflicts = FALSE) +}) + +## VIASH START +par <- list( + input = 'resources_test/task_batch_integration/cxg_immune_cell_atlas/dataset.h5ad', + output = 'output.h5ad', + dims = 30L, + k_anchor = 5L, + k_filter = 200L, + k_score = 30L +) +meta <- list( + name = "seurat_rpca" +) +## VIASH END + +cat("Read input\n") +adata <- anndata::read_h5ad(par$input) + +cat("Create Seurat object\n") +# Extract preprocessed data +counts_data <- t(adata$layers[["counts"]]) +norm_data <- t(adata$layers[["normalized"]]) +obs <- adata$obs +var <- adata$var + +# Convert to dgCMatrix if needed (Seurat v5 compatibility) +if (inherits(counts_data, "dgRMatrix")) { + dense_temp <- as.matrix(counts_data) + counts_data <- as(dense_temp, "dgCMatrix") +} + +if (inherits(norm_data, "dgRMatrix")) { + dense_temp <- as.matrix(norm_data) + norm_data <- as(dense_temp, "dgCMatrix") +} + +# Ensure proper dimnames +rownames(counts_data) <- rownames(var) +colnames(counts_data) <- rownames(obs) +rownames(norm_data) <- rownames(var) +colnames(norm_data) <- rownames(obs) + +# Create Seurat object from counts +seurat_obj <- CreateSeuratObject( + counts = counts_data, + meta.data = obs, + assay = "RNA" +) + +# Add normalized data layer +LayerData(seurat_obj, layer = "data") <- norm_data + +# Use existing HVGs from the dataset +hvg_genes <- rownames(adata$var)[adata$var$hvg] +cat("Using", length(hvg_genes), "HVGs from input dataset\n") +VariableFeatures(seurat_obj) <- hvg_genes + +# Use existing PCA from input dataset +pca_embeddings <- adata$obsm[["X_pca"]] +rownames(pca_embeddings) <- colnames(seurat_obj) +colnames(pca_embeddings) <- paste0("PC_", seq_len(ncol(pca_embeddings))) + +seurat_obj[["pca"]] <- CreateDimReducObject( + embeddings = pca_embeddings, + key = "PC_", + assay = "RNA" +) + +cat("Split object by batch\n") +# Split the object by batch for traditional integration +seurat_list <- SplitObject(seurat_obj, split.by = "batch") + +# For RPCA, we need to scale and run PCA on each dataset +cat("Scale data and run PCA for RPCA integration\n") +seurat_list <- lapply(seurat_list, function(x) { + x <- ScaleData(x, features = hvg_genes, verbose = FALSE) + x <- RunPCA(x, features = hvg_genes, npcs = par$dims, verbose = FALSE) + return(x) +}) + +cat("Perform RPCA integration\n") +# Find integration anchors using RPCA +anchors <- FindIntegrationAnchors( + object.list = seurat_list, + anchor.features = hvg_genes, + reduction = "rpca", + dims = seq_len(par$dims), + k.anchor = par$k_anchor, + k.filter = par$k_filter, + k.score = par$k_score, + verbose = FALSE +) + +# Integrate the data +integrated <- IntegrateData( + anchorset = anchors, + dims = seq_len(par$dims), + verbose = FALSE +) + +cat("Scale and run PCA on integrated data\n") +DefaultAssay(integrated) <- "integrated" +integrated <- ScaleData(integrated, verbose = FALSE) +integrated <- RunPCA(integrated, npcs = par$dims, verbose = FALSE) + +cat("Generate UMAP embedding\n") +integrated <- RunUMAP(integrated, reduction = "pca", dims = seq_len(par$dims), verbose = FALSE) + +cat("Extract embedding\n") +embedding <- Embeddings(integrated, reduction = "umap") + +cat("Store outputs\n") +output <- anndata::AnnData( + obs = adata$obs[, c()], + var = adata$var[, c()], + obsm = list( + X_emb = embedding + ), + uns = list( + dataset_id = adata$uns[["dataset_id"]], + normalization_id = adata$uns[["normalization_id"]], + method_id = meta$name + ) +) + +cat("Write output to file\n") +zzz <- output$write_h5ad(par$output, compression = "gzip") + +cat("Finished\n") From 0a32705534dcd1bb56807e7a7e3a1d7a4360f4e4 Mon Sep 17 00:00:00 2001 From: Robrecht Cannoodt Date: Wed, 20 Aug 2025 16:53:15 +0200 Subject: [PATCH 2/4] update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 041bd7299..e884bf8f0 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,6 +4,7 @@ * Added `metrics/kbet_pg` and `metrics/kbet_pg_label` components (PR #52). * Added `method/drvi` component (PR #61). +* Added `method/seurat_cca` and `method/seurat_rpca` components (PR #77). ## Minor changes From 43ba846db629119ac301679c67f8d8d863beb8f1 Mon Sep 17 00:00:00 2001 From: Robrecht Cannoodt Date: Fri, 21 Aug 2026 13:58:23 +0200 Subject: [PATCH 3/4] Rewrite seurat_cca and seurat_rpca on the Seurat v5 IntegrateLayers workflow Follow the Seurat v5 integration vignette: split the RNA assay into one layer per batch, ScaleData + RunPCA, then IntegrateLayers with CCAIntegration / RPCAIntegration and return the corrected reduction as X_emb. This replaces the v4 SplitObject / FindIntegrationAnchors / IntegrateData path and drops the UMAP that was previously exported as the embedding. Also: - convert dgRMatrix to dgCMatrix sparse-to-sparse instead of via a dense copy - stop loading the counts layer in seurat_rpca; it was unused - forward --dims/--k_anchor/--k_filter/--k_score only when set so Seurat's own defaults apply (k.filter is NA in v5, not 200) - drop the unused X_pca reduction in seurat_rpca --- src/methods/seurat_cca/config.vsh.yaml | 44 +++---- src/methods/seurat_cca/script.R | 145 ++++++++------------- src/methods/seurat_rpca/config.vsh.yaml | 46 +++---- src/methods/seurat_rpca/script.R | 165 ++++++++---------------- 4 files changed, 155 insertions(+), 245 deletions(-) diff --git a/src/methods/seurat_cca/config.vsh.yaml b/src/methods/seurat_cca/config.vsh.yaml index f689d6faa..4c9ec2680 100644 --- a/src/methods/seurat_cca/config.vsh.yaml +++ b/src/methods/seurat_cca/config.vsh.yaml @@ -1,17 +1,18 @@ __merge__: /src/api/comp_method.yaml name: seurat_cca label: Seurat CCA -summary: Integration using Seurat's anchor-based CCA integration +summary: Seurat v5 anchor-based integration using canonical correlation analysis description: | - Seurat's Canonical Correlation Analysis (CCA) integration method identifies shared - correlation structures across datasets to perform batch correction. This method is - effective for datasets with shared cell types across conditions/batches. - - The method works by: - 1. Finding highly variable features for each dataset - 2. Identifying integration anchors using CCA - 3. Using anchors to harmonize datasets - 4. Generating integrated low-dimensional embedding + Seurat's anchor-based integration identifies pairs of mutual nearest neighbours + ("anchors") between batches in a shared low-dimensional space and uses them to + correct the PCA embedding. Here the shared space is found with canonical + correlation analysis (CCA), which is recommended when cell types are largely + shared across batches but there are strong batch effects. + + This component runs the Seurat v5 workflow: the expression matrix is split + into one layer per batch, scaled and reduced with PCA, and then integrated + with `IntegrateLayers(method = CCAIntegration)`. The corrected embedding + (`integrated.cca`) is returned. references: # Stuart, T., Butler, A., Hoffman, P. et al. # Comprehensive Integration of Single-Cell Data. @@ -19,27 +20,27 @@ references: doi: 10.1016/j.cell.2019.05.031 links: repository: https://github.com/satijalab/seurat - documentation: https://satijalab.org/seurat/articles/seurat5_integration.html + documentation: https://satijalab.org/seurat/articles/seurat5_integration info: method_types: [embedding] preferred_normalization: log_cp10k arguments: - name: --dims type: integer - description: Number of dimensions to use for integration. - default: 30 + description: Number of dimensions to use for integration. Defaults to Seurat's default (30). + example: 30 - name: --k_anchor type: integer - description: Number of neighbors to use when picking anchors. - default: 5 + description: Number of neighbors to use when picking anchors. Defaults to Seurat's default (5). + example: 5 - name: --k_filter type: integer - description: Number of neighbors to use when filtering anchors. - default: 200 + description: Number of neighbors to use when filtering anchors. Defaults to Seurat's default (no filtering). + example: 200 - name: --k_score - type: integer - description: Number of neighbors to use when scoring anchors. - default: 30 + type: integer + description: Number of neighbors to use when scoring anchors. Defaults to Seurat's default (30). + example: 30 resources: - type: r_script path: script.R @@ -50,9 +51,8 @@ engines: - type: r cran: - Seurat - - SeuratObject runners: - type: executable - type: nextflow directives: - label: [lowcpu, highmem, hightime] + label: [midcpu, highmem, hightime] diff --git a/src/methods/seurat_cca/script.R b/src/methods/seurat_cca/script.R index d57df0ee1..5c1332832 100644 --- a/src/methods/seurat_cca/script.R +++ b/src/methods/seurat_cca/script.R @@ -1,112 +1,79 @@ -cat("Loading dependencies\n") +requireNamespace("anndata", quietly = TRUE) suppressPackageStartupMessages({ - requireNamespace("anndata", quietly = TRUE) - library(Matrix, warn.conflicts = FALSE) - library(Seurat, warn.conflicts = FALSE) - library(SeuratObject, warn.conflicts = FALSE) + library(Matrix) + library(SeuratObject) + library(Seurat) }) ## VIASH START par <- list( - input = 'resources_test/task_batch_integration/cxg_immune_cell_atlas/dataset.h5ad', - output = 'output.h5ad', - dims = 30L, - k_anchor = 5L, - k_filter = 200L, - k_score = 30L + input = "resources_test/task_batch_integration/cxg_immune_cell_atlas/dataset.h5ad", + output = "output.h5ad", + dims = NULL, + k_anchor = NULL, + k_filter = NULL, + k_score = NULL ) meta <- list( name = "seurat_cca" ) ## VIASH END -cat("Read input\n") -adata <- anndata::read_h5ad(par$input) - -cat("Create Seurat object using precomputed data\n") -# Extract preprocessed data -norm_data <- t(adata$layers[["normalized"]]) -obs <- adata$obs -var <- adata$var - -# Convert to dgCMatrix if needed (Seurat v5 compatibility) -if (inherits(norm_data, "dgRMatrix")) { - dense_temp <- as.matrix(norm_data) - norm_data <- as(dense_temp, "dgCMatrix") -} - -# Ensure proper dimnames for other matrix types -rownames(norm_data) <- rownames(var) -colnames(norm_data) <- rownames(obs) - -# Create Seurat object -seurat_obj <- CreateSeuratObject( - counts = norm_data, - meta.data = obs, - assay = "RNA" -) - -# In Seurat v5, we need to set the data layer for normalized data -seurat_obj[["RNA"]]$data <- norm_data - -cat("Set highly variable genes from input\n") -hvg_genes <- rownames(adata$var)[adata$var$hvg] -cat("Using", length(hvg_genes), "HVGs from input dataset\n") -VariableFeatures(seurat_obj) <- hvg_genes - -cat("Split by batch and perform CCA integration\n") -# Split the object by batch -seurat_list <- SplitObject(seurat_obj, split.by = "batch") - -# Find integration anchors using CCA -anchors <- FindIntegrationAnchors( - object.list = seurat_list, - anchor.features = hvg_genes, - dims = seq_len(par$dims), - k.anchor = par$k_anchor, - k.filter = par$k_filter, - k.score = par$k_score, - verbose = FALSE -) - -# Integrate the data -integrated <- IntegrateData( - anchorset = anchors, - dims = seq_len(par$dims), - verbose = FALSE -) - -cat("Scale integrated data and run PCA\n") -DefaultAssay(integrated) <- "integrated" -integrated <- ScaleData(integrated, verbose = FALSE) -integrated <- RunPCA(integrated, npcs = par$dims, verbose = FALSE) - -cat("Generate UMAP embedding\n") -integrated <- RunUMAP( - integrated, - reduction = "pca", - dims = seq_len(par$dims), +cat("Reading input file\n") +adata <- anndata::read_h5ad(par[["input"]]) + +cat("Create Seurat object\n") +# Seurat expects genes in rows, cells in columns, as a dgCMatrix +normalized <- Matrix::t(adata$layers[["normalized"]]) +normalized <- as(as(normalized, "CsparseMatrix"), "dgCMatrix") + +seurat_obj <- Seurat::CreateSeuratObject(counts = normalized, meta.data = adata$obs) +# The benchmark's log_cp10k normalization is Seurat's LogNormalize, so assign it to +# the "data" layer instead of calling NormalizeData(). +seurat_obj[["RNA"]]$data <- normalized +seurat_obj[["RNA"]]$counts <- NULL + +# Use the benchmark's HVGs instead of FindVariableFeatures() so that feature +# selection is the same across methods. +VariableFeatures(seurat_obj) <- rownames(adata$var)[adata$var$hvg] + +cat("Split layers by batch, scale and run PCA\n") +# Seurat v5 integration workflow, see +# https://satijalab.org/seurat/articles/seurat5_integration +seurat_obj[["RNA"]] <- split(seurat_obj[["RNA"]], f = seurat_obj$batch) +seurat_obj <- ScaleData(seurat_obj, verbose = FALSE) +pca_args <- list(object = seurat_obj, verbose = FALSE) +if (!is.null(par$dims)) pca_args$npcs <- par$dims +seurat_obj <- do.call(RunPCA, pca_args) + +cat("Run CCAIntegration\n") +# Only forward arguments that were set so Seurat's own defaults apply otherwise +integrate_args <- list( + object = seurat_obj, + method = CCAIntegration, + orig.reduction = "pca", + new.reduction = "integrated.cca", verbose = FALSE ) - -cat("Extract embedding\n") -embedding <- Embeddings(integrated, reduction = "umap") +if (!is.null(par$dims)) integrate_args$dims <- seq_len(par$dims) +if (!is.null(par$k_anchor)) integrate_args$k.anchor <- par$k_anchor +if (!is.null(par$k_filter)) integrate_args$k.filter <- par$k_filter +if (!is.null(par$k_score)) integrate_args$k.score <- par$k_score +seurat_obj <- do.call(IntegrateLayers, integrate_args) cat("Store outputs\n") output <- anndata::AnnData( - obs = adata$obs, - var = adata$var, - obsm = list( - X_emb = embedding - ), uns = list( dataset_id = adata$uns[["dataset_id"]], normalization_id = adata$uns[["normalization_id"]], method_id = meta$name + ), + obs = adata$obs, + var = adata$var, + obsm = list( + X_emb = Embeddings(seurat_obj, reduction = "integrated.cca") ) ) -cat("Write output to file\n") -zzz <- output$write_h5ad(par$output, compression = "gzip") - -cat("Finished\n") +cat("Write output AnnData to file\n") +output$write_h5ad(par[["output"]], compression = "gzip") diff --git a/src/methods/seurat_rpca/config.vsh.yaml b/src/methods/seurat_rpca/config.vsh.yaml index fa9df05c6..72486d098 100644 --- a/src/methods/seurat_rpca/config.vsh.yaml +++ b/src/methods/seurat_rpca/config.vsh.yaml @@ -1,18 +1,19 @@ __merge__: /src/api/comp_method.yaml name: seurat_rpca label: Seurat RPCA -summary: Integration using Seurat's anchor-based RPCA integration +summary: Seurat v5 anchor-based integration using reciprocal PCA description: | - Seurat's Reciprocal PCA (RPCA) integration method represents a faster and more - conservative (less correction) approach compared to CCA integration. This method - is especially useful for large datasets or datasets with similar biological - composition across batches. - - The method works by: - 1. Finding highly variable features for each dataset - 2. Running PCA on each dataset - 3. Identifying integration anchors using RPCA - 4. Using anchors to harmonize datasets with more conservative correction + Seurat's anchor-based integration identifies pairs of mutual nearest neighbours + ("anchors") between batches in a shared low-dimensional space and uses them to + correct the PCA embedding. Here the shared space is found with reciprocal PCA + (RPCA), where each batch is projected into the others' PCA space. RPCA is faster + and more conservative than CCA, and is recommended for large datasets or when + a substantial fraction of cells in one batch has no match in another. + + This component runs the Seurat v5 workflow: the expression matrix is split + into one layer per batch, scaled and reduced with PCA, and then integrated + with `IntegrateLayers(method = RPCAIntegration)`. The corrected embedding + (`integrated.rpca`) is returned. references: # Stuart, T., Butler, A., Hoffman, P. et al. # Comprehensive Integration of Single-Cell Data. @@ -20,27 +21,27 @@ references: doi: 10.1016/j.cell.2019.05.031 links: repository: https://github.com/satijalab/seurat - documentation: https://satijalab.org/seurat/articles/seurat5_integration.html + documentation: https://satijalab.org/seurat/articles/seurat5_integration info: method_types: [embedding] preferred_normalization: log_cp10k arguments: - name: --dims type: integer - description: Number of dimensions to use for integration. - default: 30 + description: Number of dimensions to use for integration. Defaults to Seurat's default (30). + example: 30 - name: --k_anchor type: integer - description: Number of neighbors to use when picking anchors. - default: 5 + description: Number of neighbors to use when picking anchors. Defaults to Seurat's default (5). + example: 5 - name: --k_filter type: integer - description: Number of neighbors to use when filtering anchors. - default: 200 + description: Number of neighbors to use when filtering anchors. Defaults to Seurat's default (no filtering). + example: 200 - name: --k_score - type: integer - description: Number of neighbors to use when scoring anchors. - default: 30 + type: integer + description: Number of neighbors to use when scoring anchors. Defaults to Seurat's default (30). + example: 30 resources: - type: r_script path: script.R @@ -51,9 +52,8 @@ engines: - type: r cran: - Seurat - - SeuratObject runners: - type: executable - type: nextflow directives: - label: [lowcpu, highmem, hightime] + label: [midcpu, midmem, midtime] diff --git a/src/methods/seurat_rpca/script.R b/src/methods/seurat_rpca/script.R index 302db7924..102eb34b1 100644 --- a/src/methods/seurat_rpca/script.R +++ b/src/methods/seurat_rpca/script.R @@ -1,136 +1,79 @@ -cat("Loading dependencies\n") +requireNamespace("anndata", quietly = TRUE) suppressPackageStartupMessages({ - requireNamespace("anndata", quietly = TRUE) - library(Matrix, warn.conflicts = FALSE) - library(Seurat, warn.conflicts = FALSE) - library(SeuratObject, warn.conflicts = FALSE) + library(Matrix) + library(SeuratObject) + library(Seurat) }) ## VIASH START par <- list( - input = 'resources_test/task_batch_integration/cxg_immune_cell_atlas/dataset.h5ad', - output = 'output.h5ad', - dims = 30L, - k_anchor = 5L, - k_filter = 200L, - k_score = 30L + input = "resources_test/task_batch_integration/cxg_immune_cell_atlas/dataset.h5ad", + output = "output.h5ad", + dims = NULL, + k_anchor = NULL, + k_filter = NULL, + k_score = NULL ) meta <- list( name = "seurat_rpca" ) ## VIASH END -cat("Read input\n") -adata <- anndata::read_h5ad(par$input) +cat("Reading input file\n") +adata <- anndata::read_h5ad(par[["input"]]) cat("Create Seurat object\n") -# Extract preprocessed data -counts_data <- t(adata$layers[["counts"]]) -norm_data <- t(adata$layers[["normalized"]]) -obs <- adata$obs -var <- adata$var - -# Convert to dgCMatrix if needed (Seurat v5 compatibility) -if (inherits(counts_data, "dgRMatrix")) { - dense_temp <- as.matrix(counts_data) - counts_data <- as(dense_temp, "dgCMatrix") -} - -if (inherits(norm_data, "dgRMatrix")) { - dense_temp <- as.matrix(norm_data) - norm_data <- as(dense_temp, "dgCMatrix") -} - -# Ensure proper dimnames -rownames(counts_data) <- rownames(var) -colnames(counts_data) <- rownames(obs) -rownames(norm_data) <- rownames(var) -colnames(norm_data) <- rownames(obs) - -# Create Seurat object from counts -seurat_obj <- CreateSeuratObject( - counts = counts_data, - meta.data = obs, - assay = "RNA" -) - -# Add normalized data layer -LayerData(seurat_obj, layer = "data") <- norm_data - -# Use existing HVGs from the dataset -hvg_genes <- rownames(adata$var)[adata$var$hvg] -cat("Using", length(hvg_genes), "HVGs from input dataset\n") -VariableFeatures(seurat_obj) <- hvg_genes - -# Use existing PCA from input dataset -pca_embeddings <- adata$obsm[["X_pca"]] -rownames(pca_embeddings) <- colnames(seurat_obj) -colnames(pca_embeddings) <- paste0("PC_", seq_len(ncol(pca_embeddings))) - -seurat_obj[["pca"]] <- CreateDimReducObject( - embeddings = pca_embeddings, - key = "PC_", - assay = "RNA" -) - -cat("Split object by batch\n") -# Split the object by batch for traditional integration -seurat_list <- SplitObject(seurat_obj, split.by = "batch") - -# For RPCA, we need to scale and run PCA on each dataset -cat("Scale data and run PCA for RPCA integration\n") -seurat_list <- lapply(seurat_list, function(x) { - x <- ScaleData(x, features = hvg_genes, verbose = FALSE) - x <- RunPCA(x, features = hvg_genes, npcs = par$dims, verbose = FALSE) - return(x) -}) - -cat("Perform RPCA integration\n") -# Find integration anchors using RPCA -anchors <- FindIntegrationAnchors( - object.list = seurat_list, - anchor.features = hvg_genes, - reduction = "rpca", - dims = seq_len(par$dims), - k.anchor = par$k_anchor, - k.filter = par$k_filter, - k.score = par$k_score, - verbose = FALSE -) - -# Integrate the data -integrated <- IntegrateData( - anchorset = anchors, - dims = seq_len(par$dims), +# Seurat expects genes in rows, cells in columns, as a dgCMatrix +normalized <- Matrix::t(adata$layers[["normalized"]]) +normalized <- as(as(normalized, "CsparseMatrix"), "dgCMatrix") + +seurat_obj <- Seurat::CreateSeuratObject(counts = normalized, meta.data = adata$obs) +# The benchmark's log_cp10k normalization is Seurat's LogNormalize, so assign it to +# the "data" layer instead of calling NormalizeData(). +seurat_obj[["RNA"]]$data <- normalized +seurat_obj[["RNA"]]$counts <- NULL + +# Use the benchmark's HVGs instead of FindVariableFeatures() so that feature +# selection is the same across methods. +VariableFeatures(seurat_obj) <- rownames(adata$var)[adata$var$hvg] + +cat("Split layers by batch, scale and run PCA\n") +# Seurat v5 integration workflow, see +# https://satijalab.org/seurat/articles/seurat5_integration +seurat_obj[["RNA"]] <- split(seurat_obj[["RNA"]], f = seurat_obj$batch) +seurat_obj <- ScaleData(seurat_obj, verbose = FALSE) +pca_args <- list(object = seurat_obj, verbose = FALSE) +if (!is.null(par$dims)) pca_args$npcs <- par$dims +seurat_obj <- do.call(RunPCA, pca_args) + +cat("Run RPCAIntegration\n") +# Only forward arguments that were set so Seurat's own defaults apply otherwise +integrate_args <- list( + object = seurat_obj, + method = RPCAIntegration, + orig.reduction = "pca", + new.reduction = "integrated.rpca", verbose = FALSE ) - -cat("Scale and run PCA on integrated data\n") -DefaultAssay(integrated) <- "integrated" -integrated <- ScaleData(integrated, verbose = FALSE) -integrated <- RunPCA(integrated, npcs = par$dims, verbose = FALSE) - -cat("Generate UMAP embedding\n") -integrated <- RunUMAP(integrated, reduction = "pca", dims = seq_len(par$dims), verbose = FALSE) - -cat("Extract embedding\n") -embedding <- Embeddings(integrated, reduction = "umap") +if (!is.null(par$dims)) integrate_args$dims <- seq_len(par$dims) +if (!is.null(par$k_anchor)) integrate_args$k.anchor <- par$k_anchor +if (!is.null(par$k_filter)) integrate_args$k.filter <- par$k_filter +if (!is.null(par$k_score)) integrate_args$k.score <- par$k_score +seurat_obj <- do.call(IntegrateLayers, integrate_args) cat("Store outputs\n") output <- anndata::AnnData( - obs = adata$obs[, c()], - var = adata$var[, c()], - obsm = list( - X_emb = embedding - ), uns = list( dataset_id = adata$uns[["dataset_id"]], normalization_id = adata$uns[["normalization_id"]], method_id = meta$name + ), + obs = adata$obs, + var = adata$var, + obsm = list( + X_emb = Embeddings(seurat_obj, reduction = "integrated.rpca") ) ) -cat("Write output to file\n") -zzz <- output$write_h5ad(par$output, compression = "gzip") - -cat("Finished\n") +cat("Write output AnnData to file\n") +output$write_h5ad(par[["output"]], compression = "gzip") From 76cc091f275e6e03006f7ade3f2b34c04e101562 Mon Sep 17 00:00:00 2001 From: Robrecht Cannoodt Date: Fri, 21 Aug 2026 14:05:34 +0200 Subject: [PATCH 4/4] Use anndataR in the seurat components read_h5ad() + as_Seurat() with layers_mapping = c(data = "normalized") builds the Seurat object directly from the h5ad, replacing the manual transpose / sparse conversion / CreateSeuratObject / drop-counts block. Also cite the Seurat v4 paper for seurat_rpca, where RPCA integration was introduced. --- src/methods/seurat_cca/script.R | 40 ++++++++++++------------- src/methods/seurat_rpca/config.vsh.yaml | 7 ++++- src/methods/seurat_rpca/script.R | 40 ++++++++++++------------- 3 files changed, 44 insertions(+), 43 deletions(-) diff --git a/src/methods/seurat_cca/script.R b/src/methods/seurat_cca/script.R index 5c1332832..41901025c 100644 --- a/src/methods/seurat_cca/script.R +++ b/src/methods/seurat_cca/script.R @@ -1,7 +1,5 @@ -requireNamespace("anndata", quietly = TRUE) suppressPackageStartupMessages({ - library(Matrix) - library(SeuratObject) + library(anndataR) library(Seurat) }) @@ -20,22 +18,22 @@ meta <- list( ## VIASH END cat("Reading input file\n") -adata <- anndata::read_h5ad(par[["input"]]) +adata <- read_h5ad(par[["input"]]) cat("Create Seurat object\n") -# Seurat expects genes in rows, cells in columns, as a dgCMatrix -normalized <- Matrix::t(adata$layers[["normalized"]]) -normalized <- as(as(normalized, "CsparseMatrix"), "dgCMatrix") - -seurat_obj <- Seurat::CreateSeuratObject(counts = normalized, meta.data = adata$obs) -# The benchmark's log_cp10k normalization is Seurat's LogNormalize, so assign it to +# The benchmark's log_cp10k normalization is Seurat's LogNormalize, so map it to # the "data" layer instead of calling NormalizeData(). -seurat_obj[["RNA"]]$data <- normalized -seurat_obj[["RNA"]]$counts <- NULL - +seurat_obj <- adata$as_Seurat( + x_mapping = NULL, + layers_mapping = c(data = "normalized"), + assay_metadata_mapping = FALSE, + reduction_mapping = FALSE, + graph_mapping = FALSE, + misc_mapping = FALSE +) # Use the benchmark's HVGs instead of FindVariableFeatures() so that feature # selection is the same across methods. -VariableFeatures(seurat_obj) <- rownames(adata$var)[adata$var$hvg] +VariableFeatures(seurat_obj) <- adata$var_names[adata$var$hvg] cat("Split layers by batch, scale and run PCA\n") # Seurat v5 integration workflow, see @@ -62,18 +60,18 @@ if (!is.null(par$k_score)) integrate_args$k.score <- par$k_score seurat_obj <- do.call(IntegrateLayers, integrate_args) cat("Store outputs\n") -output <- anndata::AnnData( +output <- AnnData( + obs = adata$obs[, character(0)], + var = adata$var[, character(0)], + obsm = list( + X_emb = Embeddings(seurat_obj, reduction = "integrated.cca") + ), uns = list( dataset_id = adata$uns[["dataset_id"]], normalization_id = adata$uns[["normalization_id"]], method_id = meta$name - ), - obs = adata$obs, - var = adata$var, - obsm = list( - X_emb = Embeddings(seurat_obj, reduction = "integrated.cca") ) ) cat("Write output AnnData to file\n") -output$write_h5ad(par[["output"]], compression = "gzip") +output$write_h5ad(par[["output"]]) diff --git a/src/methods/seurat_rpca/config.vsh.yaml b/src/methods/seurat_rpca/config.vsh.yaml index 72486d098..73ae06aca 100644 --- a/src/methods/seurat_rpca/config.vsh.yaml +++ b/src/methods/seurat_rpca/config.vsh.yaml @@ -15,10 +15,15 @@ description: | with `IntegrateLayers(method = RPCAIntegration)`. The corrected embedding (`integrated.rpca`) is returned. references: + # Hao, Y., Hao, S., Andersen-Nissen, E. et al. + # Integrated analysis of multimodal single-cell data. + # Cell 184, 3573-3587.e29 (2021). https://doi.org/10.1016/j.cell.2021.04.048 # Stuart, T., Butler, A., Hoffman, P. et al. # Comprehensive Integration of Single-Cell Data. # Cell 177, 1888-1902.e21 (2019). https://doi.org/10.1016/j.cell.2019.05.031 - doi: 10.1016/j.cell.2019.05.031 + doi: + - 10.1016/j.cell.2021.04.048 + - 10.1016/j.cell.2019.05.031 links: repository: https://github.com/satijalab/seurat documentation: https://satijalab.org/seurat/articles/seurat5_integration diff --git a/src/methods/seurat_rpca/script.R b/src/methods/seurat_rpca/script.R index 102eb34b1..ce2fe5199 100644 --- a/src/methods/seurat_rpca/script.R +++ b/src/methods/seurat_rpca/script.R @@ -1,7 +1,5 @@ -requireNamespace("anndata", quietly = TRUE) suppressPackageStartupMessages({ - library(Matrix) - library(SeuratObject) + library(anndataR) library(Seurat) }) @@ -20,22 +18,22 @@ meta <- list( ## VIASH END cat("Reading input file\n") -adata <- anndata::read_h5ad(par[["input"]]) +adata <- read_h5ad(par[["input"]]) cat("Create Seurat object\n") -# Seurat expects genes in rows, cells in columns, as a dgCMatrix -normalized <- Matrix::t(adata$layers[["normalized"]]) -normalized <- as(as(normalized, "CsparseMatrix"), "dgCMatrix") - -seurat_obj <- Seurat::CreateSeuratObject(counts = normalized, meta.data = adata$obs) -# The benchmark's log_cp10k normalization is Seurat's LogNormalize, so assign it to +# The benchmark's log_cp10k normalization is Seurat's LogNormalize, so map it to # the "data" layer instead of calling NormalizeData(). -seurat_obj[["RNA"]]$data <- normalized -seurat_obj[["RNA"]]$counts <- NULL - +seurat_obj <- adata$as_Seurat( + x_mapping = NULL, + layers_mapping = c(data = "normalized"), + assay_metadata_mapping = FALSE, + reduction_mapping = FALSE, + graph_mapping = FALSE, + misc_mapping = FALSE +) # Use the benchmark's HVGs instead of FindVariableFeatures() so that feature # selection is the same across methods. -VariableFeatures(seurat_obj) <- rownames(adata$var)[adata$var$hvg] +VariableFeatures(seurat_obj) <- adata$var_names[adata$var$hvg] cat("Split layers by batch, scale and run PCA\n") # Seurat v5 integration workflow, see @@ -62,18 +60,18 @@ if (!is.null(par$k_score)) integrate_args$k.score <- par$k_score seurat_obj <- do.call(IntegrateLayers, integrate_args) cat("Store outputs\n") -output <- anndata::AnnData( +output <- AnnData( + obs = adata$obs[, character(0)], + var = adata$var[, character(0)], + obsm = list( + X_emb = Embeddings(seurat_obj, reduction = "integrated.rpca") + ), uns = list( dataset_id = adata$uns[["dataset_id"]], normalization_id = adata$uns[["normalization_id"]], method_id = meta$name - ), - obs = adata$obs, - var = adata$var, - obsm = list( - X_emb = Embeddings(seurat_obj, reduction = "integrated.rpca") ) ) cat("Write output AnnData to file\n") -output$write_h5ad(par[["output"]], compression = "gzip") +output$write_h5ad(par[["output"]])