Add cycombine to batch

CHANGELOG.md

@@ -53,6 +53,7 @@
   * Added CytoNorm with aggregate of samples as controls (`methods/cytonorm_no_controls`).
   * Added parameters to tune CytoNorm.
 
+* Added cyCombine correction to a reference batch (PR #90).
 * Added `metrics/bras`
 
 ## MAJOR CHANGES

src/methods/cycombine_all_controls/config.vsh.yaml → src/methods/cycombine_all_controls_to_goal/config.vsh.yaml

@@ -7,19 +7,19 @@ __merge__: ../../api/comp_method.yaml
 
 # A unique identifier for your component (required).
 # Can contain only lowercase letters or underscores.
-name: cycombine_all_controls
+name: cycombine_all_controls_to_goal
 # A relatively short label, used when rendering visualisations (required)
-label: cyCombine with all controls
+label: cyCombine (all-controls, to-goal)
 # A one sentence summary of how this method works (required). Used when 
 # rendering summary tables.
-summary: "cyCombine with all control samples"
+summary: "cyCombine run with all control samples, correcting to a goal batch"
 description: |
   cyCombine perform batch integration using self-organizing maps and ComBat.
   It first uses self-organizing maps (SOM) to group similar cells into clusters, 
   then applies a ComBat-based method to correct batch effects within 
   each cluster. 
-  Here, we run cyCombine using all control samples (replicates in cyCombine terminology)
-  and square SOM grid.
+  Here, we run cyCombine using all control samples (replicates, in cyCombine terms), 
+  with a square SOM grid and correct the batches to batch 1.
   The size of the SOM grid is varied linearly between value of 6-16, with default set to 8
   following the default value in create_som function in cyCombine.
   Rlen, which is the number of times data is presented to the SOM network and can impact

src/methods/cycombine_all_controls_to_goal/script.R

@@ -0,0 +1,136 @@
+requireNamespace("cyCombine", quietly = TRUE)
+requireNamespace("anndata", 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",
+    som_grid_size = 8,
+    rlen = 10
+)
+meta <- list(name = "cycombine_all_controls_to_goal")
+## VIASH END
+
+cat("Reading input files\n")
+input_adata <- anndata::read_h5ad(par[["input"]])
+
+cat("Preparing input Anndata and df\n")
+
+markers_to_correct <- input_adata$var_names[input_adata$var$to_correct]
+
+adata_to_correct <- input_adata[, markers_to_correct]
+
+# convert adata to data.frame
+
+df_to_correct <- as.data.frame(
+    adata_to_correct$layers[["preprocessed"]]
+)
+df_to_correct$batch <- as.factor(adata_to_correct$obs$batch)
+df_to_correct$sample <- as.factor(adata_to_correct$obs$sample)
+
+# add an "anchor" column which specify which samples are the technical replicates
+# this is a bit weird as the anchor information should be a column name or character vector
+# giving all anchors the same label and every other sample a unique label.
+# so for our controls (replicate samples), i have to give it the same name,
+# while non-replicate sample should be given a unique identifier.
+# hence for non-replicate samples, we will just use the values in the "sample" column.
+# for the replicate samples, we will use the values in the "is_control" column.
+# that way, samples from the same donor will have the same label, and
+# cycombine can identify which samples are technical replicates of each other.
+# the as.character function is needed as otherwise we will get NAs for those controls..
+
+df_to_correct$anchor <- ifelse(
+    adata_to_correct$obs$is_control == 0,
+    as.character(adata_to_correct$obs$sample),
+    paste0("control_", adata_to_correct$obs$is_control)
+)
+df_to_correct$anchor <- as.factor(df_to_correct$anchor)
+
+lineage_markers <- as.vector(input_adata$var_names[
+    input_adata$var$marker_type == "lineage"
+])
+
+
+cat("Run cyCombine\n")
+
+# use the default parameters in normalize
+# do the normalisation on all markers to be corrected
+# use z-score normalisation because we are merging batches
+# from a single study.
+df_to_correct_norm <- cyCombine::normalize(
+    df = df_to_correct,
+    markers = markers_to_correct,
+    norm_method = "scale",
+    ties.method = "average"
+)
+
+# again, using default parameter values
+cluster_labels <- cyCombine::create_som(
+    df = df_to_correct_norm,
+    markers = lineage_markers,
+    rlen = par[["rlen"]],
+    seed = 42,
+    xdim = par[["som_grid_size"]],
+    ydim = par[["som_grid_size"]]
+)
+
+# Batch correct using default parameter values
+df_corrected <- cyCombine::correct_data(
+    df = df_to_correct,
+    label = cluster_labels,
+    markers = markers_to_correct,
+    method = "ComBat",
+    covar = NULL,
+    anchor = "anchor",
+    ref.batch = "1",
+    parametric = TRUE
+)
+
+cat("Preparing output Anndata\n")
+df_not_corrected <- as.data.frame(
+    input_adata[, !input_adata$var$to_correct]$layers[["preprocessed"]]
+)
+
+df_output <- cbind(
+    df_corrected[, markers_to_correct],
+    df_not_corrected
+)
+
+# reorder column to match input_adata
+df_output <- df_output[, input_adata$var_names]
+
+output <- anndata::AnnData(
+    obs = input_adata$obs[, integer(0)],
+    var = input_adata$var[colnames(df_output), integer(0)],
+    layers = list(integrated = as.matrix(df_output)),
+    uns = list(
+        dataset_id = input_adata$uns$dataset_id,
+        method_id = meta$name,
+        parameters = list(
+            "normalize" = list(
+                "markers" = markers_to_correct,
+                "norm_method" = "scale",
+                "ties.method" = "average"
+            ),
+            "create_som" = list(
+                "markers" = lineage_markers,
+                "seed" = 42,
+                "rlen" = par[["rlen"]],
+                "xdim" = par[["som_grid_size"]],
+                "ydim" = par[["som_grid_size"]]
+            ),
+            "correct_data" = list(
+                "markers" = markers_to_correct,
+                "method" = "ComBat",
+                "covar" = NULL,
+                "anchor" = "anchor",
+                "ref.batch" = "1",
+                "parametric" = TRUE
+            )
+        )
+    )
+)
+
+cat("Write output AnnData to file\n")
+
+output$write_h5ad(par[["output"]], compression = "gzip")

src/methods/cycombine_one_control/config.vsh.yaml → src/methods/cycombine_all_controls_to_mid/config.vsh.yaml

@@ -7,19 +7,19 @@ __merge__: ../../api/comp_method.yaml
 
 # A unique identifier for your component (required).
 # Can contain only lowercase letters or underscores.
-name: cycombine_one_control
+name: cycombine_all_controls_to_mid
 # A relatively short label, used when rendering visualisations (required)
-label: cyCombine with one control
+label: cyCombine (all-controls, to-middle)
 # A one sentence summary of how this method works (required). Used when 
 # rendering summary tables.
-summary: "cyCombine with one control samples"
+summary: "cyCombine run with all control samples, correcting to a midpoint"
 description: |
   cyCombine perform batch integration using self-organizing maps and ComBat.
   It first uses self-organizing maps (SOM) to group similar cells into clusters, 
   then applies a ComBat-based method to correct batch effects within 
   each cluster. 
-  Here, we run cyCombine with control samples from just one donor
-  (replicates in cyCombine terminology) using square SOM grid.
+  Here, we run cyCombine using all control samples (replicates, in cyCombine terms), 
+  with a square SOM grid and correct the batches to a midpoint derived from all batches.
   The size of the SOM grid is varied linearly between value of 6-16, with default set to 8
   following the default value in create_som function in cyCombine.
   Rlen, which is the number of times data is presented to the SOM network and can impact
@@ -53,6 +53,7 @@ argument_groups:
             upper: 20
         default: 10
         description: The number of time data is presented to SOM when clustering
+
 
 resources:
   # The script of your component (required)

src/methods/cycombine_all_controls/script.R → src/methods/cycombine_all_controls_to_mid/script.R

@@ -23,8 +23,8 @@ adata_to_correct <- input_adata[, markers_to_correct]
 df_to_correct <- as.data.frame(
     adata_to_correct$layers[["preprocessed"]]
 )
-df_to_correct$batch <- adata_to_correct$obs$batch
-df_to_correct$sample <- adata_to_correct$obs$sample
+df_to_correct$batch <- as.factor(adata_to_correct$obs$batch)
+df_to_correct$sample <- as.factor(adata_to_correct$obs$sample)
 
 # add an "anchor" column which specify which samples are the technical replicates
 # this is a bit weird as the anchor information should be a column name or character vector
@@ -42,6 +42,7 @@ df_to_correct$anchor <- ifelse(
     as.character(adata_to_correct$obs$sample),
     paste0("control_", adata_to_correct$obs$is_control)
 )
+df_to_correct$anchor <- as.factor(df_to_correct$anchor)
 
 lineage_markers <- as.vector(input_adata$var_names[
     input_adata$var$marker_type == "lineage"

src/methods/cycombine_no_controls/config.vsh.yaml → src/methods/cycombine_no_controls_to_goal/config.vsh.yaml

@@ -1,15 +1,15 @@
 __merge__: ../../api/comp_method.yaml
 
-name: cycombine_no_controls
-label: cyCombine (no-controls)
-summary: "cyCombine without control samples"
+name: cycombine_no_controls_to_goal
+label: cyCombine (no-controls, to-goal)
+summary: "cyCombine run without control samples, correcting to a goal batch"
 description: |
   cyCombine perform batch integration using self-organizing maps and ComBat.
   It first uses self-organizing maps (SOM) to group similar cells into clusters, 
   then applies a ComBat-based method to correct batch effects within 
   each cluster. 
-  Here, we run cyCombine without any control samples (replicates in cyCombine terminology) 
-  using square SOM grid.
+  Here, we run cyCombine without any control samples (replicates, in cyCombine terms), 
+  with a square SOM grid and correct the batches to batch 1.
   The size of the SOM grid is varied linearly between value of 6-16, with default set to 8
   following the default value in create_som function in cyCombine.
   Rlen, which is the number of times data is presented to the SOM network and can impact

src/methods/cycombine_no_controls_to_goal/script.R

@@ -0,0 +1,113 @@
+requireNamespace("anndata", quietly = TRUE)
+requireNamespace("cyCombine", 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",
+    som_grid_size = 8,
+    rlen = 10
+)
+meta <- list(name = "cycombine_no_controls")
+## VIASH END
+
+cat("Reading input files\n")
+input_adata <- anndata::read_h5ad(par[["input"]])
+
+cat("Preparing input Anndata and df\n")
+
+adata_to_correct <- input_adata[, input_adata$var$to_correct]
+
+markers_to_correct <- input_adata$var_names[input_adata$var$to_correct]
+
+lineage_markers <- as.vector(input_adata$var_names[
+    input_adata$var$marker_type == "lineage"
+])
+
+df_to_correct <- as.data.frame(
+    adata_to_correct$layers[["preprocessed"]]
+)
+df_to_correct$batch <- as.factor(adata_to_correct$obs$batch)
+df_to_correct$sample <- as.factor(adata_to_correct$obs$sample)
+
+cat("Run cyCombine\n")
+
+# use the default parameters in normalize
+# do the normalisation on all markers to be corrected
+df_to_correct_norm <- cyCombine::normalize(
+    df = df_to_correct,
+    markers = markers_to_correct,
+    norm_method = "scale",
+    ties.method = "average"
+)
+
+# again, using default parameter values
+cluster_labels <- cyCombine::create_som(
+    df = df_to_correct_norm,
+    markers = lineage_markers,
+    rlen = par[["rlen"]],
+    seed = 42,
+    xdim = par[["som_grid_size"]],
+    ydim = par[["som_grid_size"]]
+)
+
+# Batch correct using default parameter values
+df_corrected <- cyCombine::correct_data(
+    df = df_to_correct,
+    label = cluster_labels,
+    markers = markers_to_correct,
+    method = "ComBat",
+    covar = NULL,
+    anchor = NULL,
+    ref.batch = "1",
+    parametric = TRUE
+)
+
+cat("Preparing output Anndata\n")
+df_not_corrected <- as.data.frame(
+    input_adata[, !input_adata$var$to_correct]$layers[["preprocessed"]]
+)
+
+df_output <- cbind(
+    df_corrected[, markers_to_correct],
+    df_not_corrected
+)
+
+# reorder column to match input_adata
+df_output <- df_output[, input_adata$var_names]
+
+output <- anndata::AnnData(
+    obs = input_adata$obs[, integer(0)],
+    var = input_adata$var[colnames(df_output), integer(0)],
+    layers = list(integrated = as.matrix(df_output)),
+    uns = list(
+        dataset_id = input_adata$uns$dataset_id,
+        method_id = meta$name,
+        parameters = list(
+            "normalize" = list(
+                "markers" = markers_to_correct,
+                "norm_method" = "scale",
+                "ties.method" = "average"
+            ),
+            "create_som" = list(
+                "markers" = lineage_markers,
+                "seed" = 42,
+                "rlen" = par[["rlen"]],
+                "xdim" = par[["som_grid_size"]],
+                "ydim" = par[["som_grid_size"]]
+            ),
+            "correct_data" = list(
+                "markers" = markers_to_correct,
+                "method" = "ComBat",
+                "covar" = NULL,
+                "anchor" = NULL,
+                "ref.batch" = "1",
+                "parametric" = TRUE
+            )
+        )
+    )
+)
+
+cat("Write output AnnData to file\n")
+
+output$write_h5ad(par[["output"]], compression = "gzip")