Add Seurat components

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
 

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]

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"]])

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]

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"]])