diff --git a/R/SCP-analysis.R b/R/SCP-analysis.R
index 7a9238c1..9cc13f2a 100644
--- a/R/SCP-analysis.R
+++ b/R/SCP-analysis.R
@@ -67,7 +67,6 @@
 #' @importFrom reshape2 dcast melt
 #' @importFrom biomaRt listEnsemblArchives useMart listDatasets useDataset getBM listAttributes useEnsembl
 #' @export
-#'
 GeneConvert <- function(geneID, geneID_from_IDtype = "symbol", geneID_to_IDtype = "entrez_id",
                         species_from = "Homo_sapiens", species_to = NULL,
                         Ensembl_version = 103, biomart = NULL, mirror = NULL, max_tries = 5) {
@@ -3400,7 +3399,6 @@ PrepareDB <- function(species = c("Homo_sapiens", "Mus_musculus"),
 #' @importFrom BiocParallel bplapply bpprogressbar<- bpRNGseed<- bpworkers ipcid ipclock ipcunlock
 #' @importFrom clusterProfiler enricher simplify
 #' @export
-#'
 RunEnrichment <- function(srt = NULL, group_by = NULL, test.use = "wilcox", DE_threshold = "avg_log2FC > 0 & p_val_adj < 0.05",
                           geneID = NULL, geneID_groups = NULL, geneID_exclude = NULL, IDtype = "symbol", result_IDtype = "symbol", species = "Homo_sapiens",
                           db = "GO_BP", db_update = FALSE, db_version = "latest", db_combine = FALSE, convert_species = TRUE, Ensembl_version = 103, mirror = NULL,
@@ -3665,7 +3663,6 @@ RunEnrichment <- function(srt = NULL, group_by = NULL, test.use = "wilcox", DE_t
 #' @importFrom BiocParallel bplapply bpprogressbar<- bpRNGseed<- bpworkers ipcid ipclock ipcunlock
 #' @importFrom clusterProfiler GSEA simplify
 #' @export
-#'
 RunGSEA <- function(srt = NULL, group_by = NULL, test.use = "wilcox", DE_threshold = "p_val_adj < 0.05", scoreType = "std",
                     geneID = NULL, geneScore = NULL, geneID_groups = NULL, geneID_exclude = NULL, IDtype = "symbol", result_IDtype = "symbol", species = "Homo_sapiens",
                     db = "GO_BP", db_update = FALSE, db_version = "latest", db_combine = FALSE, convert_species = TRUE, Ensembl_version = 103, mirror = NULL,
@@ -4984,7 +4981,6 @@ py_to_r_auto <- function(x) {
 #' @importFrom Seurat GetAssayData
 #' @importFrom Matrix t
 #' @export
-#'
 srt_to_adata <- function(srt, features = NULL,
                          assay_X = "RNA", slot_X = "counts",
                          assay_layers = c("spliced", "unspliced"), slot_layers = "counts",
@@ -5418,7 +5414,6 @@ check_python_element <- function(x, depth = maxDepth(x)) {
 #' CellDimPlot(pancreas_sub, group.by = "SubCellType", reduction = "PAGAUMAP2D", paga = pancreas_sub@misc$paga)
 #'
 #' @export
-#'
 RunPAGA <- function(srt = NULL, assay_X = "RNA", slot_X = "counts", assay_layers = c("spliced", "unspliced"), slot_layers = "counts",
                     adata = NULL, group_by = NULL,
                     linear_reduction = NULL, nonlinear_reduction = NULL, basis = NULL,
@@ -5529,7 +5524,6 @@ RunPAGA <- function(srt = NULL, assay_X = "RNA", slot_X = "counts", assay_layers
 #' pancreas_sub <- RunSCVELO(srt = pancreas_sub, assay_X = "SCT", group_by = "SubCellType", linear_reduction = "Standardpca", nonlinear_reduction = "StandardTSNE2D")
 #'
 #' @export
-#'
 RunSCVELO <- function(srt = NULL, assay_X = "RNA", slot_X = "counts", assay_layers = c("spliced", "unspliced"), slot_layers = "counts",
                       adata = NULL, group_by = NULL,
                       linear_reduction = NULL, nonlinear_reduction = NULL, basis = NULL,
@@ -5636,7 +5630,6 @@ RunSCVELO <- function(srt = NULL, assay_X = "RNA", slot_X = "counts", assay_laye
 #' FeatureDimPlot(pancreas_sub, paste0(c("Alpha", "Beta", "Delta", "Epsilon"), "_diff_potential"))
 #'
 #' @export
-#'
 RunPalantir <- function(srt = NULL, assay_X = "RNA", slot_X = "counts", assay_layers = c("spliced", "unspliced"), slot_layers = "counts",
                         adata = NULL, group_by = NULL,
                         linear_reduction = NULL, nonlinear_reduction = NULL, basis = NULL,
diff --git a/R/SCP-plot.R b/R/SCP-plot.R
index cac409f9..90dca74f 100644
--- a/R/SCP-plot.R
+++ b/R/SCP-plot.R
@@ -3363,7 +3363,16 @@ FeatureDimPlot3D <- function(srt, features = NULL, reduction = NULL, dims = c(1,
 #'   group.by = "SubCellType", bg.by = "CellType", stack = TRUE, flip = TRUE
 #' ) %>% 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("Neurog3", "Rbp4", "Ins1"), group.by = "CellType", plot.by = "feature")
+#' FeatureStatPlot(pancreas_sub,
+#'   stat.by = c("Neurog3", "Rbp4", "Ins1"), group.by = "CellType", plot.by = "feature",
+#'   multiplegroup_comparisons = TRUE, sig_label = "p.format", sig_labelsize = 4
+#' )
+#' FeatureStatPlot(pancreas_sub,
+#'   stat.by = c("Neurog3", "Rbp4", "Ins1"), group.by = "CellType", plot.by = "feature",
+#'   comparisons = list(c("Neurog3", "Rbp4"), c("Rbp4", "Ins1")),
+#'   stack = TRUE
+#' )
 #' FeatureStatPlot(pancreas_sub, stat.by = c(
 #'   "Sox9", "Anxa2", "Bicc1", # Ductal
 #'   "Neurog3", "Hes6", # EPs
@@ -3372,6 +3381,13 @@ FeatureDimPlot3D <- function(srt, features = NULL, reduction = NULL, dims = c(1,
 #'   "Ins1", "Gcg", "Sst", "Ghrl" # Beta, Alpha, Delta, Epsilon
 #' ), group.by = "SubCellType", plot.by = "feature", stack = TRUE)
 #'
+#' data <- pancreas_sub@assays$RNA@data
+#' pancreas_sub@assays$RNA@scale.data <- as.matrix(data / rowMeans(data))
+#' FeatureStatPlot(pancreas_sub,
+#'   stat.by = c("Neurog3", "Rbp4", "Ins1"), group.by = "CellType",
+#'   slot = "scale.data", ylab = "FoldChange", same.y.lims = TRUE, y.max = 4
+#' )
+#'
 #' @importFrom Seurat FetchData
 #' @importFrom reshape2 melt
 #' @importFrom gtable gtable_add_cols gtable_add_rows gtable_add_grob gtable_add_padding
@@ -3379,7 +3395,7 @@ FeatureDimPlot3D <- function(srt, features = NULL, reduction = NULL, dims = c(1,
 #' @importFrom patchwork wrap_plots
 #' @export
 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,
+                            cells = NULL, slot = "data", assay = 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,
@@ -3401,7 +3417,6 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
   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)
 
@@ -3409,9 +3424,12 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
     if (length(group.by) > 1) {
       stop("The 'group.by' must have a length of 1 when 'plot.by' is set to 'feature'")
     }
+    if (!is.null(bg.by)) {
+      message("'bg.by' is invalid 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 <- FetchData(srt, vars = c(stat.by, group.by, split.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"]])
@@ -3421,7 +3439,7 @@ FeatureStatPlot <- function(srt, stat.by, group.by = NULL, split.by = NULL, bg.b
       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,
+          exp.data = exp.data, meta.data = meta.reshape, stat.by = g, group.by = "Features", split.by = split.by, bg.by = NULL, 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,
diff --git a/man/FeatureStatPlot.Rd b/man/FeatureStatPlot.Rd
index ebc6e03a..be184270 100644
--- a/man/FeatureStatPlot.Rd
+++ b/man/FeatureStatPlot.Rd
@@ -13,7 +13,7 @@ FeatureStatPlot(
   plot.by = c("group", "feature"),
   fill.by = c("group", "feature", "expression"),
   cells = NULL,
-  slot = c("data", "counts"),
+  slot = "data",
   assay = NULL,
   keep_empty = FALSE,
   individual = FALSE,
@@ -263,7 +263,16 @@ FeatureStatPlot(pancreas_sub,
   group.by = "SubCellType", bg.by = "CellType", stack = TRUE, flip = TRUE
 ) \%>\% 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("Neurog3", "Rbp4", "Ins1"), group.by = "CellType", plot.by = "feature")
+FeatureStatPlot(pancreas_sub,
+  stat.by = c("Neurog3", "Rbp4", "Ins1"), group.by = "CellType", plot.by = "feature",
+  multiplegroup_comparisons = TRUE, sig_label = "p.format", sig_labelsize = 4
+)
+FeatureStatPlot(pancreas_sub,
+  stat.by = c("Neurog3", "Rbp4", "Ins1"), group.by = "CellType", plot.by = "feature",
+  comparisons = list(c("Neurog3", "Rbp4"), c("Rbp4", "Ins1")),
+  stack = TRUE
+)
 FeatureStatPlot(pancreas_sub, stat.by = c(
   "Sox9", "Anxa2", "Bicc1", # Ductal
   "Neurog3", "Hes6", # EPs
@@ -272,4 +281,11 @@ FeatureStatPlot(pancreas_sub, stat.by = c(
   "Ins1", "Gcg", "Sst", "Ghrl" # Beta, Alpha, Delta, Epsilon
 ), group.by = "SubCellType", plot.by = "feature", stack = TRUE)
 
+data <- pancreas_sub@assays$RNA@data
+pancreas_sub@assays$RNA@scale.data <- as.matrix(data / rowMeans(data))
+FeatureStatPlot(pancreas_sub,
+  stat.by = c("Neurog3", "Rbp4", "Ins1"), group.by = "CellType",
+  slot = "scale.data", ylab = "FoldChange", same.y.lims = TRUE, y.max = 4
+)
+
 }