The assay parameter inside FeaturePlot_scCustom() does not work properly when combined with the split.by parameter

[CroixJeremy2]

Hello, it is me again,

Everything is in the title and in the following minimal example. I think something inside FeaturePlot_scCustom (and maybe other related functions?) is missing somewhere.

#Required libraries
	library(Seurat)
	library(ggplot2)
	library(patchwork)
	library(scCustomize)
	library(SeuratData)
	# SeuratData::InstallData("ifnb")

#Load data (minimal example)
	data(ifnb)
	ifnb = UpdateSeuratObject(object = ifnb)
	ifnb = Convert_Assay(seurat_object = ifnb, convert_to = "V5")

	head(ifnb)
	#                   orig.ident nCount_RNA nFeature_RNA stim seurat_annotations
	# AAACATACATTTCC.1 IMMUNE_CTRL       3017          877 CTRL          CD14 Mono
	# AAACATACCAGAAA.1 IMMUNE_CTRL       2481          713 CTRL          CD14 Mono
	# AAACATACCTCGCT.1 IMMUNE_CTRL       3420          850 CTRL          CD14 Mono
	# AAACATACCTGGTA.1 IMMUNE_CTRL       3156         1109 CTRL                pDC
	# AAACATACGATGAA.1 IMMUNE_CTRL       1868          634 CTRL       CD4 Memory T
	# AAACATACGGCATT.1 IMMUNE_CTRL       1581          557 CTRL          CD14 Mono
	# AAACATACTGCGTA.1 IMMUNE_CTRL       2747          980 CTRL        T activated
	# AAACATACTGCTGA.1 IMMUNE_CTRL       1341          581 CTRL        CD4 Naive T
	# AAACATTGAGTGTC.1 IMMUNE_CTRL       2155          880 CTRL              CD8 T
	# AAACATTGCTTCGC.1 IMMUNE_CTRL       2536          669 CTRL          CD14 Mono

	table(ifnb$seurat_annotations, ifnb$stim)
	#              CTRL STIM
	# CD14 Mono    2215 2147
	# CD4 Naive T   978 1526
	# CD4 Memory T  859  903
	# CD16 Mono     507  537
	# B             407  571
	# CD8 T         352  462
	# T activated   300  333
	# NK            298  321
	# DC            258  214
	# B Activated   185  203
	# Mk            115  121
	# pDC            51   81
	# Eryth          23   32

#Checking which meta.data column is set as the Idents()
	unique(Idents(ifnb)) #IMMUNE_CTRL IMMUNE_STIM, OK, as it corresponds to the two datasets of this study

#Quick standard pipeline in order to generate a UMAP embedding
	ifnb = NormalizeData(ifnb)
	ifnb = FindVariableFeatures(ifnb)
	ifnb = ScaleData(ifnb)
	ifnb = RunPCA(ifnb)
	ifnb = RunUMAP(ifnb, dims = 1:30)

	ifnb
	# An object of class Seurat 
	# 14053 features across 13999 samples within 1 assay 
	# Active assay: RNA (14053 features, 2000 variable features)
	#  3 layers present: counts, data, scale.data
	#  2 dimensional reductions calculated: pca, umap

#Checking the name of the default assay
	DefaultAssay(ifnb) #"RNA", as expected

#DefaultAssay behaviour, OK, it works as expected
	#Retrieving the range
		range(FetchData(object = ifnb, vars = "GAPDH", assay = "RNA")[,1], na.rm = TRUE) #0.000000 5.225659

	#Merged graph on RNA assay, OK
		FeaturePlot_scCustom(seurat_object = ifnb,
								features = "GAPDH") #color scale matches with the range of FetchData, OK

---

	#Split graph on RNA assay, OK
		FeaturePlot_scCustom(seurat_object = ifnb,
								features = "GAPDH",
								split.by = "stim") #color scale matches with the range of FetchData, OK

---

#Creating another new assay (quickly using SCT for example, but in theory it could be anything else)
	ifnb = SCTransform(object = ifnb, return.only.var.genes = FALSE)

	ifnb
	# An object of class Seurat 
	# 27515 features across 13999 samples within 2 assays 
	# Active assay: SCT (13462 features, 3000 variable features)
	#  3 layers present: counts, data, scale.data
	#  1 other assay present: RNA
	#  2 dimensional reductions calculated: pca, umap

	DefaultAssay(ifnb) # now "SCT", but let's switch it back to "RNA" (the main one in my real-non-minimal data; but for clarity in this issue I used ifnb as a minimal example)
	DefaultAssay(ifnb) = "RNA"

