diff --git a/R/calc-diversity.R b/R/calc-diversity.R
index 93aca228..32010672 100644
--- a/R/calc-diversity.R
+++ b/R/calc-diversity.R
@@ -286,6 +286,29 @@ calc_diversity <- function(input, data_col, cluster_col = NULL,
 #'   plot legend
 #' - 'none', do not display the number of cells plotted
 #'
+#' @param p_label Specification indicating how p-values should be labeled on
+#' plot, this can one of the following:
+#'
+#' - 'none', do not display p-values
+#' - 'all', show p-values for all groups
+#' - A named vector providing p-value cutoffs and labels to display,
+#'   e.g. `c('*' = 0.05, '**' = 0.01, '***' = 0.001)`. The keyword 'value' can
+#'   be used to display the p-value for those less than a certain cutoff,
+#'   e.g. `c(value = 0.05, ns = 1.1)` will show significant p-values, all others
+#'   will be labeled 'ns'.
+#'
+#' @param p_method Method to use for calculating p-values, by default when
+#' comparing two groups a t-test will be used.
+#' When comparing more than two groups the Kruskal-Wallis test will be used.
+#' p-values are adjusted for
+#' multiple testing using Bonferroni correction. Possible methods include:
+#'
+#' - 't', two sample t-test performed with `stats::t.test()`
+#' - 'wilcox', Wilcoxon rank sum test performed with `stats::wilcox.test()`
+#' - 'kruskal', Kruskal-Wallis test performed with `stats::kruskal.test()`
+#'
+#' @param p_file File path to save table containing p-values for each
+#' comparison.
 #' @param label_params Named list providing additional parameters to modify
 #' n label aesthetics, e.g. list(size = 4, color = "red")
 #' @param ... Additional arguments to pass to ggplot2, e.g. color, fill, size,
@@ -360,7 +383,8 @@ plot_diversity <- function(input, data_col, cluster_col = NULL,
                            plot_colors = NULL,
                            plot_lvls = names(plot_colors), panel_nrow = NULL,
                            panel_scales = "free", n_label = NULL,
-                           label_params = list(), ...) {
+                           p_label = "all", p_method = NULL,
+                           p_file = NULL, label_params = list(), ...) {
 
   # Check that columns are present in object
   .check_obj_cols(
@@ -378,6 +402,8 @@ plot_diversity <- function(input, data_col, cluster_col = NULL,
 
   .check_group_cols(cluster_col, group_col, input)
 
+  .check_possible_values(p_method = c("t", "wilcox", "kruskal"))
+
   if (length(method) == 1 && is.null(names(method))) {
     nm <- as.character(substitute(method))
     nm <- dplyr::last(nm)
@@ -474,8 +500,12 @@ plot_diversity <- function(input, data_col, cluster_col = NULL,
     gg_args$clst      <- cluster_col
     gg_args$grp       <- group_col
     gg_args$add_zeros <- FALSE
-    gg_args$p_label   <- "all"
-    gg_args$p_corner  <- TRUE
+
+    gg_args$p_label   <- p_label
+    gg_args$p_method  <- p_method
+    gg_args$p_file    <- p_file
+    gg_args$p_x       <- "center"
+
     gg_args$method    <- "boxplot"
     gg_args           <- append(gg_args, list(p_grp = NULL))
 
diff --git a/R/plotting-labels.R b/R/plotting-labels.R
index 0b7ab230..0e2d0bda 100644
--- a/R/plotting-labels.R
+++ b/R/plotting-labels.R
@@ -542,15 +542,31 @@
 #'
 #' @param x character vector of labels
 #' @param base base justification, increase this value to increase spacing
+#' @param exclude regular expression specifying characters/patterns to remove
+#' from x before calculating justification
+#' @param side side of label to add padding
 #' @return numeric vector of justification values
 #' @noRd
-.get_label_just <- function(x, base = 0.5) {
+.get_label_just <- function(x, base = 0.5, exclude = NULL, side = "right") {
+
+  .check_possible_values(side = c("right", "left"), .internal = TRUE)
+
+  # Remove excluded characters from x
+  if (!is.null(exclude)) x <- gsub(exclude, "", x)
 
   # Assumed height x width ratio
   char_h_w <- 1.5
 
   len <- nchar(x)
-  res <- 1 + (1 / len * (base * char_h_w))
+
+  # calculate hjust based on side
+  jst <- (1 / len * (base * char_h_w))
+
+  res <- switch(
+    side,
+    right = 1 + jst,
+    left  = 0 - jst
+  )
 
   res
 }
diff --git a/R/plotting-utils.R b/R/plotting-utils.R
index 1e733456..a5063118 100644
--- a/R/plotting-utils.R
+++ b/R/plotting-utils.R
@@ -275,8 +275,8 @@ djvdj_theme <- function(base_size = 11, base_family = "",
 #' @param p_method Method to calculate p-values
 #' @param p_grp Variable to use for grouping samples when calculating p-values.
 #' A separate p-value will be calculated for each label in p_grp.
-#' @param p_corner Should p-value be shown in the top right corner, if `FALSE`
-#' p-value will be shown above each group
+#' @param p_x Manually set x coordinate for p-label, provide x specification or
+#' 'right', 'left', or 'center'.
 #' @param p_file File path to save p-value csv
 #' @param label_params Named list with specifications to modify label
 #' aesthetics
@@ -292,7 +292,7 @@ djvdj_theme <- function(base_size = 11, base_family = "",
 .create_grouped_plot <- function(df_in, x, y, clst, grp, method = "bar",
                                  n_label = NULL, p_label = c(value = 0.05),
                                  p_y = y, p_method = NULL, p_grp = x, p_file = NULL,
-                                 p_corner = FALSE, label_params = list(),
+                                 p_x = NULL, label_params = list(),
                                  show_points = TRUE, add_zeros = TRUE,
                                  show_zeros = TRUE, ...) {
 
@@ -322,8 +322,7 @@ djvdj_theme <- function(base_size = 11, base_family = "",
   # by default only show significant p-values
   if (identical(p_label, "all")) p_label <- c(value = Inf)
 
-  add_p    <- !identical(p_label, "none")
-  p_corner <- add_p && (p_corner || is.null(p_grp))
+  add_p <- !identical(p_label, "none")
 
   if (add_p) {
     p <- .calc_pvalue(
@@ -336,8 +335,13 @@ djvdj_theme <- function(base_size = 11, base_family = "",
     p <- dplyr::group_by(p, !!!syms(c(p_grp, "p_value")))
 
     p <- dplyr::summarize(
-      p, y_min = min(!!sym(y)), y_max = max(!!sym(y)),
-      .groups = "drop"
+      p,
+      y_min = min(!!sym(y)),
+      y_max = max(!!sym(y)),
+      gap   = (max(.data$y_max) - min(.data$y_min)) * 0.05,
+      y     = .data$y_max + .data$gap,
+      n_x   = dplyr::n_distinct(!!sym(x)),  # number of x-axis groups for each
+      .groups = "drop"                      # p-value plotted, used for p_x
     )
 
     # Format p-values
@@ -352,9 +356,20 @@ djvdj_theme <- function(base_size = 11, base_family = "",
 
     if (nrow(p) == 0) add_p <- FALSE
 
-    # If only a single p-value plot in upper corner
-    if (p_corner) {
-      p <- dplyr::mutate(p, p_lab = paste0("italic(p) == ", .data$p_lab))
+    # Adjust p_x
+    if (is.null(p_x) && is.null(p_grp)) p_x <- "center"
+
+    if (!is.null(p_x)) {
+      p <- dplyr::mutate(
+        p,
+        !!sym(x) := switch(
+          as.character(p_x),  # EXPR should evaluate to character
+          right  = Inf,
+          left   = -Inf,
+          center = (n_x / 2) + 0.5,
+          p_x
+        )
+      )
     }
   }
 
@@ -395,7 +410,7 @@ djvdj_theme <- function(base_size = 11, base_family = "",
     )
 
     # Need to adjust y_max based on mean and sd to position p-value label
-    if (add_p && !p_corner) {
+    if (add_p && is.null(p_x)) {
       y_dat <- dplyr::group_by(df_in, !!sym(p_grp))
 
       y_dat <- dplyr::mutate(
@@ -423,14 +438,6 @@ djvdj_theme <- function(base_size = 11, base_family = "",
 
   # Add p-values
   if (add_p) {
-    if (!p_corner) {
-      p <- dplyr::mutate(
-        p,
-        gap = (max(.data$y_max) - min(.data$y_min)) * 0.05,
-        y   = .data$y_max + .data$gap
-      )
-    }
-
     all_sym <- !any(grepl("[a-zA-Z0-9]", names(p_label)))
 
     label_params <- .parse_label_params(label_params)$p
@@ -439,11 +446,28 @@ djvdj_theme <- function(base_size = 11, base_family = "",
       x = !!sym(x), y = .data$y, label = .data$p_lab, fill = NULL
     )
 
-    # If only one p-value calculated for plot, position in top right corner
-    if (p_corner) {
-      label_params$mapping$x <- label_params$mapping$y <- Inf
-      label_params$hjust     <- label_params$hjust %||% .get_label_just(p$p_lab)
+    # If only one p-value calculated for plot, position center of panel
+    # include 'p = ' when setting hjust for p label
+    if (!is.null(p_x)) {
+      label_params$mapping$y <- Inf
       label_params$vjust     <- label_params$vjust %||% 1.5
+
+      p_just <- ifelse(all_sym, p$p_lab, paste0("p = ", p$p_lab))
+
+      p <- dplyr::mutate(
+        p,
+        hjust = case_when(
+          orig.ident == Inf  ~ .get_label_just(p_just),
+          orig.ident == -Inf ~ .get_label_just(p_just, side = "left"),
+          TRUE               ~ 0.5
+        )
+      )
+
+      label_params$mapping$hjust <- sym("hjust")
+
+      if (!all_sym) {
+        p <- dplyr::mutate(p, p_lab = paste0("italic(p) == ", .data$p_lab))
+      }
     }
 
     label_params$data   <- p
@@ -466,7 +490,7 @@ djvdj_theme <- function(base_size = 11, base_family = "",
     res <- res +
       lift(ggplot2::geom_text)(label_params)
 
-    if (p_corner && !"corner" %in% n_label) {
+    if (!is.null(p_x) && !"corner" %in% n_label) {
       res <- res +
         ggplot2::scale_y_continuous(expand = .n_label_expansion)
     }
diff --git a/man/plot_diversity.Rd b/man/plot_diversity.Rd
index 449f5302..0bea8657 100644
--- a/man/plot_diversity.Rd
+++ b/man/plot_diversity.Rd
@@ -20,6 +20,9 @@ plot_diversity(
   panel_nrow = NULL,
   panel_scales = "free",
   n_label = NULL,
+  p_label = "all",
+  p_method = NULL,
+  p_file = NULL,
   label_params = list(),
   ...
 )
@@ -79,6 +82,32 @@ plot legend
 \item 'none', do not display the number of cells plotted
 }}
 
+\item{p_label}{Specification indicating how p-values should be labeled on
+plot, this can one of the following:
+\itemize{
+\item 'none', do not display p-values
+\item 'all', show p-values for all groups
+\item A named vector providing p-value cutoffs and labels to display,
+e.g. \code{c('*' = 0.05, '**' = 0.01, '***' = 0.001)}. The keyword 'value' can
+be used to display the p-value for those less than a certain cutoff,
+e.g. \code{c(value = 0.05, ns = 1.1)} will show significant p-values, all others
+will be labeled 'ns'.
+}}
+
+\item{p_method}{Method to use for calculating p-values, by default when
+comparing two groups a t-test will be used.
+When comparing more than two groups the Kruskal-Wallis test will be used.
+p-values are adjusted for
+multiple testing using Bonferroni correction. Possible methods include:
+\itemize{
+\item 't', two sample t-test performed with \code{stats::t.test()}
+\item 'wilcox', Wilcoxon rank sum test performed with \code{stats::wilcox.test()}
+\item 'kruskal', Kruskal-Wallis test performed with \code{stats::kruskal.test()}
+}}
+
+\item{p_file}{File path to save table containing p-values for each
+comparison.}
+
 \item{label_params}{Named list providing additional parameters to modify
 n label aesthetics, e.g. list(size = 4, color = "red")}