update Harmony code and unit test

NAMESPACE

@@ -27,83 +27,119 @@ import(Seurat)
 import(celldex)
 import(colorspace)
 import(cowplot)
-import(data.table)
 import(dplyr)
 import(gdata)
 import(ggplot2)
-import(ggpubr)
 import(ggrepel)
 import(grid)
 import(gridBase)
 import(gridExtra)
-import(harmony)
 import(httr)
 import(jsonlite)
 import(parallel)
-import(quantmod)
+import(patchwork)
 import(reshape2)
-import(rlang)
 import(scales)
 import(tidyverse)
 import(tools)
-import(utils)
 importFrom(ComplexHeatmap,pheatmap)
 importFrom(RColorBrewer,brewer.pal)
+importFrom(RColorBrewer,brewer.pal.info)
 importFrom(Seurat,AddMetaData)
+importFrom(Seurat,AddModuleScore)
+importFrom(Seurat,Assays)
 importFrom(Seurat,CreateAssayObject)
+importFrom(Seurat,DefaultAssay)
+importFrom(Seurat,DimPlot)
 importFrom(Seurat,DotPlot)
+importFrom(Seurat,FetchData)
+importFrom(Seurat,GetAssayData)
 importFrom(Seurat,Idents)
+importFrom(Seurat,Reductions)
 importFrom(Seurat,RunTSNE)
+importFrom(Seurat,RunUMAP)
+importFrom(Seurat,VariableFeatures)
 importFrom(Seurat,as.SingleCellExperiment)
+importFrom(Seurat,subset)
 importFrom(SingleR,SingleR)
 importFrom(colorspace,RGB)
 importFrom(colorspace,diverge_hcl)
 importFrom(colorspace,heat_hcl)
 importFrom(colorspace,hex)
+importFrom(cowplot,plot_grid)
+importFrom(data.table,as.factor)
 importFrom(dendextend,rotate)
 importFrom(dendsort,dendsort)
 importFrom(dplyr,across)
 importFrom(dplyr,all_of)
 importFrom(dplyr,arrange)
 importFrom(dplyr,case_when)
 importFrom(dplyr,filter)
+importFrom(dplyr,group_by)
 importFrom(dplyr,mutate)
 importFrom(dplyr,mutate_if)
 importFrom(dplyr,pull)
 importFrom(dplyr,select)
+importFrom(dplyr,summarise)
 importFrom(ggplot2,aes)
+importFrom(ggplot2,aes_string)
 importFrom(ggplot2,coord_fixed)
+importFrom(ggplot2,element_blank)
+importFrom(ggplot2,element_text)
+importFrom(ggplot2,facet_grid)
+importFrom(ggplot2,facet_wrap)
+importFrom(ggplot2,geom_boxplot)
 importFrom(ggplot2,geom_hline)
+importFrom(ggplot2,geom_jitter)
 importFrom(ggplot2,geom_point)
+importFrom(ggplot2,geom_violin)
 importFrom(ggplot2,geom_vline)
 importFrom(ggplot2,ggplot)
 importFrom(ggplot2,ggtitle)
+importFrom(ggplot2,guide_legend)
+importFrom(ggplot2,guides)
+importFrom(ggplot2,labs)
+importFrom(ggplot2,legend.text)
+importFrom(ggplot2,legend.title)
+importFrom(ggplot2,scale_color_gradient)
 importFrom(ggplot2,scale_color_identity)
+importFrom(ggplot2,scale_color_manual)
+importFrom(ggplot2,scale_y_log10)
 importFrom(ggplot2,scale_y_reverse)
 importFrom(ggplot2,theme)
 importFrom(ggplot2,theme_bw)
 importFrom(ggplot2,theme_classic)
+importFrom(ggplot2,theme_void)
 importFrom(ggplot2,xlab)
 importFrom(ggplot2,ylab)
 importFrom(ggplot2,ylim)
 importFrom(ggpubr,annotate_figure)
 importFrom(ggpubr,get_legend)
 importFrom(ggpubr,ggarrange)
+importFrom(ggpubr,text_grob)
 importFrom(grDevices,colorRampPalette)
+importFrom(grid,gpar)
 importFrom(grid,grid.draw)
+importFrom(grid,textGrob)
 importFrom(gridExtra,arrangeGrob)
+importFrom(harmony,RunHarmony)
 importFrom(htmlwidgets,saveWidget)
 importFrom(magrittr,"%>%")
 importFrom(plotly,as_widget)
 importFrom(plotly,ggplotly)
 importFrom(plotly,plot_ly)
+importFrom(purrr,compact)
+importFrom(purrr,map)
+importFrom(quantmod,paste0)
 importFrom(reshape2,melt)
