diff --git a/NAMESPACE b/NAMESPACE
index 6f15d028..8bee5a94 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -138,7 +138,6 @@ export(download)
 export(drop_data)
 export(exist_Python_pkgs)
 export(fastMNN_integrate)
-export(filter)
 export(find_conda)
 export(geom_alluvial)
 export(geom_alluvial_label)
@@ -148,7 +147,6 @@ export(geom_sankey_bump)
 export(geom_sankey_label)
 export(geom_sankey_text)
 export(get_vars)
-export(group_by)
 export(installed_Python_pkgs)
 export(invoke)
 export(isOutlier)
@@ -156,7 +154,7 @@ export(iterchunks)
 export(make_long)
 export(palette_scp)
 export(panel_fix)
-export(panel_fix_single)
+export(panel_fix_overall)
 export(scVI_integrate)
 export(segementsDf)
 export(show_palettes)
@@ -168,7 +166,6 @@ export(theme_sankey)
 export(theme_sankey_bump)
 export(theme_scp)
 export(tochunks)
-export(top_n)
 export(try_get)
 export(unnest)
 importFrom(AnnotationDbi,columns)
@@ -330,7 +327,6 @@ importFrom(dplyr,reframe)
 importFrom(dplyr,slice_head)
 importFrom(dplyr,summarise)
 importFrom(dplyr,summarise_at)
-importFrom(dplyr,top_n)
 importFrom(ggforce,geom_mark_ellipse)
 importFrom(ggforce,geom_mark_hull)
 importFrom(ggnewscale,new_scale)
