diff --git a/DESCRIPTION b/DESCRIPTION
index 0a26fcb..12fddc0 100755
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -74,6 +74,7 @@ Imports:
     Seurat (>= 4.1.1),
     gtable (>= 0.3.1),
     digest (>= 0.6.29),
+    EnhancedVolcano,
     png,
     ggExtra,
     httr,
diff --git a/NAMESPACE b/NAMESPACE
index 96b1ac4..457d856 100755
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,17 +1,21 @@
 # Generated by roxygen2: do not edit by hand
 
+export(Volcano_Plot)
 export(aggregateCounts)
 export(annotateCellTypes)
 export(appendMetadataToSeuratObject)
+export(build_modscore_plots)
 export(colorByGene)
 export(colorByMarkerTable)
 export(combineNormalize)
+export(compute_modscore_data)
 export(degGeneExpressionMarkers)
 export(dotPlotMet)
 export(dualLabeling)
 export(filterQC)
 export(filterSeuratObjectByMetadata)
 export(heatmapSC)
+export(launch_module_score_app)
 export(modScore)
 export(nameClusters)
 export(object)
@@ -46,11 +50,13 @@ import(quantmod)
 import(reshape2)
 import(rlang)
 import(scales)
+import(tidyr)
 import(tidyverse)
 import(tools)
 import(utils)
 importFrom(BiocParallel,SerialParam)
 importFrom(ComplexHeatmap,pheatmap)
+importFrom(EnhancedVolcano,EnhancedVolcano)
 importFrom(RColorBrewer,brewer.pal)
 importFrom(RColorBrewer,brewer.pal.info)
 importFrom(Seurat,AddMetaData)
@@ -76,6 +82,7 @@ importFrom(dplyr,arrange)
 importFrom(dplyr,case_when)
 importFrom(dplyr,desc)
 importFrom(dplyr,filter)
+importFrom(dplyr,group_by)
 importFrom(dplyr,if_else)
 importFrom(dplyr,mutate)
 importFrom(dplyr,mutate_if)
@@ -85,9 +92,13 @@ importFrom(dplyr,rename)
 importFrom(dplyr,row_number)
 importFrom(dplyr,select)
 importFrom(dplyr,summarise)
+importFrom(geom_point,guide_legend)
+importFrom(geom_point,guides)
+importFrom(geom_point,scale_color_gradientn)
 importFrom(ggExtra,ggMarginal)
 importFrom(ggplot2,aes)
 importFrom(ggplot2,coord_fixed)
+importFrom(ggplot2,element_blank)
 importFrom(ggplot2,geom_hline)
 importFrom(ggplot2,geom_point)
 importFrom(ggplot2,geom_vline)
@@ -104,6 +115,7 @@ importFrom(ggplot2,ylim)
 importFrom(ggpubr,annotate_figure)
 importFrom(ggpubr,get_legend)
 importFrom(ggpubr,ggarrange)
+importFrom(ggrepel,geom_text_repel)
 importFrom(grDevices,colorRampPalette)
 importFrom(grid,gTree)
 importFrom(grid,gpar)
@@ -113,7 +125,10 @@ importFrom(grid,grobHeight)
 importFrom(grid,textGrob)
 importFrom(grid,unit)
 importFrom(gridExtra,arrangeGrob)
+importFrom(gridExtra,gpar)
+importFrom(gridExtra,grid)
 importFrom(gridExtra,tableGrob)
+importFrom(gridExtra,textGrob)
 importFrom(htmlwidgets,saveWidget)
 importFrom(magrittr,"%>%")
 importFrom(plotly,as_widget)
@@ -121,6 +136,10 @@ importFrom(plotly,ggplotly)
 importFrom(plotly,plot_ly)
 importFrom(reshape2,melt)
 importFrom(scDblFinder,scDblFinder)
+importFrom(scale_y_continuous,geom_line)
+importFrom(scale_y_continuous,geom_segment)
+importFrom(scale_y_continuous,scale_x_continuous)
+importFrom(scale_y_continuous,scale_y_log10)
 importFrom(scales,rescale)
 importFrom(stats,as.hclust)
 importFrom(stats,hclust)
