'update'

feiyoung · Mar 5, 2023 · 29212d4 · 29212d4
1 parent b19266c
commit 29212d4
Show file tree

Hide file tree

Showing 57 changed files with 1,685 additions and 855 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,8 +1,8 @@
 Package: PRECAST
 Type: Package
 Title: Embedding and Clustering with Alignment for Spatial Datasets
-Version: 1.4
-Date: 2023-02-11
+Version: 1.5
+Date: 2023-03-04
 Authors@R: 
     c(person(given = "Wei",
              family = "Liu",
@@ -41,6 +41,9 @@ Imports:
     stats,
     DR.SC, 
     scales,
+    ggpubr,
+    graphics,
+    colorspace,
     Rcpp (>= 1.0.5)
 LazyData: true
 URL: https://github.com/feiyoung/PRECAST

diff --git a/NAMESPACE b/NAMESPACE
@@ -20,7 +20,8 @@ export(ICM.EM,ICM.EM_structure, selectModel.SeqK_PRECAST_Object,
        getAdj_reg, getAdj_fixedNumber)
 
 export(coordinate_rotate, firstup, SpaPlot,dimPlot,
-       plot_RGB, plot_scatter, volinPlot, boxPlot, doHeatmap, featurePlot)
+       plot_RGB, plot_scatter, volinPlot, boxPlot, doHeatmap, featurePlot,
+       drawFigs, chooseColors)
 
 importFrom(DR.SC,getAdj_auto, getAdj_manual, getAdj, RunWPCA, read10XVisium)
 # export(getAdj, getAdj_manual, getAdj_auto, RunWPCA, read10XVisium)
@@ -43,12 +44,15 @@ importFrom(GiRaF, sampler.mrf)
 importFrom(MASS, mvrnorm)
 importFrom(Matrix, t, sparseMatrix)
 importFrom(mclust, Mclust, mclustBIC)
+importFrom(ggpubr, ggarrange)
+importFrom(colorspace, adjust_transparency)
+
 #export(Mclust, mclustBIC)
 
 importFrom(Seurat, CreateDimReducObject, NormalizeData, CreateSeuratObject, FindVariableFeatures,  
            "Idents<-", Idents, DefaultAssay, "DefaultAssay<-", ScaleData)
 # export(ScaleData)
-
+importFrom("graphics", "barplot")
 importFrom(grDevices, hcl, rgb, col2rgb)
 importFrom(stats,quantile,coef, dist, rnorm, median, 
            cov, rpois, lm, sd, hclust)

diff --git a/R/SetClass.R b/R/SetClass.R
@@ -541,32 +541,45 @@ selectModel.PRECASTObj <- function(obj, criteria = 'MBIC',pen_const=1, return_pa
 # colnames(Human_HK_genes)[1:2]<- colnames(Mouse_HK_genes)[2:1] <- c('Ensembl', 'Gene')
 # usethis::use_data(Mouse_HK_genes, Human_HK_genes, overwrite = T)
 
-
-get_correct_exp <- function(XList, RfList, housekeep, q_unwanted=10){
+## reference: Removing Unwanted Variation from High Dimensional Data with Negative Controls
+get_correct_exp <- function(XList, RfList,  houseKeep, covariateList=NULL, q_unwanted=10){
 
 
   if(!all(sapply(XList, is.matrix))){
     XList <- lapply(XList, as.matrix)
   }
-  MList <- pbapply::pblapply(XList, function(x) wpca(x[,housekeep], q=q_unwanted, F)$PCs)
-  M0 <- matlist2mat(MList)
-  rm(MList)
+  XList_sub <- pbapply::pblapply(XList, function(x) x[,houseKeep])
+  M0 <- wpca(matlist2mat(XList_sub), q=q_unwanted, FALSE)$PCs 
+
+
   Rf <- matlist2mat(RfList)
   rm(RfList)
-  XList <-  lapply(XList, scale, scale=FALSE)
+  if(!is.null(covariateList)){
+    covariates <- matlist2mat(covariateList)
+    covariates <- as.matrix(covariates)
+    rm(covariateList)
+    Rf <- cbind(Rf, covariates)
+    rm(covariates)
+  }
+  ### XList <-  lapply(XList, scale, scale=FALSE)
   X0 <- matlist2mat(XList)
   rm(XList)
-  lm1 <- lm(X0~ 0+ cbind(Rf, M0))
-  hK <- ncol(Rf)
-  coefmat <- coef(lm1)[-c(1:hK),]
+  nc_M0 <- ncol(M0)
+  lm1 <- lm(X0~ 0+ cbind(M0, Rf))
+  coefmat <- coef(lm1)[c(1:nc_M0),]
   rm(lm1)
   hX <- X0 - M0 %*% coefmat
   return(hX)
 }
 
-get_correct_mean_exp <- function(XList,  hVList, hW){
+get_correct_mean_exp <- function(XList,  hVList, covariateList=NULL){
+
+  ## XList <- lapply(XList, scale, scale=FALSE)
+
+  if(!all(sapply(XList, is.matrix))){
+    XList <- lapply(XList, as.matrix)
+  }
 
-  XList <- lapply(XList, scale, scale=FALSE)
   r_max <- length(XList)
   X0 <- XList[[1]]
   hV0 <- hVList[[1]]
@@ -578,11 +591,26 @@ get_correct_mean_exp <- function(XList,  hVList, hW){
     }
   }
 
-  X0 - hV0%*% base::t(hW)
+  rm(XList)
+  if(!is.null(covariateList)){
+    covariates <- matlist2mat(covariateList)
+    covariates <- as.matrix(covariates)
+    covariates <- cbind(1, covariates)
+    rm(covariateList)
+  }else{
+    covariates <- matrix(1, nrow=nrow(hV0), ncol=1)
+  }
+
+  nc_M0 <- ncol(hV0)
+  lm1 <- lm(X0~  0+ cbind(hV0, covariates))
+  coefmat <- coef(lm1)[c(1:nc_M0),]
+  rm(lm1)
+  X0 - hV0 %*% coefmat
+
 
 }
 
-IntegrateSpaData <- function(PRECASTObj, species="Human", custom_housekeep=NULL){
+IntegrateSpaData <- function(PRECASTObj, species="Human", custom_housekeep=NULL, covariates_use=NULL){
   # suppressMessages(require(Matrix))
   # suppressMessages(require(Seurat))
 
@@ -591,14 +619,33 @@ IntegrateSpaData <- function(PRECASTObj, species="Human", custom_housekeep=NULL)
   if(is.null(PRECASTObj@seulist)) 
     stop("IntegrateSpaData: Check the argument: PRECASTObj! The slot seulist in PRECASTObj is NULL!")
 
+
   if(!tolower(species) %in% c("human", "mouse", "unknown")) 
     stop("IntegrateSpaData: Check the argument: species! it must be one of 'Human', 'Mouse' and 'Unknown'!")
 
-  XList <- lapply(PRECASTObj@seulist, function(x) Matrix::t(x[["RNA"]]@data))
-  n_r <- length(XList)
-  for(r in 1:n_r){
-    colnames(XList[[r]]) <- firstup(colnames(XList[[r]]))
+  defAssay_vec <- sapply(PRECASTObj@seulist, DefaultAssay)
+  if(any(defAssay_vec!=defAssay_vec[1])) warning("IntegrateSpaData: there are different default assays in PRECASTObj@seulist that will be used to integrating!")
+  n_r <- length(defAssay_vec)
+
+  XList <- lapply(1:n_r,  function(r) Matrix::t(PRECASTObj@seulist[[r]][[defAssay_vec[r]]]@data))
+
+  if(!is.null(covariates_use)){
+    covariateList <- lapply(PRECASTObj@seulist, function(x) x@meta.data[, covariates_use])
+  }else{
+    covariateList <- NULL
+  }
+
+  if(tolower(species) =='mouse'){
+    for(r in 1:n_r){
+      colnames(XList[[r]]) <- firstup(colnames(XList[[r]]))
+    }
+  }
+  if(tolower(species) =='human'){
+    for(r in 1:n_r){
+      colnames(XList[[r]]) <- toupper(colnames(XList[[r]]))
+    }
   }
+
   barcodes_all <- lapply(XList, row.names)
   if(any(duplicated(unlist(barcodes_all)))){
 
@@ -611,11 +658,11 @@ IntegrateSpaData <- function(PRECASTObj, species="Human", custom_housekeep=NULL)
   houseKeep <- switch (lower_species,
     human = {
       # data(Human_HK_genes)
-       intersect((genelist), Mouse_HK_genes$Gene)
+       intersect(toupper(genelist), PRECAST::Human_HK_genes$Gene)
       },
     mouse={
       #data(Mouse_HK_genes)
-      intersect((genelist), Mouse_HK_genes$Gene)
+      intersect(firstup(genelist), PRECAST::Mouse_HK_genes$Gene)
     },
     unknown={
       character()
@@ -625,10 +672,10 @@ IntegrateSpaData <- function(PRECASTObj, species="Human", custom_housekeep=NULL)
   if(length(houseKeep) < 5){
     message("Using only PRECAST results to obtain the batch corrected gene expressions since species is unknown or the genelist in PRECASTObj has less than 5 overlapp with the housekeeping genes of given species.")
     message("Users can specify the custom_housekeep by themselves to use the housekeeping genes based methods.")
-    hX <- get_correct_mean_exp(XList,PRECASTObj@resList$hV, PRECASTObj@resList$hW)
+    hX <- get_correct_mean_exp(XList,PRECASTObj@resList$hV,  covariateList=covariateList)
   }else{
     message("Using bouth housekeeping gene and PRECAST results to obtain the batch corrected gene expressions.")
-    hX <- get_correct_exp(XList, PRECASTObj@resList$Rf, houseKeep, q_unwanted=min(10, length(houseKeep)))
+    hX <- get_correct_exp(XList, PRECASTObj@resList$Rf, houseKeep=houseKeep, q_unwanted=min(10, length(houseKeep)), covariateList=covariateList)
   }
   meta_data <- data.frame(batch=factor(get_sampleID(XList)), cluster= factor(unlist(PRECASTObj@resList$cluster)))
   row.names(meta_data) <- row.names(hX)
@@ -657,7 +704,7 @@ gg_color_hue <- function(n) {
 
 SpaPlot <- function(seuInt, batch=NULL, item=NULL, point_size=2,text_size=12, 
                     cols=NULL,font_family='', border_col="gray10",
-                    fill_col='white', ncol=2, combine = TRUE, title_name="Sample"){
+                    fill_col='white', ncol=2, combine = TRUE, title_name="Sample", ...){
 
   ## Check arguments input
 
@@ -714,7 +761,7 @@ SpaPlot <- function(seuInt, batch=NULL, item=NULL, point_size=2,text_size=12,
 
       sort_id <- sort(unique(meta_data[, 'tmp_item_id']))
       p1 <- plot_scatter(embed_use, meta_data, label_name=item, 
-                         point_size=point_size, cols =cols[sort_id])
+                         point_size=point_size, cols =cols[sort_id], ...)
     }else if(item=="RGB_UMAP"){
       p1 <- plot_RGB(embed_use, seu@reductions$UMAP3@cell.embeddings, pointsize = point_size)
     }else if(item=="RGB_TSNE"){
@@ -740,7 +787,7 @@ SpaPlot <- function(seuInt, batch=NULL, item=NULL, point_size=2,text_size=12,
 }
 dimPlot <- function(seuInt, item=NULL, reduction=NULL, point_size=1,text_size=16, 
                     cols=NULL,font_family='', border_col="gray10",
-                    fill_col="white"){
+                    fill_col="white", ...){
 
 
   if(!inherits(seuInt, "Seurat")) stop("dimPlot: Check argument: seuInt! it must be a Seurat Object.")
@@ -779,7 +826,7 @@ dimPlot <- function(seuInt, item=NULL, reduction=NULL, point_size=1,text_size=16
 
   embed_use <- seuInt[[reduction]]@cell.embeddings[,c(1,2)]
   p1 <- plot_scatter(embed_use, meta_data, label_name=item, 
-                     point_size=point_size,cols =cols[sort_id])
+                     point_size=point_size,cols =cols[sort_id], ...)
 
   p1 <- p1 + mytheme_graybox(base_size = text_size, base_family = font_family, bg_fill  = fill_col,
                           border_color = border_col)
@@ -788,8 +835,6 @@ dimPlot <- function(seuInt, item=NULL, reduction=NULL, point_size=1,text_size=16
 }
 
 
-
-
 # seuList <- gendata_seulist()
 # PRECASTObj <-  CreatePRECASTObject(seuList)
 # PRECASTObj <-  AddAdjList(PRECASTObj)

diff --git a/R/Visualization.R b/R/Visualization.R
@@ -120,7 +120,7 @@ plot_scatter <- function (
     cols = NULL, 
      point_size = 0.5, point_alpha=1, 
     base_size = 12, do_points = TRUE, do_density = FALSE, border_col='gray',
-    legend_pos='right', legend_dir='vertical') {
+    legend_pos='right', legend_dir='vertical', nrow.legend=NULL) {
   # require(dplyr)
   # require(ggthemes)
   # require(ggrepel)
@@ -157,11 +157,20 @@ plot_scatter <- function (
           legend.text=element_text(size=base_size+1),
           legend.title=element_text(size=base_size+2),
           panel.background= element_rect(fill = 'white', color=border_col))+
-    guides(color = guide_legend(override.aes = list(stroke = 1, 
-                                                    alpha = 1, shape = 16, size = 4)), alpha = "none") + 
     scale_color_manual(values = cols) + scale_fill_manual(values = cols) + 
     theme(plot.title = element_text(hjust = 0.5)) + labs(x = xy_names[1], 
                                                          y = xy_names[2])
+
+  if (!is.null(nrow.legend)){
+    plt <- plt + guides(color = guide_legend(nrow = nrow.legend,override.aes = list(stroke = 1,
+                                                                                    alpha = 1, shape = 16, size = 4)))
+  }else{
+    plt <- plt + guides(color = guide_legend(override.aes = list(stroke = 1,
+                                                                 alpha = 1, shape = 16, size = 4)),
+                        alpha = "none")
+  }
+
+
   if (do_points) 
     plt <- plt + geom_point( size = point_size, alpha=point_alpha)
   if (do_density) 
@@ -412,3 +421,62 @@ mytheme <- function(legend.direction = "horizontal",
   return(th)
 }
 
+
+
+# Other plots -------------------------------------------------------------
+
+
+drawFigs <- function(pList, layout.dim = NULL, common.legend=FALSE,legend.position='right',  ...){
+  if(!is.list(pList)) stop('drawFigs: pList must be a list!')
+
+  if(is.null(layout.dim) && length(pList)>1){
+    layout.dim <- c(2, round(length(pList)/2) )
+  }
+  if(is.null(layout.dim) && length(pList) == 1){
+    layout.dim <- c(1,1)
+  }
+  ggpubr::ggarrange(plotlist = pList, ncol = layout.dim[2],
+                    nrow = layout.dim[1], common.legend = common.legend,
+                    legend = legend.position, ...)
+
+}
+
+chooseColors <- function(palettes_name= c("Nature 10", "Light 13", "Classic 20", "Blink 23", "Hue n"), n_colors = 7,
+                         alpha=1, plot_colors=FALSE){
+
+  # require(colorspace)
+  palettes_name <- match.arg(palettes_name)
+  colors <- if(palettes_name == "Classic 20"){
+    # require(ggthemes)
+    # palettes <- ggthemes_data[["tableau"]][["color-palettes"]][["regular"]]
+    pal1 <- tableau_color_pal(palettes_name)
+    pal1(n_colors)
+  }else if(palettes_name == "Nature 10"){
+    cols <- c("#E04D50", "#4374A5", "#F08A21","#2AB673", "#FCDDDE",
+              "#70B5B0", "#DFE0EE" ,"#DFCDE4", "#FACB12", "#f9decf")
+    cols[1:n_colors]
+  }else if(palettes_name == "Blink 23"){
+    cols <- c("#c10023", "#008e17", "#fb8500", "#f60000", "#FE0092", "#bc9000","#4ffc00", "#00bcac", "#0099cc",
+              "#D35400", "#00eefd", "#cf6bd6", "#99cc00", "#aa00ff", "#ff00ff", "#0053c8",
+              "#f2a287","#ffb3ff", "#800000", "#77a7b7", "#00896e", "#00cc99", "#007CC8")
+    cols[1:n_colors]
+  }else if(palettes_name == "Light 13"){
+    cols <-c( "#FD7446" ,"#709AE1", "#31A354","#9EDAE5",
+              "#DE9ED6" ,"#BCBD22", "#CE6DBD" ,"#DADAEB" ,
+              "#FF9896","#91D1C2", "#C7E9C0" ,
+              "#6B6ECF", "#7B4173" )
+    cols[1:n_colors]
+  }else if(palettes_name == "Hue n"){
+    gg_color_hue(n_colors)
+  }else{
+    stop(paste0("chooseColors: check palettes_name! Unsupported palettes_name: ", palettes_name))
+  }
+  #require(colorspace)
+  colors_new = adjust_transparency(colors,   alpha = alpha)
+  if(plot_colors){
+    barplot(rep(1, length(colors_new)), axes = FALSE, space = 0, col = colors_new)
+  }
+
+  return(colors_new)
+}
+
diff --git a/R/main.R b/R/main.R
@@ -1,8 +1,10 @@
+# library(pkgdown)
 # pkgdown::build_site()
 # pkgdown::build_reference()
-# build_article(name="PRECAST.DLPFC") # Solely compile one article for updating.
+# build_home()
+# build_article(name="PRECAST") # Solely compile one article for updating.
 # build_article(name="PRECAST.BreastCancer")
-# R CMD check --as-cran PRECAST_1.3.tar.gz
+# R CMD check --as-cran PRECAST_1.5.tar.gz
 # devtools::check_win_release()
 # iDR.SC <- function(...) UseMethod("iDR.SC")
 model_set <- function(Sigma_equal=FALSE, Sigma_diag=TRUE,mix_prop_heter=TRUE,

diff --git a/README.md b/README.md
@@ -67,11 +67,11 @@ For usage examples and guided walkthroughs, check the `vignettes` directory of t
 
 For the users that don't have set up system properly, the following setup on different systems can be referred.
 ## Setup on Windows system
-First, download [Rtools](https://cran.r-project.org/bin/windows/Rtools/); second, add the Rtools directory to the environment variable. Users can follow [here](https://helpdeskgeek.com/windows-10/add-windows-path-environment-variable/#:~:text=Go%20ahead%20and%20click%20on%20the%20Environment%20Variables,you%20have%20to%20decide%20which%20one%20to%20edit) to add Windows PATH Environment Variable.
+First, download [Rtools](https://cran.r-project.org/bin/windows/Rtools/); second, add the Rtools directory to the environment variable.
 
 
 ## Setup on MacOS system
-First, install Xcode. Installation about Xcode can be referred [here](https://stackoverflow.com/questions/8291146/xcode-installation-on-mac#:~:text=You%20get%20it%20from%20the%20Mac%20App%20Store.,find%20the%20app%2C%20and%20click%20the%20install%20button).
+First, install Xcode. Installation about Xcode can be referred [here](https://stackoverflow.com/questions/8291146/xcode-installation-on-mac).
 
 
 Second, install "gfortran" for compiling C++ and Fortran at [here](https://github.com/fxcoudert/gfortran-for-macOS).
@@ -108,11 +108,19 @@ install.packages("DR.SC")
 For an example of typical PRECAST usage, please see our [Package Website](https://feiyoung.github.io/PRECAST/index.html) for a demonstration and overview of the functions included in PRECAST.
 
 # NEWs
+
+PRECAST version 1.5 (2023-03-05)
+
+* Fix the [issue](https://github.com/feiyoung/PRECAST/issues/2) reported by anvaly. Specifically, the assay name "RNA" used in functions `IntegrateSpaData()`  is replaced  by  the default assay using `DefaultAssay` function in Seurat. Fix the typo `human = {intersect((genelist),Mouse_HK_genes$Gene)}` with replacement of `human = {intersect(toupper(genelist), PRECAST::Human_HK_genes$Gene)}`。
+
+
+* Add two functions for visualization: `chooseColors()` and `drawFigs()`.
+
 PRECAST version 1.3 (2022-10-05)
 
-* Fix the [issue](https://github.com/feiyoung/PRECAST_Analysis/issues/1) reported by Boyi Guo. Specifically, the assay name "RNA" used in functions `CreatePRECASTObject` and  `PRECAST` is replaced  by  the default assay using `DefaultAssay` function in Seurat.
+* Fix the [issue](https://github.com/feiyoung/PRECAST_Analysis/issues/1) reported by Boyi Guo. Specifically, the assay name "RNA" used in functions `CreatePRECASTObject()` and  `PRECAST()` is replaced  by  the default assay using `DefaultAssay` function in Seurat.
 
-* Provide more detailed help file for `CreatePRECASTObject` function. Users can use `?CreatePRECASTObject` in Rstudio to access the help file.
+* Provide more detailed help file for `CreatePRECASTObject()` function. Users can use `?CreatePRECASTObject` in Rstudio to access the help file.
 In detail, seuList is a list  with Seurat object as component, and each Seurat object at least includes the raw expression count matrix, and spatial coordinates in metadata for each data batch, where the spatial coordinates information must be saved in the metadata of Seurat, named "row" and "col" for each data batch. See the help file for more details.
 
 * Add the [data](https://github.com/feiyoung/PRECAST/tree/main/vignettes_data) used in [Package Website](https://feiyoung.github.io/PRECAST/index.html).