#Same graphs, but made on the second assay: PROBLEM
	#Retrieving the range
		range(FetchData(object = ifnb, vars = "GAPDH", assay = "SCT")[,1], na.rm = TRUE) #0.000000 3.433987; different than the RNA assay, OK, as expected

	#Merged graph on RNA assay, OK
		FeaturePlot_scCustom(seurat_object = ifnb,
								assay = "SCT",
								features = "GAPDH") #color scale matches with the range of FetchData for this specific assay, OK

---

	#Split graph on RNA assay, NOT OK
		FeaturePlot_scCustom(seurat_object = ifnb,
								assay = "SCT",
								features = "GAPDH",
								split.by = "stim") #the color scale does not match the SCT assay limits anymore... and goes back to the one of the RNA assay somehow...

---

#Remark: switching manually the DefaultAssay BEFORE FeaturePlot_scCustom is a workaround, but it would be way better to rather make the "assay" parameter of FeaturePlot_scCustom working properly
	DefaultAssay(ifnb) = "SCT"

	#Split graph on RNA assay, OK now
		FeaturePlot_scCustom(seurat_object = ifnb,
								# assay = "SCT", #no need to use it anymore here
								features = "GAPDH",
								split.by = "stim") #the color scale NOW does match the SCT assay limits, but it would be better to select correctly the assay inside FeaturePlot_scCustom()

---

Is there a way to make the assay parameter work when combined with split.by please?

Thanks in advance for your response, and thanks again for this great package,
Best regards,

---

sessionInfo() output

R version 4.5.3 (2026-03-11)
Platform: aarch64-apple-darwin20
Running under: macOS Tahoe 26.5.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8

time zone: Europe/Paris
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] future_1.70.0           pbmc3k.SeuratData_3.1.4 ifnb.SeuratData_3.0.0  
 [4] SeuratData_0.2.2.9002   scCustomize_3.3.0       patchwork_1.3.2        
 [7] ggplot2_4.0.3           Seurat_5.5.0            SeuratObject_5.4.0     
[10] sp_2.2-1               

loaded via a namespace (and not attached):
  [1] RColorBrewer_1.1-3          rstudioapi_0.18.0          
  [3] jsonlite_2.0.0              shape_1.4.6.1              
  [5] magrittr_2.0.5              spatstat.utils_3.2-3       
  [7] ggbeeswarm_0.7.3            farver_2.1.2               
  [9] GlobalOptions_0.1.4         vctrs_0.7.3                
 [11] ROCR_1.0-12                 DelayedMatrixStats_1.32.0  
 [13] spatstat.explore_3.8-0      paletteer_1.7.0            
 [15] janitor_2.2.1               S4Arrays_1.10.1            
 [17] htmltools_0.5.9             forcats_1.0.1              
 [19] SparseArray_1.10.10         sctransform_0.4.3          
 [21] parallelly_1.47.0           KernSmooth_2.23-26         
 [23] htmlwidgets_1.6.4           ica_1.0-3                  
 [25] plyr_1.8.9                  plotly_4.12.0              
 [27] zoo_1.8-15                  lubridate_1.9.5            
 [29] igraph_2.3.1                mime_0.13                  
 [31] lifecycle_1.0.5             pkgconfig_2.0.3            
 [33] Matrix_1.7-4                R6_2.6.1                   
 [35] fastmap_1.2.0               MatrixGenerics_1.22.0      
 [37] fitdistrplus_1.2-6          shiny_1.13.0               
 [39] snakecase_0.11.1            digest_0.6.39              
 [41] colorspace_2.1-2            rematch2_2.1.2             
 [43] S4Vectors_0.48.1            tensor_1.5.1               
 [45] RSpectra_0.16-2             irlba_2.3.7                
 [47] GenomicRanges_1.62.1        beachmat_2.26.0            
 [49] labeling_0.4.3              progressr_0.19.0           
 [51] spatstat.sparse_3.1-0       timechange_0.4.0           
 [53] httr_1.4.8                  polyclip_1.10-7            
 [55] abind_1.4-8                 compiler_4.5.3             
 [57] withr_3.0.2                 S7_0.2.2                   
 [59] fastDummies_1.7.6           MASS_7.3-65                
 [61] DelayedArray_0.36.1         rappdirs_0.3.4             
 [63] tools_4.5.3                 vipor_0.4.7                
 [65] lmtest_0.9-40               otel_0.2.0                 
 [67] beeswarm_0.4.0              httpuv_1.6.17              
 [69] future.apply_1.20.2         goftest_1.2-3              
 [71] glmGamPoi_1.22.0            glue_1.8.1                 
 [73] nlme_3.1-168                promises_1.5.0             
 [75] grid_4.5.3                  Rtsne_0.17                 
 [77] cluster_2.1.8.2             reshape2_1.4.5             
 [79] generics_0.1.4              gtable_0.3.6               
 [81] spatstat.data_3.1-9         tidyr_1.3.2                
 [83] data.table_1.18.4           XVector_0.50.0             
 [85] BiocGenerics_0.56.0         spatstat.geom_3.7-3        
 [87] RcppAnnoy_0.0.23            ggrepel_0.9.8              
 [89] RANN_2.6.2                  pillar_1.11.1              
 [91] stringr_1.6.0               spam_2.11-3                
 [93] RcppHNSW_0.6.0              ggprism_1.0.7              
 [95] later_1.4.8                 circlize_0.4.18            
 [97] splines_4.5.3               dplyr_1.2.1                
 [99] lattice_0.22-9              survival_3.8-6             