@@ -130,6 +149,7 @@ importFrom(stats,median)
 importFrom(stats,quantile)
 importFrom(stringr,str_replace_all)
 importFrom(stringr,str_sort)
+importFrom(stringr,str_split)
 importFrom(stringr,str_split_fixed)
 importFrom(stringr,str_to_title)
 importFrom(stringr,str_wrap)
@@ -138,3 +158,8 @@ importFrom(tibble,rownames_to_column)
 importFrom(tibble,tibble)
 importFrom(tidyr,fill)
 importFrom(tidyr,pivot_wider)
+importFrom(xlab,element_text)
+importFrom(xlab,geom_hline)
+importFrom(xlab,geom_violin)
+importFrom(xlab,theme_classic)
+importFrom(xlab,ylab)
diff --git a/R/Volcano_Plot.R b/R/Volcano_Plot.R
new file mode 100644
index 0000000..12b7138
--- /dev/null
+++ b/R/Volcano_Plot.R
@@ -0,0 +1,314 @@
+#' Create Volcano Plots for Differential Expression Analysis
+#'
+#' Creates enhanced volcano plots for differential expression analysis using the
+#' EnhancedVolcano package. Generates both static and interactive (plotly) versions
+#' for each comparison specified in the input data.
+#'
+#' @param DEGAnalysis Data frame containing differential expression analysis results.
+#'   Must include columns for gene/feature IDs, log fold changes, and p-values.
+#' @param label.col Character string specifying the column name containing feature
+#'   IDs (e.g., gene names). Default is "Gene".
+#' @param sig.col Character vector of column names containing significance values
+#'   (p-values or adjusted p-values). Can specify multiple comparisons.
+#' @param lfc.col Character vector of column names containing log2 fold change values.
+#'   Must have the same length as \code{sig.col}.
+#' @param pCutoff Numeric value for the significance threshold. Default is 0.001.
+#' @param FCcutoff Numeric value for the log2 fold change threshold. Default is 1.0.
+#' @param value_to_sort_the_output_dataset Character string specifying how to sort
+#'   features for labeling. Options are "p-value" or "fold-change". Default is "p-value".
+#' @param no_genes_to_label Integer specifying the number of top features to label
+#'   on the plot. Default is 30.
+#' @param use_only_addition_labels Logical indicating whether to label only the
+#'   features specified in \code{additional_labels}. Default is FALSE.
+#' @param additional_labels Character string of additional feature names to label,
+#'   separated by commas or spaces. Default is "".
+#' @param is_red Logical indicating whether to restrict labeled features to only
+#'   those passing both significance and fold change thresholds. Default is TRUE.
+#' @param labSize Numeric value for the size of feature labels. Default is 4.
+#' @param change_sig_name Character string for custom y-axis label. Default is "p-value".
+#' @param change_lfc_name Character string for custom x-axis label. Default is "log2FC".
+#' @param title Character string for the plot title. Default is "Volcano Plots".
+#' @param use_custom_lab Logical indicating whether to use custom axis labels.
+#'   Default is FALSE.
+#' @param ylim Numeric value for the maximum y-axis limit. Set to 0 for automatic
+#'   scaling. Default is 0.
+#' @param custom_xlim Character string for custom x-axis limits. Can be empty for
+#'   automatic scaling, a single number for symmetrical limits, or two comma-separated
+#'   numbers for asymmetrical limits. Default is "".
+#' @param xlim_additional Numeric value to add padding to x-axis limits. Default is 0.
+#' @param ylim_additional Numeric value to add padding to y-axis limits. Default is 0.
+#' @param axisLabSize Numeric value for the size of axis labels. Default is 24.
+#' @param pointSize Numeric value for the size of points on the plot. Default is 2.
+#'
+#' @return A named list containing:
+#'   \item{data}{Original input data frame with added rank columns for each comparison}
+#'   \item{comparison_static}{Static ggplot2 volcano plot for each comparison}
+#'   \item{comparison_interactive}{Interactive plotly volcano plot for each comparison}
+#'
+#' @examples
+#' \dontrun{
+#' # Create sample DEG data
+#' set.seed(123)
+#' deg_data <- data.frame(
+#'   Gene = paste0("Gene", 1:100),
+#'   `WT-KO_logFC` = rnorm(100, 0, 2),
+#'   `WT-KO_pval` = runif(100, 0, 0.1),
+#'   check.names = FALSE
+#' )
+#'
+#' # Generate volcano plots
+#' results <- Volcano_Plot(
+#'   DEGAnalysis = deg_data,
+#'   label.col = "Gene",
+#'   sig.col = "WT-KO_pval",
+#'   lfc.col = "WT-KO_logFC",
+#'   pCutoff = 0.05,
+#'   FCcutoff = 1.0
+#' )
+#'
+#' # Access the plots
+#' print(results$`WT-KO_static`)
+#' print(results$`WT-KO_interactive`)
+#' }
+#'
+#' @import ggplot2
+#' @import dplyr
+#' @import tidyr
+#' @importFrom stringr str_split
+#' @importFrom ggrepel geom_text_repel
+#' @importFrom EnhancedVolcano EnhancedVolcano
+#' @importFrom plotly ggplotly
+#' @importFrom grid grid.newpage
+#' @export
+Volcano_Plot <- function(DEGAnalysis,
+                         label.col = "Gene",
+                         sig.col,
+                         lfc.col,
+                         pCutoff = 0.001,
+                         FCcutoff = 1.0,
+                         value_to_sort_the_output_dataset = "p-value",
+                         no_genes_to_label = 30,
+                         use_only_addition_labels = FALSE,
+                         additional_labels = "",
+                         is_red = TRUE,
+                         labSize = 4,
+                         change_sig_name = "p-value",
+                         change_lfc_name = "log2FC",
+                         title = "Volcano Plots",
+                         use_custom_lab = FALSE,
+                         ylim = 0,
+                         custom_xlim = "",
+                         xlim_additional = 0,
+                         ylim_additional = 0,
+                         axisLabSize = 24,
+                         pointSize = 2) {
+
+  ## --------------- ##
+  ## Main Code Block ##
+  ## --------------- ##
+
+  df.orig <- DEGAnalysis
+  rank <- list()
+  plot_list <- list()
+
+  for(i in 1:length(lfc.col)){
+    lfccol <- lfc.col[i]
+    sigcol <- sig.col[i]
+    columns_of_interest <- c(label.col, lfc.col[i], sig.col[i])
+    df <- df.orig %>% 
+      dplyr::select(one_of(columns_of_interest)) %>% 
+      mutate(!!sym(lfccol) := replace_na(!!sym(lfccol), 0)) %>%
+      mutate(!!sym(sigcol) := replace_na(!!sym(sigcol), 1))
+
+    if (use_custom_lab == TRUE){
+      if (nchar(change_lfc_name) == 0){lfc_name = lfc.col[i]}
+      if (nchar(change_sig_name) == 0){sig_name = sig.col[i]}
+      colnames(df) <- c(label.col, change_lfc_name, sig_name)
+    } else {
+      lfc_name = lfc.col[i]
+      sig_name = sig.col[i]
+    }
+
+    group <- gsub("_pval|p_val_", "", sig_name)
+    rank[[i]] <- -log10(df[[sig_name]]) * sign(df[[lfc_name]]) 
+    names(rank)[i] <- paste0("C_", group, "_rank")
+
+    cat(paste0("Genes in initial dataset: ", nrow(df), "\n"))
+
+    # Select top genes by logFC or Significance
+    if (value_to_sort_the_output_dataset == "fold-change") {
+      df <- df %>% dplyr::arrange(desc(.data[[lfc_name]]))
+    } else if (value_to_sort_the_output_dataset == "p-value") {
+      df <- df %>% dplyr::arrange(.data[[sig_name]]) 
+    }
+
+    if (is_red) {
+      df_sub <- df[df[[sig_name]] <= pCutoff & abs(df[[lfc_name]]) >= FCcutoff, ]
+    } else {
+      df_sub <- df
+    }
+
+    genes_to_label <- as.character(df_sub[1:no_genes_to_label, label.col])
+
+    # Modifying Additional Labels List:
+    # Replace commas with spaces and split the string
+    split_values <- unlist(strsplit(gsub(",", " ", additional_labels), " "))
+    additional_labels_vec <- split_values[split_values != ""]
+
+    filter <- additional_labels_vec %in% df[, label.col]
+    missing_labels <- additional_labels_vec[!filter]
+    additional_labels_vec <- additional_labels_vec[filter]
+
+    if(length(missing_labels) > 0){
+      cat("Could not find:\n")
+      print(missing_labels)
+    }
+
+    if(use_only_addition_labels){
+      genes_to_label <- additional_labels_vec
+    } else {
+      genes_to_label <- unique(append(genes_to_label, additional_labels_vec))
+    }
+
+    significant = vector(length = nrow(df))
+    significant[] = "Not significant"
+    significant[which(abs(df[, 2]) > FCcutoff)] = "Fold change only"
+    significant[which(df[, 3] < pCutoff)] = "Significant only"
+    significant[which(abs(df[, 2]) > FCcutoff & df[, 3] < pCutoff)] = "Significant and fold change"
+    print(table(significant))
+
+    # Fix pvalue == 0
+    shapeCustom <- rep(19, nrow(df))
+    maxy <- max(-log10(df[[sig_name]]), na.rm = TRUE)
+    if(ylim > 0){
+      maxy <- ylim
+    }
+
+    cat(paste0("Maxy: ", maxy, "\n"))
+    if(maxy == Inf){
+      # Sometimes, pvalues == 0
+      keep <- df[[sig_name]] > 0
+      df[[sig_name]][!keep] <- min(df[[sig_name]][keep])
+      shapeCustom[!keep] <- 17
+
+      maxy <- -log10(min(df[[sig_name]][keep]))
+      cat("Some p-values equal zero. Adjusting y-limits.\n")
+      cat(paste0("Maxy adjusted: ", maxy, "\n"))
+    }
+
+    # By default, nothing will be greater than maxy. User can set this value lower
+    keep <- -log10(df[[sig_name]]) <= maxy
+    df[[sig_name]][!keep] <- maxy
+    shapeCustom[!keep] <- 17
+
+    names(shapeCustom) <- rep("Exact", length(shapeCustom))
+    names(shapeCustom)[shapeCustom == 17] <- "Adjusted"
+
+    # Remove if nothin' doin'
+    if(all(shapeCustom == 19)){
+      shapeCustom <- NULL
+    }
+
+    maxy <- ceiling(maxy)
+
+    if (grepl("log", lfc.col[i])){
+      xlab <- bquote(~Log[2]~ "fold change")
+    } else {
+      xlab <- "Fold change"
+    }
+    if (grepl("adj", sig.col[i])){
+      ylab <- bquote(~-Log[10]~ "FDR")
+    } else {
+      ylab <- bquote(~-Log[10]~ "p-value")
+    }
+    if(use_custom_lab){
+      if(lfc_name != lfc.col[i]){
+        xlab <- gsub("_", " ", lfc_name)
+      }
+      if (sig_name != sig.col[i]){ 
+        ylab <- gsub("_", " ", sig_name)
+      }
+    }
+
+    # X-axis custom range change:
+    if (custom_xlim == "") {
+      xlim = c(floor(min(df[, lfc_name])) - xlim_additional, 
+               ceiling(max(df[, lfc_name])) + xlim_additional)
+    } else if (grepl(",", custom_xlim) == FALSE) {
+      xlim = c(-1 * as.numeric(trimws(custom_xlim)), 
+               as.numeric(trimws(custom_xlim)))
+    } else {
+      split_values <- strsplit(custom_xlim, ",")[[1]]
+      x_min <- as.numeric(trimws(split_values[1]))
+      x_max <- as.numeric(trimws(split_values[2]))
+      xlim <- c(x_min, x_max)
+    }
+
+    # Create static plot
+    p <- EnhancedVolcano(df, x = lfc_name, y = sig_name,
+                         lab = df[, label.col],
+                         selectLab = genes_to_label,
+                         title = title,
+                         subtitle = group,
+                         xlab = xlab,
+                         ylab = ylab,
+                         xlim = xlim,
+                         ylim = c(0, maxy + ylim_additional),
+                         pCutoff = pCutoff,
+                         FCcutoff = FCcutoff,
+                         axisLabSize = axisLabSize,
+                         labSize = labSize,
+                         pointSize = pointSize,
+                         shapeCustom = shapeCustom)
+
+    # Create interactive plot with no labels
+    p_empty <- EnhancedVolcano(df, x = lfc_name, y = sig_name,
+                               lab = rep("", nrow(df)),
+                               selectLab = NULL,
+                               title = title,
+                               subtitle = group,
+                               xlab = xlab,
+                               ylab = ylab,
+                               xlim = xlim,
+                               ylim = c(0, maxy + ylim_additional),
+                               pCutoff = pCutoff,
+                               FCcutoff = FCcutoff,
+                               axisLabSize = axisLabSize,
+                               labSize = labSize,
+                               pointSize = pointSize,
+                               shapeCustom = shapeCustom)
+
+    # Extract the data used for plotting
+    plot_data <- ggplot_build(p_empty)$data[[1]]
+
+    pxx <- p_empty +
+      xlab("Fold Change") +
+      ylab("Significance") +
+      theme_minimal() +
+      geom_point(aes(
+        text = paste("Gene:", df[[label.col]], 
+                     "<br>Log2FC:", df[[lfc_name]], 
+                     "<br>P-value:", df[[sig_name]]),
+        colour = as.character(plot_data$colour),
+        fill = as.character(plot_data$colour)
+      ), 
+      shape = 21,
+      size = 2,
+      stroke = 0.1) + 
+      scale_fill_identity()
+
+    # Add interactive hover labels for the gene names
+    interactive_plot <- ggplotly(pxx, tooltip = c("text"))
+
+    # Store plots in the list with descriptive names
+    plot_list[[paste0(group, "_static")]] <- p
+    plot_list[[paste0(group, "_interactive")]] <- interactive_plot
+  }
+
+  # Combine original data with rank columns
+  df.final <- cbind(df.orig, do.call(cbind, rank))
+
+  # Return list with data first, then all plots
+  result <- c(list(data = df.final), plot_list)
+  return(result)
+}
\ No newline at end of file
diff --git a/man/Volcano_Plot.Rd b/man/Volcano_Plot.Rd
new file mode 100644
index 0000000..889f40c
--- /dev/null
+++ b/man/Volcano_Plot.Rd
@@ -0,0 +1,127 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Volcano_Plot.R
+\name{Volcano_Plot}
+\alias{Volcano_Plot}
+\title{Create Volcano Plots for Differential Expression Analysis}
+\usage{
+Volcano_Plot(
+  DEGAnalysis,
+  label.col = "Gene",
+  sig.col,
+  lfc.col,
+  pCutoff = 0.001,
+  FCcutoff = 1,
+  value_to_sort_the_output_dataset = "p-value",
+  no_genes_to_label = 30,
+  use_only_addition_labels = FALSE,
+  additional_labels = "",
+  is_red = TRUE,
+  labSize = 4,
+  change_sig_name = "p-value",
+  change_lfc_name = "log2FC",
+  title = "Volcano Plots",
+  use_custom_lab = FALSE,
+  ylim = 0,
+  custom_xlim = "",
+  xlim_additional = 0,
+  ylim_additional = 0,
+  axisLabSize = 24,
+  pointSize = 2
+)
+}
+\arguments{
+\item{DEGAnalysis}{Data frame containing differential expression analysis results.
+Must include columns for gene/feature IDs, log fold changes, and p-values.}
+
+\item{label.col}{Character string specifying the column name containing feature
+IDs (e.g., gene names). Default is "Gene".}
+
+\item{sig.col}{Character vector of column names containing significance values
+(p-values or adjusted p-values). Can specify multiple comparisons.}
+
+\item{lfc.col}{Character vector of column names containing log2 fold change values.
+Must have the same length as \code{sig.col}.}
+
+\item{pCutoff}{Numeric value for the significance threshold. Default is 0.001.}
+
+\item{FCcutoff}{Numeric value for the log2 fold change threshold. Default is 1.0.}
+
+\item{value_to_sort_the_output_dataset}{Character string specifying how to sort
+features for labeling. Options are "p-value" or "fold-change". Default is "p-value".}
+
+\item{no_genes_to_label}{Integer specifying the number of top features to label
+on the plot. Default is 30.}
+
+\item{use_only_addition_labels}{Logical indicating whether to label only the
+features specified in \code{additional_labels}. Default is FALSE.}
+
+\item{additional_labels}{Character string of additional feature names to label,
+separated by commas or spaces. Default is "".}
+
+\item{is_red}{Logical indicating whether to restrict labeled features to only
+those passing both significance and fold change thresholds. Default is TRUE.}
+
+\item{labSize}{Numeric value for the size of feature labels. Default is 4.}
+
+\item{change_sig_name}{Character string for custom y-axis label. Default is "p-value".}
+
+\item{change_lfc_name}{Character string for custom x-axis label. Default is "log2FC".}
+
+\item{title}{Character string for the plot title. Default is "Volcano Plots".}
+
+\item{use_custom_lab}{Logical indicating whether to use custom axis labels.
+Default is FALSE.}
+
+\item{ylim}{Numeric value for the maximum y-axis limit. Set to 0 for automatic
+scaling. Default is 0.}
+
+\item{custom_xlim}{Character string for custom x-axis limits. Can be empty for
+automatic scaling, a single number for symmetrical limits, or two comma-separated
+numbers for asymmetrical limits. Default is "".}
+
+\item{xlim_additional}{Numeric value to add padding to x-axis limits. Default is 0.}
+
+\item{ylim_additional}{Numeric value to add padding to y-axis limits. Default is 0.}
+
+\item{axisLabSize}{Numeric value for the size of axis labels. Default is 24.}
+
+\item{pointSize}{Numeric value for the size of points on the plot. Default is 2.}
+}
+\value{
+A named list containing:
+\item{data}{Original input data frame with added rank columns for each comparison}
+\item{comparison_static}{Static ggplot2 volcano plot for each comparison}
+\item{comparison_interactive}{Interactive plotly volcano plot for each comparison}
+}
+\description{
+Creates enhanced volcano plots for differential expression analysis using the
+EnhancedVolcano package. Generates both static and interactive (plotly) versions
+for each comparison specified in the input data.
+}
+\examples{
+\dontrun{
+# Create sample DEG data
+set.seed(123)
+deg_data <- data.frame(
+  Gene = paste0("Gene", 1:100),
+  `WT-KO_logFC` = rnorm(100, 0, 2),
+  `WT-KO_pval` = runif(100, 0, 0.1),
+  check.names = FALSE
+)
+
+# Generate volcano plots
+results <- Volcano_Plot(
+  DEGAnalysis = deg_data,
+  label.col = "Gene",
+  sig.col = "WT-KO_pval",
+  lfc.col = "WT-KO_logFC",
+  pCutoff = 0.05,
+  FCcutoff = 1.0
+)
+
+# Access the plots
+print(results$`WT-KO_static`)
+print(results$`WT-KO_interactive`)
+}
+
+}
diff --git a/tests/testthat/test-Volcano_Plot.R b/tests/testthat/test-Volcano_Plot.R
new file mode 100644
index 0000000..b5b98cd
--- /dev/null
+++ b/tests/testthat/test-Volcano_Plot.R
@@ -0,0 +1,139 @@
+# Tests for Volcano_Plot function
+
+library(testthat)
+library(dplyr)
+
+# Helper function to create test data
+create_test_deg_data <- function(n_genes = 100) {
+  set.seed(123)
+  data.frame(
+    Gene = paste0("Gene", 1:n_genes),
+    `WT-KO_logFC` = rnorm(n_genes, 0, 2),
+    `WT-KO_pval` = runif(n_genes, 0, 0.1),
+    `Treated-Control_logFC` = rnorm(n_genes, 0, 1.5),
+    `Treated-Control_pval` = runif(n_genes, 0, 0.05),
+    check.names = FALSE
+  )
+}
+
+test_that("Volcano_Plot runs without errors on valid input", {
+  deg_data <- create_test_deg_data()
+  
+  expect_no_error(
+    result <- Volcano_Plot(
+      DEGAnalysis = deg_data,
+      label.col = "Gene",
+      sig.col = "WT-KO_pval",
+      lfc.col = "WT-KO_logFC",
+      pCutoff = 0.05,
+      FCcutoff = 1.0
+    )
+  )
+})
+
+test_that("Volcano_Plot returns correct structure", {
+  deg_data <- create_test_deg_data()
+  
+  result <- Volcano_Plot(
+    DEGAnalysis = deg_data,
+    label.col = "Gene",
+    sig.col = "WT-KO_pval",
+    lfc.col = "WT-KO_logFC",
+    pCutoff = 0.05,
+    FCcutoff = 1.0
+  )
+  
+  # Check that result is a list
+  expect_type(result, "list")
+  
+  # Check that 'data' element exists
+  expect_true("data" %in% names(result))
+  
+  # Check that data is a data frame
+  expect_s3_class(result$data, "data.frame")
+  
+  # Check that rank column was added
+  expect_true(any(grepl("_rank$", colnames(result$data))))
+  
+  # Check that static plot exists
+  expect_true(any(grepl("_static$", names(result))))
+  
+  # Check that interactive plot exists
+  expect_true(any(grepl("_interactive$", names(result))))
+  
+  # Check that static plot is a ggplot object
+  static_plot_name <- names(result)[grepl("_static$", names(result))][1]
+  expect_s3_class(result[[static_plot_name]], "ggplot")
+  
+  # Check that interactive plot is a plotly object
+  interactive_plot_name <- names(result)[grepl("_interactive$", names(result))][1]
+  expect_s3_class(result[[interactive_plot_name]], "plotly")
+})
+
+test_that("Volcano_Plot handles multiple comparisons correctly", {
+  deg_data <- create_test_deg_data()
+  
+  result <- Volcano_Plot(
+    DEGAnalysis = deg_data,
+    label.col = "Gene",
+    sig.col = c("WT-KO_pval", "Treated-Control_pval"),
+    lfc.col = c("WT-KO_logFC", "Treated-Control_logFC"),
+    pCutoff = 0.05,
+    FCcutoff = 1.0
+  )
+  
+  # Should have data + 4 plots (2 comparisons × 2 plot types)
+  expect_equal(length(result), 5)
+  
+  # Check that both comparison names appear
+  expect_true(any(grepl("WT-KO", names(result))))
+  expect_true(any(grepl("Treated-Control", names(result))))
+  
+  # Check that we have both static and interactive for each comparison
+  expect_equal(sum(grepl("_static$", names(result))), 2)
+  expect_equal(sum(grepl("_interactive$", names(result))), 2)
+})
+
+test_that("Volcano_Plot handles edge cases", {
+  # Test with p-values = 0
+  deg_data <- create_test_deg_data(n_genes = 50)
+  deg_data$`WT-KO_pval`[1:5] <- 0  # Set some p-values to exactly 0
+  
+  expect_no_error(
+    result <- Volcano_Plot(
+      DEGAnalysis = deg_data,
+      label.col = "Gene",
+      sig.col = "WT-KO_pval",
+      lfc.col = "WT-KO_logFC",
+      pCutoff = 0.05,
+      FCcutoff = 1.0
+    )
+  )
+  
+  # Check that the function completed and returned data
+  expect_true("data" %in% names(result))
+  expect_equal(nrow(result$data), 50)
+  
+  # Test with NA values
+  deg_data_na <- create_test_deg_data(n_genes = 50)
+  deg_data_na$`WT-KO_logFC`[1:3] <- NA
+  deg_data_na$`WT-KO_pval`[4:6] <- NA
+  
+  expect_no_error(
+    result_na <- Volcano_Plot(
+      DEGAnalysis = deg_data_na,
+      label.col = "Gene",
+      sig.col = "WT-KO_pval",
+      lfc.col = "WT-KO_logFC",
+      pCutoff = 0.05,
+      FCcutoff = 1.0
+    )
+  )
+  
+  # Check that the function completed and returned valid data
+  # (NA values are handled internally for plotting but preserved in returned data)
+  expect_true("data" %in% names(result_na))
+  expect_equal(nrow(result_na$data), 50)
+  expect_s3_class(result_na$`WT-KO_static`, "ggplot")
+  expect_s3_class(result_na$`WT-KO_interactive`, "plotly")
+})
\ No newline at end of file