Delete duplicate layer in AnnData objects

bin/add_celltypes_to_sce.R

@@ -191,13 +191,13 @@ if (!is.null(opt$cellassign_predictions)) {
   if (file.size(opt$cellassign_predictions) > 0) {
     predictions <- readr::read_tsv(opt$cellassign_predictions)
   } else {
-    # if it's empty, then sce could not be converted to anndata and cell assign was not run
-    sce$cellassign_celltype_annotation <- "Not run"
+    predictions <- NULL
   }
 
-  # if the only column is the barcode column then CellAssign didn't complete successfully
+  # if the only column is the barcode column or if the predictions file was empty
+  # then CellAssign didn't complete successfully
   # otherwise add in cell type annotations and metadata to SCE
-  if (all(colnames(predictions) == "barcode")) {
+  if (is.null(predictions) || all(colnames(predictions) == "barcode")) {
     # if failed then note that in the cell type column
     sce$cellassign_celltype_annotation <- "Not run"
   } else {

bin/merge_sces.R

@@ -247,7 +247,7 @@ adt_present_columns <- sce_list |>
 
 # ensure that there are indeed no "adt" altExps if adt_present_columns is empty
 adt_altexps <- sce_list |>
-  purrr::map(\(sce) "adt" %in% altExpNames(sce))
+  purrr::map_lgl(\(sce) "adt" %in% altExpNames(sce))
 if (is.null(adt_present_columns) && sum(adt_altexps) > 0) {
   stop("Error in determining which adt altExp columns should be retained.")
 }

bin/move_counts_anndata.py

@@ -49,6 +49,7 @@
     # move logcounts to X and rename
     object.X = object.layers["logcounts"]
     object.uns["X_name"] = "logcounts"
+    del object.layers["logcounts"]
 
     # export object
     object.write_h5ad(args.anndata_file, compression="gzip" if args.compress else None)

bin/sce_qc_report.R

@@ -285,17 +285,20 @@ if (opt$celltype_report_file != "") {
     stop("Supplemental cell types report template not found.")
   }
 
-  # render report
-  rmarkdown::render(
-    input = opt$celltype_report_template,
-    output_file = basename(opt$celltype_report_file),
-    output_dir = dirname(opt$celltype_report_file),
-    intermediates_dir = tempdir(),
-    knit_root_dir = tempdir(),
-    envir = new.env(),
-    params = list(
-      library = metadata_list$library_id,
-      processed_sce = processed_sce
+  # only render supplemental report if there's more than one cell
+  if (ncol(processed_sce) > 1) {
+    # render report
+    rmarkdown::render(
+      input = opt$celltype_report_template,
+      output_file = basename(opt$celltype_report_file),
+      output_dir = dirname(opt$celltype_report_file),
+      intermediates_dir = tempdir(),
+      knit_root_dir = tempdir(),
+      envir = new.env(),
+      params = list(
+        library = metadata_list$library_id,
+        processed_sce = processed_sce
+      )
     )
-  )
+  }
 }

merge.nf

@@ -178,15 +178,26 @@ workflow {
         log.warn("Not merging ${it.project_id} because it contains multiplexed libraries.")
       }
 
+    // print out warning message for any libraries not included in merging
+    filtered_libraries_ch.single_sample
+      .map{[
+        it.library_id,
+        file("${params.results_dir}/${it.project_id}/${it.sample_id}/${it.library_id}_processed.rds")
+      ]}
+    .filter{!(it[1].exists() && it[1].size() > 0)}
+    .subscribe{
+      log.warn("Processed files do not exist for ${it[0]}. This library will not be included in the merged object.")
+    }
+
     grouped_libraries_ch = filtered_libraries_ch.single_sample
       // create tuple of [project id, library_id, processed_sce_file]
       .map{[
         it.project_id,
         it.library_id,
         file("${params.results_dir}/${it.project_id}/${it.sample_id}/${it.library_id}_processed.rds")
       ]}
-      // only include libraries that have been processed through scpca-nf
-      .filter{file(it[2]).exists()}
+      // only include libraries that have been processed through scpca-nf and aren't empty
+      .filter{it[2].exists() && it[2].size() > 0}
       // only one row per library ID, this removes all the duplicates that may be present due to CITE/hashing
       .unique()
       // group tuple by project id: [project_id, [library_id1, library_id2, ...], [sce_file1, sce_file2, ...]]

templates/qc_report/celltypes_qc.rmd

@@ -1,5 +1,10 @@
 # Cell type Annotation Summary
 
+<!--
+This file is meant to be run as a child report within either `main_qc_report.rmd` or `celltypes_supplemental_report.rmd`. 
+-->
+
+
 ```{r}
 ## function definitions ##
 
@@ -86,14 +91,14 @@ lump_wrap_celltypes <- function(df, n_celltypes = 7, wrap = 35) {
 #' @param color_variable Column in data frame to color by, not a string.
 #' @param legend_title Title for legend.
 #' @param legend_nrow Number of rows in legend. Default is 2.
-#' @param point_size Point size
+#' @param point_size Point size. Default is 1
 #'
 #' @return UMAP plot as a ggplot2 object
 plot_umap <- function(
     umap_df,
     color_variable,
     legend_title,
-    point_size = umap_point_size,
+    point_size = 1,
     legend_nrow = 2) {
   ggplot(umap_df) +
     aes(
@@ -131,14 +136,14 @@ plot_umap <- function(
 #' @param umap_df Data frame with UMAP1 and UMAP2 columns
 #' @param n_celltypes The number of cell types (facets) displayed in the plot
 #' @param annotation_column Column containing cell type annotations
-#' @param point_size Point size
+#' @param point_size Point size. Default is 1
 #'
 #' @return ggplot object containing a faceted UMAP where each cell type is a facet.
 #'   In each panel, the cell type of interest is colored and all other cells are grey.
 faceted_umap <- function(umap_df,
                          n_celltypes,
                          annotation_column,
-                         point_size = umap_facet_point_size) {
+                         point_size = 1) {
   # Determine legend y-coordinate based on n_celltypes
   if (n_celltypes %in% 7:8) {
     legend_y <- 0.33
@@ -283,6 +288,11 @@ glue::glue("
 ```{r, warning = FALSE}
 # Create data frame of cell types
 celltype_df <- create_celltype_df(processed_sce)
+
+# determine UMAP point sizing
+umap_points_sizes <- determine_umap_point_size(ncol(processed_sce))
+umap_point_size <- umap_points_sizes[1]
+umap_facet_point_size <- umap_points_sizes[2]
 ```
 
 
@@ -370,7 +380,7 @@ create_celltype_n_table(celltype_df, cellassign_celltype_annotation) |>
 ```
 
 
-```{r, eval = has_umap}
+```{r, eval = has_umap && has_clusters}
 knitr::asis_output(glue::glue("
 ## UMAPs
 
@@ -388,7 +398,7 @@ umap_df <- lump_wrap_celltypes(celltype_df)
 
 <!-- First UMAP: clusters -->
 
-```{r, eval = has_umap && has_multiplex, results='asis'}
+```{r, eval = has_umap && has_multiplex && has_clusters, results='asis'}
 glue::glue("
   <div class=\"alert alert-info\">
     This library contains multiple samples that have not been batch-corrected, which may confound clustering assignments.
@@ -398,11 +408,12 @@ glue::glue("
 ```
 
 
-```{r eval = has_umap, message=FALSE, warning=FALSE}
+```{r eval = has_umap && has_clusters, message=FALSE, warning=FALSE}
 clusters_plot <- plot_umap(
   umap_df,
   cluster,
-  "Cluster"
+  "Cluster",
+  point_size = umap_point_size
 ) +
   ggtitle("UMAP colored by cluster identity")
 
@@ -418,7 +429,7 @@ if (length(levels(umap_df$cluster)) <= 8) {
 ```
 
 
-```{r, eval = has_umap}
+```{r, eval = has_umap && has_celltypes}
 knitr::asis_output(
   'Next, we show UMAPs colored by cell types.
 For each cell typing method, we show a separate faceted UMAP.
@@ -447,7 +458,8 @@ if (has_submitter & has_umap) {
 faceted_umap(
   umap_df,
   submitter_n_celltypes,
-  submitter_celltype_annotation_lumped
+  submitter_celltype_annotation_lumped,
+  point_size = umap_facet_point_size
 ) +
   ggtitle("UMAP colored by submitter-provided annotations")
 ```
@@ -469,7 +481,8 @@ if (has_singler & has_umap) {
 faceted_umap(
   umap_df,
   singler_n_celltypes,
-  singler_celltype_annotation_lumped
+  singler_celltype_annotation_lumped,
+  point_size = umap_facet_point_size
 ) +
   ggtitle("UMAP colored by SingleR annotations")
 ```
@@ -490,7 +503,8 @@ if (has_cellassign & has_umap) {
 faceted_umap(
   umap_df,
   cellassign_n_celltypes,
-  cellassign_celltype_annotation_lumped
+  cellassign_celltype_annotation_lumped,
+  point_size = umap_facet_point_size
 ) +
   ggtitle("UMAP colored by CellAssign annotations")
 ```

templates/qc_report/celltypes_supplemental_report.rmd

@@ -264,8 +264,12 @@ has_cellassign <- "cellassign" %in% metadata(processed_sce)$celltype_methods
 has_submitter <- "submitter" %in% metadata(processed_sce)$celltype_methods &&
   !all(is.na(processed_sce$submitter_celltype_annotation)) # make sure they aren't all NA
 
-# check for umap
+# If at least 1 is present, we have cell type annotations.
+has_celltypes <- any(has_singler, has_cellassign, has_submitter)
+
+# check for umap and clusters
 has_umap <- "UMAP" %in% reducedDimNames(processed_sce)
+has_clusters <- "cluster" %in% names(colData(processed_sce))
 
 # what celltypes are available?
 available_celltypes <- c(
@@ -296,11 +300,6 @@ plot_height <- 1
 #  sample_id should be defined with length > 1
 sample_id <- metadata(processed_sce)$sample_id
 has_multiplex <- length(sample_id) > 1
-
-# determine UMAP point sizing
-umap_points_sizes <- determine_umap_point_size(ncol(processed_sce))
-umap_point_size <- umap_points_sizes[1]
-umap_facet_point_size <- umap_points_sizes[2]
 ```
 
 <!-- If multiplexed, open with warning  --> 
@@ -416,7 +415,7 @@ glue::glue("
 ```
 
 <!-- If not multiplexed, show the header, text, and heatmap --> 
-```{r, eval = !has_multiplex, results='asis'}
+```{r, eval = !has_multiplex && has_clusters, results='asis'}
 glue::glue("
   ## Unsupervised clustering
 
@@ -448,7 +447,7 @@ plot_height <- calculate_plot_height(
 ```
 
 
-```{r, eval = !has_multiplex, fig.height = plot_height, fig.width = 8.5, warning = FALSE}
+```{r, eval = !has_multiplex && has_clusters, fig.height = plot_height, fig.width = 8.5, warning = FALSE}
 jaccard_cluster_matrices |>
   create_heatmap_list(
     column_title = "Clusters",

templates/qc_report/main_qc_report.rmd

@@ -117,7 +117,7 @@ if (has_cellhash) {
 # check for umap and celltypes, but need to be sure that processed_sce exists first
 if (has_processed) {
   has_umap <- "UMAP" %in% reducedDimNames(processed_sce)
-
+  has_clusters <- "cluster" %in% names(colData(processed_sce))
   has_singler <- "singler" %in% metadata(processed_sce)$celltype_methods
   has_cellassign <- "cellassign" %in% metadata(processed_sce)$celltype_methods
   has_submitter <- "submitter" %in% metadata(processed_sce)$celltype_methods &&
@@ -129,6 +129,7 @@ if (has_processed) {
   is_supplemental <- FALSE # this is not the celltype supp report
 } else {
   has_umap <- FALSE
+  has_clusters <- FALSE
   has_singler <- FALSE
   has_cellassign <- FALSE
   has_submitter <- FALSE
@@ -143,11 +144,6 @@ if ((has_singler | has_cellassign) & is.null(params$celltype_report)) {
 # check if we have multiplex
 has_multiplex <- length(sample_id) > 1
 sample_types <- metadata(unfiltered_sce)$sample_type
-
-# determine UMAP point sizing
-umap_points_sizes <- determine_umap_point_size(ncol(processed_sce))
-umap_point_size <- umap_points_sizes[1]
-umap_facet_point_size <- umap_points_sizes[2]
 ```
 
 

templates/qc_report/umap_qc.rmd

@@ -2,6 +2,14 @@
 
 The below plot shows the UMAP (Uniform Manifold Approximation and Projection) embeddings for each cell, coloring each cell by the total number of genes detected per cell.
 
+```{r}
+# determine UMAP point sizing
+umap_points_sizes <- determine_umap_point_size(ncol(processed_sce))
+umap_point_size <- umap_points_sizes[1]
+umap_facet_point_size <- umap_points_sizes[2]
+```
+
+
 ```{r message=FALSE}
 # create UMAP colored by number of genes detected
 scater::plotUMAP(