---
license: mit
repo: https://github.com/kahkengwong/CEP-IP_Framework
task_categories:
- other
tags:
- single-cell
- scRNA-seq
- prostate-cancer
- glioblastoma
- brain
- seurat
- transcriptomics
- bioinformatics
- TRPM4
- Ribosome
- GAM
- REML
- TPRS
- cancer-research
- gene-expression
language:
- en
data_files:
- "GSE185344_Seurat_processed.RData"
- "AllenMTG_Seurat_processed.RData"
- "NeftelGBM_SS2_AdultMalignant_Seurat_processed.RData"
size_categories:
- 10K<n<100K
arxiv: 2509.12073
doi: 10.1016/j.cmpb.2026.109372
configs:
- config_name: default
  data_files:
  - split: prostate_cancer
    path: Dataset_Overview.csv
---

# CEP-IP: An Explainable Framework for Cell Subpopulation Identification in Single-cell Transcriptomics (by Kah Keng Wong) ([Published in *Computer Methods and Programs in Biomedicine*](https://doi.org/10.1016/j.cmpb.2026.109372))

## 🧬 Abstract

**Background and objective:** Single-cell RNA sequencing (scRNA-seq) frameworks lack explainable approaches for identifying cell subpopulations harboring strong pairwise monotonic gene-module relationships between a gene of interest (GOI) and its co-expressed genes. In this study, CEP-IP is introduced as a novel explainable machine learning framework to address this gap.

**Methods:** Prostate cancer (PCa) scRNA-seq dataset was used as the initial dataset, whereby _TRPM4_ served as the GOI and its co-expressed ribosomal genes (Ribo) were identified via Spearman-Kendall dual-filter (i.e., dual-filtered genes, DFG). Next, generalized additive modeling quantified the strength of _TRPM4_-Ribo relationship, represented by deviance explained (DE). _TRPM4_-Ribo’s DE was then assigned to individual cells via cell explanatory power (CEP) classification, identifying cells harboring the _TRPM4_-Ribo module [i.e., top-ranked explanatory power (TREP) cells]. _TRPM4_-Ribo transcriptional space was then stratified into pre-IP and post-IP regions using inflection point (IP) analysis, producing four distinct cell subpopulations per patient for pathway analysis. Validation was performed in the Allen middle temporal gyrus (MTG) and Neftel glioblastoma multiforme (GBM) transcriptomically heterogeneous datasets.

**Results:** _TRPM4_-Ribo modeling outperformed alternative gene set modules (FDR<0.05). In each PCa patient, CEP-IP yielded four cell subpopulations, where pre-IP TREP cells showed enrichment of immune-related processes, and post-IP TREP cells were enriched for ribosomal, translation, and cell adhesion pathways. In the MTG validation dataset (_CARM1P1_-DFG module), post-IP TREP cells showed enrichment of neuron projection ontologies. In the GBM dataset, _FOXM1_ was the sole GOI yielding mesenchymal-state DFG, with _FOXM1_-DFG post-IP TREP cells enriched for cell division and microtubule pathways; 3D trajectory analysis demonstrated continuous trajectories of TREP cells that were obscured in 2D embeddings.

**Conclusions:** CEP-IP identifies biologically distinct cell subpopulations in three independent scRNA-seq datasets. The framework may generalize to other pairwise GOI-DFG module in single-cell transcriptomics beyond the datasets investigated in this study.

> Wong KK (2026). *Computer Methods and Programs in Biomedicine*, 282, 109372. https://doi.org/10.1016/j.cmpb.2026.109372

## 📊 Dataset Overview
This repository contains processed single-cell RNA-seq Seurat objects for three independent datasets used in the CEP-IP framework study: (1) a prostate cancer (PCa) primary dataset, and two validation datasets comprising (2) the Neftel glioblastoma multiforme (GBM) dataset and (3) the Allen Human Middle Temporal Gyrus (MTG) dataset. All datasets are quality-controlled and normalized using Seurat, with clustering applied where appropriate (see Technical Details).

**Contents:**
- Processed Seurat object for the primary PCa dataset: `GSE185344_Seurat_processed.RData` (9.52 GB)
- Processed Seurat object for the Allen MTG validation dataset: `AllenMTG_Seurat_processed.RData`
- Processed Seurat object for the Neftel GBM validation dataset: `NeftelGBM_SS2_AdultMalignant_Seurat_processed.RData`
- Quality-controlled expression data; batch-corrected UMAP embeddings and cluster annotations where applicable (PCa and Allen MTG)
- Cluster marker genes and metadata
- Ready for subsequent generalized additive model (GAM) analysis with the `mgcv` package


