diff --git a/GGPLOT2_V4_KNOWN_ISSUES.md b/GGPLOT2_V4_KNOWN_ISSUES.md
new file mode 100644
index 000000000..42f2297ac
--- /dev/null
+++ b/GGPLOT2_V4_KNOWN_ISSUES.md
@@ -0,0 +1,90 @@
+# ggplot2 v4.0.0 Known Issues in Seurat
+
+## Overview
+ggplot2 version 4.0.0 (released September 2025) introduced breaking changes by migrating from S3 to S7 object system. This document tracks known compatibility issues and their status.
+
+## Fixed Issues
+
+### 1. DarkTheme() - Deprecated `size` parameter ✅
+**Status**: FIXED  
+**Files**: `R/visualization.R` (lines 6476-6491)  
+**Changes**: 
+- Replaced `element_rect(fill = 'black', size = 0)` with `linewidth = 0`
+- Replaced `element_line(colour = 'white', size = 1)` with `linewidth = 1`
+- Replaced `element_line(size = 0)` with `linewidth = 0`
+
+### 2. LabelClusters() - S7 object color retrieval ✅
+**Status**: FIXED  
+**Issues**: #10156, #10176  
+**Files**: `R/visualization.R` (lines 6134-6227)  
+**Problem**: With S7 objects, `pb$data` returns numeric coordinates instead of valid colors  
+**Solution**: Extract colors directly from plot scales using `pb$plot$scales$get_scales()`
+**Changes**:
+- Added color scale extraction from plot scales
+- Created group_colors mapping with validation
+- Fixed data.medians$color assignment to use explicit group_colors mapping
+
+## Outstanding Issues (Require ggplot2 Package Fix)
+
+### 3. VlnPlot - S4SXP Error 🔴
+**Status**: NEEDS GGPLOT2 FIX  
+**Issues**: #10101, #10188, #10160  
+**Files**: Not Seurat-specific code  
+**Error**: 
+```
+Error in deparse(substitute(e2, env = caller_env(2))) :
+'S4SXP': should not happen - please report
+```
+**Root Cause**: Deep incompatibility between ggplot2 v4 S7 objects and R's `substitute()` function calls within ggplot2 internals  
+**Workaround**: Downgrade to ggplot2 3.5.2
+```r
+remotes::install_version("ggplot2", version = "3.5.2", repos = "https://cran.r-project.org")
+```
+**Note**: Some packages (e.g., clusterProfiler >= 4.12.0) require ggplot2 >= 4.0.0, creating dependency conflicts
+
+## Potential Issues (To Monitor)
+
+### 4. ggplot_build()$plot$data Access Pattern ⚠️
+**Status**: MONITORING  
+**Files**: `R/visualization.R` (lines 5948, 7795)  
+**Functions**: `GGpointToPlotly()`, `GGpointToPlotlyBuild()`  
+**Concern**: Direct `$plot$data` access may behave differently with S7 objects  
+**Impact**: Unknown - needs testing with real-world usage
+
+### 5. plot$layers Access Pattern ⚠️
+**Status**: MONITORING  
+**Files**: `R/visualization.R` (multiple locations)  
+**Concern**: `plot$layers[[i]]$data` and `plot$layers[[i]]$mapping` access may need S7 updates  
+**Impact**: Unknown - backward compatibility should handle this, but monitor for issues
+
+## Migration Notes
+
+### What Changed in ggplot2 v4
+- **S7 Object System**: ggplot objects now use S7 instead of S3
+- **Deprecated Parameters**: 
+  - `size` → `linewidth` in `element_rect()` and `element_line()`
+  - Keep `size` in `element_text()` (NOT deprecated)
+- **S7 Backward Compatibility**: `$data` and `$layers` access still works, but internal data structures changed
+- **New Features**: 
+  - Theme improvements (ink/paper/accent, palette.*, theme_sub_*)
+  - `geom_label()` new aesthetics (linewidth, border.colour, text.colour)
+
+### Testing Recommendations
+1. Test all visualization functions with ggplot2 v4
+2. Verify label.box functionality in DimPlot variants
+3. Check VlnPlot alternatives if S4SXP persists
+4. Monitor Plotly integration functions for S7 issues
+
+### Dependencies
+- **Minimum ggplot2 version**: 3.5.2 (for compatibility)
+- **Target ggplot2 version**: 4.0.0+ (once S4SXP issue resolved upstream)
+
+## References
+- ggplot2 v4 blog post: https://www.tidyverse.org/blog/2025/09/ggplot2-4-0-0/
+- Issue #10101: VlnPlot S4SXP error
+- Issue #10156: DimPlot label.box color error
+- Issue #10176: label.box returning numeric values
+- Issue #10188: Temporary S4SXP solution thread
+
+## Last Updated
+2025-11-27
diff --git a/R/integration.R b/R/integration.R
index 55569c943..3baca56ca 100644
--- a/R/integration.R
+++ b/R/integration.R
@@ -823,10 +823,27 @@ FindTransferAnchors <- function(
           features = features,
           verbose = FALSE
           ))
