diff --git a/CHANGELOG.md b/CHANGELOG.md index 3054e439..9355c063 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -28,6 +28,8 @@ We propose adding this metric to substitute iLISI. * Added `method/limma_removebatcheffect` component (PR #79). +* Added `methods/seurat_cca` and `methods/seurat_rpca` components (PR #77). + - Seurat v5 anchor-based integration (`IntegrateLayers` with `CCAIntegration` / `RPCAIntegration`). ## Minor changes diff --git a/src/methods/seurat_cca/config.vsh.yaml b/src/methods/seurat_cca/config.vsh.yaml new file mode 100644 index 00000000..4c9ec268 --- /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: Seurat v5 anchor-based integration using canonical correlation analysis +description: | + 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. + # 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 +info: + method_types: [embedding] + preferred_normalization: log_cp10k +arguments: + - name: --dims + type: integer + 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. Defaults to Seurat's default (5). + example: 5 + - name: --k_filter + type: integer + 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. Defaults to Seurat's default (30). + example: 30 +resources: + - type: r_script + path: script.R +engines: + - type: docker + image: openproblems/base_r:1 + setup: + - type: r + cran: + - Seurat +runners: + - type: executable + - type: nextflow + directives: + label: [midcpu, highmem, hightime] diff --git a/src/methods/seurat_cca/script.R b/src/methods/seurat_cca/script.R new file mode 100644 index 00000000..41901025 --- /dev/null +++ b/src/methods/seurat_cca/script.R @@ -0,0 +1,77 @@ +suppressPackageStartupMessages({ + library(anndataR) + library(Seurat) +}) + +## VIASH START +par <- list( + 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("Reading input file\n") +adata <- read_h5ad(par[["input"]]) + +cat("Create Seurat object\n") +# The benchmark's log_cp10k normalization is Seurat's LogNormalize, so map it to +# the "data" layer instead of calling NormalizeData(). +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) <- adata$var_names[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 +) +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( + 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 + ) +) + +cat("Write output AnnData to file\n") +output$write_h5ad(par[["output"]]) diff --git a/src/methods/seurat_rpca/config.vsh.yaml b/src/methods/seurat_rpca/config.vsh.yaml new file mode 100644 index 00000000..73ae06ac --- /dev/null +++ b/src/methods/seurat_rpca/config.vsh.yaml @@ -0,0 +1,64 @@ +__merge__: /src/api/comp_method.yaml +name: seurat_rpca +label: Seurat RPCA +summary: Seurat v5 anchor-based integration using reciprocal PCA +description: | + 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: + # 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.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 +info: + method_types: [embedding] + preferred_normalization: log_cp10k +arguments: + - name: --dims + type: integer + 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. Defaults to Seurat's default (5). + example: 5 + - name: --k_filter + type: integer + 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. Defaults to Seurat's default (30). + example: 30 +resources: + - type: r_script + path: script.R +engines: + - type: docker + image: openproblems/base_r:1 + setup: + - type: r + cran: + - Seurat +runners: + - type: executable + - type: nextflow + directives: + label: [midcpu, midmem, midtime] diff --git a/src/methods/seurat_rpca/script.R b/src/methods/seurat_rpca/script.R new file mode 100644 index 00000000..ce2fe519 --- /dev/null +++ b/src/methods/seurat_rpca/script.R @@ -0,0 +1,77 @@ +suppressPackageStartupMessages({ + library(anndataR) + library(Seurat) +}) + +## VIASH START +par <- list( + 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("Reading input file\n") +adata <- read_h5ad(par[["input"]]) + +cat("Create Seurat object\n") +# The benchmark's log_cp10k normalization is Seurat's LogNormalize, so map it to +# the "data" layer instead of calling NormalizeData(). +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) <- adata$var_names[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 +) +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( + 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 + ) +) + +cat("Write output AnnData to file\n") +output$write_h5ad(par[["output"]])