## 📚 Original Data Sources
This repository contains processed versions of the following publicly available datasets:

**Primary dataset (PCa):**
Wong HY, Sheng Q, Hesterberg AB, Croessmann S et al. Single cell analysis of cribriform prostate cancer reveals cell intrinsic and tumor microenvironmental pathways of aggressive disease. Nat Commun 2022;13(1):6036. https://doi.org/10.1038/s41467-022-33780-1

**Validation dataset 1 (Allen MTG):**
Hodge RD, Bakken TE, Miller JA, Smith KA et al. Conserved cell types with divergent features in human versus mouse cortex. Nature 2019;573:61-68. https://doi.org/10.1038/s41586-019-1506-7

**Validation dataset 2 (GBM):**
Neftel C, Laffy J, Filbin MG, Hara T et al. An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. Cell 2019;178:835-849. https://doi.org/10.1016/j.cell.2019.06.024


## 🛠️ Technical Details

### Processing Strategies
Each dataset was processed at the level appropriate for its intended use in the CEP-IP framework. The PCa and Allen MTG datasets were raw count matrices processed through a complete Seurat pipeline (QC → normalization → integration → clustering). The Neftel GBM dataset was published with annotated cell state derived from the authors' own TPM-based scoring method; accordingly, it was prepared for CEP-IP analysis without re-clustering, using the original Neftel et al. cell state labels (AC, MES, NPC, OPC) directly.

### Dataset-Specific Processing Summary

| Parameter | PCa (GSE185344) | Allen MTG | Neftel GBM (SS2) |
|---|---|---|---|
| Total Cells | 30,932 (tumor) / 12,205 (benign) | 13,616 | 4,916 |
| Total Genes | 846 (tumor) / 779 (benign) | 50,281 (1,328 integrated) | 23,686 |
| Samples | 7 tumor / 7 benign | 8 donors | 20 tumors |
| Cell Groups | 14 clusters each | 17 clusters | 4 states (AC/MES/NPC/OPC) |
| Seurat Version | 5.1.0 | 5.3.0 | 5.3.0 |
| Normalization | LogNormalize + SCTransform | LogNormalize + SCTransform | LogNormalize |
| Integration | CCA | CCA | N/A (published annotations) |
| Dimensionality Reduction | PCA + UMAP | PCA + UMAP | N/A |
| Clustering Resolution | 0.5 | 0.5 | N/A |
| Variable Features | 2,000 (VST) | 1,328 (integration features) | 2,000 (VST) |
| PCA Dimensions | 1:15 | 1:8 | N/A |
| UMAP Metric | Cosine | Cosine | N/A |
| Cell Cycle Regression | Yes (S.Score + G2M.Score) | Yes (S.Score + G2M.Score) | Scored only (G1S + G2M) |
| Doublet Removal | Yes (scDblFinder) | Yes (scDblFinder) | N/A (published data) |
| Batch Correction | Yes (SCTransform + CCA) | Yes (SCTransform + CCA) | N/A |
| Min Features/Cell | >500 | >500 | N/A (Neftel QC) |
| Ribo Filter | 90th percentile | 90th percentile | N/A |
| MT Filter | 90th percentile | 90th percentile | N/A |

### Key R Packages Used
- `future` / `parallel` (parallel processing)
- `scDblFinder` (doublet detection)
- `Seurat` (single-cell analysis)
- `SingleCellExperiment` (data conversion)

### Key Processing Steps

**PCa and Allen MTG datasets (full pipeline):**
- Quality control filtering (>500 features per cell)
- Ribosomal gene filtering (cells above 90th percentile removed)
- Mitochondrial gene filtering (cells above 90th percentile removed)
- Cell cycle scoring and regression (S.Score + G2M.Score)
- Doublet removal using scDblFinder
- Batch effect correction via SCTransform + CCA integration
- Dimensionality reduction: PCA (elbow plot assessment) + UMAP (cosine metric)
- Clustering at resolution 0.5

**Neftel GBM dataset (annotation-based):**
- Published Smart-seq2 data with pre-assigned malignant cell state labels (AC, MES, NPC, OPC) from Neftel et al. 2019
- Pure malignant adult GBM cells (n=4,916) extracted directly from the Neftel et al. processed dataset, where malignant identity was established by the original authors
- Variable feature selection (2,000 VST) and scaling for GAM expression matrix preparation
- Cell cycle scoring retained from original metadata (G1S, G2M scores)
- No re-clustering performed; original Neftel state assignments used directly for CEP-IP analysis