diff --git a/R/SCP-analysis.R b/R/SCP-analysis.R
index 4e5c2731..7a9238c1 100644
--- a/R/SCP-analysis.R
+++ b/R/SCP-analysis.R
@@ -1244,6 +1244,7 @@ FindConservedMarkers2 <- function(object, grouping.var, ident.1, ident.2 = NULL,
 #'
 #' @seealso \code{\link{RunEnrichment}} \code{\link{RunGSEA}} \code{\link{GroupHeatmap}}
 #' @examples
+#' library(dplyr)
 #' data("pancreas_sub")
 #' pancreas_sub <- RunDEtest(pancreas_sub, group_by = "SubCellType")
 #' AllMarkers <- filter(pancreas_sub@tools$DEtest_SubCellType$AllMarkers_wilcox, p_val_adj < 0.05 & avg_log2FC > 1)
@@ -4051,7 +4052,7 @@ RunMonocle2 <- function(srt, assay = NULL, slot = "counts", expressionFamily = "
                         root_state = NULL, seed = 11) {
   set.seed(seed)
   check_R(c("monocle", "DDRTree", "BiocGenerics", "Biobase", "VGAM"))
-  library(monocle)
+  attachNamespace("DDRTree")
 
   assay <- assay %||% DefaultAssay(srt)
   expr_matrix <- as.sparse(GetAssayData(srt, assay = assay, slot = slot))
diff --git a/R/SCP-app.R b/R/SCP-app.R
index 55bbc7ca..b69f7e80 100644
--- a/R/SCP-app.R
+++ b/R/SCP-app.R
@@ -538,7 +538,8 @@ CreateSeuratObject2 <- function(counts, project = "SeuratProject", assay = "RNA"
 #'   initial_dataset = "mouse_pancreas",
 #'   initial_group = "CellType",
 #'   initial_feature = "Neurog3",
-#'   workers = 2, overwrite = TRUE
+#'   session_workers = 2,
+#'   overwrite = TRUE
 #' )
 #' list.files("./SCExplorer") # This directory can be used as site directory for Shiny Server.
 #'
@@ -1274,6 +1275,28 @@ ui <- fluidPage(
             height = "100px",
             placeholder = paste(sample(all_features, 4), collapse = "\\n")
           ),
+          fluidRow(
+            column(
+              width = 6, align = "center",
+              radioButtons(
+                inputId = "plotby4",
+                label = "Plot by",
+                choices = c("group", "feature"),
+                selected = "group",
+                inline = TRUE
+              )
+            ),
+            column(
+              width = 6, align = "center",
+              radioButtons(
+                inputId = "fillby4",
+                label = "Fill by",
+                choices = c("group", "feature", "expression"),
+                selected = "group",
+                inline = TRUE
+              )
+            )
+          ),
           fluidRow(
             column(
               width = 4, align = "center",
@@ -1338,12 +1361,6 @@ ui <- fluidPage(
               )
             )
           ),
-          selectInput(
-            inputId = "fillby4",
-            label = "Fill by",
-            choices = c("group", "feature", "expression"),
-            selected = "group"
-          ),
           selectInput(
             inputId = "palette4",
             label = "Select a palette",
@@ -2038,13 +2055,14 @@ server <- function(input, output, session) {
     plottype4 <- input$plottype4
     features4 <- input$features4
     feature_area4 <- input$feature_area4
+    plotby4 <- input$plotby4
+    fillby4 <- input$fillby4
     coExp4 <- input$coExp4 == "Yes"
     stack4 <- input$stack4 == "Yes"
     flip4 <- input$flip4 == "Yes"
     addbox4 <- input$addbox4 == "Yes"
     addpoint4 <- input$addpoint4 == "Yes"
     addtrend4 <- input$addtrend4 == "Yes"
-    fillby4 <- input$fillby4
     palette4 <- input$palette4
     theme4 <- input$theme4
     sameylims4 <- input$sameylims4 == "Yes"
@@ -2104,7 +2122,7 @@ server <- function(input, output, session) {
           srt = srt_tmp, stat.by = features4, group.by = group4, split.by = split4, cells = cells, slot = "data", plot_type = plottype4,
           calculate_coexp = coExp4, stack = stack4, flip = flip4,
           add_box = addbox4, add_point = addpoint4, add_trend = addtrend4,
-          fill.by = fillby4, palette = palette4, theme_use = theme4, same.y.lims = sameylims4,
+          plot.by = plotby4, fill.by = fillby4, palette = palette4, theme_use = theme4, same.y.lims = sameylims4,
           aspect.ratio = as.numeric(aspect.ratio), # must be class of numeric instead of integer
           ncol = ncol4, byrow = byrow4, force = TRUE
         )
diff --git a/R/SCP-plot.R b/R/SCP-plot.R
index 87738235..cac409f9 100644
--- a/R/SCP-plot.R
+++ b/R/SCP-plot.R
@@ -13,7 +13,6 @@
 #' p + theme_scp()
 #' @importFrom ggplot2 theme element_blank element_text element_rect margin
 #' @export
-#'
 theme_scp <- function(aspect.ratio = NULL, base_size = 12, ...) {
   text_size_scale <- base_size / 12
   args1 <- list(
@@ -62,6 +61,9 @@ theme_scp <- function(aspect.ratio = NULL, base_size = 12, ...) {
 
 #' Blank theme
 #'
+#' This function creates a theme with all elements blank except for axis lines and labels.
+#' It can optionally add coordinate axes in the plot.
+#'
 #' @param add_coord Whether to add coordinate arrows. Default is \code{TRUE}.
 #' @param xlen_npc The length of the x-axis arrow in "npc".
 #' @param ylen_npc The length of the y-axis arrow in "npc".
@@ -489,7 +491,7 @@ panel_fix <- function(x = NULL, panel_index = NULL, respect = NULL,
         subgrob <- gtable$grobs[[i]][["children"]][[j]][["children"]][[1]][["children"]][[1]]
         # print(subgrob$grobs[[1]][["children"]])
         if (length(subgrob$grobs[[1]][["children"]]) > 0 && all(sapply(subgrob$grobs[[1]][["children"]], function(x) inherits(x, "recordedGrob")))) {
-          subgrob <- panel_fix_single(subgrob$grobs[[1]][["children"]], width = width, height = height, margin = padding, units = units, raster = raster, dpi = dpi, return_grob = TRUE)
+          subgrob <- panel_fix_overall(subgrob$grobs[[1]][["children"]], width = width, height = height, margin = padding, units = units, raster = raster, dpi = dpi, return_grob = TRUE)
         } else {
           subgrob <- panel_fix(subgrob, width = width, height = height, margin = padding, units = units, raster = raster, dpi = dpi, return_grob = TRUE, verbose = verbose, depth = depth + 1)
         }
@@ -498,7 +500,7 @@ panel_fix <- function(x = NULL, panel_index = NULL, respect = NULL,
       }
       sum_width <- convertWidth(sum(subgrob[["widths"]]), unitTo = units, valueOnly = TRUE) / as.numeric(gtable$grobs[[i]][["children"]][[j]]$vp$width)
       sum_height <- convertHeight(sum(subgrob[["heights"]]), unitTo = units, valueOnly = TRUE) / as.numeric(gtable$grobs[[i]][["children"]][[j]]$vp$height)
-      gtable <- panel_fix_single(gtable, panel_index = i, width = sum_width, height = sum_height, margin = ifelse(depth == 1, margin, 0), units = units, raster = FALSE, return_grob = TRUE)
+      gtable <- panel_fix_overall(gtable, panel_index = i, width = sum_width, height = sum_height, margin = ifelse(depth == 1, margin, 0), units = units, raster = FALSE, return_grob = TRUE)
     } else if (gtable$grobs[[i]]$name == "layout" || inherits(x, "patchwork")) {
       if (isTRUE(verbose)) {
         message("panel ", i, " is detected as generated by patchwork.")
@@ -508,7 +510,7 @@ panel_fix <- function(x = NULL, panel_index = NULL, respect = NULL,
       # }
       subgrob <- gtable$grobs[[i]]
       if (length(subgrob[["children"]]) > 0 && all(sapply(subgrob[["children"]], function(x) inherits(x, "recordedGrob")))) {
-        subgrob <- panel_fix_single(subgrob[["children"]], width = width, height = height, margin = 0, units = units, raster = raster, dpi = dpi, return_grob = TRUE)
+        subgrob <- panel_fix_overall(subgrob[["children"]], width = width, height = height, margin = 0, units = units, raster = raster, dpi = dpi, return_grob = TRUE)
       } else {
         subgrob <- panel_fix(subgrob, width = width, height = height, margin = 0, units = units, raster = raster, dpi = dpi, return_grob = TRUE, verbose = verbose, depth = depth + 1)
       }
@@ -533,13 +535,13 @@ panel_fix <- function(x = NULL, panel_index = NULL, respect = NULL,
       sum_width <- convertWidth(sum(subgrob[["widths"]]), unitTo = units, valueOnly = TRUE)
       sum_height <- convertHeight(sum(subgrob[["heights"]]), unitTo = units, valueOnly = TRUE)
 
-      gtable <- panel_fix_single(gtable, panel_index = i, width = sum_width, height = sum_height, margin = ifelse(depth == 1 & add_margin, margin, 0), units = units, raster = FALSE, respect = TRUE, return_grob = TRUE)
+      gtable <- panel_fix_overall(gtable, panel_index = i, width = sum_width, height = sum_height, margin = ifelse(depth == 1 & add_margin, margin, 0), units = units, raster = FALSE, respect = TRUE, return_grob = TRUE)
       if (depth == 1 & add_margin) {
         add_margin <- FALSE
       }
     } else {
       # print("fix the gtable")
-      gtable <- panel_fix_single(gtable, panel_index = i, width = width, height = height, margin = margin, units = units, raster = raster, dpi = dpi, return_grob = TRUE)
+      gtable <- panel_fix_overall(gtable, panel_index = i, width = width, height = height, margin = margin, units = units, raster = raster, dpi = dpi, return_grob = TRUE)
     }
   }
 
@@ -583,10 +585,10 @@ panel_fix <- function(x = NULL, panel_index = NULL, respect = NULL,
 #' @importFrom patchwork wrap_plots
 #' @importFrom grDevices dev.off
 #' @export
-panel_fix_single <- function(x, panel_index = NULL, respect = NULL,
-                             width = NULL, height = NULL, margin = 1, units = "in",
-                             raster = FALSE, dpi = 300, BPPARAM = BiocParallel::SerialParam(),
-                             return_grob = FALSE, bg_color = "white", save = NULL, verbose = TRUE) {
+panel_fix_overall <- function(x, panel_index = NULL, respect = NULL,
+                              width = NULL, height = NULL, margin = 1, units = "in",
+                              raster = FALSE, dpi = 300, BPPARAM = BiocParallel::SerialParam(),
+                              return_grob = FALSE, bg_color = "white", save = NULL, verbose = TRUE) {
   if (!inherits(x, "gtable")) {
     if (inherits(x, "gTree")) {
       x <- x[["children"]]
@@ -3252,6 +3254,7 @@ FeatureDimPlot3D <- function(srt, features = NULL, reduction = NULL, dims = c(1,
 #' @param stat.by A character vector specifying the features to plot.
 #' @param group.by A character vector specifying the groups to group by. Default is NULL.
 #' @param split.by A character vector specifying the variable to split the plot by. Default is NULL.
+#' @param plot.by A character vector specifying how to plot the data, by group or feature. Possible values are "group", "feature". Default is "group".
 #' @param bg.by A character vector specifying the variable to use as the background color. Default is NULL.
 #' @param fill.by A string specifying what to fill the plot by. Possible values are "group", "feature", or "expression". Default is "group".
 #' @param cells A character vector specifying the cells to include in the plot. Default is NULL.
@@ -3358,15 +3361,24 @@ FeatureDimPlot3D <- function(srt, features = NULL, reduction = NULL, dims = c(1,
 #'   ),
 #'   fill.by = "feature", plot_type = "box",
 #'   group.by = "SubCellType", bg.by = "CellType", stack = TRUE, flip = TRUE
-#' ) %>% panel_fix_single(width = 8, height = 5) # Because the plot is made by combining, we want to adjust the overall height and width
-#' @importFrom Seurat DefaultAssay GetAssayData
+#' ) %>% panel_fix_overall(width = 8, height = 5) # As the plot is created by combining, we can adjust the overall height and width directly.
+#'
+#' FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "CellType", plot.by = "feature")
+#' FeatureStatPlot(pancreas_sub, stat.by = c(
+#'   "Sox9", "Anxa2", "Bicc1", # Ductal
+#'   "Neurog3", "Hes6", # EPs
+#'   "Fev", "Neurod1", # Pre-endocrine
+#'   "Rbp4", "Pyy", # Endocrine
+#'   "Ins1", "Gcg", "Sst", "Ghrl" # Beta, Alpha, Delta, Epsilon
+#' ), group.by = "SubCellType", plot.by = "feature", stack = TRUE)
+#'
+#' @importFrom Seurat FetchData
+#' @importFrom reshape2 melt
 #' @importFrom gtable gtable_add_cols gtable_add_rows gtable_add_grob gtable_add_padding
-#' @importFrom ggplot2 geom_blank geom_violin geom_rect geom_boxplot geom_count geom_col geom_vline geom_hline layer_data layer_scales position_jitterdodge position_dodge stat_summary scale_x_discrete element_line element_text element_blank annotate mean_sdl after_stat scale_shape_identity
 #' @importFrom grid grobHeight grobWidth
 #' @importFrom patchwork wrap_plots
-#' @importFrom Matrix rowSums
 #' @export
-FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.by = NULL, fill.by = c("group", "feature", "expression"),
+FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.by = NULL, plot.by = c("group", "feature"), fill.by = c("group", "feature", "expression"),
                             cells = NULL, slot = c("data", "counts"), assay = NULL, keep_empty = FALSE, individual = FALSE,
                             plot_type = c("violin", "box", "bar", "dot", "col"),
                             palette = "Paired", palcolor = NULL, alpha = 1,
@@ -3386,17 +3398,207 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
                             legend.position = "right", legend.direction = "vertical",
                             theme_use = "theme_scp", theme_args = list(),
                             combine = TRUE, nrow = NULL, ncol = NULL, byrow = TRUE, force = FALSE, seed = 11) {
+  meta.data <- srt@meta.data
+  meta.data[["cells"]] <- rownames(meta.data)
+  assay <- assay %||% DefaultAssay(srt)
+  slot <- match.arg(slot)
+  exp.data <- slot(srt@assays[[assay]], slot)
+  plot.by <- match.arg(plot.by)
+
+  if (plot.by == "feature") {
+    if (length(group.by) > 1) {
+      stop("The 'group.by' must have a length of 1 when 'plot.by' is set to 'feature'")
+    }
+    message("Setting 'group.by' to 'Features' as 'plot.by' is set to 'feature'")
+    srt@assays[setdiff(names(srt@assays), assay)] <- NULL
+    meta.reshape <- FetchData(srt, vars = c(stat.by, group.by, split.by, bg.by), cells = cells %||% rownames(meta.data), slot = slot)
+    meta.reshape[["cells"]] <- rownames(meta.reshape)
+    meta.reshape <- melt(meta.reshape, measure.vars = stat.by, variable.name = "Features", value.name = "Stat.by")
+    rownames(meta.reshape) <- paste0(meta.reshape[["cells"]], "-", meta.reshape[["Features"]])
+    exp.data <- matrix(0, nrow = 1, ncol = nrow(meta.reshape), dimnames = list("Stat.by", rownames(meta.reshape)))
+    plist <- list()
+    for (g in unique(meta.reshape[[group.by]])) {
+      if (length(rownames(meta.reshape)[meta.reshape[[group.by]] == g]) > 0) {
+        meta.reshape[[g]] <- meta.reshape[["Stat.by"]]
+        p <- ExpressionStatPlot(
+          exp.data = exp.data, meta.data = meta.reshape, stat.by = g, group.by = "Features", split.by = split.by, bg.by = bg.by, plot.by = "group", fill.by = fill.by,
+          cells = rownames(meta.reshape)[meta.reshape[[group.by]] == g], keep_empty = keep_empty, individual = individual,
+          plot_type = plot_type,
+          palette = palette, palcolor = palcolor, alpha = alpha,
+          bg_palette = bg_palette, bg_palcolor = bg_palcolor, bg_alpha = bg_alpha,
+          add_box = add_box, box_color = box_color, box_width = box_width, box_ptsize = box_ptsize,
+          add_point = add_point, pt.color = pt.color, pt.size = pt.size, pt.alpha = pt.alpha, jitter.width = jitter.width,
+          add_trend = add_trend, trend_color = trend_color, trend_linewidth = trend_linewidth, trend_ptsize = trend_ptsize,
+          add_stat = add_stat, stat_color = stat_color, stat_size = stat_size,
+          cells.highlight = cells.highlight, cols.highlight = cols.highlight, sizes.highlight = sizes.highlight, alpha.highlight = alpha.highlight,
+          calculate_coexp = calculate_coexp,
+          same.y.lims = same.y.lims, y.min = y.min, y.max = y.max, y.trans = y.trans, y.nbreaks = y.nbreaks,
+          sort = sort, stack = stack, flip = flip,
+          comparisons = comparisons, ref_group = ref_group, pairwise_method = pairwise_method,
+          multiplegroup_comparisons = multiplegroup_comparisons, multiple_method = multiple_method,
+          sig_label = sig_label, sig_labelsize = sig_labelsize,
+          aspect.ratio = aspect.ratio, title = title, subtitle = subtitle, xlab = xlab, ylab = ylab,
+          legend.position = legend.position, legend.direction = legend.direction,
+          theme_use = theme_use, theme_args = theme_args,
+          force = force, seed = seed
+        )
+        meta.reshape[[g]] <- NULL
+        plist <- append(plist, p)
+      }
+    }
+    group.by <- "Features"
+  } else {
+    plist <- ExpressionStatPlot(
+      exp.data = exp.data, meta.data = meta.data, stat.by = stat.by, group.by = group.by, split.by = split.by, bg.by = bg.by, plot.by = "group", fill.by = fill.by,
+      cells = cells, keep_empty = keep_empty, individual = individual,
+      plot_type = plot_type,
+      palette = palette, palcolor = palcolor, alpha = alpha,
+      bg_palette = bg_palette, bg_palcolor = bg_palcolor, bg_alpha = bg_alpha,
+      add_box = add_box, box_color = box_color, box_width = box_width, box_ptsize = box_ptsize,
+      add_point = add_point, pt.color = pt.color, pt.size = pt.size, pt.alpha = pt.alpha, jitter.width = jitter.width,
+      add_trend = add_trend, trend_color = trend_color, trend_linewidth = trend_linewidth, trend_ptsize = trend_ptsize,
+      add_stat = add_stat, stat_color = stat_color, stat_size = stat_size,
+      cells.highlight = cells.highlight, cols.highlight = cols.highlight, sizes.highlight = sizes.highlight, alpha.highlight = alpha.highlight,
+      calculate_coexp = calculate_coexp,
+      same.y.lims = same.y.lims, y.min = y.min, y.max = y.max, y.trans = y.trans, y.nbreaks = y.nbreaks,
+      sort = sort, stack = stack, flip = flip,
+      comparisons = comparisons, ref_group = ref_group, pairwise_method = pairwise_method,
+      multiplegroup_comparisons = multiplegroup_comparisons, multiple_method = multiple_method,
+      sig_label = sig_label, sig_labelsize = sig_labelsize,
+      aspect.ratio = aspect.ratio, title = title, subtitle = subtitle, xlab = xlab, ylab = ylab,
+      legend.position = legend.position, legend.direction = legend.direction,
+      theme_use = theme_use, theme_args = theme_args,
+      force = force, seed = seed
+    )
+  }
+
+  plist_stack <- list()
+  if (isTRUE(stack) && length(stat.by) > 1 && isFALSE(individual)) {
+    for (g in group.by) {
+      plist_g <- plist[sapply(strsplit(names(plist), ":"), function(x) x[2]) == g]
+      legend <- get_legend(plist_g[[1]])
+      if (isTRUE(flip)) {
+        lab <- textGrob(label = ifelse(is.null(ylab), "Expression level", ylab), hjust = 0.5)
+        plist_g <- lapply(seq_along(plist_g), FUN = function(i) {
+          p <- plist_g[[i]]
+          if (i != 1) {
+            suppressWarnings(p <- p + theme(
+              legend.position = "none",
+              panel.grid = element_blank(),
+              plot.title = element_blank(),
+              plot.subtitle = element_blank(),
+              axis.title = element_blank(),
+              axis.text.y = element_blank(),
+              axis.text.x = element_text(vjust = c(1, 0)),
+              axis.ticks.length.y = unit(0, "pt"),
+              plot.margin = unit(c(0, -0.5, 0, 0), "mm")
+            ))
+          } else {
+            suppressWarnings(p <- p + theme(
+              legend.position = "none",
+              panel.grid = element_blank(),
+              axis.title.x = element_blank(),
+              axis.text.x = element_text(vjust = c(1, 0)),
+              axis.ticks.length.y = unit(0, "pt"),
+              plot.margin = unit(c(0, -0.5, 0, 0), "mm")
+            ))
+          }
+          return(as_grob(p))
+        })
+        gtable <- do.call(cbind, plist_g)
+        gtable <- add_grob(gtable, lab, "bottom", clip = "off")
+        gtable <- add_grob(gtable, legend, legend.position)
+      } else {
+        lab <- textGrob(label = ifelse(is.null(ylab), "Expression level", ylab), rot = 90, hjust = 0.5)
+        plist_g <- lapply(seq_along(plist_g), FUN = function(i) {
+          p <- plist_g[[i]]
+          if (i != length(plist_g)) {
+            suppressWarnings(p <- p + theme(
+              legend.position = "none",
+              panel.grid = element_blank(),
+              axis.title = element_blank(),
+              axis.text.x = element_blank(),
+              axis.text.y = element_text(vjust = c(0, 1)),
+              axis.ticks.length.x = unit(0, "pt"),
+              plot.margin = unit(c(-0.5, 0, 0, 0), "mm")
+            ))
+            if (i == 1) {
+              p <- p + theme(plot.title = element_blank(), plot.subtitle = element_blank())
+            }
+          } else {
+            suppressWarnings(p <- p + theme(
+              legend.position = "none",
+              panel.grid = element_blank(),
+              axis.title.y = element_blank(),
+              axis.text.y = element_text(vjust = c(0, 1)),
+              axis.ticks.length.x = unit(0, "pt"),
+              plot.margin = unit(c(-0.5, 0, 0, 0), "mm")
+            ))
+          }
+          return(as_grob(p))
+        })
+        gtable <- do.call(rbind, plist_g)
+        gtable <- add_grob(gtable, lab, "left", clip = "off")
+        gtable <- add_grob(gtable, legend, legend.position)
+      }
+      gtable <- gtable_add_padding(gtable, unit(c(1, 1, 1, 1), units = "cm"))
+      plot <- wrap_plots(gtable)
+      plist_stack[[g]] <- plot
+    }
+  }
+
+  if (length(plist_stack) > 0) {
+    plist <- plist_stack
+  }
+  if (isTRUE(combine)) {
+    if (length(plist) > 1) {
+      plot <- wrap_plots(plotlist = plist, nrow = nrow, ncol = ncol, byrow = byrow)
+    } else {
+      plot <- plist[[1]]
+    }
+    return(plot)
+  } else {
+    return(plist)
+  }
+}
+
+#' @importFrom Seurat DefaultAssay GetAssayData
+#' @importFrom ggplot2 geom_blank geom_violin geom_rect geom_boxplot geom_count geom_col geom_vline geom_hline layer_data layer_scales position_jitterdodge position_dodge stat_summary scale_x_discrete element_line element_text element_blank annotate mean_sdl after_stat scale_shape_identity
+#' @importFrom Matrix rowSums
+ExpressionStatPlot <- function(exp.data, meta.data, stat.by, group.by = NULL, split.by = NULL, bg.by = NULL, plot.by = c("group", "feature"), fill.by = c("group", "feature", "expression"),
+                               cells = NULL, keep_empty = FALSE, individual = FALSE,
+                               plot_type = c("violin", "box", "bar", "dot", "col"),
+                               palette = "Paired", palcolor = NULL, alpha = 1,
+                               bg_palette = "Paired", bg_palcolor = NULL, bg_alpha = 0.2,
+                               add_box = FALSE, box_color = "black", box_width = 0.1, box_ptsize = 2,
+                               add_point = FALSE, pt.color = "grey30", pt.size = NULL, pt.alpha = 1, jitter.width = 0.5,
+                               add_trend = FALSE, trend_color = "black", trend_linewidth = 1, trend_ptsize = 2,
+                               add_stat = c("none", "mean", "median"), stat_color = "black", stat_size = 1,
+                               cells.highlight = NULL, cols.highlight = "red", sizes.highlight = 1, alpha.highlight = 1,
+                               calculate_coexp = FALSE,
+                               same.y.lims = FALSE, y.min = NULL, y.max = NULL, y.trans = "identity", y.nbreaks = 5,
+                               sort = FALSE, stack = FALSE, flip = FALSE,
+                               comparisons = NULL, ref_group = NULL, pairwise_method = "wilcox.test",
+                               multiplegroup_comparisons = FALSE, multiple_method = "kruskal.test",
+                               sig_label = c("p.signif", "p.format"), sig_labelsize = 3.5,
+                               aspect.ratio = NULL, title = NULL, subtitle = NULL, xlab = NULL, ylab = "Expression level",
+                               legend.position = "right", legend.direction = "vertical",
+                               theme_use = "theme_scp", theme_args = list(),
+                               force = FALSE, seed = 11) {
   set.seed(seed)
 
+  plot.by <- match.arg(plot.by)
   plot_type <- match.arg(plot_type)
   fill.by <- match.arg(fill.by)
-  slot <- match.arg(slot)
   sig_label <- match.arg(sig_label)
   add_stat <- match.arg(add_stat)
 
-  meta.data <- srt@meta.data
-  assay <- assay %||% DefaultAssay(srt)
-  exp.data <- slot(srt@assays[[assay]], slot)
+  if (missing(exp.data)) {
+    exp.data <- matrix(0, nrow = 1, ncol = nrow(meta.data), dimnames = list("", rownames(meta.data)))
+  }
+
+  allfeatures <- rownames(exp.data)
+  allcells <- rownames(meta.data)
 
   if (plot_type == "col") {
     if (isTRUE(add_box) || isTRUE(add_point) || isTRUE(add_trend) || isTRUE(add_stat != "none")) {
@@ -3427,7 +3629,7 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
   }
   for (i in unique(c(group.by, split.by, bg.by))) {
     if (!i %in% colnames(meta.data)) {
-      stop(paste0(i, " is not in the meta.data of srt object."))
+      stop(paste0(i, " is not in the meta.data."))
     }
     if (!is.factor(meta.data[[i]])) {
       meta.data[[i]] <- factor(meta.data[[i]], levels = unique(meta.data[[i]]))
@@ -3449,16 +3651,16 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
     }
   }
   if (!is.null(cells.highlight) && !isTRUE(cells.highlight)) {
-    if (!any(cells.highlight %in% colnames(srt@assays[[1]]))) {
+    if (!any(cells.highlight %in% allcells)) {
       stop("No cells in 'cells.highlight' found.")
     }
-    if (!all(cells.highlight %in% colnames(srt@assays[[1]]))) {
+    if (!all(cells.highlight %in% allcells)) {
       warning("Some cells in 'cells.highlight' not found.", immediate. = TRUE)
     }
-    cells.highlight <- intersect(cells.highlight, colnames(srt@assays[[1]]))
+    cells.highlight <- intersect(cells.highlight, allcells)
   }
   if (isTRUE(cells.highlight)) {
-    cells.highlight <- colnames(srt@assays[[1]])
+    cells.highlight <- allcells
   }
   if (!is.null(cells.highlight) && isFALSE(add_point)) {
     warning("'cells.highlight' is valid only when add_point=TRUE.", immediate. = TRUE)
@@ -3478,7 +3680,7 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
   stat.by <- unique(stat.by)
   features_drop <- stat.by[!stat.by %in% c(rownames(exp.data), colnames(meta.data))]
   if (length(features_drop) > 0) {
-    warning(paste0(features_drop, collapse = ","), " are not in the features of srt.", immediate. = TRUE)
+    warning(paste0(features_drop, collapse = ","), " are not found.", immediate. = TRUE)
     stat.by <- stat.by[!stat.by %in% features_drop]
   }
 
@@ -3490,7 +3692,7 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
       warning(paste(features_meta, collapse = ","), "is not used when calculating co-expression", immediate. = TRUE)
     }
     status <- check_DataType(data = exp.data)
-    message("Data type detected in ", slot, " slot: ", status)
+    message("Data type: ", status)
     if (status %in% c("raw_counts", "raw_normalized_counts")) {
       meta.data[["CoExp"]] <- apply(exp.data[features_gene, , drop = FALSE], 2, function(x) exp(mean(log(x))))
     } else if (status == "log_normalized_counts") {
@@ -3502,18 +3704,18 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
     features_meta <- c(features_meta, "CoExp")
   }
   if (length(features_gene) > 0) {
-    if (all(rownames(srt@assays[[assay]]) %in% features_gene)) {
+    if (all(allfeatures %in% features_gene)) {
       dat_gene <- t(exp.data)
     } else {
       dat_gene <- t(exp.data[features_gene, , drop = FALSE])
     }
   } else {
-    dat_gene <- matrix(nrow = ncol(srt@assays[[1]]), ncol = 0)
+    dat_gene <- matrix(nrow = length(allcells), ncol = 0)
   }
   if (length(features_meta) > 0) {
     dat_meta <- as.matrix(meta.data[, features_meta, drop = FALSE])
   } else {
-    dat_meta <- matrix(nrow = ncol(srt@assays[[1]]), ncol = 0)
+    dat_meta <- matrix(nrow = length(allcells), ncol = 0)
   }
   dat_exp <- cbind(dat_gene, dat_meta)
   stat.by <- unique(stat.by[stat.by %in% c(features_gene, features_meta)])
@@ -3521,12 +3723,15 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
   if (!is.numeric(dat_exp) && !inherits(dat_exp, "Matrix")) {
     stop("'stat.by' must be type of numeric variable.")
   }
-  dat_group <- meta.data[, unique(c(group.by, bg.by, split.by)), drop = FALSE]
+  dat_group <- meta.data[, unique(c("cells", group.by, bg.by, split.by)), drop = FALSE]
   dat_use <- cbind(dat_group, dat_exp[row.names(dat_group), , drop = FALSE])
   if (!is.null(cells)) {
     dat_group <- dat_group[intersect(rownames(dat_group), cells), , drop = FALSE]
     dat_use <- dat_use[intersect(rownames(dat_use), cells), , drop = FALSE]
   }
+  if (nrow(dat_group) == 0) {
+    stop("No specified cells found.")
+  }
 
   if (is.null(pt.size)) {
     pt.size <- min(3000 / nrow(dat_group), 0.5)
@@ -3559,7 +3764,6 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
   }
 
   plist <- list()
-  plist_stack <- list()
 
   comb_list <- list()
   comb <- expand.grid(group_name = group.by, stat_name = stat.by, stringsAsFactors = FALSE)
@@ -3804,7 +4008,7 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
       if (isTRUE(comparisons)) {
         group_use <- names(which(rowSums(table(dat[["group.by"]], dat[["split.by"]]) >= 2) >= 2))
         if (any(rowSums(table(dat[["group.by"]], dat[["split.by"]]) >= 2) >= 3)) {
-          warning("Detected more than 2 groups. Use multiple_method for comparison", immediate. = TRUE)
+          message("Detected more than 2 groups. Use multiple_method for comparison")
           method <- multiple_method
         } else {
           method <- pairwise_method
@@ -4017,98 +4221,11 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
         fill = guide_colorbar(frame.colour = "black", ticks.colour = "black", barheight = 4, barwidth = 1, title.hjust = 0, order = 1)
       )
     }
-    plist[[paste0(f, ":", g, ":", paste0(single_group, collapse = ","), ":", paste0(sp, collapse = ","))]] <- p
+    # plist[[paste0(f, ":", g, ":", paste0(single_group, collapse = ","), ":", paste0(sp, collapse = ","))]] <- p
   })
 
-  if (isTRUE(stack) && length(stat.by) > 1 && isFALSE(individual)) {
-    for (g in group.by) {
-      plist_g <- plist[sapply(strsplit(names(plist), ":"), function(x) x[2]) == g]
-      legend <- get_legend(plist_g[[1]])
-      if (isTRUE(flip)) {
-        lab <- textGrob(label = ifelse(is.null(ylab), "Expression level", ylab), hjust = 0.5)
-        plist_g <- lapply(seq_along(plist_g), FUN = function(i) {
-          p <- plist_g[[i]]
-          if (i != 1) {
-            suppressWarnings(p <- p + theme(
-              legend.position = "none",
-              panel.grid = element_blank(),
-              plot.title = element_blank(),
-              plot.subtitle = element_blank(),
-              axis.title = element_blank(),
-              axis.text.y = element_blank(),
-              axis.text.x = element_text(vjust = c(1, 0)),
-              axis.ticks.length.y = unit(0, "pt"),
-              plot.margin = unit(c(0, -0.5, 0, 0), "mm")
-            ))
-          } else {
-            suppressWarnings(p <- p + theme(
-              legend.position = "none",
-              panel.grid = element_blank(),
-              axis.title.x = element_blank(),
-              axis.text.x = element_text(vjust = c(1, 0)),
-              axis.ticks.length.y = unit(0, "pt"),
-              plot.margin = unit(c(0, -0.5, 0, 0), "mm")
-            ))
-          }
-          return(as_grob(p))
-        })
-        gtable <- do.call(cbind, plist_g)
-        gtable <- add_grob(gtable, lab, "bottom", clip = "off")
-        gtable <- add_grob(gtable, legend, legend.position)
-      } else {
-        lab <- textGrob(label = ifelse(is.null(ylab), "Expression level", ylab), rot = 90, hjust = 0.5)
-        plist_g <- lapply(seq_along(plist_g), FUN = function(i) {
-          p <- plist_g[[i]]
-          if (i != length(plist_g)) {
-            suppressWarnings(p <- p + theme(
-              legend.position = "none",
-              panel.grid = element_blank(),
-              axis.title = element_blank(),
-              axis.text.x = element_blank(),
-              axis.text.y = element_text(vjust = c(0, 1)),
-              axis.ticks.length.x = unit(0, "pt"),
-              plot.margin = unit(c(-0.5, 0, 0, 0), "mm")
-            ))
-            if (i == 1) {
-              p <- p + theme(plot.title = element_blank(), plot.subtitle = element_blank())
-            }
-          } else {
-            suppressWarnings(p <- p + theme(
-              legend.position = "none",
-              panel.grid = element_blank(),
-              axis.title.y = element_blank(),
-              axis.text.y = element_text(vjust = c(0, 1)),
-              axis.ticks.length.x = unit(0, "pt"),
-              plot.margin = unit(c(-0.5, 0, 0, 0), "mm")
-            ))
-          }
-          return(as_grob(p))
-        })
-        gtable <- do.call(rbind, plist_g)
-        gtable <- add_grob(gtable, lab, "left", clip = "off")
-        gtable <- add_grob(gtable, legend, legend.position)
-      }
-      gtable <- gtable_add_padding(gtable, unit(c(1, 1, 1, 1), units = "cm"))
-      plot <- wrap_plots(gtable)
-      plist_stack[[g]] <- plot
-    }
-  }
-
-  if (length(plist_stack) > 0) {
-    plist <- plist_stack
-  }
-  if (isTRUE(combine)) {
-    if (length(plist) > 1) {
-      plot <- wrap_plots(plotlist = plist, nrow = nrow, ncol = ncol, byrow = byrow)
-    } else {
-      plot <- plist[[1]]
-    }
-    return(plot)
-  } else {
-    return(plist)
-  }
+  return(plist)
 }
-
 #' Statistical plot of cells
 #'
 #' @inheritParams StatPlot
@@ -4988,11 +5105,11 @@ FeatureCorPlot <- function(srt, features, group.by = NULL, split.by = NULL, cell
   if (slot == "data" && status != "log_normalized_counts") {
     if (status == "raw_counts") {
       warning("Data in the 'data' slot is raw counts. Perform NormalizeData(LogNormalize) on the data.", immediate. = TRUE)
-      srt <- suppressWarnings(NormalizeData(object = srt, assay = assay, normalization.method = "LogNormalize", verbose = FALSE))
+      srt <- NormalizeData(object = srt, assay = assay, normalization.method = "LogNormalize", verbose = FALSE)
     }
     if (status == "raw_normalized_counts") {
       warning("Data in the 'data' slot is raw_normalized_counts. Perform NormalizeData(LogNormalize) on the data.", immediate. = TRUE)
-      srt <- suppressWarnings(NormalizeData(object = srt, assay = assay, normalization.method = "LogNormalize", verbose = FALSE))
+      srt <- NormalizeData(object = srt, assay = assay, normalization.method = "LogNormalize", verbose = FALSE)
     }
     if (status == "unknown") {
       stop("Data in the 'data' slot is unknown. Please check the data type.")
@@ -7493,6 +7610,7 @@ mestimate <- function(data) {
 #'   }
 #'
 #' @examples
+#' library(dplyr)
 #' data("pancreas_sub")
 #' ht1 <- GroupHeatmap(pancreas_sub,
 #'   features = c(
@@ -8702,7 +8820,7 @@ GroupHeatmap <- function(srt, features = NULL, group.by = NULL, split.by = NULL,
   }
 
   if (isTRUE(fix)) {
-    p <- panel_fix_single(gTree, width = as.numeric(ht_width), height = as.numeric(ht_height), units = units)
+    p <- panel_fix_overall(gTree, width = as.numeric(ht_width), height = as.numeric(ht_height), units = units)
   } else {
     p <- wrap_plots(gTree)
   }
@@ -8726,6 +8844,7 @@ GroupHeatmap <- function(srt, features = NULL, group.by = NULL, split.by = NULL,
 #' @seealso \code{\link{RunDEtest}}
 #'
 #' @examples
+#' library(dplyr)
 #' data("pancreas_sub")
 #' pancreas_sub <- RunDEtest(pancreas_sub, group_by = "CellType")
 #' de_filter <- filter(pancreas_sub@tools$DEtest_CellType$AllMarkers_wilcox, p_val_adj < 0.05 & avg_log2FC > 1)
@@ -9645,7 +9764,7 @@ FeatureHeatmap <- function(srt, features = NULL, cells = NULL, group.by = NULL,
   }
 
   if (isTRUE(fix)) {
-    p <- panel_fix_single(gTree, width = as.numeric(ht_width), height = as.numeric(ht_height), units = units)
+    p <- panel_fix_overall(gTree, width = as.numeric(ht_width), height = as.numeric(ht_height), units = units)
   } else {
     p <- wrap_plots(gTree)
   }
@@ -10547,7 +10666,7 @@ CellCorHeatmap <- function(srt_query, srt_ref = NULL, bulk_ref = NULL,
   }
 
   if (isTRUE(fix)) {
-    p <- panel_fix_single(gTree, width = as.numeric(ht_width), height = as.numeric(ht_height), units = units)
+    p <- panel_fix_overall(gTree, width = as.numeric(ht_width), height = as.numeric(ht_height), units = units)
   } else {
     p <- wrap_plots(gTree)
   }
@@ -11681,7 +11800,7 @@ DynamicHeatmap <- function(srt, lineages, features = NULL, use_fitted = FALSE, b
   }
 
   if (isTRUE(fix)) {
-    p <- panel_fix_single(gTree, width = as.numeric(ht_width), height = as.numeric(ht_height), units = units)
+    p <- panel_fix_overall(gTree, width = as.numeric(ht_width), height = as.numeric(ht_height), units = units)
   } else {
     p <- wrap_plots(gTree)
   }
@@ -13204,7 +13323,7 @@ adjustlayout <- function(graph, layout, width, height = 2, scale = 100, iter = 1
 #' GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", group_use = "Ductal")
 #' GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", group_use = "Ductal", id_use = "GO:0006412")
 #' GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", group_use = "Endocrine", id_use = c("GO:0046903", "GO:0015031", "GO:0007600")) %>%
-#'   panel_fix_single(height = 6) # Because the plot is made by combining, we need to adjust the overall height and width
+#'   panel_fix_overall(height = 6) # As the plot is created by combining, we can adjust the overall height and width directly.
 #'
 #' GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", topTerm = 3, plot_type = "comparison")
 #' GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", topTerm = 3, plot_type = "comparison", direction = "neg")
@@ -13663,8 +13782,8 @@ GSEAPlot <- function(srt, db = "GO_BP", group_by = NULL, test.use = "wilcox", re
         plotlist <- lapply(plotlist[subplots], as_grob)
         rel_heights <- rel_heights[subplots]
         for (i in seq_along(plotlist)) {
-          plotlist[[i]] <- panel_fix_single(plotlist[[i]], height = rel_heights[i], units = "null", margin = 0, respect = TRUE, return_grob = TRUE)
-          plotlist[[i]] <- panel_fix_single(plotlist[[i]], width = rel_width, units = "null", margin = 0, respect = TRUE, return_grob = TRUE)
+          plotlist[[i]] <- panel_fix_overall(plotlist[[i]], height = rel_heights[i], units = "null", margin = 0, respect = TRUE, return_grob = TRUE)
+          plotlist[[i]] <- panel_fix_overall(plotlist[[i]], width = rel_width, units = "null", margin = 0, respect = TRUE, return_grob = TRUE)
         }
         p_out <- do.call(rbind, c(plotlist, size = "first"))
 
diff --git a/R/SCP-projection.R b/R/SCP-projection.R
index 79d2115d..865226a2 100644
--- a/R/SCP-projection.R
+++ b/R/SCP-projection.R
@@ -529,7 +529,7 @@ RunCSSMap <- function(srt_query, srt_ref, query_assay = NULL, ref_assay = srt_re
   CSSmodel <- srt_ref[[ref_css]]@misc$model
   raw_assay <- DefaultAssay(srt_query)
   DefaultAssay(srt_query) <- query_assay
-  srt_query <- simspec::css_project(object = srt_query, model = CSSmodel)
+  srt_query <- invoke(.fn = get("css_project", envir = getNamespace("simspec")), .args = list(object = srt_query, model = CSSmodel))
   DefaultAssay(srt_query) <- raw_assay
 
   message("Run UMAP projection")
diff --git a/R/SCP-workflow.R b/R/SCP-workflow.R
index d4ba65df..c90761f9 100644
--- a/R/SCP-workflow.R
+++ b/R/SCP-workflow.R
@@ -187,11 +187,11 @@ check_srtList <- function(srtList, batch, assay = NULL,
     if (isTRUE(do_normalization)) {
       if (normalization_method == "LogNormalize") {
         cat("Perform NormalizeData(LogNormalize) on the data ", i, "/", length(srtList), " of the srtList...\n", sep = "")
-        srtList[[i]] <- suppressWarnings(NormalizeData(object = srtList[[i]], assay = assay, normalization.method = "LogNormalize", verbose = FALSE))
+        srtList[[i]] <- NormalizeData(object = srtList[[i]], assay = assay, normalization.method = "LogNormalize", verbose = FALSE)
       }
       if (normalization_method == "TFIDF") {
         cat("Perform RunTFIDF on the data ", i, "/", length(srtList), " of the srtList...\n", sep = "")
-        srtList[[i]] <- suppressWarnings(RunTFIDF(object = srtList[[i]], assay = assay, verbose = FALSE))
+        srtList[[i]] <- RunTFIDF(object = srtList[[i]], assay = assay, verbose = FALSE)
       }
     } else if (is.null(do_normalization)) {
       status <- check_DataType(srtList[[i]], slot = "data", assay = assay)
@@ -201,11 +201,11 @@ check_srtList <- function(srtList, batch, assay = NULL,
       if (status %in% c("raw_counts", "raw_normalized_counts")) {
         if (normalization_method == "LogNormalize") {
           cat("Data ", i, "/", length(srtList), " of the srtList is ", status, ". Perform NormalizeData(LogNormalize) on the data ...\n", sep = "")
-          srtList[[i]] <- suppressWarnings(NormalizeData(object = srtList[[i]], assay = assay, normalization.method = "LogNormalize", verbose = FALSE))
+          srtList[[i]] <- NormalizeData(object = srtList[[i]], assay = assay, normalization.method = "LogNormalize", verbose = FALSE)
         }
         if (normalization_method == "TFIDF") {
           cat("Data ", i, "/", length(srtList), " of the srtList is ", status, ". Perform RunTFIDF on the data ...\n", sep = "")
-          srtList[[i]] <- suppressWarnings(RunTFIDF(object = srtList[[i]], assay = assay, verbose = FALSE))
+          srtList[[i]] <- RunTFIDF(object = srtList[[i]], assay = assay, verbose = FALSE)
         }
       }
       if (status == "unknown") {
@@ -216,11 +216,11 @@ check_srtList <- function(srtList, batch, assay = NULL,
       if (isTRUE(do_HVF_finding) || is.null(do_HVF_finding) || length(VariableFeatures(srtList[[i]], assay = assay)) == 0) {
         # if (type == "RNA") {
         cat("Perform FindVariableFeatures on the data ", i, "/", length(srtList), " of the srtList...\n", sep = "")
-        srtList[[i]] <- suppressWarnings(FindVariableFeatures(srtList[[i]], assay = assay, nfeatures = nHVF, selection.method = HVF_method, verbose = FALSE))
+        srtList[[i]] <- FindVariableFeatures(srtList[[i]], assay = assay, nfeatures = nHVF, selection.method = HVF_method, verbose = FALSE)
         # }
         # if (type == "Chromatin") {
         #   cat("Perform FindTopFeatures on the data ", i, "/", length(srtList), " of the srtList...\n", sep = "")
-        #   srtList[[i]] <- suppressWarnings(FindTopFeatures(srtList[[i]], assay = assay, min.cutoff = HVF_min_cutoff, verbose = FALSE))
+        #   srtList[[i]] <- FindTopFeatures(srtList[[i]], assay = assay, min.cutoff = HVF_min_cutoff, verbose = FALSE)
         # }
       }
     }
@@ -1313,7 +1313,7 @@ Seurat_integrate <- function(srtMerge = NULL, batch = NULL, append = TRUE, srtLi
     if (is.null(FindIntegrationAnchors_params[["dims"]])) {
       FindIntegrationAnchors_params[["dims"]] <- 2:min(linear_reduction_dims, 30)
     }
-    srtMerge <- suppressWarnings(RunTFIDF(object = srtMerge, assay = DefaultAssay(srtMerge), verbose = FALSE))
+    srtMerge <- RunTFIDF(object = srtMerge, assay = DefaultAssay(srtMerge), verbose = FALSE)
     srtMerge <- RunDimReduction(
       srt = srtMerge, prefix = "", features = HVF, assay = DefaultAssay(srtMerge),
       linear_reduction = "svd", linear_reduction_dims = linear_reduction_dims, linear_reduction_params = linear_reduction_params, force_linear_reduction = force_linear_reduction,
@@ -1410,9 +1410,7 @@ Seurat_integrate <- function(srtMerge = NULL, batch = NULL, append = TRUE, srtLi
     for (nm in names(FindIntegrationAnchors_params)) {
       params1[[nm]] <- FindIntegrationAnchors_params[[nm]]
     }
-    suppressWarnings({
-      srt_anchors <- invoke(.fn = FindIntegrationAnchors, .args = params1)
-    })
+    srt_anchors <- invoke(.fn = FindIntegrationAnchors, .args = params1)
 
     cat(paste0("[", Sys.time(), "]", " Perform integration(Seurat) on the data...\n"))
     params2 <- list(
@@ -2820,7 +2818,7 @@ CSS_integrate <- function(srtMerge = NULL, batch = NULL, append = TRUE, srtList
   for (nm in names(CSS_params)) {
     params[[nm]] <- CSS_params[[nm]]
   }
-  srtIntegrated <- invoke(.fn = simspec::cluster_sim_spectrum, .args = params)
+  srtIntegrated <- invoke(.fn = get("cluster_sim_spectrum", envir = getNamespace("simspec")), .args = params)
 
   if (any(is.na(srtIntegrated@reductions[["CSS"]]@cell.embeddings))) {
     stop("NA detected in the CSS embeddings. You can try to use a lower resolution value in the CSS_param.")
diff --git a/R/reexports.R b/R/reexports.R
index c71d14b1..eef1a972 100644
--- a/R/reexports.R
+++ b/R/reexports.R
@@ -11,18 +11,3 @@ rlang::`%||%`
 #' @export
 #'
 dplyr::`%>%`
-
-#' @importFrom dplyr filter
-#' @export
-#'
-dplyr::filter
-
-#' @importFrom dplyr group_by
-#' @export
-#'
-dplyr::group_by
-
-#' @importFrom dplyr top_n
-#' @export
-#'
-dplyr::top_n
diff --git a/man/FeatureHeatmap.Rd b/man/FeatureHeatmap.Rd
index cf76f6d9..fa5d7063 100644
--- a/man/FeatureHeatmap.Rd
+++ b/man/FeatureHeatmap.Rd
@@ -315,6 +315,7 @@ FeatureHeatmap(
 FeatureHeatmap
 }
 \examples{
+library(dplyr)
 data("pancreas_sub")
 pancreas_sub <- RunDEtest(pancreas_sub, group_by = "CellType")
 de_filter <- filter(pancreas_sub@tools$DEtest_CellType$AllMarkers_wilcox, p_val_adj < 0.05 & avg_log2FC > 1)
diff --git a/man/FeatureStatPlot.Rd b/man/FeatureStatPlot.Rd
index febf0cfb..ebc6e03a 100644
--- a/man/FeatureStatPlot.Rd
+++ b/man/FeatureStatPlot.Rd
@@ -10,6 +10,7 @@ FeatureStatPlot(
   group.by = NULL,
   split.by = NULL,
   bg.by = NULL,
+  plot.by = c("group", "feature"),
   fill.by = c("group", "feature", "expression"),
   cells = NULL,
   slot = c("data", "counts"),
@@ -87,6 +88,8 @@ FeatureStatPlot(
 
 \item{bg.by}{A character vector specifying the variable to use as the background color. Default is NULL.}
 
+\item{plot.by}{A character vector specifying how to plot the data, by group or feature. Possible values are "group", "feature". Default is "group".}
+
 \item{fill.by}{A string specifying what to fill the plot by. Possible values are "group", "feature", or "expression". Default is "group".}
 
 \item{cells}{A character vector specifying the cells to include in the plot. Default is NULL.}
@@ -258,5 +261,15 @@ FeatureStatPlot(pancreas_sub,
   ),
   fill.by = "feature", plot_type = "box",
   group.by = "SubCellType", bg.by = "CellType", stack = TRUE, flip = TRUE
-) \%>\% panel_fix_single(width = 8, height = 5) # Because the plot is made by combining, we want to adjust the overall height and width
+) \%>\% panel_fix_overall(width = 8, height = 5) # As the plot is created by combining, we can adjust the overall height and width directly.
+
+FeatureStatPlot(pancreas_sub, stat.by = c("G2M_score", "Fev"), group.by = "CellType", plot.by = "feature")
+FeatureStatPlot(pancreas_sub, stat.by = c(
+  "Sox9", "Anxa2", "Bicc1", # Ductal
+  "Neurog3", "Hes6", # EPs
+  "Fev", "Neurod1", # Pre-endocrine
+  "Rbp4", "Pyy", # Endocrine
+  "Ins1", "Gcg", "Sst", "Ghrl" # Beta, Alpha, Delta, Epsilon
+), group.by = "SubCellType", plot.by = "feature", stack = TRUE)
+
 }
diff --git a/man/GSEAPlot.Rd b/man/GSEAPlot.Rd
index 365b53df..086277ed 100644
--- a/man/GSEAPlot.Rd
+++ b/man/GSEAPlot.Rd
@@ -185,7 +185,7 @@ pancreas_sub <- RunGSEA(pancreas_sub, group_by = "CellType", db = "GO_BP", speci
 GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", group_use = "Ductal")
 GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", group_use = "Ductal", id_use = "GO:0006412")
 GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", group_use = "Endocrine", id_use = c("GO:0046903", "GO:0015031", "GO:0007600")) \%>\%
-  panel_fix_single(height = 6) # Because the plot is made by combining, we need to adjust the overall height and width
+  panel_fix_overall(height = 6) # As the plot is created by combining, we can adjust the overall height and width directly.
 
 GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", topTerm = 3, plot_type = "comparison")
 GSEAPlot(pancreas_sub, db = "GO_BP", group_by = "CellType", topTerm = 3, plot_type = "comparison", direction = "neg")
diff --git a/man/GroupHeatmap.Rd b/man/GroupHeatmap.Rd
index 76bc3194..17ec0d2c 100644
--- a/man/GroupHeatmap.Rd
+++ b/man/GroupHeatmap.Rd
@@ -363,6 +363,7 @@ A list with the following elements:
 GroupHeatmap
 }
 \examples{
+library(dplyr)
 data("pancreas_sub")
 ht1 <- GroupHeatmap(pancreas_sub,
   features = c(
diff --git a/man/RunDEtest.Rd b/man/RunDEtest.Rd
index 862802b9..f84df542 100644
--- a/man/RunDEtest.Rd
+++ b/man/RunDEtest.Rd
@@ -151,6 +151,7 @@ This function utilizes the Seurat package to perform a differential expression (
 Users have the flexibility to specify custom cell groups, marker types, and various options for DE analysis.
 }
 \examples{
+library(dplyr)
 data("pancreas_sub")
 pancreas_sub <- RunDEtest(pancreas_sub, group_by = "SubCellType")
 AllMarkers <- filter(pancreas_sub@tools$DEtest_SubCellType$AllMarkers_wilcox, p_val_adj < 0.05 & avg_log2FC > 1)
diff --git a/man/RunSCExplorer.Rd b/man/RunSCExplorer.Rd
index 9e5058c8..ca2401ea 100644
--- a/man/RunSCExplorer.Rd
+++ b/man/RunSCExplorer.Rd
@@ -95,7 +95,8 @@ app <- RunSCExplorer(
   initial_dataset = "mouse_pancreas",
   initial_group = "CellType",
   initial_feature = "Neurog3",
-  workers = 2, overwrite = TRUE
+  session_workers = 2,
+  overwrite = TRUE
 )
 list.files("./SCExplorer") # This directory can be used as site directory for Shiny Server.
 
diff --git a/man/panel_fix.Rd b/man/panel_fix.Rd
index 3c72fab9..1a28b865 100644
--- a/man/panel_fix.Rd
+++ b/man/panel_fix.Rd
@@ -2,7 +2,7 @@
 % Please edit documentation in R/SCP-plot.R
 \name{panel_fix}
 \alias{panel_fix}
-\alias{panel_fix_single}
+\alias{panel_fix_overall}
 \title{Set the panel width/height of a plot object to a fixed value.}
 \usage{
 panel_fix(
@@ -24,7 +24,7 @@ panel_fix(
   ...
 )
 
-panel_fix_single(
+panel_fix_overall(
   x,
   panel_index = NULL,
   respect = NULL,
diff --git a/man/reexports.Rd b/man/reexports.Rd
index b8cb05b0..2ef063f3 100644
--- a/man/reexports.Rd
+++ b/man/reexports.Rd
@@ -5,9 +5,6 @@
 \alias{reexports}
 \alias{\%||\%}
 \alias{\%>\%}
-\alias{filter}
-\alias{group_by}
-\alias{top_n}
 \title{Objects exported from other packages}
 \keyword{internal}
 \description{
@@ -15,7 +12,7 @@ These objects are imported from other packages. Follow the links
 below to see their documentation.
 
 \describe{
-  \item{dplyr}{\code{\link[dplyr:reexports]{\%>\%}}, \code{\link[dplyr]{filter}}, \code{\link[dplyr]{group_by}}, \code{\link[dplyr]{top_n}}}
+  \item{dplyr}{\code{\link[dplyr:reexports]{\%>\%}}}
 
   \item{rlang}{\code{\link[rlang:op-null-default]{\%||\%}}}
 }}
diff --git a/man/theme_blank.Rd b/man/theme_blank.Rd
index 8af376de..712c7c79 100644
--- a/man/theme_blank.Rd
+++ b/man/theme_blank.Rd
@@ -30,7 +30,8 @@ theme_blank(
 \item{...}{Arguments passed to the \code{\link[ggplot2]{theme}}.}
 }
 \description{
-Blank theme
+This function creates a theme with all elements blank except for axis lines and labels.
+It can optionally add coordinate axes in the plot.
 }
 \examples{
 library(ggplot2)