Add rpca

CHANGELOG.md

@@ -56,7 +56,9 @@
 
 * Added CytoNorm correction to a goal batch (PR #92).
 * Added cyCombine correction to a reference batch (PR #90).
-* Added `metrics/bras`
+* Added `metrics/bras` (PR #91).
+
+* Added Seurat rPCA (PR #95).
 
 
 ## MAJOR CHANGES

src/methods/rpca_to_goal/config.vsh.yaml

@@ -0,0 +1,104 @@
+# The API specifies which type of component this is.
+# It contains specifications for:
+#   - The input/output files
+#   - Common parameters
+#   - A unit test
+__merge__: ../../api/comp_method.yaml
+
+# A unique identifier for your component (required).
+# Can contain only lowercase letters or underscores.
+name: rpca_to_goal
+# A relatively short label, used when rendering visualisations (required)
+label: Seurat RPCA (to-goal)
+
+# A one sentence summary of how this method works (required). Used when 
+# rendering summary tables.
+summary: "Batch integrate data to a goal batch using mutual nearest neighbors identified via Seurat reciprocal PCA."
+
+# A multi-line description of how this component works (required). Used
+# when rendering reference documentation.
+description: |
+  Seurat RPCA performs batch integration by projecting each query dataset into the PCA space 
+  of a goal batch, and identifying anchors using reciprocal PCA (RPCA). 
+  RPCA identifies mutual nearest neighbors between the goal batch and the remaining batches 
+  in their shared low-dimensional space (PCA), which are used to align 
+  and integrate the batches into the goal batch. 
+
+  We ran Seurat RPCA implemented in Seurat v4.4.0, available on their github.
+  This is because subsequent version (>= v5) does not support getting corrected
+  count matrix, only corrected PC space.
+  See: https://github.com/satijalab/seurat/issues/8551.
+
+  We varied the number of PCs and nearest neighbours considered when running RPCA.
+
+references:
+  doi: 
+    - 10.1016/j.cell.2021.04.048
+  bibtex:
+    - |
+      @article{hao2021integrated,
+        title={Integrated analysis of multimodal single-cell data},
+        author={Hao, Yuhan and Hao, Stephanie and Andersen-Nissen, Erica and Mauck, William M and Zheng, Shiwei and Butler, Andrew and Lee, Madison J and Wilk, Aaron J and Darby, Charlotte and Zager, Michael and others},
+        journal={Cell},
+        volume={184},
+        number={13},
+        pages={3573--3587},
+        year={2021},
+        publisher={Elsevier}
+      }
+
+links:
+  # URL to the documentation for this method (required).
+  documentation: https://satijalab.org/seurat/articles/integration_rpca.html
+  # URL to the code repository for this method (required).
+  repository: https://github.com/satijalab/seurat
+
+argument_groups:
+  - name: Parameters
+    arguments:
+      - type: integer
+        name: --npcs
+        info:
+          optimize:
+            type: linear
+            lower: 10
+            upper: 20
+        default: 10
+        description: The number of principal component (PC) dimensions to use when computing anchors between batches.
+      - type: integer
+        name: --n_neighbours
+        info:
+          optimize:
+            type: linear
+            lower: 5
+            upper: 50
+        default: 5
+        description: The number of mutual nearest neighbors that are used when identifying anchors between batches
+
+# Resources required to run the component
+resources:
+  # The script of your component (required)
+  - type: r_script
+    path: script.R
+
+engines:
+  # Specifications for the Docker image for this component.
+  - type: docker
+    image: openproblems/base_r:1
+    # Add custom dependencies here (optional). For more information, see
+    # https://viash.io/reference/config/engines/docker/#setup .
+    setup:
+      - type: r
+        packages: [ anndata ]
+        url: [
+          https://cran.r-project.org/src/contrib/Archive/Seurat/Seurat_4.4.0.tar.gz,
+          https://cran.r-project.org/src/contrib/Archive/SeuratObject/SeuratObject_4.1.4.tar.gz
+        ]
+
+runners:
+  # This platform allows running the component natively
+  - type: executable
+  # Allows turning the component into a Nextflow module / pipeline.
+  - type: nextflow
+    directives:
+      label: [midtime,midmem,midcpu]

src/methods/rpca_to_goal/script.R

@@ -0,0 +1,147 @@
+requireNamespace("anndata", quietly = TRUE)
+requireNamespace("Seurat", quietly = TRUE)
+
+## VIASH START
+par <- list(
+    input = "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/unintegrated_censored.h5ad",
+    output = "resources_test/task_cyto_batch_integration/mouse_spleen_flow_cytometry_subset/output.h5ad",
+    npcs = 10,
+    n_neighbours = 50
+)
+meta <- list(
+    name = "rpca_to_goal"
+)
+## VIASH END
+
+cat("Reading input files\n")
+input_adata <- anndata::read_h5ad(par[["input"]])
+
+cat("Preparing input Anndata\n")
+input_adata$obs$batch <- as.factor(input_adata$obs$batch)
+
+adata_to_correct <- input_adata[, input_adata$var$to_correct]
+markers_to_correct <- input_adata$var_names[input_adata$var$to_correct]
+
+cat("Creating Seurat object and preprocess\n")
+
+# create one seurat object per batch
+batches <- unique(input_adata$obs$batch)
+
+seurat_objs <- lapply(batches, function(batch) {
+
+    cat(paste("Processing batch", batch))
+
+    adata_batch <- input_adata[
+        input_adata$obs$batch == batch,
+        input_adata$var$to_correct
+    ]
+    # batch <- batches[1]
+    mat <- Matrix::as.matrix(adata_batch$layers["preprocessed"])
+
+    # have to transpose so cells are columns..
+    mat <- Matrix::t(mat)
+
+    # convert to sparse matrix
+    mat <- Matrix::Matrix(mat, sparse = TRUE)
+
+    seurat_obj <- Seurat::CreateSeuratObject(
+        counts = mat,
+        data = mat,
+        assay = "cyto",
+        meta.data = adata_batch$obs
+    )
+
+    # save RAM
+    rm(mat)
+
+    # scale all features/markers
+    seurat_obj <- Seurat::ScaleData(
+        object = seurat_obj,
+        features = markers_to_correct,
+        assay = "cyto",
+        verbose = FALSE
+    )
+
+    # run pca. mandatory
+    # if num pcs is more than number of markers, it'll be capped at
+    # the number of markers
+    # not using approximate pca as we don't have many markers
+    seurat_obj <- Seurat::RunPCA(
+        object = seurat_obj,
+        features = markers_to_correct,
+        assay = "cyto",
+        npcs = par[["npcs"]],
+        approx = FALSE,
+        verbose = FALSE
+    )
+
+    return(seurat_obj)
+})
+
+names(seurat_objs) <- batches
+
+cat("Finding anchors\n")
+
+# get how many PCs we have calculated
+# if the number of PCs is more than how many markers
+# seurat set that to how many markers.
+# hence we can't just use par[["npcs"]] below.
+npcs_computed <- dim(seurat_objs[[1]][["pca"]])[2]
+
+anchors <- Seurat::FindIntegrationAnchors(
+    object.list = seurat_objs,
+    anchor.features = markers_to_correct,
+    dims = seq(npcs_computed),
+    k.anchor = par[["n_neighbours"]],
+    reduction = "rpca",
+    verbose = FALSE,
+    reference = which(names(seurat_objs) == "1")
+)
+
+cat("Batch correct\n")
+
+# Warning will say Layer counts isn't present in the assay object; returning NULL
+# Even though the original assay has counts layer.
+# Not sure why. But the object has data layer.
+batch_corrected_seurat_obj <- Seurat::IntegrateData(
+    anchorset = anchors,
+    features = markers_to_correct,
+    features.to.integrate = markers_to_correct,
+    dims = seq(npcs_computed),
+    verbose = FALSE
+)
+# just to be sure!
+Seurat::DefaultAssay(batch_corrected_seurat_obj) <- "integrated"
+
+cat("Creating output AnnData\n")
+
+batch_corrected_mat <- Matrix::t(
+    Matrix::as.matrix(Seurat::GetAssayData(batch_corrected_seurat_obj))
+)
+# cbind corrected matrix to matrix containing markers not corrected
+batch_corrected_mat <- cbind(
+    batch_corrected_mat,
+    input_adata[, !input_adata$var$to_correct]$layers[["preprocessed"]]
+)
+
+# make sure the row and column orders are matching
+# between input adata and the batch corrected matrix
+batch_corrected_mat <- batch_corrected_mat[
+    input_adata$obs_names, input_adata$var_names
+]
+
+cat("Write output AnnData to file\n")
+output <- anndata::AnnData(
+    obs = input_adata$obs[, integer(0)],
+    var = input_adata$var[colnames(batch_corrected_mat), integer(0)],
+    layers = list(integrated = batch_corrected_mat),
+    uns = list(
+        dataset_id = input_adata$uns$dataset_id,
+        method_id = meta$name,
+        parameters = list(
+            "npcs" = par[["npcs"]],
+            "n_neighbours" = par[["n_neighbours"]]
+        )
+    )
+)
+output$write_h5ad(par[["output"]], compression = "gzip")

src/methods/rpca_to_mid/config.vsh.yaml

@@ -0,0 +1,106 @@
+# The API specifies which type of component this is.
+# It contains specifications for:
+#   - The input/output files
+#   - Common parameters
+#   - A unit test
+__merge__: ../../api/comp_method.yaml
+
+# A unique identifier for your component (required).
+# Can contain only lowercase letters or underscores.
+name: rpca_to_mid
+# A relatively short label, used when rendering visualisations (required)
+label: Seurat RPCA (to-middle)
+
+# A one sentence summary of how this method works (required). Used when 
+# rendering summary tables.
+summary: "Batch integrate data to a midpoint using mutual nearest neighbors identified via Seurat reciprocal PCA."
+
+# A multi-line description of how this component works (required). Used
+# when rendering reference documentation.
+description: |
+  Seurat RPCA performs batch integration by identifying mutual nearest neighbors (anchors) 
+  between all batches using reciprocal PCA (RPCA). 
+  Instead of merging batches into a single PCA space, each batch is projected into the 
+  PCA space of the others, and anchors are found where cells are mutual nearest neighbors 
+  across these projections. 
+  These anchors are then used to compute batch correction vectors, 
+  aligning the batches into a shared corrected space. 
+
+  We ran Seurat RPCA implemented in Seurat v4.4.0, available on their github.
+  This is because subsequent version (>= v5) does not support getting corrected
+  count matrix, only corrected PC space.
+  See: https://github.com/satijalab/seurat/issues/8551
+
+  We varied the number of PCs and nearest neighbours considered when running RPCA.
+
+references:
+  doi: 
+    - 10.1016/j.cell.2021.04.048
+  bibtex:
+    - |
+      @article{hao2021integrated,
+        title={Integrated analysis of multimodal single-cell data},
+        author={Hao, Yuhan and Hao, Stephanie and Andersen-Nissen, Erica and Mauck, William M and Zheng, Shiwei and Butler, Andrew and Lee, Madison J and Wilk, Aaron J and Darby, Charlotte and Zager, Michael and others},
+        journal={Cell},
+        volume={184},
+        number={13},
+        pages={3573--3587},
+        year={2021},
+        publisher={Elsevier}
+      }
+
+links:
+  # URL to the documentation for this method (required).
+  documentation: https://satijalab.org/seurat/articles/integration_rpca.html
+  # URL to the code repository for this method (required).
+  repository: https://github.com/satijalab/seurat
+
+argument_groups:
+  - name: Parameters
+    arguments:
+      - type: integer
+        name: --npcs
+        info:
+          optimize:
+            type: linear
+            lower: 10
+            upper: 20
+        default: 10
+        description: The number of principal component (PC) dimensions to use when computing anchors between batches.
+      - type: integer
+        name: --n_neighbours
+        info:
+          optimize:
+            type: linear
+            lower: 5
+            upper: 50
+        default: 5
+        description: The number of mutual nearest neighbors that are used when identifying anchors between batches
+
+# Resources required to run the component
+resources:
+  # The script of your component (required)
+  - type: r_script
+    path: script.R
+
+engines:
+  # Specifications for the Docker image for this component.
+  - type: docker
+    image: openproblems/base_r:1
+    # Add custom dependencies here (optional). For more information, see
+    # https://viash.io/reference/config/engines/docker/#setup .
+    setup:
+      - type: r
+        packages: [ anndata ]
+        url: [
+          https://cran.r-project.org/src/contrib/Archive/Seurat/Seurat_4.4.0.tar.gz,
+          https://cran.r-project.org/src/contrib/Archive/SeuratObject/SeuratObject_4.1.4.tar.gz
+        ]
+
+runners:
+  # This platform allows running the component natively
+  - type: executable
+  # Allows turning the component into a Nextflow module / pipeline.
+  - type: nextflow
+    directives:
+      label: [midtime,midmem,midcpu]