Update Visualization Templates for use with Harmony Assay

R/Color_by_Gene.R

@@ -112,6 +112,7 @@ colorByGene <- function(object,
           p1 <- DimPlot(object.sub,
                         reduction = "tsne",
                         group.by = "ident")
+          colnames(p1$data) <- gsub("tsne_","tSNE_",colnames(p1$data))
           clus.mat = data.frame(
             umap1 = p1$data$tSNE_1,
             umap2 = p1$data$tSNE_2,
@@ -122,6 +123,7 @@ colorByGene <- function(object,
           p1 <- DimPlot(object.sub,
                         reduction = "umap",
                         group.by = "ident")
+          colnames(p1$data) <- gsub("umap_","UMAP_",colnames(p1$data))
           clus.mat = data.frame(
             umap1 = p1$data$UMAP_1,
             umap2 = p1$data$UMAP_2,
@@ -131,6 +133,7 @@ colorByGene <- function(object,
           p1 <- DimPlot(object.sub,
                         reduction = "pca",
                         group.by = "ident")
+          colnames(p1$data) <- gsub("pc_","PC_",colnames(p1$data))
           clus.mat = data.frame(
             umap1 = p1$data$PC_1,
             umap2 = p1$data$PC_2,

R/Color_by_Genes_Automatic.R

@@ -73,6 +73,7 @@ colorByMarkerTable <- function(object,
       if (!(cite.seq)) {
         if (reduction.type == "tsne") {
           p1 <- DimPlot(object.sub, reduction = "tsne", group.by = "ident")
+          colnames(p1$data) <- gsub("tsne_","tSNE_",colnames(p1$data))
           clusmat = data.frame(
             umap1 = p1$data$tSNE_1,
             umap2 = p1$data$tSNE_2,
@@ -82,6 +83,7 @@ colorByMarkerTable <- function(object,
         }
         else if (reduction.type == "umap") {
           p1 <- DimPlot(object.sub, reduction = "umap", group.by = "ident")
+          colnames(p1$data) <- gsub("umap_","UMAP_",colnames(p1$data))
           clusmat = data.frame(
             umap1 = p1$data$UMAP_1,
             umap2 = p1$data$UMAP_2,
@@ -91,6 +93,7 @@ colorByMarkerTable <- function(object,
         }
         else{
           p1 <- DimPlot(object.sub, reduction = "pca", group.by = "ident")
+          colnames(p1$data) <- gsub("pc_","PC_",colnames(p1$data))
           clusmat = data.frame(
             umap1 = p1$data$PC_1,
             umap2 = p1$data$PC_2,

R/Dotplot_by_Metadata.R

@@ -29,6 +29,7 @@ dotPlotMet <- function(object,
                        metadata,
                        cells,
                        markers,
+                       use_assay = "SCT",
                        plot.reverse = FALSE,
                        cell.reverse.sort = FALSE,
                        dot.color = "darkblue") {
@@ -113,7 +114,7 @@ dotPlotMet <- function(object,
 
   #Run Seurat Dotplot function
   dp <- DotPlot(object,
-                assay = "SCT",
+                assay = use_assay,
                 features = markers,
                 dot.scale = 4,
                 cols = c("lightgrey", dot.color)

R/Dual_Labeling.R

@@ -89,6 +89,8 @@ dualLabeling <- function (object,
                           display.unscaled.values = FALSE) 
 {
 
+  # Feb28 2025 ver.
+
   #### Error Messages ####
 
   #Errors for genes not available in dataset/slot

R/Harmony.R

@@ -31,11 +31,20 @@
 #' @return A list: adj.object with harmony-adjusted gene expression (SCT slot) 
 #'                 adj.tsne: harmonized tSNE plot
 
-harmonyBatchCorrect <- function(object, 
-                                nvar = 2000, 
+object = readRDS('tests/testthat/fixtures/BRCA/BRCA_Combine_and_Renormalize_SO_downsample.rds')
+
+harmonyBatchCorrect <- function(object,
+                                nvar = 200,
                                 genes.to.add = c(),
                                 group.by.var,
-                                npc = 20) {
+                                return_lognorm = T,
+                                npc = 30) {
+
+  library(patchwork)  
+  library(harmony)
+  library(Seurat)
+  library(ggplot2)
+  library(RColorBrewer)
 
   # Error and Warning Messages
   if(is.null(genes.to.add)){
@@ -59,6 +68,7 @@ harmonyBatchCorrect <- function(object,
   sdat.tsne.orig <- data.frame(as.vector(object@reductions$tsne@cell.embeddings[,1]),
                                as.vector(object@reductions$tsne@cell.embeddings[,2]),
                                object@meta.data[eval(parse(text = "group.by.var"))])
+
   names(sdat.tsne.orig) <- c("TSNE1","TSNE2","ident")
 
   sdat.umap.orig <- data.frame(as.vector(object@reductions$umap@cell.embeddings[,1]),
@@ -130,18 +140,40 @@ harmonyBatchCorrect <- function(object,
   object@reductions$pca@feature.loadings <- ppldngs
   object@reductions$pca@stdev <- pppca$d
 
+  # Store original log-normalized data and scaling parameters for back-calculation
+  if (return_lognorm) {
+    library(Matrix)
+    # Get log-normalized data for the variable features
+    lognorm_data <- object@assays$SCT@data[mvf, , drop = FALSE]
+    print(str(object))
+    print("hello")
+    print(class(lognorm_data))
+    print(dim(lognorm_data))
+
+    # Calculate scaling parameters from the original scaled data
+    #scale_center <- Matrix::rowMeans(lognorm_data)
+    scale_center <- Matrix::rowMeans(as.matrix(lognorm_data))
+    scale_scale <- apply(lognorm_data, 1, sd)
+
+    # Store these for later reconstruction
+    scaling_params <- list(
+      center = scale_center,
+      scale = scale_scale,
+      genes = mvf
+    )
+  }
+
   # By default, Harmony corrects pca embeddings. 
   # Set do_pca to FALSE to use your own pca embeddings. 
   # Stores adjusted embeddings in harmony reduction slot
-
   object <- RunHarmony(object, 
                        group.by.var,
                        do_pca=FALSE,
                        assay.use = "SCT",
                        plot_convergence = FALSE)
 
-  object <- RunUMAP(object, reduction = "harmony", dims=1:npc)
-  object <- RunTSNE(object, reduction = "harmony", dims=1:npc)
+  object <- RunUMAP(object, reduction = "harmony", dims = 1:npc)
+  object <- RunTSNE(object, reduction = "harmony", dims = 1:npc)
 
   # Plot harmony embeddings annotated by variable to batch correct for
   sdat.tsne <- data.frame(as.vector(object@reductions$tsne@cell.embeddings[,1]),
@@ -181,16 +213,36 @@ harmonyBatchCorrect <- function(object,
     annotate("text", x = Inf, y = -Inf, label = "Harmonized UMAP", hjust = 1.1, vjust = -1, size = 5)
 
   print((orig.tsne + harm.tsne) + plot_layout(ncol = 2))
-
   print((orig.umap + harm.umap) + plot_layout(ncol = 2))
 
   # Calculate adjusted gene expression from embeddings
   harm.embeds <- object@reductions$harmony@cell.embeddings
-  harm.lvl.backcalc <- harm.embeds %*% t(ppldngs)
-
+  harm.lvl.backcalc.scaled <- harm.embeds %*% t(ppldngs)
+
+  # Store batch-corrected scaled data in Harmony assay
+
+  if (return_lognorm) {
+    # Fast conversion back to log-normalized space
+    # Direct vectorized operations on the transposed matrix
+    harm.lvl.backcalc.lognorm <- t(harm.lvl.backcalc.scaled) * scaling_params$scale[mvf] + scaling_params$center[mvf]
+
+    print("Batch-corrected data stored in 'Harmony' assay:")
+    print("- Log-normalized data: object@assays$Harmony@data")
+    print("- Scaled data: object@assays$Harmony@scale.data")
+  } else {
+    print("Batch-corrected scaled data stored in object@assays$Harmony@scale.data")
+  }
 
   # Insert back-calculated data into seurat
-  object[["Harmony"]] <- CreateAssayObject(data = t(harm.lvl.backcalc))
+  object[["Harmony"]] <- CreateAssayObject(data = harm.lvl.backcalc.lognorm)
+  object@assays$Harmony@scale.data <- t(harm.lvl.backcalc.scaled)
+
+  object <- ScaleData(object, assay = "Harmony", verbose = FALSE)
+
+  # re-run PCA on harmony embeddings using top variable genes (mvf)
+  object <- RunPCA(object, assay = "Harmony", verbose = FALSE, features = rownames(object))
+
+  object <- FindNeighbors(object, reduction = "harmony", dims = 1:10, assay = "Harmony")
 
   return(object)
 }

R/Heatmap.R

@@ -51,6 +51,7 @@ heatmapSC <- function(object,
                       sample.names,
                       metadata,
                       transcripts,
+                      use_assay = 'SCT',
                       proteins = NULL,
                       heatmap.color = "Bu Yl Rd",
                       plot.title = "Heatmap",
@@ -309,16 +310,27 @@ heatmapSC <- function(object,
   }
 
 
-  #collect transcript expression data from SCT slot
+  #collect transcript expression data from SCT / Harmony slot
   df.mat1 = NULL
   if (length(transcripts) > 0) {
     if (length(transcripts) == 1) {
-      df.mat1 <-
-        vector(mode = "numeric",
-               length = length(object$SCT@scale.data[transcripts,]))
-      df.mat1 <- object$SCT@scale.data[transcripts,]
+      if(use_assay == 'SCT'){
+        df.mat1 <-
+          vector(mode = "numeric",
+                 length = length(object$SCT@scale.data[transcripts,]))
+        df.mat1 <- object$SCT@scale.data[transcripts,]
+      } else if (use_assay == 'Harmony'){
+        df.mat1 <-
+          vector(mode = "numeric",
+                 length = length(object$Harmony@scale.data[transcripts,]))
+        df.mat1 <- object$Harmony@scale.data[transcripts,]
+        }
     } else {
-      df.mat1 <- as.matrix(object$SCT@scale.data[transcripts,])
+      if(use_assay == 'SCT'){
+        df.mat1 <- as.matrix(object$SCT@scale.data[transcripts,])
+      } else if (use_assay == 'Harmony'){
+        df.mat1 <- as.matrix(object$Harmony@scale.data[transcripts,])
+        }
     }
   }
 

R/ModuleScore.R

@@ -60,7 +60,7 @@
 #'         distribution of cell marker gene, Seurat Object with cell 
 #'         classification metadata
 
-modScore <- function(object, marker.table, ms.threshold, 
+modScore <- function(object, marker.table, ms.threshold, use_assay = "SCT",
                      general.class, lvl.vec = c(), reduction = "tsne", 
                      nbins = 10, gradient.ft.size = 6, 
                      violin.ft.size = 6, step.size = 0.1) 
@@ -155,7 +155,8 @@ modScore <- function(object, marker.table, ms.threshold,
   # Calculate MS, make density plots
   for (celltype_name in names(marker.list)) {
     object = AddModuleScore(object, marker.list[celltype_name], 
-                            name = celltype_name, nbin = nbins, assay = "SCT")
+                            name = celltype_name, nbin = nbins, 
+                            assay = use_assay)
     m = paste0(celltype_name, "1")
     object@meta.data[[m]] <- scales::rescale(object@meta.data[[m]], 
                                              to = c(0, 1))