[101] deldir_2.0-4                tidyselect_1.2.1           
[103] miniUI_0.1.2                pbapply_1.7-4              
[105] gridExtra_2.3               Seqinfo_1.0.0              
[107] IRanges_2.44.0              SummarizedExperiment_1.40.0
[109] scattermore_1.2             stats4_4.5.3               
[111] Biobase_2.70.0              matrixStats_1.5.0          
[113] stringi_1.8.7               lazyeval_0.2.3             
[115] codetools_0.2-20            tibble_3.3.1               
[117] cli_3.6.6                   uwot_0.2.4                 
[119] xtable_1.8-8                reticulate_1.46.0          
[121] Rcpp_1.1.1-1.1              globals_0.19.1             
[123] spatstat.random_3.4-5       mcprogress_0.1.1           
[125] png_0.1-9                   ggrastr_1.0.2              
[127] spatstat.univar_3.2-0       parallel_4.5.3             
[129] dotCall64_1.2               sparseMatrixStats_1.22.0   
[131] listenv_0.10.1              viridisLite_0.4.3          
[133] scales_1.4.0                ggridges_0.5.7             
[135] crayon_1.5.3                purrr_1.2.2                
[137] rlang_1.2.0                 cowplot_1.2.0              

[samuel-marsh]

Hi @CroixJeremy2,

Thanks so much for bug report and never worry about posting them too much :). I almost certainly know what's causing this and will work on fix. Issue is that assay is new parameter for FeaturePlot and works on update to FetchData that also adds assay parameter. So while scCustomize passes assay in its FeaturePlot calls it doesn't for FetchData. I will set up conditional solution that depends on users Seurat version. In meantime workaround is just to set default assay before plotting so that scCustomize extracts data correctly. Will update when fix is live.

Thanks!
Sam

[CroixJeremy2]

Workaround/patch by AI... (look for "✅ FIX" comments). I don't know if the solution/logic is OK or not though:

# Patched FeaturePlot_scCustom ----
# Fixes: assay parameter (passed via ...) was not forwarded to FetchData()
# in the split.by code paths, causing wrong color scale limits.
# Original issue: https://github.com/samuel-marsh/scCustomize/issues/267