## 🧮 Code Availability
The complete analysis pipeline of this CEP-IP framework project, framework diagrams, and graphical abstract are available on GitHub:
[GitHub repository](https://github.com/kahkengwong/CEP-IP_Framework).

**Primary dataset (PCa):** After downloading `GSE185344_Seurat_processed.RData`, run the code starting from `Part_2_UMAP_Heatmap_Spearman-Kendall's-matrix.r` until `Part_3.15_CEP-IP_in_Monocle3_Trajectory.r`.

**Validation datasets (Allen MTG and Neftel GBM):** After downloading `AllenMTG_Seurat_processed.RData` and `NeftelGBM_SS2_AdultMalignant_Seurat_processed.RData`, run the corresponding validation scripts available on GitHub.


## 🎯 Citation
If you use any of these processed datasets, please cite:

Wong KK (2026). CEP-IP: An explainable framework for cell subpopulation identification in single-cell transcriptomics. Computer Methods and Programs in Biomedicine, 282, 109372. https://doi.org/10.1016/j.cmpb.2026.109372

Please also cite the respective source datasets:

Wong HY, Sheng Q, Hesterberg AB, Croessmann S et al (2022). Single cell analysis of cribriform prostate cancer reveals cell intrinsic and tumor microenvironmental pathways of aggressive disease. Nat Commun 13(1):6036. https://doi.org/10.1038/s41467-022-33780-1

Hodge RD, Bakken TE, Miller JA, Smith KA et al (2019). Conserved cell types with divergent features in human versus mouse cortex. Nature 573:61-68. https://doi.org/10.1038/s41586-019-1506-7

Neftel C, Laffy J, Filbin MG, Hara T et al (2019). An Integrative Model of Cellular States, Plasticity, and Genetics for Glioblastoma. Cell 178:835-849. https://doi.org/10.1016/j.cell.2019.06.024


## 📋 License
This dataset is licensed under the MIT License.

<details>
<summary><strong>📝 Click to view complete MIT License</strong></summary>

```text
MIT License

Copyright (c) 2025 Kah Keng Wong

Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
in the Software without restriction, including without limitation the rights
to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
SOFTWARE.
```
</details>


## 🔬 Detailed Analysis Pipeline
The complete PCa and Allen MTG processing pipeline includes: data loading → quality control → cell cycle regression → doublet removal → batch correction → clustering → UMAP visualization. The Neftel GBM dataset followed an annotation-based preparation workflow; see Key Processing Steps above.
<details>
<summary><strong>📝 Click to view complete processing code</strong></summary>
  
```r
##########################################
# A. Dataset Description
##########################################
This dataset contains scRNA-seq data processed using Seurat v5.1.0, and the dataset was obtained from: 
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE185344

The original data was published by Wong et al. Nat Commun 2022;13(1):6036 (doi: 10.1038/s41467-022-33780-1)
available from: https://pubmed.ncbi.nlm.nih.gov/36229464/ 

Note: The code below covers the PCa primary dataset only. GBM and Allen MTG processing scripts are available at https://github.com/kahkengwong/CEP-IP_Framework

###############################################
# B. Processing Information with Seurat in R
###############################################
setwd("C:/...") # Set the local directory
library(dplyr)  
library(future)         
library(ggplot2)      
library(parallel)      
library(scales)        
library(scDblFinder) 
library(Seurat)          
library(SingleCellExperiment) 
library(viridis)  

# =========================================
# 1. scRNA-seq Dataset Pre-processing
# =========================================
# Load scRNA-seq dataset
loaded_df <- readRDS("C:/...directory.../GSE185344_PH_scRNA.final.rds")

# Extract Seurat object, the core data structure
seurat_obj <- loaded_df$obj

# Define sample names (tumor vs benign)
prostate_ca_samples <- c("HYW_4701_Tumor", "HYW_4847_Tumor", "HYW_4880_Tumor", 
                         "HYW_4881_Tumor", "HYW_5386_Tumor", "HYW_5742_Tumor", 
                         "HYW_5755_Tumor")
non_cancerous_samples <- c("HYW_4701_Benign", "HYW_4847_Benign", "HYW_4880_Benign", 
                           "HYW_4881_Benign", "HYW_5386_Benign", "HYW_5742_Benign", 
                           "HYW_5755_Benign")

# Subset and process Seurat object for normalization and feature selection
process_seurat <- function(seurat_obj, sample_names, project_name) {
    subset_obj <- subset(seurat_obj, subset = orig.ident %in% sample_names)  # Subset by sample ID
    subset_obj <- NormalizeData(subset_obj, normalization.method = "LogNormalize", scale.factor = 10000)  # Log-normalize (expression scaling)
    subset_obj <- FindVariableFeatures(subset_obj, selection.method = "vst", nfeatures = 2000)  # Top 2000 variable features (VST method)
    return(subset_obj)
}

# Process samples to split into tumor and benign
prostate_ca_seurat <- process_seurat(seurat_obj, prostate_ca_samples, "prostate-ca")
non_cancerous_seurat <- process_seurat(seurat_obj, non_cancerous_samples, "NonCancerous")

# Filter by feature and count thresholds to remove low-quality cells
prostate_ca_seurat <- subset(prostate_ca_seurat, subset = nFeature_RNA > 500 & nCount_RNA > 0)
non_cancerous_seurat <- subset(non_cancerous_seurat, subset = nFeature_RNA > 500 & nCount_RNA > 0)

# Filter high ribosomal content to mitigate bias from cells with highest expression of ribosomal genes (top 10th percentile)
filter_ribosomal <- function(seurat_obj, method = "fixed", cutoff = 10) {
    rp_genes <- grep("^RP[SL]|^MRP[SL]", rownames(seurat_obj), value = TRUE)  # Ribosomal genes (RP/MRP prefixes)
    seurat_obj[["percent.ribo"]] <- PercentageFeatureSet(seurat_obj, features = rp_genes)  # % ribosomal expression (per cell)
    threshold <- if (method == "percentile") quantile(seurat_obj$percent.ribo, probs = cutoff) else cutoff  # Threshold (percentile or fixed)
    plot <- ggplot(seurat_obj@meta.data, aes(x = percent.ribo)) +  # Plot distribution (with threshold line)
        geom_histogram(bins = 100) +
        geom_vline(xintercept = threshold, color = "red", linetype = "dashed") +
        ggtitle("Distribution of Ribosomal Gene Percentage")
    print(plot)
    genes_before <- nrow(seurat_obj); cells_before <- ncol(seurat_obj)  # Pre-filter counts (genes, cells)
    seurat_obj_filtered <- subset(seurat_obj, subset = percent.ribo < threshold)  # Apply filter (below threshold)
    genes_after <- nrow(seurat_obj_filtered); cells_after <- ncol(seurat_obj_filtered)  # Post-filter counts
    cat("Threshold:", threshold, "\nCells before:", cells_before, "\nCells after:", cells_after, 
        "\nRemoved:", round((cells_before - cells_after) / cells_before * 100, 2), "%\n",
        "Genes before:", genes_before, "\nGenes after:", genes_after, "\n")
    return(seurat_obj_filtered)
}

# Apply ribosomal filter (90th percentile cutoff)
cat("Filtering prostate cancer samples (ribosomal)\n")
prostate_ca_seurat <- filter_ribosomal(prostate_ca_seurat, method = "percentile", cutoff = 0.90)
cat("Filtering non-cancerous samples (ribosomal)\n")
non_cancerous_seurat <- filter_ribosomal(non_cancerous_seurat, method = "percentile", cutoff = 0.90)

# Filter high mitochondrial content to mitigate dying cells
filter_mitochondrial <- function(seurat_obj, method = "fixed", cutoff = 10) {
    mt_genes <- grep("^MT-", rownames(seurat_obj), value = TRUE)  # Mitochondrial genes (MT- prefix)
    seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, features = mt_genes)  # % mitochondrial expression (per cell)
    threshold <- if (method == "percentile") quantile(seurat_obj$percent.mt, probs = cutoff) else cutoff  # Threshold (percentile or fixed)
    plot <- ggplot(seurat_obj@meta.data, aes(x = percent.mt)) +  # Plot distribution (with threshold line)
        geom_histogram(bins = 100) +
        geom_vline(xintercept = threshold, color = "red", linetype = "dashed") +
        ggtitle("Distribution of Mitochondrial Gene Percentage")
    print(plot)
    genes_before <- nrow(seurat_obj); cells_before <- ncol(seurat_obj)  # Pre-filter counts (genes, cells)
    seurat_obj_filtered <- subset(seurat_obj, subset = percent.mt < threshold)  # Apply filter (below threshold)
    genes_after <- nrow(seurat_obj_filtered); cells_after <- ncol(seurat_obj_filtered)  # Post-filter counts
    cat("Threshold:", threshold, "\nCells before:", cells_before, "\nCells after:", cells_after, 
        "\nRemoved:", round((cells_before - cells_after) / cells_before * 100, 2), "%\n",
        "Genes before:", genes_before, "\nGenes after:", genes_after, "\n")
    return(seurat_obj_filtered)
}

# Apply mitochondrial filter (90th percentile cutoff)
cat("Filtering prostate cancer samples (mitochondrial)\n")
prostate_ca_seurat <- filter_mitochondrial(prostate_ca_seurat, method = "percentile", cutoff = 0.90)
cat("Filtering non-cancerous samples (mitochondrial)\n")
non_cancerous_seurat <- filter_mitochondrial(non_cancerous_seurat, method = "percentile", cutoff = 0.90)

# =========================================
# 2. Cell Cycle Regression
# =========================================
# Check pre-regression cell count (after the previous steps)
cat("Cells before cell cycle regression (prostate cancer):", ncol(prostate_ca_seurat), "\n")

# Default cell cycle genes (common S and G2M phase markers)
s_genes_default <- c("MCM5", "PCNA", "TYMS", "FEN1", "MCM2", "MCM4", "RRM1", "UNG", "GINS2", "MCM6", "CDCA7", "DTL", "PRIM1", "UHRF1", "MLF1IP", "HELLS", "RFC2", "RPA2", "NASP", "RAD51AP1", "GMNN", "WDR76", "SLBP", "CCNE2", "UBR7", "POLD3", "MSH2", "ATAD2", "RAD51", "RRM2", "CDC45", "CDC6", "EXO1", "TIPIN", "DSCC1", "BLM", "CASP8AP2", "USP1", "CLSPN", "POLA1", "CHAF1B", "BRIP1", "E2F8")
g2m_genes_default <- c("HMGB2", "CDK1", "NUSAP1", "UBE2C", "BIRC5", "TPX2", "TOP2A", "NDC80", "CKS2", "NUF2", "CKS1B", "MKI67", "TMPO", "CENPF", "TACC3", "FAM64A", "SMC4", "CCNB2", "CKAP2L", "CKAP2", "AURKB", "BUB1", "KIF11", "ANP32E", "TUBB4B", "GTSE1", "KIF20B", "HJURP", "CDCA3", "HN1", "CDC20", "TTK", "CDC25C", "KIF2C", "RANGAP1", "NCAPD2", "DLGAP5", "CDCA2", "CDCA8", "ECT2", "KIF23", "HMMR", "AURKA", "PSRC1", "ANLN", "LBR", "CKAP5", "CENPE", "CTCF", "NEK2", "G2E3", "GAS2L3", "CBX5", "CENPA")

# Score and regress out cell cycle to remove phase effects (prostate cancer)
prostate_ca_seurat <- CellCycleScoring(prostate_ca_seurat, s.features = s_genes_default, g2m.features = g2m_genes_default, set.ident = TRUE)
prostate_ca_seurat <- ScaleData(prostate_ca_seurat, vars.to.regress = c("S.Score", "G2M.Score"))

# Post-regression cell count to verify no cell loss
cat("Cells after cell cycle regression (prostate cancer):", ncol(prostate_ca_seurat), "\n") # 22796 (no changes)

# Score and regress out cell cycle (non-cancerous)
cat("Cells before cell cycle regression (non-cancerous):", ncol(non_cancerous_seurat), "\n")
non_cancerous_seurat <- CellCycleScoring(non_cancerous_seurat, s.features = s_genes_default, g2m.features = g2m_genes_default, set.ident = TRUE)
non_cancerous_seurat <- ScaleData(non_cancerous_seurat, vars.to.regress = c("S.Score", "G2M.Score"))

# =========================================
# 3. Doublets removal
# =========================================
# Set up parallel processing to speed up doublet detection
plan(multisession, workers = availableCores())

# Remove doublets using scDblFinder
remove_doublets <- function(seurat_obj) {
    sce_obj <- as.SingleCellExperiment(seurat_obj)  # Convert to SCE for scDblFinder
    samples <- seurat_obj@meta.data$orig.ident  # Batch info by sample IDs
    doublet_scores <- scDblFinder(sce_obj, samples = samples, k = 30, nfeatures = 2000)  # Run scDblFinder 
    batch_thresholds <- tapply(doublet_scores$scDblFinder.score, samples, function(x) quantile(x, probs = 0.95))  # Batch-specific thresholds (95th percentile)
    cat("Batch thresholds:\n"); print(batch_thresholds)
    doublet_cells <- colnames(sce_obj)[mapply(function(x, y) x > batch_thresholds[y], doublet_scores$scDblFinder.score, samples)]  # Identify doublets (above threshold)
    cat("Doublets:", length(doublet_cells), "\n")
    seurat_obj$doublet <- colnames(seurat_obj) %in% doublet_cells  # Mark doublets as TRUE/FALSE
    seurat_obj <- subset(seurat_obj, subset = doublet == FALSE)  # Remove doublets
    cat("Cells after removal:", ncol(seurat_obj), "\n")
    seurat_obj$filtered <- "filtered"  # Update metadata and flag filtered cells
    return(seurat_obj)
}

# Apply doublet removal to prostate ca and benign cases
cat("Processing prostate cancer samples\n")
prostate_ca_seurat <- remove_doublets(prostate_ca_seurat)
cat("Processing non-cancerous samples\n")
non_cancerous_seurat <- remove_doublets(non_cancerous_seurat)

# =========================================
# 4. Batch effects correction
# =========================================
# Disable parallel processing to avoid integration issues
plan(sequential)

# Correct batch effects by integrating across samples
correct_batch_effects <- function(seurat_obj) {
    cat("Metadata columns:\n"); print(colnames(seurat_obj@meta.data))
    cat("Unique orig.ident:\n"); print(unique(seurat_obj@meta.data$orig.ident))
    seurat_obj_before_integration <- FindVariableFeatures(seurat_obj, selection.method = "vst", nfeatures = 500)  # Pre-integration features (500 VST)
    seurat_obj_before_integration <- RunPCA(seurat_obj_before_integration, verbose = FALSE)  # PCA (dimensionality reduction)
    seurat_obj_before_integration <- RunUMAP(seurat_obj_before_integration, dims = 1:8, verbose = FALSE)  # UMAP (pre-integration visualization)
    plasma_colors <- viridis(n = length(unique(seurat_obj_before_integration$orig.ident)), option = "plasma")
    p1 <- DimPlot(seurat_obj_before_integration, group.by = "orig.ident", pt.size = 0.5, label = FALSE, repel = TRUE, cols = plasma_colors) + 
        ggtitle("UMAP Before Integration") + theme(legend.position = "right")  # Pre-integration UMAP (batch-colored)
    print(p1)
    sample_list <- SplitObject(seurat_obj_before_integration, split.by = "orig.ident")  # Split by batch (sample IDs)
    sample_list <- lapply(sample_list, function(x) {  # SCTransform and clean NAs (per sample)
        x <- SCTransform(x, verbose = FALSE, variable.features.n = 500, vst.flavor = "v2")
        x@meta.data <- x@meta.data[complete.cases(x@meta.data), ]
        x
    })
    anchors <- FindIntegrationAnchors(object.list = sample_list, dims = 1:5, verbose = FALSE)  # Find anchors (dims 1-5)
    seurat_obj_integrated <- IntegrateData(anchorset = anchors, dims = 1:5, verbose = FALSE)  # Integrate (batch-corrected)
    seurat_obj_integrated <- ScaleData(seurat_obj_integrated, verbose = FALSE)  # Scale (post-integration)
    seurat_obj_integrated <- RunPCA(seurat_obj_integrated, verbose = FALSE)  # PCA (integrated)
    seurat_obj_integrated <- FindNeighbors(seurat_obj_integrated, dims = 1:8)  # Neighbors (for clustering)
    seurat_obj_integrated <- FindClusters(seurat_obj_integrated, resolution = 0.5)  # Clusters (res 0.5)
    seurat_obj_integrated <- RunUMAP(seurat_obj_integrated, dims = 1:8, verbose = FALSE, umap.method = "uwot", metric = "cosine")  # UMAP (post-integration)
    plasma_colors <- viridis(n = length(unique(seurat_obj_integrated$orig.ident)), option = "plasma")
    p3 <- DimPlot(seurat_obj_integrated, group.by = "orig.ident", pt.size = 0.5, label = FALSE, repel = TRUE, cols = plasma_colors) + 
        ggtitle("UMAP After Integration") + theme(legend.position = "right")  # Post-integration UMAP (batch-colored)
    print(p3)
    p4 <- DimPlot(seurat_obj_integrated, group.by = "orig.ident", pt.size = 0.5, label = TRUE, repel = TRUE) + 
        ggtitle("UMAP After Integration") + theme(legend.position = "right")  # Labeled UMAP by batch IDs
    print(p4)
    return(seurat_obj_integrated)
}

# Remove parallelization limits to ensure stability
options(future.globals.maxSize = Inf)

# Apply batch correction for prostate ca and benign cases
cat("Processing prostate cancer samples (batch effects)\n")
prostate_ca_seurat_integrated <- correct_batch_effects(prostate_ca_seurat)
cat("Processing non-cancerous samples (batch effects)\n")
non_cancerous_seurat_integrated <- correct_batch_effects(non_cancerous_seurat)

# =========================================
# 5. UMAP Clusters
# =========================================
# Generate elbow plots to assess PCA dimensionality reduction
generate_elbow_plot <- function(seurat_obj_integrated, output_prefix) {
    seurat_obj_integrated <- RunPCA(seurat_obj_integrated, verbose = FALSE)  # PCA (dimensionality reduction)
    elbow_plot <- ElbowPlot(seurat_obj_integrated, ndims = 50) +  # Elbow plot
        labs(title = paste("Elbow Plot for", output_prefix), x = "Principal Components", y = "Standard Deviation") +
        theme(plot.title = element_text(size = 14, face = "bold"), axis.title.x = element_text(size = 12), 
              axis.title.y = element_text(size = 12), axis.text.x = element_text(size = 10), axis.text.y = element_text(size = 10))
    ggsave(paste0(output_prefix, "_ElbowPlot.pdf"), elbow_plot, width = 6.83, height = 6.41)
    print(elbow_plot)
}

# Generate elbow plots (prostate ca and benign)
generate_elbow_plot(prostate_ca_seurat_integrated, "prostate_ca")
generate_elbow_plot(non_cancerous_seurat_integrated, "non_cancerous")

# Downstream analyses with UMAP (clustering and markers)
downstream_analyses <- function(seurat_obj_integrated, gene_of_interest, output_prefix, dims = 15) {
    set.seed(10)
    DefaultAssay(seurat_obj_integrated) <- "RNA"
    seurat_obj_integrated <- FindVariableFeatures(seurat_obj_integrated, selection.method = "vst", nfeatures = 2000)  # Variable features (2000 VST)
    seurat_obj_integrated <- ScaleData(seurat_obj_integrated, verbose = FALSE)  # Scale (center and normalize)
    seurat_obj_integrated <- RunPCA(seurat_obj_integrated, verbose = FALSE)  # PCA (dims reduction)
    seurat_obj_integrated <- RunUMAP(seurat_obj_integrated, dims = 1:dims, verbose = FALSE)  # UMAP (dims 1-15)
    set.seed(11); seurat_obj_integrated <- FindNeighbors(seurat_obj_integrated, dims = 1:dims)  # Neighbors (kNN graph)
    set.seed(12); seurat_obj_integrated <- FindClusters(seurat_obj_integrated, resolution = 0.5)  # Clusters (Louvain, res 0.5)
    set.seed(13); cluster_markers <- FindAllMarkers(seurat_obj_integrated, only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.25)  # Marker genes (positive, logFC > 0.25)
    print(paste("Cluster markers:", nrow(cluster_markers)))
    if (nrow(cluster_markers) == 0) {
        print("Cluster levels:"); print(levels(Idents(seurat_obj_integrated)))
        print("Cells per cluster:"); print(table(Idents(seurat_obj_integrated)))
    }
    umap_data <- as.data.frame(Embeddings(seurat_obj_integrated, "umap")); umap_data$cluster_id <- Idents(seurat_obj_integrated)  # UMAP data (coords + clusters)
    umap_data_mean <- aggregate(. ~ cluster_id, data = umap_data, FUN = mean)  # Mean coords (per cluster)
    plasma_func <- colorRampPalette(viridis::viridis(100, direction = -1, option = "plasma")); portion <- 0.8  # Colors (plasma palette)
    n_colors <- round(length(unique(umap_data$cluster_id)) / portion); plasma_colors <- plasma_func(n_colors)
    set.seed(14); umap_plot_with_labels <- ggplot(umap_data, aes(x = umap_1, y = umap_2, color = as.factor(cluster_id))) +  # Labeled UMAP (cluster IDs)
        geom_point(size = 0.3, alpha = 0.5) + scale_color_manual(values = plasma_colors) +
        geom_text(data = umap_data_mean, aes(label = cluster_id, x = umap_1, y = umap_2), color = "black", size = 3, fontface = "bold", check_overlap = TRUE) +
        theme(panel.border = element_rect(fill = NA, color = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
              axis.text.x = element_text(color = "black"), axis.ticks.x = element_line(color = "black"), axis.text.y = element_text(color = "black"),
              axis.ticks.y = element_line(color = "black"), panel.background = element_rect(fill = "white")) +
        labs(title = "UMAP plot colored by cluster (with labels)", x = "umap_1", y = "umap_2", color = "Cluster") +
        guides(color = guide_legend(override.aes = list(size = 3)))
    print(umap_plot_with_labels)
    set.seed(15); umap_plot_no_labels <- ggplot(umap_data, aes(x = umap_1, y = umap_2, color = as.factor(cluster_id))) +  # Unlabeled UMAP (clusters only)
        geom_point(size = 0.3, alpha = 0.5) + scale_color_manual(values = plasma_colors) +
        theme(panel.border = element_rect(fill = NA, color = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
              axis.text.x = element_text(color = "black"), axis.ticks.x = element_line(color = "black"), axis.text.y = element_text(color = "black"),
              axis.ticks.y = element_line(color = "black"), panel.background = element_rect(fill = "white")) +
        labs(title = "UMAP plot colored by cluster (without labels)", x = "umap_1", y = "umap_2", color = "Cluster") +
        guides(color = guide_legend(override.aes = list(size = 3)))
    print(umap_plot_no_labels)
    if (gene_of_interest %in% rownames(seurat_obj_integrated)) {  # Gene expression UMAP (if gene exists)
        gene_colors_alpha <- c(scales::alpha("lightgray", 0.85), scales::alpha("lightpink", 0.85), scales::alpha("#FF6666", 0.85), 
                               scales::alpha("#BC2727", 0.85), scales::alpha("#660000", 0.85))
        set.seed(16); feature_plot <- FeaturePlot(seurat_obj_integrated, features = gene_of_interest, min.cutoff = 'q10', max.cutoff = 'q90',
                                                  pt.size = 0.2, cols = gene_colors_alpha) +
            theme(panel.border = element_rect(fill = NA, color = "black"), panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
                  axis.text.x = element_text(color = "black"), axis.ticks.x = element_line(color = "black"), axis.text.y = element_text(color = "black"),
                  axis.ticks.y = element_line(color = "black"), panel.background = element_rect(fill = "white")) +
            labs(title = paste("UMAP plot colored by", gene_of_interest, "expression"), x = "umap_1", y = "umap_2")
        print(feature_plot)
    } else {
        cat(paste("Warning: Gene", gene_of_interest, "not found.\nAvailable genes:\n"))
        print(head(rownames(seurat_obj_integrated), 20))
    }
    if (nrow(cluster_markers) > 0) {  # Top markers (50 per cluster)
        top_markers <- cluster_markers %>% group_by(cluster) %>% top_n(n = 50, wt = avg_log2FC)
    } else {
        top_markers <- data.frame()
        warning("No cluster markers found.")
    }
    write.table(top_markers, file = paste0(output_prefix, "_top_markers_for_each_cluster_vRibo.tsv"), sep = "\t", col.names = TRUE, row.names = TRUE, quote = FALSE)  # Save markers (TSV)
    return(list(seurat_obj = seurat_obj_integrated, cluster_markers = cluster_markers, top_markers = top_markers))
}

# Apply downstream analyses (TRPM4 focus)
set.seed(42)
prostate_results <- downstream_analyses(prostate_ca_seurat_integrated, "TRPM4", "prostate_ca", dims = 15)
non_cancerous_results <- downstream_analyses(non_cancerous_seurat_integrated, "TRPM4", "non_cancerous", dims = 15)

# Save workspace
save.image(file = "GSE185344_Seurat_processed.RData")

# For subsequent analysis, load the saved file
load("GSE185344_Seurat_processed.RData")


#########################################################
# C. Subsequent GAM-REML-TPRS Analysis Code Availability
#########################################################
The subsequent analysis pipeline of this project is available on GitHub: https://github.com/kahkengwong/CEP-IP_Framework
```