+        # Check if scale.data exists and has features
+        if (is.null(reference[[reference.assay]]$scale.data) || 
+            nrow(reference[[reference.assay]]$scale.data) == 0) {
+          stop(
+            "No scale.data found in the ", reference.assay, " assay of the reference object. ",
+            "For SCT normalization, please run GetResidual() on the reference object before using FindTransferAnchors. ",
+            "ScaleData() is not required for SCT assays.",
+            call. = FALSE
+          )
+        }
         features <- intersect(
           x = features,
           y = rownames(reference[[reference.assay]]$scale.data)
         )
+        if (length(features) == 0) {
+          stop(
+            "No features remaining after intersecting with scale.data in ", reference.assay, " assay. ",
+            "For SCT normalization, please ensure GetResidual() has been run on the reference object with appropriate features.",
+            call. = FALSE
+          )
+        }
         VariableFeatures(reference) <- features
       }
       if (IsSCT(assay = query[[query.assay]])) {
@@ -3617,7 +3634,33 @@ TransferData <- function(
          prediction.scores <- (prediction.scores + bridge.prediction.scores)/2
          prediction.scores <- as.matrix(x = prediction.scores)
       }
-      prediction.ids <- possible.ids[apply(X = prediction.scores, MARGIN = 1, FUN = which.max)]
+      
+      # Check for NaN values in prediction scores and handle them
+      if (any(is.nan(prediction.scores))) {
+        warning(
+          "NaN values detected in prediction scores. This may indicate issues with ",
+          "normalization, k.weight parameter, or data quality. NaN scores will be replaced with 0.",
+          call. = FALSE,
+          immediate. = TRUE
+        )
+        prediction.scores[is.nan(prediction.scores)] <- 0
+      }
+      
+      # Use a safer version of which.max that handles edge cases
+      # A safer version of which.max that handles the all-NA case.
+      # Negative and zero values are handled by which.max as usual.
+      safe_which_max <- function(x) {
+        if (all(is.na(x))) {
+          return(1L)  # Default to first class if all are NA
+        }
+        result <- which.max(x)
+        if (length(result) == 0) {
+          return(1L)  # Default to first class if which.max returns empty
+        }
+        return(result)
+      }
+      
+      prediction.ids <- possible.ids[apply(X = prediction.scores, MARGIN = 1, FUN = safe_which_max)]
       prediction.ids <- as.character(prediction.ids)
       prediction.max <- apply(X = prediction.scores, MARGIN = 1, FUN = max)
       if (is.null(x = query)) {
diff --git a/R/preprocessing.R b/R/preprocessing.R
index fa95b7e2e..0f1937357 100644
--- a/R/preprocessing.R
+++ b/R/preprocessing.R
@@ -2846,7 +2846,8 @@ ReadXenium <- function(
         col.use = c(
           x_location = letters[24+flip.xy],
           y_location = letters[25-flip.xy],
-          feature_name = 'gene'
+          feature_name = 'gene',
+          qv = 'qv'
         )
 
         for(option in Filter(function(x) x$req, list(
@@ -2879,6 +2880,13 @@ ReadXenium <- function(
         colnames(transcripts) <- col.use
 
         transcripts$gene <- binary_to_string(transcripts$gene)
+        
+        # Apply QV threshold filtering
+        if (!is.null(mols.qv.threshold) && 'qv' %in% colnames(transcripts)) {
+          transcripts <- transcripts[transcripts$qv >= mols.qv.threshold, ]
+          # Remove qv column after filtering (not needed in output)
+          transcripts <- transcripts[, setdiff(colnames(transcripts), 'qv')]
+        }
 
         pmicrons(type = 'finish')
 
diff --git a/R/utilities.R b/R/utilities.R
index 59668765d..39fb1d785 100644
--- a/R/utilities.R
+++ b/R/utilities.R
@@ -2625,12 +2625,21 @@ Online <- function(url, strict = FALSE, seconds = 5L) {
 # Parenting(parent.find = 'Seurat', features = features[features > 7])
 #
 Parenting <- function(parent.find = 'Seurat', ...) {
-  calls <- as.character(x = sys.calls())
-  calls <- lapply(
-    X = strsplit(x = calls, split = '(', fixed = TRUE),
-    FUN = '[',
-    1
-  )
+  # Extract function names from call stack without serializing entire objects
+  # This avoids the performance issue when large objects are passed via do.call
+  calls <- vapply(sys.calls(), function(call) {
+    if (is.call(call) && length(call) > 0) {
+      # Get just the function name, not the full call with arguments
+      func <- call[[1]]
+      if (is.name(func)) {
+        return(as.character(func))
+      } else if (is.call(func)) {
+        # Handle cases like pkg::function
+        return(paste(as.character(func), collapse = "::"))
+      }
+    }
+    return("")
+  }, character(1))
   parent.index <- grep(pattern = parent.find, x = calls)
   if (length(x = parent.index) != 1) {
     warning(
diff --git a/R/visualization.R b/R/visualization.R
index 52033b3a1..a70e9e7e5 100644
--- a/R/visualization.R
+++ b/R/visualization.R
@@ -2179,7 +2179,24 @@ VariableFeaturePlot <- function(
   )
   status.col <- colnames(hvf.info)[grepl("variable", colnames(hvf.info))][[1]]
   var.status <- c('no', 'yes')[unlist(hvf.info[[status.col]]) + 1]
-  if (colnames(x = hvf.info)[3] == 'dispersion.scaled') {
+  
+  # Handle cases where hvf.info may have different numbers of columns
+  if (ncol(hvf.info) < 3) {
+    warning(
+      "HVFInfo returned fewer than 3 columns. ",
+      "This may occur when using different normalization methods on the same object. ",
+      "Using available columns for plotting.",
+      call. = FALSE,
+      immediate. = TRUE
+    )
+    hvf.info <- hvf.info[, 1:min(2, ncol(hvf.info)), drop = FALSE]
+    if (ncol(hvf.info) < 2) {
+      stop(
+        "HVFInfo returned fewer than 2 columns. ",
+        "Cannot generate plot. Please check your normalization and feature selection methods."
+      )
+    }
+  } else if (colnames(x = hvf.info)[3] == 'dispersion.scaled') {
     hvf.info <- hvf.info[, c(1, 2)]
   } else if (colnames(x = hvf.info)[3] == 'variance.expected') {
     hvf.info <- hvf.info[, c(1, 4)]
@@ -6132,12 +6149,44 @@ LabelClusters <- function(
     }
   }
 
-  # Retrieve colour from built data
-  col_choice <- intersect(c("colour", "color"), names(pb$data[[1]]))
-  if (length(col_choice) > 0) {
-    data <- cbind(data, color = pb$data[[1]][[col_choice[1]]])
-  } else {
-    data <- cbind(data, color = NA_character_)
+  # Retrieve colour mapping for each group
+  # Extract colors from the plot's scales
+  color_scale <- NULL
+  if (!is.null(pb$plot$scales$get_scales("colour"))) {
+    color_scale <- pb$plot$scales$get_scales("colour")
+  } else if (!is.null(pb$plot$scales$get_scales("color"))) {
+    color_scale <- pb$plot$scales$get_scales("color")
+  }
+  
+  # Create color mapping for groups
+  group_colors <- setNames(rep(NA_character_, length(groups)), groups)
+  if (!is.null(color_scale) && !is.null(color_scale$palette)) {
+    if (is.function(color_scale$map)) {
+      # Use the scale's map function to get the correct color for each group
+      mapped_colors <- color_scale$map(groups)
+      group_colors <- setNames(mapped_colors, groups)
+    }
+  }
+  
+  # If no colors from scale, try to get from built data
+  if (all(is.na(group_colors))) {
+    col_choice <- intersect(c("colour", "color"), names(pb$data[[1]]))
+    if (length(col_choice) > 0) {
+      data <- cbind(data, color = pb$data[[1]][[col_choice[1]]])
+      # Try to map colors from data
+      for (group in groups) {
+        group_data <- if (inherits(data, "sf")) data[data[[id]] == group, , drop = FALSE] else data[data[, id] == group, , drop = FALSE]
+        if (nrow(group_data) > 0) {
+          group_color <- group_data$color[1]
+          # Check if it's a valid color
+          if (!is.na(group_color) && !is.numeric(group_color) && grepl("^#[0-9A-Fa-f]{6,8}$", group_color)) {
+            group_colors[group] <- group_color
+          }
+        }
+      }
+    } else {
+      data <- cbind(data, color = NA_character_)
+    }
   }
 
   labels.loc <- lapply(
@@ -6189,7 +6238,8 @@ LabelClusters <- function(
         )))
       }
       data.medians[, id] <- group
-      data.medians$color <- data.use$color[1]
+      # Assign color from group_colors mapping
+      data.medians$color <- group_colors[as.character(group)]
       return(data.medians)
     }
   )
@@ -6430,7 +6480,7 @@ CenterTitle <- function(...) {
 DarkTheme <- function(...) {
   #   Some constants for easier changing in the future
   black.background <- element_rect(fill = 'black')
-  black.background.no.border <- element_rect(fill = 'black', size = 0)
+  black.background.no.border <- element_rect(fill = 'black', linewidth = 0)
   font.margin <- 4
   white.text <- element_text(
     colour = 'white',
@@ -6441,8 +6491,8 @@ DarkTheme <- function(...) {
       l = font.margin
     )
   )
-  white.line <- element_line(colour = 'white', size = 1)
-  no.line <- element_line(size = 0)
+  white.line <- element_line(colour = 'white', linewidth = 1)
+  no.line <- element_line(linewidth = 0)
   #   Create the dark theme
   dark.theme <- theme(
     #   Set background colors