Merge branch 'imputation_qmd' of https://github.com/bcbio/singlecell-reports into imputation_qmd

Author: lpantano

06_imputation/correlation_workflow.qmd

@@ -60,8 +60,6 @@ library(pheatmap)
 library(gridExtra)
 library(RColorBrewer)
 library(viridis)
-
-library(CSCORE)
 library(SAVER)
 library(Rmagic)
 
@@ -108,7 +106,7 @@ min_perc <- 0.2
 
 This report intends to both calculate correlations of gene of interests and also compares how imputation changes those values. 
 
-Throughout this report, we assess two main ways of calculating gene correlations scores (**`Spearman` and `CS-CORE`**). 
+Throughout this report, we assess Spearman correlation estimates.
 
 We further assessed two ways of imputating gene expression values: MAGIC and SAVER.
 
@@ -218,36 +216,45 @@ We compare three alternative methods of estimating expression levels to log norm
 # Store output so we don't have to re-run imputation each time
 filename <- glue("{path_outs}/imputed.RDS")
 if (!file.exists(filename)) {
-  # Get raw counts
-  raw_rna <- LayerData(seurat, assay = "RNA", layer = "counts")
+  # Get raw counts (genes x cells matrix)
+  raw_rna <- as.matrix(GetAssayData(seurat[["RNA"]], slot = "counts"))
 
   # SCT
   # Re-run SCT on subset data
   seurat <- SCTransform(seurat, return.only.var.genes = FALSE, min_cells = 1)
 
   # Creating new seurat object for genes of interest only
-  data_raw <- FetchData(seurat, assay = "RNA", layer = "counts", vars = corr_genes)
-  data_rna <- FetchData(seurat, assay = "RNA", layer = "data", vars = corr_genes)
-  data_sct <- FetchData(seurat, assay = "SCT", layer = "data", vars = corr_genes)
+  # Use GetAssayData to reliably extract counts/data (features x cells)
+  data_raw <- as.matrix(GetAssayData(seurat[["RNA"]], slot = "counts")[corr_genes, , drop = FALSE])
+  data_rna <- as.matrix(GetAssayData(seurat[["RNA"]], slot = "data")[corr_genes, , drop = FALSE])
+  data_sct <- as.matrix(GetAssayData(seurat[["SCT"]], slot = "data")[corr_genes, , drop = FALSE])
 
   seurat_imputed <- CreateSeuratObject(
-    counts = t(data_raw),
-    data = t(data_rna),
+    counts = data_raw,
+    data = data_rna,
     meta.data = seurat@meta.data
   )
-  seurat_imputed[["SCT"]] <- CreateAssayObject(data = t(data_sct))
+  seurat_imputed[["SCT"]] <- CreateAssayObject(data = data_sct)
   seurat_imputed[["RAW"]] <- CreateAssayObject(counts = raw_rna)
 
   # Delete the original seurat object to save memory
   rm(seurat)
 
   data_magic <- magic(t(raw_rna), genes = corr_genes)$result
-  seurat_imputed[["MAGIC"]] <- CreateAssayObject(data = t(data_magic))
+  # magic may return a data.frame/tibble or matrix (cells x genes). Force to matrix and transpose to features x cells
+  if (!is.matrix(data_magic)) data_magic <- as.matrix(data_magic)
+  data_magic_t <- t(data_magic)
+  # ensure dimnames align: rows = genes, cols = cells
+  if (is.null(rownames(data_magic_t))) rownames(data_magic_t) <- corr_genes
+  if (is.null(colnames(data_magic_t))) colnames(data_magic_t) <- colnames(raw_rna)
+  seurat_imputed[["MAGIC"]] <- CreateAssayObject(data = data_magic_t)
 
   # SAVER
   # Generate SAVER predictions for those genes
   genes.ind <- which(rownames(raw_rna) %in% corr_genes)
-  data_saver <- saver(raw_rna,
+  saverD <- raw_rna
+  attr(saverD, "class") <- "matrix"
+  data_saver <- saver(saverD,
     pred.genes = genes.ind,
     pred.genes.only = TRUE,
     estimates.only = TRUE,
@@ -293,8 +300,7 @@ We have a few different ways to compute correlation scores with their associated
   - `SCTransform` counts -> spearman correlation matrix
   - `MAGIC` imputed -> spearman correlation matrix
   - `SAVER` imputed -> spearman correlation matrix
-2. `CS-CORE`
-    - Raw RNA counts -> co-expression matrix
+2. (removed) `CS-CORE` (this report no longer runs CS-CORE)
 
 ```{r correlations}
 # Store output so we don't have to re-run correlation each time
@@ -321,8 +327,15 @@ if (!file.exists(filename)) {
     gene_2 <- genes_comb[idx, 2]
 
     for (assay_ in assays) {
-      gene_exp <- t(seurat_imputed[[assay_]]$data[c(gene_1, gene_2), ]) %>%
-        as.data.frame()
+      # extract assay data safely (features x cells) and subset to the two genes
+      assay_mat <- tryCatch(as.matrix(GetAssayData(seurat_imputed[[assay_]], slot = "data")), error = function(e) NULL)
+      if (is.null(assay_mat)) {
+        gene_exp <- data.frame()
+      } else {
+        sub_mat <- assay_mat[c(gene_1, gene_2), , drop = FALSE]
+        # transpose to cells x genes for cor.test
+        gene_exp <- as.data.frame(t(sub_mat))
+      }
 
       if (all(gene_exp[[gene_1]] == 0) | all(gene_exp[[gene_2]] == 0)) {
         corr_val <- 0.0
@@ -342,15 +355,7 @@ if (!file.exists(filename)) {
     }
   }
 
-  # Run CS-CORE
-  DefaultAssay(seurat_imputed) <- "RAW"
-  CSCORE_result <- CSCORE(seurat_imputed, genes = corr_genes)
-
-  # Store CS-CORE results
-  tmp <- reshape2::melt(as.matrix(CSCORE_result$est)) %>% rename(CSCORE = value)
-  df_corr <- left_join(df_corr, tmp)
-  tmp <- reshape2::melt(as.matrix(CSCORE_result$p_value)) %>% rename(CSCORE = value)
-  df_p_val <- left_join(df_p_val, tmp)
+  # CS-CORE removed: no additional co-expression estimates are appended here
 
   # Save output
   write.csv(df_corr, filename)
@@ -367,7 +372,7 @@ Showing the patterns of correlation for each method. The x-axis and y-axis are t
 
 ```{r visualize-cors}
 #| fig-width: 7
-methods <- c("RNA", "SCT", "MAGIC", "SAVER", "CSCORE")
+methods <- c("RNA", "SCT", "MAGIC", "SAVER")
 
 cor_List <- purrr::map(methods, \(method){
   corr <- df_corr[c("Var1", "Var2", method)]
@@ -381,9 +386,6 @@ cor_List <- purrr::map(methods, \(method){
 
   breaks <- seq(-1, 1, by = 0.1)
   show_legend <- F
-  if (method == "CSCORE") {
-    show_legend <- T
-  }
   p <- pheatmap(mtx,
     color = inferno(10),
     show_rownames = FALSE,
@@ -403,14 +405,14 @@ plot(cor_comb)
 
 # Compare correlation estimates across methods 
 
-Comparing the correlation scores for each gene pair for MAGIC, SAVER, and CS-CORE.
+Comparing the correlation scores for each gene pair for MAGIC and SAVER.
 
 In these scatterplots, the gene-pairs that are colored red have different results for significance.
 
 ```{r cor-compare}
 #| fig-width: 3
 #| fig-height: 8
-methods <- c("MAGIC", "SAVER", "CSCORE")
+methods <- c("MAGIC", "SAVER")
 methods_comb <- data.frame(t(combn(methods, 2)))
 plot_list <- list()
 
@@ -447,8 +449,6 @@ grid.arrange(grobs = plot_list, nrow = 3)
 
 # Method description 
 
-[`CS-CORE`](https://github.com/ChangSuBiostats/CS-CORE) is a R package for cell-type-specific co-expression inference from single cell RNA-sequencing data. It provides an implementation for the statistical method CS-CORE proposed in this [paper](https://www.nature.com/articles/s41467-023-40503-7). 
-
 [`MAGIC`](https://github.com/KrishnaswamyLab/MAGIC) is an algorithm for denoising high-dimensional data most commonly applied to single-cell RNA sequencing data. MAGIC learns the manifold data, using the resultant graph to smooth the features and restore the structure of the data.
 
 [`SAVER`](https://github.com/mohuangx/SAVER) implements a regularized regression prediction and empirical Bayes method to recover the true gene expression profile in noisy and sparse single-cell RNA-seq data.
@@ -460,3 +460,4 @@ List and version of tools used for the QC report generation.
 ```{r}
 sessionInfo()
 ```
+