+importFrom(rlang,"%||%")
 importFrom(scDblFinder,scDblFinder)
 importFrom(scales,rescale)
 importFrom(stats,as.hclust)
 importFrom(stats,hclust)
 importFrom(stats,kmeans)
 importFrom(stats,quantile)
+importFrom(stringr,str_detect)
 importFrom(stringr,str_replace_all)
 importFrom(stringr,str_sort)
 importFrom(stringr,str_to_title)
@@ -112,3 +148,5 @@ importFrom(tibble,deframe)
 importFrom(tibble,rownames_to_column)
 importFrom(tibble,tibble)
 importFrom(tidyr,pivot_wider)
+importFrom(tidyverse,arrange)
+importFrom(tidyverse,mutate)

R/Harmony.R

@@ -12,14 +12,15 @@
 #'                     expressed genes
 #' @param group.by.var Which variable should be accounted for when running 
 #'                     batch correction
-#' @param npc Number of principal components to use when running Harmony
-#'            (Default: 20)
-
-#' @import Seurat 
-#' @import harmony
-#' @import gridExtra
-#' @import RColorBrewer
-#' @import ggplot2
+#' @param npc Number of principal components to use when generating tSNE from
+#'            harmony embeddings (Default: 20)
+#'    
+#' @import patchwork    
+#' @importFrom harmony RunHarmony
+#' @importFrom Seurat VariableFeatures RunUMAP RunTSNE CreateAssayObject
+#' @importFrom ggplot2 ggplot theme_bw theme geom_point scale_color_manual 
+#'  guides guide_legend element_blank element_text
+#' @importFrom RColorBrewer brewer.pal brewer.pal.info
 #'   
 #' @export
 #' @example Do not run: harmonyBatchCorrect(object = seurat,
@@ -94,17 +95,21 @@ harmonyBatchCorrect <- function(object,
   object <- RunHarmony(object, group.by.var,
                    do_pca=FALSE,
                    assay.use = "SCT",
-                   plot_convergence = TRUE,
                    return_object=TRUE)
 
   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 <- data.frame(as.vector(object@reductions$tsne@cell.embeddings[,1]),
+  sdat.tsne <- 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) <- c("TSNE1","TSNE2","ident")
+  names(sdat.tsne) <- c("TSNE1","TSNE2","ident")
+
+  sdat.umap <- data.frame(as.vector(object@reductions$umap@cell.embeddings[,1]),
+                          as.vector(object@reductions$umap@cell.embeddings[,2]),
+                          object@meta.data[eval(parse(text = "group.by.var"))])
+  names(sdat.umap) <- c("UMAP1","UMAP2","ident")
 
   n <- 60
   qual.col.pals = brewer.pal.info[brewer.pal.info$category == 'qual',]
@@ -113,7 +118,21 @@ harmonyBatchCorrect <- function(object,
                        rownames(qual.col.pals)))
 
   set.seed(10)
-  g1 <- ggplot(sdat, aes(x=TSNE1, y=TSNE2)) +
+  g1 <- ggplot(sdat.tsne, aes(x=TSNE1, y=TSNE2)) +
+    theme_bw() +
+    theme(legend.title=element_blank()) +
+    geom_point(aes(colour=ident),size=1) +
+    scale_color_manual(values=cols) +
+    theme(panel.grid.major = element_blank(), 
+          panel.grid.minor = element_blank(), 
+          legend.position="top",
+          panel.background = element_blank(),
+          legend.text=element_text(size=rel(0.5))) +
+    guides(colour = guide_legend(override.aes = list(size=5, alpha = 1))
+    ) 
+
+  set.seed(11)
+  g2 <- ggplot(sdat.umap, aes(x=UMAP1, y=UMAP2)) +
     theme_bw() +
     theme(legend.title=element_blank()) +
     geom_point(aes(colour=ident),size=1) +
@@ -127,15 +146,13 @@ harmonyBatchCorrect <- function(object,
     ) 
 
   # Calculate adjusted gene expression from embeddings
-  seur.SCT <- t(object@assays$SCT@scale.data)
   harm.embeds <- object@reductions$harmony@cell.embeddings
   harm.lvl.backcalc <- harm.embeds %*% t(ppldngs)
 
-  # Insert back-calculated data into seurat
-  harmSCT <- CreateAssayObject(data = t(harm.lvl.backcalc))
-  object[["harmSCT"]] <- harmSCT
+  # Replace SCT scale.data expression with backcalculated data
+  object$SCT@scale.data <- t(harm.lvl.backcalc)
 
-  harmony.res <- list(adj.object = object,
-                      adj.tsne = g1)
+  harm_res <- list(harm.object = object, harm.figures = g1 + g2)
 
+  return(harm_res)
 }