FeaturePlot_scCustom_MANUAL_UPDATE_2 = function(
	seurat_object,
	features,
	colors_use = viridis_plasma_dark_high,
	na_color = "lightgray",
	order = TRUE,
	pt.size = NULL,
	reduction = NULL,
	na_cutoff = 1e-09,
	raster = NULL,
	raster.dpi = c(512, 512),
	split.by = NULL,
	split_collect = NULL,
	aspect_ratio = NULL,
	figure_plot = FALSE,
	num_columns = NULL,
	layer = "data",
	alpha_exp = NULL,
	alpha_na_exp = NULL,
	label = FALSE,
	label_feature_yaxis = FALSE,
	max.cutoff = NA,
	min.cutoff = NA,
	combine = TRUE,
	...
) {
	Is_Seurat(seurat_object = seurat_object)

	# ✅ FIX #1: capture assay from ... so we can forward it to FetchData() later.
	# If the user did not pass assay, this is NULL and FetchData() will use DefaultAssay.
	dots = list(...)
	assay_use = dots[["assay"]]

	if (!is.null(x = split.by)) {
		split.by = Meta_Present(
			object = seurat_object,
			meta_col_names = split.by,
			print_msg = FALSE,
			omit_warn = FALSE
		)[[1]]
	}

	if (is.null(x = split_collect)) {
		if (length(x = features) == 1) {
			split_collect = TRUE
		} else {
			split_collect = FALSE
		}
	}

	if (!is.null(x = split_collect)) {
		if (length(x = features) > 1 && isTRUE(x = split_collect)) {
			cli::cli_abort(message = "{.code split_collect} cannot be set to {.field TRUE} if the number of features is greater than 1.")
		}
	}

	all_found_features = Feature_PreCheck(
		object = seurat_object,
		features = features
	)

	if (is.null(x = split.by) && isTRUE(x = label_feature_yaxis)) {
		cli::cli_abort(message = "Setting {.code label_feature_yaxis = TRUE} is only supported when also setting {.code split.by}.")
	}

	if (!is.null(x = split.by)) {
		split.by_length = length(x = unique(x = seurat_object@meta.data[[split.by]]))
		if (!is.null(x = num_columns) && isTRUE(x = label_feature_yaxis)) {
			cli::cli_warn(message = c(
				"Setting number of columns is not permitted if {.code label_feature_yaxis = TRUE}",
				i = "Number of columns be automatically set to number of levels in `split.by` ({.field {split.by_length}})."
			))
			num_columns = split.by_length
		}
		if (is.null(x = num_columns)) {
			num_columns = split.by_length
		}
		num_rows = ceiling(split.by_length / num_columns)
		if (num_columns > split.by_length) {
			cli::cli_abort(message = c(
				"The number of columns specified is greater than the number of meta data variables.",
				`*` = "{.val {split.by}} only contains {.field {split.by_length}} variables.",
				i = "Please adjust {.code num_columns} to be less than or equal to {.field {split.by_length}}."
			))
		}
	}

	if (any(all_found_features %in% colnames(x = seurat_object@meta.data))) {
		cli::cli_warn(message = c(
			"Some of the plotted features are from meta.data slot.",
			`*` = "Please check that {.code na_cutoff} param is being set appropriately for those features."
		))
	}

	raster = raster %||% (length(x = Cells(x = seurat_object)) > 2e+05)

	if (is.null(x = pt.size)) {
		if (!is.null(x = split.by)) {
			cells_by_meta = data.frame(table(seurat_object@meta.data[, split.by]))
			max_cells = max(cells_by_meta$Freq)
			pt.size = AutoPointSize_scCustom(data = max_cells, raster = raster)
		} else {
			cells_total = nrow(x = seurat_object@meta.data)
			pt.size = AutoPointSize_scCustom(data = cells_total, raster = raster)
		}
	}

	if (is.null(x = na_cutoff)) {
		na_cutoff = NA
	}

	reduction = reduction %||% DefaultDimReduc(object = seurat_object)
	seurat_version = packageVersion("Seurat")

	if (!is.null(x = alpha_exp) && seurat_version < "5") {
		colors_use = scales::alpha(colors_use, alpha_exp)
	}
	if (!is.null(x = alpha_na_exp) && seurat_version < "5") {
		na_color = scales::alpha(na_color, alpha_na_exp)
	}
	if (!is.null(x = alpha_na_exp) && seurat_version >= "5") {
		cli::cli_warn(message = "{.code alpha_na_exp} is not currently supported for Seurat v5+")
	}
	if (is.null(x = alpha_exp) && seurat_version >= "5") {
		alpha_exp = 1
	}

	if (length(x = min.cutoff) > 1 && length(x = min.cutoff) != length(x = all_found_features)) {
		cli::cli_abort(message = "The length of {.code min.cutoff} ({.field {length(x = min.cutoff)}}) must be single value or equal in length to number of features being plotted ({.field {length(x = all_found_features)}}).")
	}
	if (length(x = max.cutoff) > 1 && length(x = max.cutoff) != length(x = all_found_features)) {
		cli::cli_abort(message = "The length of {.code max.cutoff} ({.field {length(x = max.cutoff)}}) must be single value or equal in length to number of features being plotted ({.field {length(x = all_found_features)}}).")
	}

	# --- Branch 1: no split.by, combine = TRUE ---
	if (is.null(x = split.by) && isTRUE(x = combine)) {
		if (seurat_version >= "5") {
			plot = suppressMessages(
				FeaturePlot(
					object = seurat_object,
					features = all_found_features,
					order = order,
					pt.size = pt.size,
					reduction = reduction,
					raster = raster,
					split.by = split.by,
					ncol = num_columns,
					combine = combine,
					raster.dpi = raster.dpi,
					label = label,
					alpha = alpha_exp,
					min.cutoff = min.cutoff,
					max.cutoff = max.cutoff,
					...
				) & scale_color_gradientn(
					colors = colors_use,
					limits = c(na_cutoff, NA),
					na.value = na_color
				)
			)
		} else {
			plot = suppressMessages(
				FeaturePlot(
					object = seurat_object,
					features = all_found_features,
					order = order,
					pt.size = pt.size,
					reduction = reduction,
					raster = raster,
					split.by = split.by,
					ncol = num_columns,
					combine = combine,
					raster.dpi = raster.dpi,
					label = label,
					min.cutoff = min.cutoff,
					max.cutoff = max.cutoff,
					...
				) & scale_color_gradientn(
					colors = colors_use,
					limits = c(na_cutoff, NA),
					na.value = na_color
				)
			)
		}
	}

	# --- Branch 2: no split.by, combine = FALSE ---
	if (is.null(x = split.by) && isFALSE(x = combine)) {
		if (seurat_version >= "5") {
			plot_list = suppressMessages(FeaturePlot(
				object = seurat_object,
				features = all_found_features,
				order = order,
				pt.size = pt.size,
				reduction = reduction,
				raster = raster,
				split.by = split.by,
				ncol = num_columns,
				combine = combine,
				raster.dpi = raster.dpi,
				label = label,
				alpha = alpha_exp,
				min.cutoff = min.cutoff,
				max.cutoff = max.cutoff,
				...
			))
			plot = lapply(1:length(x = plot_list), function(i) {
				plot_list[[i]] = suppressMessages(
					plot_list[[i]] + scale_color_gradientn(
						colors = colors_use,
						limits = c(na_cutoff, NA),
						na.value = na_color
					)
				)
			})
		} else {
			plot_list = suppressMessages(FeaturePlot(
				object = seurat_object,
				features = all_found_features,
				order = order,
				pt.size = pt.size,
				reduction = reduction,
				raster = raster,
				split.by = split.by,
				ncol = num_columns,
				combine = combine,
				raster.dpi = raster.dpi,
				label = label,
				min.cutoff = min.cutoff,
				max.cutoff = max.cutoff,
				...
			))
			plot = lapply(1:length(x = plot_list), function(i) {
				plot_list[[i]] = suppressMessages(
					plot_list[[i]] + scale_color_gradientn(
						colors = colors_use,
						limits = c(na_cutoff, NA),
						na.value = na_color
					)
				)
			})
		}
	}

	# --- Branch 3: split.by + single feature ---
	if (!is.null(x = split.by) && length(x = all_found_features) == 1) {
		# ✅ FIX #2: pass assay_use to FetchData so expression range is
		# computed from the correct assay (not necessarily DefaultAssay).
		feature_data = FetchData(
			object = seurat_object,
			vars = all_found_features,
			layer = layer,
			assay = assay_use # ✅ FIX #2
		)

		min.cutoff = mapply(
			FUN = function(cutoff, feature) {
				return(ifelse(
					test = is.na(x = cutoff),
					yes = min(feature_data[, feature]),
					no = cutoff
				))
			},
			cutoff = min.cutoff,
			feature = all_found_features
		)
		max.cutoff = mapply(
			FUN = function(cutoff, feature) {
				return(ifelse(
					test = is.na(x = cutoff),
					yes = max(feature_data[, feature]),
					no = cutoff
				))
			},
			cutoff = max.cutoff,
			feature = all_found_features
		)

		check.lengths = unique(x = vapply(
			X = list(all_found_features, min.cutoff, max.cutoff),
			FUN = length,
			FUN.VALUE = numeric(length = 1)
		))
		if (length(x = check.lengths) != 1) {
			cli::cli_abort(message = "There must be the same number of minimum and maximum cuttoffs as there are features")
		}

		names(x = min.cutoff) = names(x = max.cutoff) = all_found_features
		max_exp_value = max.cutoff
		min_exp_value = min.cutoff

		if (seurat_version >= "5") {
			plot = suppressMessages(
				FeaturePlot(
					object = seurat_object,
					features = all_found_features,
					order = order,
					pt.size = pt.size,
					reduction = reduction,
					raster = raster,
					split.by = split.by,
					raster.dpi = raster.dpi,
					label = label,
					alpha = alpha_exp,
					min.cutoff = min.cutoff,
					max.cutoff = max.cutoff,
					...
				) & scale_color_gradientn(
					colors = colors_use,
					limits = c(na_cutoff, max_exp_value),
					na.value = na_color,
					name = all_found_features
				)
			) & RestoreLegend() & theme(axis.title.y.right = element_blank())
		} else {
			plot = suppressMessages(
				FeaturePlot(
					object = seurat_object,
					features = all_found_features,
					order = order,
					pt.size = pt.size,
					reduction = reduction,
					raster = raster,
					split.by = split.by,
					raster.dpi = raster.dpi,
					label = label,
					min.cutoff = min.cutoff,
					max.cutoff = max.cutoff,
					...
				) & scale_color_gradientn(
					colors = colors_use,
					limits = c(na_cutoff, max_exp_value),
					na.value = na_color,
					name = all_found_features
				)
			) & RestoreLegend() & theme(axis.title.y.right = element_blank())
		}

		if (isTRUE(x = label_feature_yaxis)) {
			plot = plot + plot_layout(nrow = num_rows, ncol = num_columns)
			plot = plot & theme(legend.title = element_blank())
			plot = suppressMessages(
				plot + scale_y_continuous(sec.axis = dup_axis(name = all_found_features))
			) + No_Right()
		} else {
			if (isTRUE(x = split_collect)) {
				if (hasArg("keep.scale")) {
					cli::cli_abort(message = "The parameter {.code keep.scale} cannot be set different from default if {.code split_collect - TRUE}.")
				}
				plot = plot + plot_layout(nrow = num_rows, ncol = num_columns, guides = "collect")
			} else {
				plot = plot + plot_layout(nrow = num_rows, ncol = num_columns)
			}
		}
	}

	# --- Branch 4: split.by + multiple features ---
	if (!is.null(x = split.by) && length(x = all_found_features) > 1) {
		if (length(x = min.cutoff) == 1) {
			min.cutoff = rep(x = min.cutoff, length(x = all_found_features))
		}
		if (length(x = max.cutoff) == 1) {
			max.cutoff = rep(x = max.cutoff, length(x = all_found_features))
		}

		plot_list = lapply(1:length(x = all_found_features), function(i) {
			# ✅ FIX #3: pass assay_use to FetchData (same fix as #2, inside lapply).
			feature_data = FetchData(
				object = seurat_object,
				vars = all_found_features[i],
				layer = layer,
				assay = assay_use # ✅ FIX #3
			)

			min.cutoff = mapply(
				FUN = function(cutoff, feature) {
					return(ifelse(
						test = is.na(x = cutoff),
						yes = min(feature_data[, feature]),
						no = cutoff
					))
				},
				cutoff = min.cutoff[i],
				feature = all_found_features[i]
			)
			max.cutoff = mapply(
				FUN = function(cutoff, feature) {
					return(ifelse(
						test = is.na(x = cutoff),
						yes = max(feature_data[, feature]),
						no = cutoff
					))
				},
				cutoff = max.cutoff[i],
				feature = all_found_features[i]
			)

			check.lengths = unique(x = vapply(
				X = list(all_found_features[i], min.cutoff, max.cutoff),
				FUN = length,
				FUN.VALUE = numeric(length = 1)
			))
			if (length(x = check.lengths) != 1) {
				cli::cli_abort(message = "There must be the same number of minimum and maximum cuttoffs as there are features")
			}

			names(x = min.cutoff) = names(x = max.cutoff) = all_found_features[i]
			max_exp_value = max.cutoff
			min_exp_value = min.cutoff

			if (seurat_version >= "5") {
				single_plot = suppressMessages(
					FeaturePlot(
						object = seurat_object,
						features = all_found_features[i],
						order = order,
						pt.size = pt.size,
						reduction = reduction,
						raster = raster,
						split.by = split.by,
						raster.dpi = raster.dpi,
						label = label,
						alpha = alpha_exp,
						min.cutoff = min.cutoff,
						max.cutoff = max.cutoff,
						...
					) & scale_color_gradientn(
						colors = colors_use,
						limits = c(na_cutoff, max_exp_value),
						na.value = na_color,
						name = all_found_features[i]
					)
				) & RestoreLegend() & theme(axis.title.y.right = element_blank())
			} else {
				single_plot = suppressMessages(
					FeaturePlot(
						object = seurat_object,
						features = all_found_features[i],
						order = order,
						pt.size = pt.size,
						reduction = reduction,
						raster = raster,
						split.by = split.by,
						raster.dpi = raster.dpi,
						label = label,
						min.cutoff = min.cutoff,
						max.cutoff = max.cutoff,
						...
					) & scale_color_gradientn(
						colors = colors_use,
						limits = c(na_cutoff, max_exp_value),
						na.value = na_color,
						name = all_found_features[i] # ✅ FIX #4: was features[i], should be all_found_features[i] for consistency
					)
				) & RestoreLegend() & theme(axis.title.y.right = element_blank())
			}

			if (isTRUE(x = label_feature_yaxis)) {
				single_plot = single_plot + plot_layout(nrow = num_rows, ncol = num_columns)
				single_plot = single_plot & theme(legend.title = element_blank())
				single_plot = suppressMessages(
					single_plot + scale_y_continuous(sec.axis = dup_axis(name = all_found_features[i]))
				) + No_Right()
			} else {
				single_plot = single_plot + plot_layout(nrow = num_rows, ncol = num_columns)
			}
		})

		plot = wrap_plots(plot_list) + plot_layout(ncol = 1)
	}

	if (getOption(x = "scCustomize_warn_na_cutoff", default = TRUE) &&
		!is.na(x = na_cutoff) && na_cutoff == 1e-09) {
		cli::cli_inform(message = c(
			"",
			"NOTE: {.field FeaturePlot_scCustom} uses a specified {.code na_cutoff} when plotting to",
			"color cells with no expression as background color separate from color scale.",
			"Please ensure {.code na_cutoff} value is appropriate for feature being plotted.",
			"Default setting is appropriate for use when plotting from 'RNA' assay.\n",
			"When {.code na_cutoff} not appropriate (e.g., module scores) set to NULL to \n",
			"plot all cells in gradient color palette.",
			"",
			"-----This message will be shown once per session.-----"
		))
		options(scCustomize_warn_na_cutoff = FALSE)
	}

	if (getOption(x = "scCustomize_warn_zero_na_cutoff", default = TRUE) &&
		!is.na(x = na_cutoff) && na_cutoff == 0) {
		cli::cli_inform(message = c(
			"",
			"NOTE: Specified {.code na_cutoff} is set to {.field zero (0)}. This means that only cells/nuclei",
			"with expression less than zero will be plotted with {.code na_color}: {.val {na_color}}.",
			"To plot cells with expression values of zero using {.code na_color} leave",
			"default {.code na_cutoff} value. If you want to plot full spectrum without",
			"{.code na_cutoff} (e.g., for module scores) then set {.code na_cutoff = NULL`}.",
			"",
			"-----This message will be shown once per session.-----"
		))
		options(scCustomize_warn_na_cutoff = FALSE)
	}

	if (!is.null(x = aspect_ratio)) {
		if (!is.numeric(x = aspect_ratio)) {
			cli::cli_abort(message = "{.code aspect_ratio} must be a {.field numeric} value.")
		}
		plot = plot & theme(aspect.ratio = aspect_ratio)
	}

	if (isTRUE(x = figure_plot)) {
		if (length(x = all_found_features) == 1) {
			plot = Figure_Plot(plot = plot)
		} else {
			plot_list = lapply(1:length(x = all_found_features), function(j) {
				fig_plot = Figure_Plot(plot = plot[[j]])
			})
			plot = wrap_plots(plot_list, ncol = num_columns)
		}
	}

	return(plot)
}

# [PATCH] Bind the function to scCustomize's namespace so it can access
# internal (non-exported) helpers: Is_Seurat, Meta_Present, Feature_PreCheck,
# AutoPointSize_scCustom, Figure_Plot, No_Right, etc.
environment(FeaturePlot_scCustom_MANUAL_UPDATE_2) = asNamespace("scCustomize")
# [END PATCH]

After that, the expected behaviour is happening:

DefaultAssay(ifnb) #RNA
FeaturePlot_scCustom_MANUAL_UPDATE_2(seurat_object = ifnb,
								assay = "SCT",
								features = "GAPDH",
								split.by = "stim")