Cluster-wise Stability Assessment of Glioblastoma scRNA-seq Data
Andi Rosilala
2026-02-21
# Kelompok Bioconductor lengkap — install sekaligus
BiocManager::install(
c("AnnotationDbi",
"clusterProfiler",
"enrichplot",
"DOSE",
"msigdbr"
),
type = "binary",
update = FALSE,
ask = FALSE
)pkgs <- c("GenomeInfoDbData", "AnnotationDbi", "clusterProfiler",
"org.Hs.eg.db", "enrichplot", "msigdbr")
for (pkg in pkgs) {
ok <- file.exists(system.file(package = pkg))
cat(ifelse(ok, "✅", "❌"), pkg, "\n")
}## ✅ GenomeInfoDbData
## ✅ AnnotationDbi
## ✅ clusterProfiler
## ✅ org.Hs.eg.db
## ✅ enrichplot
## ✅ msigdbr# Perpanjang timeout dari 60 detik jadi 300 detik
options(timeout = 300)
options(pkgType = "binary")
options(install.packages.compile.from.source = "never")
# Install satu per satu agar kalau gagal mudah diulang
pkgs_to_install <- c(
"SeuratObject",
"Seurat",
"tidyverse",
"patchwork",
"ggplot2",
"fpc",
"cluster",
"igraph",
"ggridges",
"viridis",
"RColorBrewer",
"pheatmap",
"cowplot",
"Matrix",
"data.table",
"future",
"future.apply",
"kableExtra",
"knit",
"rmarkdown"
)
for (pkg in pkgs_to_install) {
if (!file.exists(system.file(package = pkg))) {
cat("Installing:", pkg, "\n")
tryCatch(
install.packages(pkg, type = "binary"),
error = function(e) cat("❌ GAGAL:", pkg, "-", e$message, "\n")
)
} else {
cat("✅ Sudah ada:", pkg, "\n")
}
}## ✅ Sudah ada: SeuratObject
## ✅ Sudah ada: Seurat
## ✅ Sudah ada: tidyverse
## ✅ Sudah ada: patchwork
## ✅ Sudah ada: ggplot2
## ✅ Sudah ada: fpc
## ✅ Sudah ada: cluster
## ✅ Sudah ada: igraph
## ✅ Sudah ada: ggridges
## ✅ Sudah ada: viridis
## ✅ Sudah ada: RColorBrewer
## ✅ Sudah ada: pheatmap
## ✅ Sudah ada: cowplot
## ✅ Sudah ada: Matrix
## ✅ Sudah ada: data.table
## ✅ Sudah ada: future
## ✅ Sudah ada: future.apply
## ✅ Sudah ada: kableExtra
## Installing: knit## ❌ GAGAL: knit - trying to use CRAN without setting a mirror
## ✅ Sudah ada: rmarkdown# Install AnnotationDbi dulu
BiocManager::install(
"AnnotationDbi",
type = "binary",
ask = FALSE,
update = FALSE
)
# Install GenomeInfoDbData dari file lokal
install.packages(
"C:/Users/Lenovo/Downloads/GenomeInfoDbData_1.2.11.tar.gz",
repos = NULL,
type = "source"
)
# Install dari file lokal — tidak butuh internet
install.packages(
"C:/Users/Lenovo/Downloads/org.Hs.eg.db_3.18.0.tar.gz",
repos = NULL, # tidak pakai internet
type = "source" # install dari file tar.gz
)pkgs <- c(
"SeuratObject", "Seurat", "tidyverse", "patchwork",
"fpc", "cluster", "igraph", "ggridges", "viridis",
"RColorBrewer", "pheatmap", "cowplot", "Matrix",
"data.table", "future", "future.apply",
"kableExtra", "knitr", "rmarkdown"
)
for (pkg in pkgs) {
ok <- file.exists(system.file(package = pkg))
cat(ifelse(ok, "✅", "❌"), pkg, "\n")
}## ✅ SeuratObject
## ✅ Seurat
## ✅ tidyverse
## ✅ patchwork
## ✅ fpc
## ✅ cluster
## ✅ igraph
## ✅ ggridges
## ✅ viridis
## ✅ RColorBrewer
## ✅ pheatmap
## ✅ cowplot
## ✅ Matrix
## ✅ data.table
## ✅ future
## ✅ future.apply
## ✅ kableExtra
## ✅ knitr
## ✅ rmarkdownOverview
This R Markdown script performs comprehensive cluster stability assessment on glioblastoma single-cell RNA-seq data using bootstrap resampling and Jaccard coefficient, following the methodology of Hennig (2007).
Key Steps:
- Data loading and preprocessing
- Quality control
- Normalization and feature selection
- Dimensionality reduction (PCA, UMAP)
- Graph-based clustering
- Cluster-wise stability assessment using bootstrap
- Biological characterization
- Visualization and reporting
1. Load Required Libraries
# Core Seurat and tidyverse
library(SeuratObject)
library(Seurat)
library(tidyverse)
library(patchwork)
# For clustering and stability
library(fpc)
library(cluster)
library(igraph)
# For visualization
library(ggplot2)
library(ggridges)
library(viridis)
library(RColorBrewer)
library(pheatmap)
library(cowplot)
# For pathway analysis
library(clusterProfiler)
library(org.Hs.eg.db)
library(enrichplot)
library(msigdbr)
# For data manipulation
library(Matrix)
library(data.table)
# For parallel processing
library(future)
library(future.apply)
# Set future global options for parallel processing
plan("multicore", workers = 8) # Adjust based on your CPU cores
options(future.globals.maxSize = 8000 * 1024^2) # 8GB
# Set seed for reproducibility
set.seed(123)
# Custom theme for plots
theme_publication <- theme_bw() +
theme(
axis.text = element_text(size = 12),
axis.title = element_text(size = 14, face = "bold"),
plot.title = element_text(size = 16, face = "bold", hjust = 0.5),
legend.text = element_text(size = 11),
legend.title = element_text(size = 12, face = "bold"),
strip.text = element_text(size = 12, face = "bold")
)2. Load and Prepare Data
setwd("D:/GBM_Cluster")
data_path <- "./GSE131928/"
expr_10x <- read.delim(
paste0(data_path, "GSM3828673_10X_GBM_IDHwt_processed_TPM.tsv"),
row.names = 1,
check.names = FALSE
)
expr_10x <- as.matrix(expr_10x)
expr_ss2 <- read.delim(
paste0(data_path, "GSM3828672_Smartseq2_GBM_IDHwt_processed_TPM.tsv"),
row.names = 1,
check.names = FALSE
)
expr_ss2 <- as.matrix(expr_ss2)
meta_10x <- data.frame(
cell = colnames(expr_10x),
platform = "10X",
disease = "GBM",
IDH_status = "IDHwt"
)
rownames(meta_10x) <- meta_10x$cell
meta_ss2 <- data.frame(
cell = colnames(expr_ss2),
platform = "SmartSeq2",
disease = "GBM",
IDH_status = "IDHwt"
)
rownames(meta_ss2) <- meta_ss2$cell
gbm_10x <- CreateSeuratObject(
counts = expr_10x,
project = "GBM_10X",
min.cells = 3,
min.features = 50,
meta.data = meta_10x
)
gbm_ss2 <- CreateSeuratObject(
counts = expr_ss2,
project = "GBM_SmartSeq2",
min.cells = 3,
min.features = 50,
meta.data = meta_ss2
)
gbm_10x@assays$RNA@misc$data_type <- "TPM"
gbm_ss2@assays$RNA@misc$data_type <- "TPM"
gbm_10x## An object of class Seurat
## 22185 features across 16201 samples within 1 assay
## Active assay: RNA (22185 features, 0 variable features)
## 1 layer present: counts## An object of class Seurat
## 22197 features across 7930 samples within 1 assay
## Active assay: RNA (22197 features, 0 variable features)
## 1 layer present: counts##
## 10X
## 16201##
## SmartSeq2
## 79303. Quality Control
DefaultAssay(gbm_10x) <- "RNA"
# Hitung persentase gen mitokondria
gbm_10x[["percent.mt"]] <- PercentageFeatureSet(
gbm_10x,
pattern = "^MT-"
)
# Hitung persentase gen ribosomal (opsional)
gbm_10x[["percent.ribo"]] <- PercentageFeatureSet(
gbm_10x,
pattern = "^RP[SL]"
)
VlnPlot(
gbm_10x,
features = c("nFeature_RNA", "percent.mt"),
ncol = 2,
pt.size = 0.1
)
# Scatter: kompleksitas vs mitokondria
FeatureScatter(
gbm_10x,
feature1 = "nFeature_RNA",
feature2 = "percent.mt"
) +
geom_hline(yintercept = 25, linetype = "dashed", color = "red") +
geom_vline(xintercept = 200, linetype = "dashed", color = "red")
##
## === QC Summary (TPM-based, 10X) ===## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 154 854 1798 2184 3069 9204## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.525 2.924 3.554 4.639 78.571DefaultAssay(gbm_ss2) <- "RNA"
# Persentase gen mitokondria
gbm_ss2[["percent.mt"]] <- PercentageFeatureSet(
gbm_ss2,
pattern = "^MT-"
)
# Persentase gen ribosomal
gbm_ss2[["percent.ribo"]] <- PercentageFeatureSet(
gbm_ss2,
pattern = "^RP[SL]"
)
# Violin plot QC
VlnPlot(
gbm_ss2,
features = c("nFeature_RNA", "percent.mt"),
ncol = 2,
pt.size = 0.1
)
# Scatter: kompleksitas vs mitokondria
FeatureScatter(
gbm_ss2,
feature1 = "nFeature_RNA",
feature2 = "percent.mt"
) +
geom_hline(yintercept = 20, linetype = "dashed", color = "red") +
geom_vline(xintercept = 200, linetype = "dashed", color = "red")
##
## === QC Summary (TPM-based, Smart-seq2) ===## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 3003 4639 5660 5735 6752 11890## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0 0 0 0 0 0# Apply QC filters
gbm_10x <- subset(
gbm_10x,
subset =
nFeature_RNA > 200 & # remove low-complexity cells
percent.mt < 25 # remove stressed / dying cells
)
cat("\n=== After QC Filtering (10X, TPM-based) ===\n")##
## === After QC Filtering (10X, TPM-based) ===## Retained 15072 cells
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 201 1046 1914 2334 3211 9204## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 1.596 2.952 3.491 4.612 24.9414. Normalization and Feature Selection
# Normalize data using LogNormalize
gbm_10x <- NormalizeData(
object = gbm_10x,
normalization.method = "LogNormalize",
scale.factor = 10000,
verbose = TRUE
)
cat("✔ Normalization completed (10X)\n")## ✔ Normalization completed (10X)## Data slot after normalization:## [1] "layers" "cells" "features" "default" "assay.orig"
## [6] "meta.data" "misc" "key"# Identify highly variable features
gbm_10x <- FindVariableFeatures(
object = gbm_10x,
selection.method = "vst",
nfeatures = 2000,
verbose = T
)
# Check number of HVGs
cat("Number of variable features:",
length(VariableFeatures(gbm_10x)), "\n")## Number of variable features: 2000# Top 10 variable genes
top10 <- head(VariableFeatures(gbm_10x), 10)
# Plot variable features
hvg_plot <- VariableFeaturePlot(gbm_10x) +
theme_publication +
ggtitle("Highly Variable Genes (10X, TPM-based)")
hvg_plot_labeled <- LabelPoints(
plot = hvg_plot,
points = top10,
repel = TRUE
)
print(hvg_plot_labeled)
##
## Top 10 highly variable genes:## [1] "HBB" "HBA2" "HBA1" "TPSAB1" "SAA1" "SLPI" "PLA2G2A"
## [8] "SAA2" "PI3" "MT1H"5. Dimensionality Reduction
gbm_10x <- ScaleData(
object = gbm_10x,
features = VariableFeatures(gbm_10x),
verbose = T
)
cat("ScaleData completed (10X, TPM-based)\n")## ScaleData completed (10X, TPM-based)graphics.off()
# Run PCA on 10X data
gbm_10x <- RunPCA(
object = gbm_10x,
features = VariableFeatures(gbm_10x),
npcs = 50,
verbose = T
)
cat("PCA completed (10X, TPM-based)\n")## PCA completed (10X, TPM-based)## PC_ 1
## Positive: C1orf61, APLP1, CNP, FABP7, S100B, MT3, CLDND1, CLU, PLP1, SEC61G
## PPP1R14A, CRYAB, OLIG1, SEPT4, SOX4, QDPR, CA2, TUBB4A, AMER2, ABCA2
## Negative: SRGN, C1QC, C1QA, CD74, C1QB, HLA-DRB1, HLA-DRA, AIF1, RGS1, HLA-DPB1
## HLA-DPA1, ALOX5AP, CD14, VSIG4, HLA-DMA, APOC1, HLA-DRB5, FTL, HLA-DMB, PLAUR
## PC_ 2
## Positive: SEC61G, CLU, S100A16, FABP7, TIMP1, OCIAD2, CSRP2, HOPX, LGALS3, CHI3L1
## SPOCD1, MT2A, IGFBP3, COL6A2, SPARCL1, IGFBP7, METTL7B, C1orf61, FKBP10, MT1X
## Negative: MAG, ERMN, MOBP, MOG, NKX6-2, EDIL3, CNDP1, TMEM144, CAPN3, ENPP2
## PPP1R14A, TF, CLDN11, CARNS1, PLP1, MAL, CNTN2, KLK6, HAPLN2, SPOCK3
## PC_ 3
## Positive: CHI3L1, TIMP1, MT2A, CHI3L2, CLU, MGST1, C1R, SLPI, LGALS3, LTF
## SOD2, SERPING1, MT1X, S100A6, CRYAB, C8orf4, SAA1, C1S, EFEMP1, MGP
## Negative: SOX4, SOX11, STMN2, NNAT, STMN4, DLX5, HES6, GADD45G, DLL3, DLX6-AS1
## B4GALNT1, DLX2, HIST1H4C, CCT2, RAB3IP, DCTN2, DTX3, KIF5A, NUP107, YEATS4
## PC_ 4
## Positive: RPS17, UQCR11, STMN2, SLPI, LTF, SAA1, DLX5, STMN4, NNAT, DLX6-AS1
## CHI3L2, AC009501.4, PSMA2, ATXN8OS, DLX2, MTRNR2L12, SAA2, MGP, MTRNR2L1, CRYAB
## Negative: CST3, FTH1, HOPX, IGFBP7, PTTG1, POSTN, APOE, CENPF, BIRC5, UBE2C
## NUSAP1, TOP2A, FABP7, AURKA, CDK1, CDKN3, CDC20, COL6A2, IGFBP3, SGOL2
## PC_ 5
## Positive: AQP1, FAM107A, NTRK2, GJA1, OS9, DTX3, GFAP, ITM2C, SLC14A1, APOE
## TSPAN31, RP11-620J15.3, C2orf40, MDM2, AGAP2-AS1, NUP107, AQP4, AEBP1, BCAN, ATP1A2
## Negative: UBE2C, BIRC5, CDKN3, CENPF, TOP2A, NUSAP1, AURKA, PTTG1, CDC20, PBK
## CDK1, POSTN, CCNB1, SEC61G, HJURP, SGOL1, IGFBP3, UBE2T, SGOL2, COL6A2pca_plot1 <- DimPlot(
gbm_10x,
reduction = "pca"
) +
theme_publication +
ggtitle("PCA Visualization (10X, TPM-based)")
print(pca_plot1)
pca_plot2 <- VizDimLoadings(
gbm_10x,
dims = 1:2,
reduction = "pca",
ncol = 2
) +
theme_publication
print(pca_plot2)# Elbow plot to determine number of PCs
elbow <- ElbowPlot(
object = gbm_10x,
ndims = 50
) +
theme_publication +
ggtitle("Elbow Plot for PC Selection (10X, TPM-based)") +
geom_vline(xintercept = 30, color = "red", linetype = "dashed")
print(elbow)
# Use top PCs based on elbow plot
n_pcs <- 30
cat(
"\nUsing", n_pcs,
"principal components for downstream analysis (10X)\n"
)##
## Using 30 principal components for downstream analysis (10X)# Run UMAP
graphics.off()
gbm_10x <- RunUMAP(
object = gbm_10x,
dims = 1:n_pcs,
n.neighbors = 50,
min.dist = 0.3,
metric = "cosine",
verbose = T
)
cat("✔ UMAP completed\n")## ✔ UMAP completedpdf("UMAP_clusters_10X.pdf", width = 8, height = 6)
DimPlot(
object = gbm_10x,
reduction = "umap",
label = TRUE,
repel = TRUE,
pt.size = 0.3
)
dev.off()## pdf
## 2pdf("UMAP_QC_10X.pdf", width = 10, height = 4)
FeaturePlot(gbm_10x, features = "nFeature_RNA", reduction = "umap")
FeaturePlot(gbm_10x, features = "nCount_RNA", reduction = "umap")
FeaturePlot(gbm_10x, features = "percent.mt", reduction = "umap")
dev.off()## pdf
## 26. Graph-based Clustering
# Build KNN graph
gbm_10x <- FindNeighbors(
gbm_10x,
dims = 1:n_pcs,
k.param = 50,
verbose = T
)
cat("KNN graph constructed\n")## KNN graph constructed# Test multiple resolutions
resolutions <- c(0.1, 0.15, 0.2, 0.25, 0.3)
for (res in resolutions) {
gbm_10x <- FindClusters(
gbm_10x,
resolution = res,
algorithm = 1, # Louvain
verbose = T
)
}## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 15072
## Number of edges: 1340579
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9778
## Number of communities: 14
## Elapsed time: 3 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 15072
## Number of edges: 1340579
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9720
## Number of communities: 15
## Elapsed time: 3 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 15072
## Number of edges: 1340579
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9667
## Number of communities: 17
## Elapsed time: 3 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 15072
## Number of edges: 1340579
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9624
## Number of communities: 17
## Elapsed time: 5 seconds
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 15072
## Number of edges: 1340579
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9582
## Number of communities: 17
## Elapsed time: 8 seconds## ✔ Clustering completed for resolutions: 0.1, 0.15, 0.2, 0.25, 0.3# Visualize clustering at different resolutions
pdf("UMAP_clustering_resolutions_10X.pdf", width = 12, height = 8)
for (res in resolutions) {
col_name <- paste0("RNA_snn_res.", res)
p <- DimPlot(
object = gbm_10x,
reduction = "umap",
group.by = col_name,
label = TRUE
) +
ggtitle(paste("Resolution =", res)) +
theme_publication +
NoLegend()
print(p)
}
dev.off()## png
## 2# Select optimal resolution (0.1 based on biological interpretation)
optimal_resolution <- 0.1
Idents(gbm_10x) <- paste0("RNA_snn_res.", optimal_resolution)
# Simpan identitas cluster final ke kolom standar
gbm_10x$seurat_clusters <- Idents(gbm_10x)
cat("\nSelected resolution:", optimal_resolution, "\n")##
## Selected resolution: 0.1## Number of clusters: 14##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 3997 2176 1881 1590 1548 1228 749 641 420 252 207 148 125 110# Visualize final clustering
pdf(
file = "UMAP_final_clustering_res0.1_10X.pdf",
width = 8,
height = 6
)
DimPlot(
object = gbm_10x,
reduction = "umap",
label = TRUE,
label.size = 5,
pt.size = 0.4
) +
theme_publication +
ggtitle(
paste("Final Clustering (Resolution =", optimal_resolution, ")")
)
dev.off()## png
## 27. Find Cluster Markers
# Find markers for all clusters
DefaultAssay(gbm_10x) <- "RNA"
# Install presto dulu
# devtools::install_github("immunogenomics/presto")
library(presto)
Idents(gbm_10x) <- "seurat_clusters"
table(Idents(gbm_10x))##
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 3997 2176 1881 1590 1548 1228 749 641 420 252 207 148 125 110expr_mat <- GetAssayData(
gbm_10x,
assay = "RNA",
layer = "data"
)
# lebih cepat!
cluster_labels <- Idents(gbm_10x)
all_markers <- wilcoxauc(
X = expr_mat,
y = cluster_labels
)
tail(all_markers)## [1] "feature" "group" "avgExpr" "logFC" "statistic" "auc"
## [7] "pval" "padj" "pct_in" "pct_out"sig_markers <- all_markers %>%
filter(
padj < 0.05, # adjusted p-value
logFC > 0.5, # effect size
pct_in > 0.25, # prevalensi dalam cluster
auc > 0.6 # separabilitas
) %>%
arrange(group, desc(auc))
top10_markers <- sig_markers %>%
group_by(group) %>%
slice_head(n = 10)
print(top10_markers)## # A tibble: 140 × 10
## # Groups: group [14]
## feature group avgExpr logFC statistic auc pval padj pct_in pct_out
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 RNASET2 0 2.78 2.52 43298713 0.978 0 0 97.8 29.6
## 2 CYBA 0 2.92 2.71 43270426 0.977 0 0 98.0 16.9
## 3 GPX1 0 3.32 2.01 41730086. 0.943 0 0 99.3 74.0
## 4 LAPTM5 0 2.66 2.27 41361636 0.934 0 0 97.6 22.1
## 5 TYROBP 0 3.30 2.73 40490626. 0.915 0 0 99.1 26.2
## 6 SRGN 0 3.31 2.66 40487207 0.915 0 0 98.0 30.8
## 7 C1QC 0 3.60 2.95 40383725 0.912 0 0 95.1 32.7
## 8 CD74 0 4.55 3.06 40352537 0.912 0 0 99.5 65.7
## 9 HLA-DRB1 0 3.38 2.62 40300513 0.910 0 0 97.4 35.5
## 10 RGS1 0 2.84 2.38 40179986. 0.908 0 0 94.4 24.1
## # ℹ 130 more rows# Save markers
write.csv(
sig_markers,
file = "cluster_markers_all_res0.1_presto.csv",
row.names = FALSE
)
write.csv(
top10_markers,
file = "cluster_markers_top10_res0.1_presto.csv",
row.names = FALSE
)# Heatmap of top markers
library(dplyr)
top5_markers <- sig_markers %>%
group_by(group) %>%
slice_head(n = 5)
# Pastikan gen unik
marker_genes <- unique(top5_markers$feature)
# Gunakan assay RNA (data sudah dinormalisasi)
DefaultAssay(gbm_10x) <- "RNA"
heatmap_plot <- DoHeatmap(
object = gbm_10x,
features = marker_genes,
group.by = "seurat_clusters",
size = 3,
angle = 0
) +
scale_fill_viridis_c(option = "magma") +
theme(
axis.text.y = element_text(size = 6),
axis.text.x = element_text(size = 8)
) +
ggtitle("Top 5 Discriminative Markers per Cluster (presto)")
print(heatmap_plot)
# Dot plot of selected markers
dotplot <- DotPlot(
object = gbm_10x,
features = marker_genes,
group.by = "seurat_clusters",
assay = "RNA"
) +
theme_publication +
theme(
axis.text.x = element_text(
angle = 90,
hjust = 1,
vjust = 0.5,
size = 8
),
axis.text.y = element_text(size = 8)
) +
ggtitle("Expression of Top Discriminative Markers per Cluster")
print(dotplot)
8. Annotate Cell States
# Define cellular state signatures based on literature (Neftel et al. 2019)
# MES-like signature
mes_signature <- c(
"CHI3L1", "LGALS1", "CD44", "VIM", "FN1", "COL1A1", "COL1A2",
"COL3A1", "COL5A1", "MMP2", "MMP9", "IL6", "TNF", "CXCL8",
"SERPINE1", "TIMP1"
)
# OPC-like signature
opc_signature <- c(
"PDGFRA", "OLIG1", "OLIG2", "SOX10", "OMG", "MBP", "PLP1",
"MAG", "MOG", "CSPG4"
)
# AC-like signature
ac_signature <- c(
"AQP4", "GFAP", "S100B", "ALDH1L1", "SLC1A2", "SLC1A3",
"GLUL", "GJA1", "CD44"
)
# NPC-like signature
npc_signature <- c(
"SOX2", "NES", "DCX", "NEUROD1", "ASCL1", "HES1", "ID1", "ID3"
)
# Cycling signature
cycling_signature <- c(
"MKI67", "TOP2A", "PCNA", "CCNB1", "CDK1", "UBE2C",
"CENPE", "AURKA", "AURKB"
)
# Hypoxic signature
hypoxic_signature <- c(
"HIF1A", "VEGFA", "BNIP3", "LDHA", "PGK1", "ENO1",
"PDK1", "CA9", "SLC2A1"
)# Calculate module scores for GBM cellular states
DefaultAssay(gbm_10x) <- "RNA"
# MES-like
gbm_10x <- AddModuleScore(
object = gbm_10x,
features = list(mes_signature),
name = "MES"
)
# OPC-like
gbm_10x <- AddModuleScore(
object = gbm_10x,
features = list(opc_signature),
name = "OPC"
)
# AC-like
gbm_10x <- AddModuleScore(
object = gbm_10x,
features = list(ac_signature),
name = "AC"
)
# NPC-like
gbm_10x <- AddModuleScore(
object = gbm_10x,
features = list(npc_signature),
name = "NPC"
)
# Cycling
gbm_10x <- AddModuleScore(
object = gbm_10x,
features = list(cycling_signature),
name = "Cycling"
)
# Hypoxic
gbm_10x <- AddModuleScore(
object = gbm_10x,
features = list(hypoxic_signature),
name = "Hypoxic"
)
# Visualize module scores
score_features <- c(
"MES1", "OPC1", "AC1",
"NPC1", "Cycling1", "Hypoxic1"
)
score_umap <- FeaturePlot(
gbm_10x,
features = score_features,
reduction = "umap",
ncol = 3,
order = TRUE
) &
scale_color_viridis_c(option = "magma") &
theme_publication
print(score_umap)
# Violin plots of module scores by cluster
score_violin <- VlnPlot(
gbm_10x,
features = score_features,
group.by = "seurat_clusters",
ncol = 3,
pt.size = 0
) &
theme_publication
print(score_violin)
# Annotate clusters based on marker expression and module scores
# This is a manual annotation based on the top markers and scores
# Adjust based on your actual data
# ── Anotasi terkoreksi berdasarkan marker genes aktual (menggantikan anotasi lama) ──
# Justifikasi lengkap ada di Section 18
cluster_annotations <- c(
"0" = "Microglia_Macrophage",
"1" = "MES_AC_Transitional",
"2" = "NPC-like_mature",
"3" = "MES-like",
"4" = "NPC-like",
"5" = "Macrophage_complement",
"6" = "AC-like",
"7" = "Cycling_Tumor",
"8" = "Oligodendrocyte",
"9" = "Endothelial_Astrocyte",
"10" = "T_cell",
"11" = "Erythrocyte",
"12" = "Macrophage_LAM",
"13" = "Pericyte_Stromal"
)
# Apply annotations — via @meta.data untuk menghindari masalah factor
cluster_vec <- as.character(gbm_10x@meta.data$seurat_clusters)
gbm_10x@meta.data$cell_state <- cluster_annotations[cluster_vec]
table(gbm_10x$cell_state, useNA = "ifany")##
## AC-like Cycling_Tumor Endothelial_Astrocyte
## 749 641 252
## Erythrocyte Macrophage_complement Macrophage_LAM
## 148 1228 125
## MES-like MES_AC_Transitional Microglia_Macrophage
## 1590 2176 3997
## NPC-like NPC-like_mature Oligodendrocyte
## 1548 1881 420
## Pericyte_Stromal T_cell
## 110 207# Visualize annotated clusters
annotated_plot <- DimPlot(
gbm_10x,
reduction = "umap",
group.by = "cell_state",
label = TRUE,
label.size = 5,
pt.size = 0.5,
repel = TRUE
) +
theme_publication +
ggtitle("Annotated Cell States")
print(annotated_plot)
##
## AC-like Cycling_Tumor Endothelial_Astrocyte
## 749 641 252
## Erythrocyte Macrophage_complement Macrophage_LAM
## 148 1228 125
## MES-like MES_AC_Transitional Microglia_Macrophage
## 1590 2176 3997
## NPC-like NPC-like_mature Oligodendrocyte
## 1548 1881 420
## Pericyte_Stromal T_cell
## 110 2079. Cluster Stability Assessment - Core Functions
##
## ========================================## RE-DEFINING STABILITY FUNCTIONS## ========================================# Load required libraries upfront
library(FNN)
library(igraph)
# -------------------------------------------------------
# Fungsi Louvain clustering yang benar untuk bootstrap
# Menerima data matrix (cells x features), sudah scaled/PCA
# k disesuaikan secara otomatis agar tidak melebihi n_cells-1
# -------------------------------------------------------
louvain_cluster_seurat <- function(data_matrix, k = 20, resolution = 0.1, dims = 1:20) {
n_cells_tmp <- nrow(data_matrix)
# Adaptif k: pastikan k < jumlah sel
k_actual <- min(k, n_cells_tmp - 1)
# Batasi dims agar tidak melebihi kolom tersedia
max_dims <- min(max(dims), ncol(data_matrix) - 1)
dims_actual <- 1:min(max(dims), max_dims)
# Run PCA (center=FALSE karena sudah scaled)
pca_result <- prcomp(data_matrix, center = FALSE, scale. = FALSE,
rank. = max(dims_actual))
# Extract PC embeddings
pca_embeddings <- pca_result$x[, dims_actual, drop = FALSE]
# Get k nearest neighbors
knn_result <- FNN::get.knn(pca_embeddings, k = k_actual)
# Build graph
edges <- data.frame(
from = rep(1:n_cells_tmp, each = k_actual),
to = as.vector(t(knn_result$nn.index))
)
# Hapus self-loops dan duplikat
edges <- edges[edges$from != edges$to, ]
g <- igraph::graph_from_data_frame(edges, directed = FALSE,
vertices = data.frame(name = 1:n_cells_tmp))
g <- igraph::simplify(g)
# Run Louvain clustering
louvain_result <- igraph::cluster_louvain(g, resolution = resolution)
return(as.integer(igraph::membership(louvain_result)))
}
# -------------------------------------------------------
# Fungsi Jaccard coefficient
# -------------------------------------------------------
compute_jaccard <- function(set1, set2) {
intersection <- sum(set1 & set2)
union <- sum(set1 | set2)
if (union == 0) return(NA_real_)
return(intersection / union)
}
# -------------------------------------------------------
# Fungsi bootstrap stability (Hennig 2007) yang BENAR
#
# Pendekatan:
# 1. Ambil bootstrap sample (dengan replacement)
# 2. Ambil sel UNIK dari sample (untuk menghindari duplikat)
# 3. Clustering pada sel unik tersebut
# 4. Hitung Jaccard antara label asli vs label bootstrap
# hanya pada irisan sel yang ada di kedua set
# -------------------------------------------------------
compute_bootstrap_stability_once <- function(
pca_embeddings, # matrix: n_cells x n_pcs (dari data asli)
original_labels, # integer vector panjang n_cells
unique_clusters, # cluster IDs yang ada
k = 20,
resolution = 0.1
) {
n_cells_orig <- nrow(pca_embeddings)
# 1. Bootstrap sample
boot_idx <- sample(n_cells_orig, n_cells_orig, replace = TRUE)
# 2. Ambil sel unik saja (menghapus duplikat untuk clustering yang bersih)
unique_boot_idx <- unique(boot_idx) # indeks sel unik dalam sampel
boot_data <- pca_embeddings[unique_boot_idx, , drop = FALSE]
n_boot <- nrow(boot_data)
if (n_boot < 10) return(rep(NA_real_, length(unique_clusters)))
k_actual <- min(k, n_boot - 1)
# 3. KNN graph pada bootstrap data
knn_result <- tryCatch(
FNN::get.knn(boot_data, k = k_actual),
error = function(e) NULL
)
if (is.null(knn_result)) return(rep(NA_real_, length(unique_clusters)))
edges <- data.frame(
from = rep(1:n_boot, each = k_actual),
to = as.vector(t(knn_result$nn.index))
)
edges <- edges[edges$from != edges$to, ]
g <- igraph::graph_from_data_frame(edges, directed = FALSE,
vertices = data.frame(name = 1:n_boot))
g <- igraph::simplify(g)
louvain_result <- tryCatch(
igraph::cluster_louvain(g, resolution = resolution),
error = function(e) NULL
)
if (is.null(louvain_result)) return(rep(NA_real_, length(unique_clusters)))
boot_labels <- as.integer(igraph::membership(louvain_result))
# 4. Label matching: Hungarian-style maximum overlap assignment
# Untuk setiap original cluster, cari boot cluster yang paling overlap
orig_labels_boot <- original_labels[unique_boot_idx] # label asli untuk sel di bootstrap
# Buat contingency table
tab <- table(orig_labels_boot, boot_labels)
# Untuk setiap original cluster, cari boot cluster dengan overlap terbesar (greedy)
# dan hitung Jaccard
n_orig_clusters <- length(unique_clusters)
jaccard_vec <- numeric(n_orig_clusters)
for (ci in seq_along(unique_clusters)) {
cl <- unique_clusters[ci]
# Sel dari cluster ini yang ada di bootstrap (unique_boot_idx)
in_orig <- (orig_labels_boot == cl)
if (sum(in_orig) == 0) {
jaccard_vec[ci] <- NA_real_
next
}
# Distribusi boot labels untuk sel ini
boot_labels_for_cl <- boot_labels[in_orig]
# Boot cluster dengan frekuensi terbanyak (best match)
best_boot_cl <- as.integer(names(which.max(table(boot_labels_for_cl))))
# Sel yang bootstrap-assigned ke best_boot_cl
in_boot <- (boot_labels == best_boot_cl)
# Jaccard: intersection / union
intersection <- sum(in_orig & in_boot)
union <- sum(in_orig | in_boot)
jaccard_vec[ci] <- if (union == 0) NA_real_ else intersection / union
}
return(jaccard_vec)
}
cat("✓ Fungsi stability (Hennig 2007 corrected) didefinisikan\n")## ✓ Fungsi stability (Hennig 2007 corrected) didefinisikan##
## Menjalankan test fungsi...# Gunakan PCA embeddings dari objek Seurat yang sudah ada
pca_embed_test <- Embeddings(gbm_10x, reduction = "pca")[, 1:20]
orig_labs_test <- as.integer(as.character(Idents(gbm_10x)))
unique_cl_test <- sort(unique(orig_labs_test))
tryCatch({
test_result <- compute_bootstrap_stability_once(
pca_embeddings = pca_embed_test,
original_labels = orig_labs_test,
unique_clusters = unique_cl_test,
k = 20,
resolution = 0.1
)
cat("✓ Test berhasil!\n")
cat(" - Jaccard per cluster:", round(test_result, 3), "\n")
cat(" - Range:", round(range(test_result, na.rm=TRUE), 3), "\n")
cat("\n✓ Siap untuk bootstrap penuh!\n")
}, error = function(e) {
cat("ERROR:", e$message, "\n")
stop("Function test failed.")
})## ✓ Test berhasil!
## - Jaccard per cluster: 0.986 0.996 0.424 0.985 0.349 0.828 0.176 0.698 0.992 0.301 0.145 0.029 0.974 0.017
## - Range: 0.017 0.996
##
## ✓ Siap untuk bootstrap penuh!# ========================================
# 9. Prepare Data for Bootstrap
# Gunakan PCA embeddings yang sudah ada di objek Seurat
# (lebih efisien dan konsisten dengan clustering asli)
# ========================================
cat("\n========================================\n")##
## ========================================## PREPARING DATA FOR BOOTSTRAP## ========================================# Gunakan PCA embeddings dari Seurat (sudah ada dari RunPCA sebelumnya)
# Ini jauh lebih efisien daripada re-komputasi scaled data
n_pcs_bootstrap <- 20 # jumlah PC yang digunakan untuk bootstrap
if (!"pca" %in% names(gbm_10x@reductions)) {
stop("PCA reduction tidak ditemukan! Pastikan RunPCA sudah dijalankan.")
}
# Ambil PCA embeddings (cells x PCs)
pca_embeddings_full <- Embeddings(gbm_10x, reduction = "pca")[, 1:n_pcs_bootstrap]
cat("PCA embeddings untuk bootstrap:\n")## PCA embeddings untuk bootstrap:## - Dimensi: 15072 20## - PCs digunakan: 1 - 20cat(" - Range nilai: [", round(min(pca_embeddings_full), 3), ",",
round(max(pca_embeddings_full), 3), "]\n")## - Range nilai: [ -52.935 , 44.267 ]# Get original clustering
original_clusters <- as.integer(as.character(Idents(gbm_10x)))
n_cells <- length(original_clusters)
unique_clusters <- sort(unique(original_clusters))
cat("\nOriginal clustering info:\n")##
## Original clustering info:## - Total cells: 15072## - Number of clusters: 14## - Cluster IDs: 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13## - Cluster sizes:## original_clusters
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13
## 3997 2176 1881 1590 1548 1228 749 641 420 252 207 148 125 110# Verifikasi dimensi
if (nrow(pca_embeddings_full) != n_cells) {
stop("Dimension mismatch! PCA rows != number of cells")
}
cat("\n✓ Data preparation complete\n")##
## ✓ Data preparation complete# ========================================
# RUN Bootstrap Analysis (100 iterasi)
# Implementasi Hennig (2007) yang benar
# ========================================
cat("\n========================================\n")##
## ========================================## BOOTSTRAP STABILITY ANALYSIS## ========================================# Bootstrap parameters
n_bootstrap <- 100
optimal_resolution <- 0.1
k_param <- 20 # lebih kecil dari sebelumnya (50) agar stabil di subsample
n_pcs_boot <- n_pcs_bootstrap # dari chunk sebelumnya
set.seed(123)
cat("Bootstrap parameters:\n")## Bootstrap parameters:## - Iterasi: 100## - Resolution: 0.1## - k.param: 20## - PCs: 20## - Jumlah cluster: 14## Metode: Hennig (2007) - unique-cell Jaccard## Estimasi waktu: 5-15 menit## Mulai bootstrap...# Initialize storage
bootstrap_results <- vector("list", n_bootstrap)
# Progress tracking
pb <- txtProgressBar(min = 0, max = n_bootstrap, style = 3)## | | | 0%start_time <- Sys.time()
error_count <- 0
# Run bootstrap
for (b in 1:n_bootstrap) {
result <- tryCatch({
compute_bootstrap_stability_once(
pca_embeddings = pca_embeddings_full,
original_labels = original_clusters,
unique_clusters = unique_clusters,
k = k_param,
resolution = optimal_resolution
)
}, error = function(e) {
error_count <<- error_count + 1
if (error_count <= 5) {
cat(sprintf("\nError iterasi %d: %s\n", b, e$message))
}
rep(NA_real_, length(unique_clusters))
})
bootstrap_results[[b]] <- result
# GC setiap 10 iterasi
if (b %% 10 == 0) gc(verbose = FALSE)
# Update progress
setTxtProgressBar(pb, b)
# Estimasi waktu setiap 20 iterasi
if (b %% 20 == 0) {
elapsed <- as.numeric(difftime(Sys.time(), start_time, units = "mins"))
est_remaining <- (elapsed / b) * (n_bootstrap - b)
cat(sprintf("\n %d/%d | Elapsed: %.1f min | Remaining: ~%.1f min | Errors: %d\n",
b, n_bootstrap, elapsed, est_remaining, error_count))
}
}## | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20%
## 20/100 | Elapsed: 0.6 min | Remaining: ~2.2 min | Errors: 0
## | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40%
## 40/100 | Elapsed: 1.0 min | Remaining: ~1.6 min | Errors: 0
## | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60%
## 60/100 | Elapsed: 1.6 min | Remaining: ~1.1 min | Errors: 0
## | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80%
## 80/100 | Elapsed: 2.1 min | Remaining: ~0.5 min | Errors: 0
## | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
## 100/100 | Elapsed: 2.9 min | Remaining: ~0.0 min | Errors: 0## used (Mb) gc trigger (Mb) max used (Mb)
## Ncells 8988803 480.1 14541112 776.6 14541112 776.6
## Vcells 979420428 7472.4 2012804046 15356.5 2012800200 15356.5end_time <- Sys.time()
total_time <- difftime(end_time, start_time, units = "mins")
cat("\n\n✓ Bootstrap selesai!\n")##
##
## ✓ Bootstrap selesai!cat(sprintf("Waktu total: %.1f menit | Errors: %d/%d\n",
as.numeric(total_time), error_count, n_bootstrap))## Waktu total: 2.9 menit | Errors: 0/100# ========================================
# Process Results
# ========================================
cat("\n========================================\n")##
## ========================================## PROCESSING RESULTS## ========================================# Filter valid results
valid_results <- bootstrap_results[sapply(bootstrap_results, function(x) {
!is.null(x) && !all(is.na(x))
})]
cat("Iterasi valid:", length(valid_results), "/", n_bootstrap, "\n\n")## Iterasi valid: 100 / 100if (length(valid_results) == 0) {
stop("ERROR: Tidak ada hasil valid!")
}
# Convert ke matrix
stability_matrix <- do.call(rbind, valid_results)
colnames(stability_matrix) <- paste0("Cluster_", unique_clusters)
stability_matrix <- stability_matrix[complete.cases(stability_matrix), ]
cat("Baris complete (tanpa NA):", nrow(stability_matrix), "\n")## Baris complete (tanpa NA): 100##
## Distribusi Jaccard per cluster:## Cluster_0 Cluster_1 Cluster_2 Cluster_3 Cluster_4 Cluster_5 Cluster_6
## Min. 0.931 0.857 0.248 0.432 0.339 0.412 0.158
## 1st Qu. 0.959 0.994 0.320 0.573 0.354 0.832 0.172
## Median 0.983 0.996 0.418 0.976 0.444 0.965 0.203
## Mean 0.974 0.989 0.383 0.822 0.476 0.891 0.292
## 3rd Qu. 0.987 0.997 0.429 0.982 0.613 0.973 0.372
## Max. 0.991 0.999 0.518 0.988 0.646 0.988 0.978
## Cluster_7 Cluster_8 Cluster_9 Cluster_10 Cluster_11 Cluster_12
## Min. 0.691 0.638 0.117 0.125 0.025 0.024
## 1st Qu. 0.712 0.989 0.267 0.150 0.030 0.030
## Median 0.719 0.989 0.278 0.978 0.038 0.974
## Mean 0.728 0.983 0.276 0.699 0.042 0.602
## 3rd Qu. 0.729 0.993 0.285 0.992 0.052 0.987
## Max. 0.995 1.000 0.357 1.000 0.137 1.000
## Cluster_13
## Min. 0.014
## 1st Qu. 0.020
## Median 0.027
## Mean 0.034
## 3rd Qu. 0.035
## Max. 0.603# Summary statistics
stability_summary <- data.frame(
cluster = unique_clusters,
cell_state = unname(cluster_annotations[as.character(unique_clusters)]),
n_cells = as.numeric(table(original_clusters)),
mean_stability = colMeans(stability_matrix, na.rm = TRUE),
sd_stability = apply(stability_matrix, 2, sd, na.rm = TRUE),
median_stability = apply(stability_matrix, 2, median, na.rm = TRUE),
q25_stability = apply(stability_matrix, 2, quantile, probs = 0.25, na.rm = TRUE),
q75_stability = apply(stability_matrix, 2, quantile, probs = 0.75, na.rm = TRUE),
ci_lower = apply(stability_matrix, 2, quantile, probs = 0.025, na.rm = TRUE),
ci_upper = apply(stability_matrix, 2, quantile, probs = 0.975, na.rm = TRUE),
dissolution_freq = apply(stability_matrix, 2, function(x) mean(x < 0.5, na.rm = TRUE)),
stringsAsFactors = FALSE
)
stability_summary$stability_category <- cut(
stability_summary$mean_stability,
breaks = c(-Inf, 0.6, 0.75, 0.85, Inf),
labels = c("Unstable", "Moderately Stable", "Stable", "Highly Stable")
)
stability_summary <- stability_summary %>% arrange(desc(mean_stability))
cat("\nSTABILITY SUMMARY:\n")##
## STABILITY SUMMARY:print(stability_summary %>% dplyr::select(cluster, cell_state, n_cells, mean_stability, stability_category))## cluster cell_state n_cells mean_stability
## Cluster_1 1 MES_AC_Transitional 2176 0.98911648
## Cluster_8 8 Oligodendrocyte 420 0.98349564
## Cluster_0 0 Microglia_Macrophage 3997 0.97356710
## Cluster_5 5 Macrophage_complement 1228 0.89096841
## Cluster_3 3 MES-like 1590 0.82165614
## Cluster_7 7 Cycling_Tumor 641 0.72776413
## Cluster_10 10 T_cell 207 0.69932756
## Cluster_12 12 Macrophage_LAM 125 0.60233116
## Cluster_4 4 NPC-like 1548 0.47614106
## Cluster_2 2 NPC-like_mature 1881 0.38264315
## Cluster_6 6 AC-like 749 0.29181945
## Cluster_9 9 Endothelial_Astrocyte 252 0.27643262
## Cluster_11 11 Erythrocyte 148 0.04222826
## Cluster_13 13 Pericyte_Stromal 110 0.03430782
## stability_category
## Cluster_1 Highly Stable
## Cluster_8 Highly Stable
## Cluster_0 Highly Stable
## Cluster_5 Highly Stable
## Cluster_3 Stable
## Cluster_7 Moderately Stable
## Cluster_10 Moderately Stable
## Cluster_12 Moderately Stable
## Cluster_4 Unstable
## Cluster_2 Unstable
## Cluster_6 Unstable
## Cluster_9 Unstable
## Cluster_11 Unstable
## Cluster_13 Unstable# Save
write.csv(stability_summary, "cluster_stability_summary.csv", row.names = FALSE)
write.csv(stability_matrix, "bootstrap_stability_matrix.csv", row.names = FALSE)
cat("\n✓ File tersimpan:\n")##
## ✓ File tersimpan:## - cluster_stability_summary.csv## - bootstrap_stability_matrix.csv##
## Stability distribution:##
## Unstable Moderately Stable Stable Highly Stable
## 6 3 1 4##
## ========================================## ANALYSIS COMPLETE!## ========================================# ========================================
# DEBUG: Verifikasi hasil bootstrap
# ========================================
cat("\n========================================\n")##
## ========================================## VERIFIKASI BOOTSTRAP RESULTS## ========================================# Ringkasan status hasil
n_null <- sum(sapply(bootstrap_results, is.null))
n_all_na <- sum(sapply(bootstrap_results, function(x) !is.null(x) && all(is.na(x))))
n_valid <- sum(sapply(bootstrap_results, function(x) !is.null(x) && !all(is.na(x))))
cat("Status hasil bootstrap:\n")## Status hasil bootstrap:## - NULL results: 0## - All-NA results: 0## - Valid results: 100## Cek 5 iterasi pertama:for (i in 1:min(5, length(bootstrap_results))) {
cat(sprintf(" Iterasi %d: %s\n", i,
paste(round(bootstrap_results[[i]], 3), collapse = ", ")))
}## Iterasi 1: 0.963, 0.996, 0.31, 0.567, 0.525, 0.979, 0.703, 0.729, 0.992, 0.266, 0.993, 0.052, 0.029, 0.029
## Iterasi 2: 0.984, 0.995, 0.428, 0.982, 0.642, 0.846, 0.355, 0.724, 0.989, 0.273, 0.139, 0.055, 0.974, 0.043
## Iterasi 3: 0.984, 0.998, 0.411, 0.983, 0.602, 0.966, 0.403, 0.722, 0.993, 0.275, 0.984, 0.048, 0.975, 0.033
## Iterasi 4: 0.986, 0.996, 0.25, 0.449, 0.632, 0.963, 0.217, 0.723, 0.993, 0.273, 0.984, 0.039, 0.987, 0.022
## Iterasi 5: 0.988, 0.997, 0.259, 0.571, 0.416, 0.412, 0.203, 0.703, 0.993, 0.293, 1, 0.037, 1, 0.022##
## Distribusi Jaccard mean per cluster:if (exists("stability_summary")) {
print(stability_summary %>%
dplyr::select(cluster, cell_state, mean_stability, stability_category) %>%
arrange(desc(mean_stability)))
}## cluster cell_state mean_stability stability_category
## Cluster_1 1 MES_AC_Transitional 0.98911648 Highly Stable
## Cluster_8 8 Oligodendrocyte 0.98349564 Highly Stable
## Cluster_0 0 Microglia_Macrophage 0.97356710 Highly Stable
## Cluster_5 5 Macrophage_complement 0.89096841 Highly Stable
## Cluster_3 3 MES-like 0.82165614 Stable
## Cluster_7 7 Cycling_Tumor 0.72776413 Moderately Stable
## Cluster_10 10 T_cell 0.69932756 Moderately Stable
## Cluster_12 12 Macrophage_LAM 0.60233116 Moderately Stable
## Cluster_4 4 NPC-like 0.47614106 Unstable
## Cluster_2 2 NPC-like_mature 0.38264315 Unstable
## Cluster_6 6 AC-like 0.29181945 Unstable
## Cluster_9 9 Endothelial_Astrocyte 0.27643262 Unstable
## Cluster_11 11 Erythrocyte 0.04222826 Unstable
## Cluster_13 13 Pericyte_Stromal 0.03430782 Unstable# ========================================
# 12. Ringkasan Akhir Stability Results
# (stability_summary sudah dibuat di chunk bootstrap)
# ========================================
cat("\n========================================\n")##
## ========================================## RINGKASAN STABILITAS CLUSTER## ========================================## Valid bootstrap iterations: 100 dari 100## Stability Summary (terurut dari tertinggi):print(stability_summary %>%
dplyr::select(cluster, cell_state, n_cells, mean_stability, sd_stability,
median_stability, dissolution_freq, stability_category))## cluster cell_state n_cells mean_stability sd_stability
## Cluster_1 1 MES_AC_Transitional 2176 0.98911648 0.01849675
## Cluster_8 8 Oligodendrocyte 420 0.98349564 0.04969743
## Cluster_0 0 Microglia_Macrophage 3997 0.97356710 0.01486751
## Cluster_5 5 Macrophage_complement 1228 0.89096841 0.13406342
## Cluster_3 3 MES-like 1590 0.82165614 0.22174131
## Cluster_7 7 Cycling_Tumor 641 0.72776413 0.04862425
## Cluster_10 10 T_cell 207 0.69932756 0.40134647
## Cluster_12 12 Macrophage_LAM 125 0.60233116 0.47048976
## Cluster_4 4 NPC-like 1548 0.47614106 0.12397384
## Cluster_2 2 NPC-like_mature 1881 0.38264315 0.06899437
## Cluster_6 6 AC-like 749 0.29181945 0.20692969
## Cluster_9 9 Endothelial_Astrocyte 252 0.27643262 0.02312202
## Cluster_11 11 Erythrocyte 148 0.04222826 0.01677463
## Cluster_13 13 Pericyte_Stromal 110 0.03430782 0.05817144
## median_stability dissolution_freq stability_category
## Cluster_1 0.99596969 0.00 Highly Stable
## Cluster_8 0.98926651 0.00 Highly Stable
## Cluster_0 0.98327261 0.00 Highly Stable
## Cluster_5 0.96515244 0.06 Highly Stable
## Cluster_3 0.97630779 0.14 Stable
## Cluster_7 0.71869450 0.00 Moderately Stable
## Cluster_10 0.97793640 0.34 Moderately Stable
## Cluster_12 0.97402159 0.40 Moderately Stable
## Cluster_4 0.44362624 0.54 Unstable
## Cluster_2 0.41802858 0.98 Unstable
## Cluster_6 0.20269057 0.91 Unstable
## Cluster_9 0.27846750 1.00 Unstable
## Cluster_11 0.03777670 1.00 Unstable
## Cluster_13 0.02715031 0.99 Unstable##
## Distribusi kategori stabilitas:##
## Unstable Moderately Stable Stable Highly Stable
## 6 3 1 4##
## ========================================## STABILITY ANALYSIS COMPLETE!## ========================================## Next step: Run visualization scripts11. Visualize Stability Results
# Bar plot of mean stability
stability_barplot <- ggplot(
stability_summary,
aes(x = reorder(cell_state, mean_stability), y = mean_stability, fill = stability_category)
) +
geom_bar(stat = "identity", width = 0.7) +
geom_errorbar(
aes(ymin = ci_lower, ymax = ci_upper),
width = 0.2
) +
geom_hline(yintercept = c(0.6, 0.75, 0.85), linetype = "dashed", color = "gray40") +
scale_fill_manual(
values = c(
"Unstable" = "#d73027",
"Moderately Stable" = "#fee090",
"Stable" = "#91bfdb",
"Highly Stable" = "#4575b4"
)
) +
coord_flip() +
labs(
title = "Cluster Stability Assessment",
x = "Cell State",
y = "Mean Jaccard Coefficient",
fill = "Stability Category"
) +
theme_publication +
theme(legend.position = "bottom")
print(stability_barplot)
# Convert stability matrix to long format for visualization
stability_long <- stability_matrix %>%
as.data.frame() %>%
mutate(iteration = 1:n()) %>%
pivot_longer(
cols = -iteration,
names_to = "cluster_col",
values_to = "jaccard"
) %>%
mutate(cluster = as.integer(sub("Cluster_", "", cluster_col))) %>%
left_join(
stability_summary %>% dplyr::select(cluster, cell_state),
by = "cluster"
)
# Raincloud plot (violin + boxplot + points)
raincloud_plot <- ggplot(
stability_long,
aes(x = reorder(cell_state, jaccard, FUN = median), y = jaccard, fill = cell_state)
) +
geom_violin(alpha = 0.6, scale = "width") +
geom_boxplot(width = 0.2, alpha = 0.8, outlier.shape = NA) +
geom_hline(yintercept = c(0.6, 0.75, 0.85), linetype = "dashed", color = "gray40") +
coord_flip() +
scale_fill_brewer(palette = "Set3") +
labs(
title = "Distribution of Jaccard Coefficients Across Bootstrap Iterations",
x = "Cell State",
y = "Jaccard Coefficient"
) +
theme_publication +
theme(legend.position = "none")
print(raincloud_plot)
# Ridge plot for distributions
ridge_plot <- ggplot(
stability_long,
aes(x = jaccard, y = reorder(cell_state, jaccard, FUN = median), fill = cell_state)
) +
geom_density_ridges(alpha = 0.7, scale = 2) +
geom_vline(xintercept = c(0.6, 0.75, 0.85), linetype = "dashed", color = "gray40") +
scale_fill_brewer(palette = "Set3") +
labs(
title = "Stability Distribution per Cell State",
x = "Jaccard Coefficient",
y = "Cell State"
) +
theme_publication +
theme(legend.position = "none")
print(ridge_plot)
# Heatmap of bootstrap iterations
# Sample 100 iterations for visualization
sample_iterations <- sample(1:nrow(stability_matrix), min(100, nrow(stability_matrix)))
stability_heatmap_data <- stability_matrix[sample_iterations, ]
# Reorder columns by mean stability
col_order <- stability_summary %>%
arrange(desc(mean_stability)) %>%
pull(cluster)
stability_heatmap_data <- stability_heatmap_data[, paste0("Cluster_", col_order)]
colnames(stability_heatmap_data) <- cluster_annotations[as.character(col_order)]
# Create heatmap
pheatmap(
t(stability_heatmap_data),
color = viridis(100),
cluster_rows = FALSE,
cluster_cols = FALSE,
show_colnames = FALSE,
main = "Stability Across Bootstrap Iterations (100 samples)",
fontsize = 10,
border_color = NA,
filename = "stability_heatmap.png",
width = 12,
height = 6
)# Plot dissolution frequency
dissolution_plot <- ggplot(
stability_summary,
aes(x = reorder(cell_state, dissolution_freq), y = dissolution_freq * 100, fill = cell_state)
) +
geom_bar(stat = "identity", width = 0.7) +
coord_flip() +
scale_fill_brewer(palette = "Set3") +
labs(
title = "Cluster Dissolution Frequency",
x = "Cell State",
y = "Dissolution Frequency (%)",
caption = "Proportion of bootstrap iterations with Jaccard < 0.5"
) +
theme_publication +
theme(legend.position = "none")
print(dissolution_plot)
12. Biological Characterization of Stable vs. Unstable Clusters
# Calculate marker specificity for each cluster
# sig_markers dari presto menggunakan kolom: feature, group, logFC, pct_in, pct_out, auc, padj
marker_specificity <- sig_markers %>%
group_by(group) %>%
summarize(
mean_logFC = mean(logFC, na.rm = TRUE),
max_logFC = max(logFC, na.rm = TRUE),
n_markers = n(),
specificity_score = mean(pct_in / (pct_out + 1), na.rm = TRUE) # pct_in/pct_out dari presto
) %>%
mutate(cluster = as.integer(as.character(group))) %>%
left_join(
stability_summary %>% dplyr::select(cluster, cell_state, mean_stability),
by = "cluster"
)
print(marker_specificity)## # A tibble: 14 × 8
## group mean_logFC max_logFC n_markers specificity_score cluster cell_state
## <chr> <dbl> <dbl> <int> <dbl> <int> <chr>
## 1 0 1.01 3.06 279 4.24 0 Microglia_Mac…
## 2 1 0.740 2.90 644 3.44 1 MES_AC_Transi…
## 3 10 0.931 2.63 212 6.51 10 T_cell
## 4 11 1.74 6.42 16 8.03 11 Erythrocyte
## 5 12 0.980 3.65 340 4.23 12 Macrophage_LAM
## 6 13 0.812 2.37 121 3.90 13 Pericyte_Stro…
## 7 2 0.669 1.51 120 1.85 2 NPC-like_matu…
## 8 3 1.23 3.57 41 2.47 3 MES-like
## 9 4 0.742 2.13 210 1.98 4 NPC-like
## 10 5 1.10 2.72 80 2.18 5 Macrophage_co…
## 11 6 0.786 2.06 179 2.13 6 AC-like
## 12 7 0.878 3.27 350 2.95 7 Cycling_Tumor
## 13 8 0.992 4.91 349 7.29 8 Oligodendrocy…
## 14 9 0.901 3.37 447 3.86 9 Endothelial_A…
## # ℹ 1 more variable: mean_stability <dbl># Correlation between specificity and stability
specificity_vs_stability <- ggplot(
marker_specificity,
aes(x = specificity_score, y = mean_stability, label = cell_state)
) +
geom_point(size = 4, aes(color = cell_state)) +
geom_smooth(method = "lm", se = TRUE, color = "gray50", linetype = "dashed") +
geom_text(vjust = -1, hjust = 0.5, size = 3) +
scale_color_brewer(palette = "Set3") +
labs(
title = "Marker Specificity vs. Cluster Stability",
x = "Marker Specificity Score",
y = "Mean Jaccard Coefficient",
subtitle = sprintf(
"Spearman correlation: %.2f",
cor(marker_specificity$specificity_score, marker_specificity$mean_stability, method = "spearman")
)
) +
theme_publication +
theme(legend.position = "none")
print(specificity_vs_stability)
# Calculate within-cluster transcriptomic homogeneity
# This is computationally intensive - sample cells if dataset is large
calculate_cluster_homogeneity <- function(seurat_obj, cluster_id, max_cells = 500) {
cells_in_cluster <- WhichCells(seurat_obj, idents = cluster_id)
# Sample if too many cells
if (length(cells_in_cluster) > max_cells) {
cells_in_cluster <- sample(cells_in_cluster, max_cells)
}
# Get expression matrix - compatible Seurat v4 & v5
expr_matrix <- tryCatch({
# Seurat v5
sd_layer <- LayerData(seurat_obj, assay = "RNA", layer = "scale.data")
t(as.matrix(sd_layer[VariableFeatures(seurat_obj), cells_in_cluster]))
}, error = function(e) {
# Seurat v4 fallback
t(as.matrix(seurat_obj@assays$RNA@scale.data[VariableFeatures(seurat_obj), cells_in_cluster]))
})
# Calculate pairwise correlations
cor_matrix <- cor(t(expr_matrix))
# Return mean correlation (excluding diagonal)
mean_cor <- mean(cor_matrix[upper.tri(cor_matrix)])
return(mean_cor)
}
cat("Calculating within-cluster homogeneity...\n")## Calculating within-cluster homogeneity...# Alias: gbm_seurat = gbm_10x (konsistensi penamaan di seluruh script)
gbm_seurat <- gbm_10x
homogeneity_results <- data.frame(
cluster = unique_clusters,
mean_correlation = sapply(unique_clusters, function(cid) {
calculate_cluster_homogeneity(gbm_seurat, cid)
})
) %>%
left_join(
stability_summary %>% dplyr::select(cluster, cell_state, mean_stability),
by = "cluster"
)
print(homogeneity_results)## cluster mean_correlation cell_state mean_stability
## 1 0 0.08931513 Microglia_Macrophage 0.97356710
## 2 1 0.17543694 MES_AC_Transitional 0.98911648
## 3 2 0.04180014 NPC-like_mature 0.38264315
## 4 3 0.06578758 MES-like 0.82165614
## 5 4 0.06624653 NPC-like 0.47614106
## 6 5 0.08914908 Macrophage_complement 0.89096841
## 7 6 0.10594823 AC-like 0.29181945
## 8 7 0.20116611 Cycling_Tumor 0.72776413
## 9 8 0.33225613 Oligodendrocyte 0.98349564
## 10 9 0.27601443 Endothelial_Astrocyte 0.27643262
## 11 10 0.29145627 T_cell 0.69932756
## 12 11 0.21913672 Erythrocyte 0.04222826
## 13 12 0.49147993 Macrophage_LAM 0.60233116
## 14 13 0.12488089 Pericyte_Stromal 0.03430782# Plot homogeneity vs stability
homogeneity_vs_stability <- ggplot(
homogeneity_results,
aes(x = mean_correlation, y = mean_stability, label = cell_state)
) +
geom_point(size = 4, aes(color = cell_state)) +
geom_smooth(method = "lm", se = TRUE, color = "gray50", linetype = "dashed") +
geom_text(vjust = -1, hjust = 0.5, size = 3) +
scale_color_brewer(palette = "Set3") +
labs(
title = "Within-cluster Homogeneity vs. Cluster Stability",
x = "Mean Pairwise Correlation",
y = "Mean Jaccard Coefficient",
subtitle = sprintf(
"Spearman correlation: %.2f",
cor(homogeneity_results$mean_correlation, homogeneity_results$mean_stability, method = "spearman")
)
) +
theme_publication +
theme(legend.position = "none")
print(homogeneity_vs_stability)
# Calculate silhouette scores
n_pcs <- n_pcs_bootstrap # gunakan n_pcs yang sama dengan bootstrap
dist_matrix <- dist(Embeddings(gbm_seurat, "pca")[, 1:n_pcs])
sil_scores <- silhouette(as.numeric(Idents(gbm_seurat)), dist_matrix)
# Summarize by cluster
silhouette_summary <- data.frame(
cluster = unique_clusters,
mean_silhouette = sapply(unique_clusters, function(cid) {
mean(sil_scores[Idents(gbm_seurat) == cid, "sil_width"])
})
) %>%
left_join(
stability_summary %>% dplyr::select(cluster, cell_state, mean_stability),
by = "cluster"
)
print(silhouette_summary)## cluster mean_silhouette cell_state mean_stability
## 1 0 0.08191495 Microglia_Macrophage 0.97356710
## 2 1 0.43894833 MES_AC_Transitional 0.98911648
## 3 2 -0.06384047 NPC-like_mature 0.38264315
## 4 3 0.25668648 MES-like 0.82165614
## 5 4 0.22694009 NPC-like 0.47614106
## 6 5 0.29924473 Macrophage_complement 0.89096841
## 7 6 0.29959318 AC-like 0.29181945
## 8 7 0.44016356 Cycling_Tumor 0.72776413
## 9 8 0.46316610 Oligodendrocyte 0.98349564
## 10 9 0.59015106 Endothelial_Astrocyte 0.27643262
## 11 10 0.53679959 T_cell 0.69932756
## 12 11 0.47475159 Erythrocyte 0.04222826
## 13 12 0.71679661 Macrophage_LAM 0.60233116
## 14 13 0.19159496 Pericyte_Stromal 0.03430782# Plot silhouette vs stability
silhouette_vs_stability <- ggplot(
silhouette_summary,
aes(x = mean_silhouette, y = mean_stability, label = cell_state)
) +
geom_point(size = 4, aes(color = cell_state)) +
geom_smooth(method = "lm", se = TRUE, color = "gray50", linetype = "dashed") +
geom_text(vjust = -1, hjust = 0.5, size = 3) +
scale_color_brewer(palette = "Set3") +
labs(
title = "Silhouette Score vs. Cluster Stability",
x = "Mean Silhouette Score",
y = "Mean Jaccard Coefficient",
subtitle = sprintf(
"Spearman correlation: %.2f",
cor(silhouette_summary$mean_silhouette, silhouette_summary$mean_stability, method = "spearman")
)
) +
theme_publication +
theme(legend.position = "none")
print(silhouette_vs_stability)
# Combine all metrics
combined_metrics <- stability_summary %>%
left_join(marker_specificity %>% dplyr::select(cluster, specificity_score), by = "cluster") %>%
left_join(homogeneity_results %>% dplyr::select(cluster, mean_correlation), by = "cluster") %>%
left_join(silhouette_summary %>% dplyr::select(cluster, mean_silhouette), by = "cluster")
print(combined_metrics)## cluster cell_state n_cells mean_stability sd_stability
## 1 1 MES_AC_Transitional 2176 0.98911648 0.01849675
## 2 8 Oligodendrocyte 420 0.98349564 0.04969743
## 3 0 Microglia_Macrophage 3997 0.97356710 0.01486751
## 4 5 Macrophage_complement 1228 0.89096841 0.13406342
## 5 3 MES-like 1590 0.82165614 0.22174131
## 6 7 Cycling_Tumor 641 0.72776413 0.04862425
## 7 10 T_cell 207 0.69932756 0.40134647
## 8 12 Macrophage_LAM 125 0.60233116 0.47048976
## 9 4 NPC-like 1548 0.47614106 0.12397384
## 10 2 NPC-like_mature 1881 0.38264315 0.06899437
## 11 6 AC-like 749 0.29181945 0.20692969
## 12 9 Endothelial_Astrocyte 252 0.27643262 0.02312202
## 13 11 Erythrocyte 148 0.04222826 0.01677463
## 14 13 Pericyte_Stromal 110 0.03430782 0.05817144
## median_stability q25_stability q75_stability ci_lower ci_upper
## 1 0.99596969 0.99427746 0.99707227 0.95704355 0.99822370
## 2 0.98926651 0.98858227 0.99259943 0.98461220 0.99633060
## 3 0.98327261 0.95943967 0.98650759 0.94617455 0.98937837
## 4 0.96515244 0.83227163 0.97322217 0.43301558 0.98361636
## 5 0.97630779 0.57289390 0.98203833 0.43800378 0.98746705
## 6 0.71869450 0.71198569 0.72916492 0.69715083 0.87627373
## 7 0.97793640 0.14965842 0.99212598 0.13065673 1.00000000
## 8 0.97402159 0.02997149 0.98722006 0.02605096 1.00000000
## 9 0.44362624 0.35392493 0.61335229 0.34350975 0.64112817
## 10 0.41802858 0.31979761 0.42940061 0.25309842 0.44447255
## 11 0.20269057 0.17199583 0.37171713 0.16312184 0.97403419
## 12 0.27846750 0.26678378 0.28471518 0.25056154 0.30383738
## 13 0.03777670 0.03006813 0.05199217 0.02698853 0.07690434
## 14 0.02715031 0.02045780 0.03544552 0.01660308 0.05361316
## dissolution_freq stability_category specificity_score mean_correlation
## 1 0.00 Highly Stable 3.440992 0.17543694
## 2 0.00 Highly Stable 7.289669 0.33225613
## 3 0.00 Highly Stable 4.237465 0.08931513
## 4 0.06 Highly Stable 2.176000 0.08914908
## 5 0.14 Stable 2.467094 0.06578758
## 6 0.00 Moderately Stable 2.951864 0.20116611
## 7 0.34 Moderately Stable 6.505862 0.29145627
## 8 0.40 Moderately Stable 4.227366 0.49147993
## 9 0.54 Unstable 1.977918 0.06624653
## 10 0.98 Unstable 1.847073 0.04180014
## 11 0.91 Unstable 2.132470 0.10594823
## 12 1.00 Unstable 3.858486 0.27601443
## 13 1.00 Unstable 8.033775 0.21913672
## 14 0.99 Unstable 3.900804 0.12488089
## mean_silhouette
## 1 0.43894833
## 2 0.46316610
## 3 0.08191495
## 4 0.29924473
## 5 0.25668648
## 6 0.44016356
## 7 0.53679959
## 8 0.71679661
## 9 0.22694009
## 10 -0.06384047
## 11 0.29959318
## 12 0.59015106
## 13 0.47475159
## 14 0.19159496write.csv(combined_metrics, "combined_cluster_metrics.csv", row.names = FALSE)
# Correlation matrix of metrics
metric_columns <- c("mean_stability", "specificity_score", "mean_correlation", "mean_silhouette")
cor_matrix <- cor(combined_metrics[, metric_columns], use = "complete.obs", method = "spearman")
# Visualize correlation matrix
pheatmap(
cor_matrix,
display_numbers = TRUE,
number_format = "%.2f",
color = colorRampPalette(c("blue", "white", "red"))(100),
breaks = seq(-1, 1, length.out = 101),
main = "Correlation Between Cluster Quality Metrics",
fontsize = 12,
filename = "metric_correlation_heatmap.png",
width = 8,
height = 6
)13. Pathway Enrichment Analysis
# Perform GO enrichment for each cluster
# This requires converting gene symbols to Entrez IDs
perform_enrichment <- function(cluster_markers, cluster_id) {
# Get significant markers for this cluster
cluster_genes <- cluster_markers %>%
filter(group == as.character(cluster_id), padj < 0.05) %>%
pull(feature)
# Convert to Entrez IDs
entrez_ids <- bitr(
cluster_genes,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db
)
if (nrow(entrez_ids) < 5) {
return(NULL)
}
# GO enrichment
ego <- enrichGO(
gene = entrez_ids$ENTREZID,
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05,
qvalueCutoff = 0.2,
readable = TRUE
)
return(ego)
}
# Run enrichment for all clusters
cat("Performing pathway enrichment analysis...\n")## Performing pathway enrichment analysis...enrichment_results <- list()
for (cid in unique_clusters) {
cat(paste("Processing cluster", cid, "...\n"))
enrichment_results[[as.character(cid)]] <- perform_enrichment(sig_markers, cid)
}## Processing cluster 0 ...## Processing cluster 1 ...## Processing cluster 2 ...## Processing cluster 3 ...## Processing cluster 4 ...## Processing cluster 5 ...## Processing cluster 6 ...## Processing cluster 7 ...## Processing cluster 8 ...## Processing cluster 9 ...## Processing cluster 10 ...## Processing cluster 11 ...## Processing cluster 12 ...## Processing cluster 13 ...# Filter out NULL results
enrichment_results <- enrichment_results[!sapply(enrichment_results, is.null)]# Visualize enrichment for selected clusters
# Focus on highly stable and unstable clusters
if (length(enrichment_results) > 0) {
# Plot for highly stable cluster (e.g., MES-like)
stable_cluster <- stability_summary %>%
filter(stability_category == "Highly Stable") %>%
slice_head(n = 1) %>%
pull(cluster)
if (as.character(stable_cluster) %in% names(enrichment_results)) {
ego_stable <- enrichment_results[[as.character(stable_cluster)]]
if (!is.null(ego_stable) && nrow(ego_stable) > 0) {
p1 <- dotplot(ego_stable, showCategory = 15) +
ggtitle(paste("GO Enrichment:", cluster_annotations[as.character(stable_cluster)])) +
theme_publication
print(p1)
ggsave(
paste0("enrichment_stable_cluster_", stable_cluster, ".png"),
p1,
width = 12,
height = 8,
dpi = 300
)
}
}
# Plot for unstable cluster (e.g., NPC-like)
unstable_cluster <- stability_summary %>%
filter(stability_category == "Unstable") %>%
slice_head(n = 1) %>%
pull(cluster)
if (as.character(unstable_cluster) %in% names(enrichment_results)) {
ego_unstable <- enrichment_results[[as.character(unstable_cluster)]]
if (!is.null(ego_unstable) && nrow(ego_unstable) > 0) {
p2 <- dotplot(ego_unstable, showCategory = 15) +
ggtitle(paste("GO Enrichment:", cluster_annotations[as.character(unstable_cluster)])) +
theme_publication
print(p2)
ggsave(
paste0("enrichment_unstable_cluster_", unstable_cluster, ".png"),
p2,
width = 12,
height = 8,
dpi = 300
)
}
}
}
# GSEA using Hallmark gene sets
# Get Hallmark gene sets
hallmark_sets <- msigdbr(species = "Homo sapiens", category = "H")
perform_gsea <- function(cluster_markers, cluster_id, gene_sets) {
# Get all genes and their log2FC for this cluster
cluster_genes <- cluster_markers %>%
filter(group == as.character(cluster_id)) %>%
arrange(desc(logFC))
# Create ranked gene list
gene_list <- cluster_genes$logFC
names(gene_list) <- cluster_genes$feature
# Format gene sets
gene_sets_list <- split(gene_sets$gene_symbol, gene_sets$gs_name)
# Run GSEA
gsea_result <- GSEA(
geneList = gene_list,
TERM2GENE = gene_sets[, c("gs_name", "gene_symbol")],
pvalueCutoff = 0.05,
pAdjustMethod = "BH"
)
return(gsea_result)
}
# Run GSEA for selected clusters
cat("Performing GSEA with Hallmark gene sets...\n")## Performing GSEA with Hallmark gene sets...gsea_results <- list()
for (cid in unique_clusters) {
cat(paste("Processing cluster", cid, "...\n"))
tryCatch({
gsea_results[[as.character(cid)]] <- perform_gsea(
all_markers,
cid,
hallmark_sets
)
}, error = function(e) {
cat(paste("Error in cluster", cid, ":", e$message, "\n"))
})
}## Processing cluster 0 ...## Processing cluster 1 ...## Processing cluster 2 ...## Processing cluster 3 ...## Processing cluster 4 ...## Processing cluster 5 ...## Processing cluster 6 ...## Processing cluster 7 ...## Processing cluster 8 ...## Processing cluster 9 ...## Processing cluster 10 ...## Processing cluster 11 ...## Processing cluster 12 ...## Processing cluster 13 ...# Visualize GSEA results
if (length(gsea_results) > 0) {
# For highly stable cluster
if (as.character(stable_cluster) %in% names(gsea_results)) {
gsea_stable <- gsea_results[[as.character(stable_cluster)]]
if (!is.null(gsea_stable) && nrow(gsea_stable) > 0) {
p3 <- dotplot(gsea_stable, showCategory = 15, split = ".sign") +
facet_grid(~.sign) +
ggtitle(paste("GSEA Hallmark:", cluster_annotations[as.character(stable_cluster)])) +
theme_publication
print(p3)
ggsave(
paste0("gsea_stable_cluster_", stable_cluster, ".png"),
p3,
width = 14,
height = 8,
dpi = 300
)
}
}
}
14. Alternative Resampling Methods Comparison
# Subsampling (without replacement) untuk perbandingan
# Menggunakan pendekatan Jaccard yang benar: hanya pada sel yang di-subsample
n_subsample <- 100 # dikurangi agar tidak terlalu lama
subsample_prop <- 0.8
cat(paste("Running", n_subsample, "subsampling iterations...\n"))## Running 100 subsampling iterations...## | | | 0%for (b in 1:n_subsample) {
tryCatch({
# Subsample tanpa replacement
subsample_size <- floor(n_cells * subsample_prop)
subsample_indices <- sample(n_cells, subsample_size, replace = FALSE)
subsample_embed <- pca_embeddings_full[subsample_indices, , drop = FALSE]
orig_labs_sub <- original_clusters[subsample_indices]
k_sub <- min(k_param, nrow(subsample_embed) - 1)
# KNN + Louvain
knn_sub <- FNN::get.knn(subsample_embed, k = k_sub)
edges_sub <- data.frame(
from = rep(1:nrow(subsample_embed), each = k_sub),
to = as.vector(t(knn_sub$nn.index))
)
edges_sub <- edges_sub[edges_sub$from != edges_sub$to, ]
g_sub <- igraph::simplify(
igraph::graph_from_data_frame(edges_sub, directed = FALSE,
vertices = data.frame(name = 1:nrow(subsample_embed)))
)
sub_clusters <- as.integer(igraph::membership(
igraph::cluster_louvain(g_sub, resolution = optimal_resolution)
))
# Jaccard: label matching
jac_vec <- sapply(unique_clusters, function(cl) {
in_orig <- (orig_labs_sub == cl)
if (sum(in_orig) == 0) return(NA_real_)
best_boot_cl <- as.integer(names(which.max(table(sub_clusters[in_orig]))))
in_boot <- (sub_clusters == best_boot_cl)
intersection <- sum(in_orig & in_boot)
union <- sum(in_orig | in_boot)
if (union == 0) NA_real_ else intersection / union
})
subsample_results[[b]] <- jac_vec
}, error = function(e) {
subsample_results[[b]] <- rep(NA_real_, length(unique_clusters))
})
setTxtProgressBar(pb, b)
}## | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%# Process subsample results
subsample_matrix <- do.call(rbind, subsample_results)
colnames(subsample_matrix) <- paste0("Cluster_", unique_clusters)
subsample_matrix <- subsample_matrix[complete.cases(subsample_matrix), ]
subsample_summary <- data.frame(
cluster = unique_clusters,
cell_state = unname(cluster_annotations[as.character(unique_clusters)]),
mean_stability_subsample = colMeans(subsample_matrix, na.rm = TRUE)
)
print(subsample_summary)## cluster cell_state mean_stability_subsample
## Cluster_0 0 Microglia_Macrophage 0.97032259
## Cluster_1 1 MES_AC_Transitional 0.99386861
## Cluster_2 2 NPC-like_mature 0.41803241
## Cluster_3 3 MES-like 0.88355256
## Cluster_4 4 NPC-like 0.48336639
## Cluster_5 5 Macrophage_complement 0.87221898
## Cluster_6 6 AC-like 0.41858927
## Cluster_7 7 Cycling_Tumor 0.71971277
## Cluster_8 8 Oligodendrocyte 0.99060631
## Cluster_9 9 Endothelial_Astrocyte 0.27551556
## Cluster_10 10 T_cell 0.64882014
## Cluster_11 11 Erythrocyte 0.04113362
## Cluster_12 12 Macrophage_LAM 0.68801379
## Cluster_13 13 Pericyte_Stromal 0.10236947# Noise injection untuk perbandingan
# Tambahkan noise Gaussian ke PCA embeddings, lalu re-cluster
n_noise <- 100
noise_sd <- 0.5 # disesuaikan untuk PCA space (bukan gene expression space)
cat(paste("Running", n_noise, "noise injection iterations...\n"))## Running 100 noise injection iterations...## | | | 0%for (b in 1:n_noise) {
tryCatch({
# Tambahkan Gaussian noise ke PCA embeddings
noise_embed <- pca_embeddings_full +
matrix(rnorm(n_cells * ncol(pca_embeddings_full), 0, noise_sd),
nrow = n_cells)
k_actual <- min(k_param, n_cells - 1)
# KNN + Louvain
knn_noise <- FNN::get.knn(noise_embed, k = k_actual)
edges_noise <- data.frame(
from = rep(1:n_cells, each = k_actual),
to = as.vector(t(knn_noise$nn.index))
)
edges_noise <- edges_noise[edges_noise$from != edges_noise$to, ]
g_noise <- igraph::simplify(
igraph::graph_from_data_frame(edges_noise, directed = FALSE,
vertices = data.frame(name = 1:n_cells))
)
noise_clusters <- as.integer(igraph::membership(
igraph::cluster_louvain(g_noise, resolution = optimal_resolution)
))
# Jaccard dengan label matching
jac_vec <- sapply(unique_clusters, function(cl) {
in_orig <- (original_clusters == cl)
if (sum(in_orig) == 0) return(NA_real_)
best_boot_cl <- as.integer(names(which.max(table(noise_clusters[in_orig]))))
in_noise_cl <- (noise_clusters == best_boot_cl)
intersection <- sum(in_orig & in_noise_cl)
union <- sum(in_orig | in_noise_cl)
if (union == 0) NA_real_ else intersection / union
})
noise_results[[b]] <- jac_vec
}, error = function(e) {
noise_results[[b]] <- rep(NA_real_, length(unique_clusters))
})
setTxtProgressBar(pb, b)
}## | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%# Process noise results
noise_matrix <- do.call(rbind, noise_results)
colnames(noise_matrix) <- paste0("Cluster_", unique_clusters)
noise_matrix <- noise_matrix[complete.cases(noise_matrix), ]
noise_summary <- data.frame(
cluster = unique_clusters,
cell_state = unname(cluster_annotations[as.character(unique_clusters)]),
mean_stability_noise = colMeans(noise_matrix, na.rm = TRUE)
)
print(noise_summary)## cluster cell_state mean_stability_noise
## Cluster_0 0 Microglia_Macrophage 0.96860385
## Cluster_1 1 MES_AC_Transitional 0.99516831
## Cluster_2 2 NPC-like_mature 0.36287551
## Cluster_3 3 MES-like 0.70574136
## Cluster_4 4 NPC-like 0.54106981
## Cluster_5 5 Macrophage_complement 0.83897119
## Cluster_6 6 AC-like 0.37392569
## Cluster_7 7 Cycling_Tumor 0.71901714
## Cluster_8 8 Oligodendrocyte 0.99173847
## Cluster_9 9 Endothelial_Astrocyte 0.27615235
## Cluster_10 10 T_cell 0.74076172
## Cluster_11 11 Erythrocyte 0.04636565
## Cluster_12 12 Macrophage_LAM 0.58222749
## Cluster_13 13 Pericyte_Stromal 0.29241785# Compare all three methods
method_comparison <- data.frame(
cluster = stability_summary$cluster,
cell_state = stability_summary$cell_state,
bootstrap = stability_summary$mean_stability,
subsampling = subsample_summary$mean_stability_subsample[
match(stability_summary$cluster, subsample_summary$cluster)],
noise_injection = noise_summary$mean_stability_noise[
match(stability_summary$cluster, noise_summary$cluster)]
)
print(method_comparison)## cluster cell_state bootstrap subsampling noise_injection
## 1 1 MES_AC_Transitional 0.98911648 0.99386861 0.99516831
## 2 8 Oligodendrocyte 0.98349564 0.99060631 0.99173847
## 3 0 Microglia_Macrophage 0.97356710 0.97032259 0.96860385
## 4 5 Macrophage_complement 0.89096841 0.87221898 0.83897119
## 5 3 MES-like 0.82165614 0.88355256 0.70574136
## 6 7 Cycling_Tumor 0.72776413 0.71971277 0.71901714
## 7 10 T_cell 0.69932756 0.64882014 0.74076172
## 8 12 Macrophage_LAM 0.60233116 0.68801379 0.58222749
## 9 4 NPC-like 0.47614106 0.48336639 0.54106981
## 10 2 NPC-like_mature 0.38264315 0.41803241 0.36287551
## 11 6 AC-like 0.29181945 0.41858927 0.37392569
## 12 9 Endothelial_Astrocyte 0.27643262 0.27551556 0.27615235
## 13 11 Erythrocyte 0.04222826 0.04113362 0.04636565
## 14 13 Pericyte_Stromal 0.03430782 0.10236947 0.29241785write.csv(method_comparison, "method_comparison.csv", row.names = FALSE)
# Visualize comparison
comparison_long <- method_comparison %>%
pivot_longer(
cols = c(bootstrap, subsampling, noise_injection),
names_to = "method",
values_to = "stability"
)
comparison_plot <- ggplot(
comparison_long,
aes(x = reorder(cell_state, stability, FUN = mean), y = stability, fill = method)
) +
geom_bar(stat = "identity", position = "dodge", width = 0.7) +
coord_flip() +
scale_fill_brewer(palette = "Set2") +
labs(
title = "Comparison of Stability Across Resampling Methods",
x = "Cell State",
y = "Mean Stability",
fill = "Method"
) +
theme_publication +
theme(legend.position = "bottom")
print(comparison_plot)
ggsave("method_comparison.png", comparison_plot, width = 12, height = 8, dpi = 300)
# Correlation between methods
cor_bootstrap_subsample <- cor(
method_comparison$bootstrap,
method_comparison$subsampling,
method = "spearman"
)
cor_bootstrap_noise <- cor(
method_comparison$bootstrap,
method_comparison$noise_injection,
method = "spearman"
)
cat(paste("\nSpearman correlation (Bootstrap vs Subsampling):", round(cor_bootstrap_subsample, 3)))##
## Spearman correlation (Bootstrap vs Subsampling): 0.982##
## Spearman correlation (Bootstrap vs Noise): 0.96515. Create Comprehensive Summary Report
# Create comprehensive summary table
summary_table <- combined_metrics %>%
dplyr::select(
cluster,
cell_state,
n_cells,
mean_stability,
stability_category,
dissolution_freq,
specificity_score,
mean_correlation,
mean_silhouette
) %>%
arrange(desc(mean_stability))
# Format for presentation
summary_table_formatted <- summary_table %>%
mutate(
n_cells = as.integer(n_cells),
mean_stability = round(mean_stability, 3),
dissolution_freq = round(dissolution_freq * 100, 1),
specificity_score = round(specificity_score, 2),
mean_correlation = round(mean_correlation, 3),
mean_silhouette = round(mean_silhouette, 3)
)
print(summary_table_formatted)## cluster cell_state n_cells mean_stability stability_category
## 1 1 MES_AC_Transitional 2176 0.989 Highly Stable
## 2 8 Oligodendrocyte 420 0.983 Highly Stable
## 3 0 Microglia_Macrophage 3997 0.974 Highly Stable
## 4 5 Macrophage_complement 1228 0.891 Highly Stable
## 5 3 MES-like 1590 0.822 Stable
## 6 7 Cycling_Tumor 641 0.728 Moderately Stable
## 7 10 T_cell 207 0.699 Moderately Stable
## 8 12 Macrophage_LAM 125 0.602 Moderately Stable
## 9 4 NPC-like 1548 0.476 Unstable
## 10 2 NPC-like_mature 1881 0.383 Unstable
## 11 6 AC-like 749 0.292 Unstable
## 12 9 Endothelial_Astrocyte 252 0.276 Unstable
## 13 11 Erythrocyte 148 0.042 Unstable
## 14 13 Pericyte_Stromal 110 0.034 Unstable
## dissolution_freq specificity_score mean_correlation mean_silhouette
## 1 0 3.44 0.175 0.439
## 2 0 7.29 0.332 0.463
## 3 0 4.24 0.089 0.082
## 4 6 2.18 0.089 0.299
## 5 14 2.47 0.066 0.257
## 6 0 2.95 0.201 0.440
## 7 34 6.51 0.291 0.537
## 8 40 4.23 0.491 0.717
## 9 54 1.98 0.066 0.227
## 10 98 1.85 0.042 -0.064
## 11 91 2.13 0.106 0.300
## 12 100 3.86 0.276 0.590
## 13 100 8.03 0.219 0.475
## 14 99 3.90 0.125 0.192# Save as formatted table
write.csv(summary_table_formatted, "final_summary_table.csv", row.names = FALSE)
# Create a nice HTML table
library(kableExtra)
summary_table_formatted %>%
kbl(
caption = "Comprehensive Cluster Stability and Quality Metrics",
col.names = c(
"Cluster",
"Cell State",
"N Cells",
"Mean Stability",
"Category",
"Dissolution (%)",
"Specificity",
"Homogeneity",
"Silhouette"
)
) %>%
kable_styling(
bootstrap_options = c("striped", "hover", "condensed"),
full_width = FALSE
) %>%
column_spec(5, color = "white", background = spec_color(
as.numeric(factor(summary_table_formatted$stability_category)),
palette = c("#d73027", "#fee090", "#91bfdb", "#4575b4")
))| Cluster | Cell State | N Cells | Mean Stability | Category | Dissolution (%) | Specificity | Homogeneity | Silhouette |
|---|---|---|---|---|---|---|---|---|
| 1 | MES_AC_Transitional | 2176 | 0.989 | Highly Stable | 0 | 3.44 | 0.175 | 0.439 |
| 8 | Oligodendrocyte | 420 | 0.983 | Highly Stable | 0 | 7.29 | 0.332 | 0.463 |
| 0 | Microglia_Macrophage | 3997 | 0.974 | Highly Stable | 0 | 4.24 | 0.089 | 0.082 |
| 5 | Macrophage_complement | 1228 | 0.891 | Highly Stable | 6 | 2.18 | 0.089 | 0.299 |
| 3 | MES-like | 1590 | 0.822 | Stable | 14 | 2.47 | 0.066 | 0.257 |
| 7 | Cycling_Tumor | 641 | 0.728 | Moderately Stable | 0 | 2.95 | 0.201 | 0.440 |
| 10 | T_cell | 207 | 0.699 | Moderately Stable | 34 | 6.51 | 0.291 | 0.537 |
| 12 | Macrophage_LAM | 125 | 0.602 | Moderately Stable | 40 | 4.23 | 0.491 | 0.717 |
| 4 | NPC-like | 1548 | 0.476 | Unstable | 54 | 1.98 | 0.066 | 0.227 |
| 2 | NPC-like_mature | 1881 | 0.383 | Unstable | 98 | 1.85 | 0.042 | -0.064 |
| 6 | AC-like | 749 | 0.292 | Unstable | 91 | 2.13 | 0.106 | 0.300 |
| 9 | Endothelial_Astrocyte | 252 | 0.276 | Unstable | 100 | 3.86 | 0.276 | 0.590 |
| 11 | Erythrocyte | 148 | 0.042 | Unstable | 100 | 8.03 | 0.219 | 0.475 |
| 13 | Pericyte_Stromal | 110 | 0.034 | Unstable | 99 | 3.90 | 0.125 | 0.192 |
# Create final annotated UMAP with stability information
gbm_seurat$stability <- stability_summary$mean_stability[match(
as.character(Idents(gbm_seurat)),
stability_summary$cluster
)]
# UMAP colored by stability
umap_stability <- FeaturePlot(
gbm_seurat,
features = "stability",
reduction = "umap",
pt.size = 0.5
) +
scale_color_viridis(option = "plasma") +
labs(
title = "UMAP Colored by Cluster Stability",
color = "Jaccard\nCoefficient"
) +
theme_publication
print(umap_stability)
ggsave("umap_colored_by_stability.png", umap_stability, width = 10, height = 8, dpi = 300)
# Side-by-side: clusters and stability
final_comparison <- (
DimPlot(
gbm_seurat,
reduction = "umap",
group.by = "cell_state",
label = TRUE,
repel = TRUE
) +
ggtitle("Cell States") +
theme_publication
) | (
FeaturePlot(
gbm_seurat,
features = "stability",
reduction = "umap"
) +
scale_color_viridis(option = "plasma") +
ggtitle("Cluster Stability") +
theme_publication
)
print(final_comparison)
ggsave("final_comparison_states_stability.png", final_comparison, width = 16, height = 7, dpi = 300)16. Tumor vs. Non-Tumor Separation
Cluster yang mengandung sel non-tumor (TME, sel normal, kontaminasi) harus dipisahkan sebelum Neftel state annotation dilakukan. Ini adalah praktek standar dalam publikasi GBM single-cell (lihat: Neftel et al. 2019, Couturier et al. 2020).
# ── Definisi kategori sel ──────────────────────────────────────────────────────
tumor_clusters <- c("0" = FALSE, # Microglia — TME
"1" = TRUE, # MES-AC Transitional — Tumor
"2" = TRUE, # NPC-like mature — Tumor
"3" = TRUE, # MES-like — Tumor
"4" = TRUE, # NPC-like — Tumor
"5" = FALSE, # Macrophage — TME
"6" = TRUE, # AC-like — Tumor
"7" = TRUE, # Cycling Tumor — Tumor
"8" = FALSE, # Oligodendrocyte — Normal brain
"9" = FALSE, # Endothelial/Astrocyte — TME/Normal
"10" = FALSE, # T cell — TME
"11" = FALSE, # Erythrocyte — Kontaminasi
"12" = FALSE, # Macrophage LAM — TME
"13" = FALSE) # Pericyte/Stromal — TME
# Konversi factor → character dulu (sama seperti corrected_annotation)
cluster_vec <- as.character(gbm_10x@meta.data$seurat_clusters)
gbm_10x@meta.data$is_tumor <- tumor_clusters[cluster_vec]
# Kategori broad — gunakan cluster_vec agar tidak ada masalah factor
gbm_10x@meta.data$cell_category <- case_when(
cluster_vec %in% c("0", "5", "12") ~ "Microglia_Macrophage",
cluster_vec %in% c("10") ~ "T_cell",
cluster_vec %in% c("8") ~ "Oligodendrocyte",
cluster_vec %in% c("9", "13") ~ "Vascular_Stromal",
cluster_vec %in% c("11") ~ "Erythrocyte",
gbm_10x@meta.data$is_tumor == TRUE ~ "Tumor_cell",
TRUE ~ "Unclassified"
)
# Sync ke gbm_seurat
gbm_seurat@meta.data$is_tumor <- gbm_10x@meta.data$is_tumor
gbm_seurat@meta.data$cell_category <- gbm_10x@meta.data$cell_category
cat("Komposisi sel berdasarkan kategori:\n")## Komposisi sel berdasarkan kategori:##
## Erythrocyte Microglia_Macrophage Oligodendrocyte
## 148 5350 420
## T_cell Tumor_cell Vascular_Stromal
## 207 8585 362##
## Jumlah tumor cells: 8585## Jumlah non-tumor cells: 6487cat("Proporsi tumor cells:",
round(mean(gbm_10x@meta.data$is_tumor == TRUE, na.rm = TRUE) * 100, 1), "%\n")## Proporsi tumor cells: 57 %# ── Plot 1: Tumor vs non-tumor ─────────────────────────────────────────────────
p_tumor <- DimPlot(
gbm_10x,
group.by = "is_tumor",
reduction = "umap",
cols = c("TRUE" = "#E63946", "FALSE" = "#457B9D"),
pt.size = 0.3
) +
ggtitle("Tumor vs. Non-Tumor Cells") +
labs(color = "Is Tumor") +
theme_minimal(base_size = 12)
# ── Plot 2: Cell category ──────────────────────────────────────────────────────
category_colors <- c(
"Tumor_cell" = "#E63946",
"Microglia_Macrophage" = "#2A9D8F",
"T_cell" = "#E9C46A",
"Oligodendrocyte" = "#264653",
"Vascular_Stromal" = "#F4A261",
"Erythrocyte" = "#A8DADC",
"Unclassified" = "#CCCCCC"
)
p_category <- DimPlot(
gbm_10x,
group.by = "cell_category",
reduction = "umap",
cols = category_colors,
label = TRUE,
repel = TRUE,
label.size = 3,
pt.size = 0.3
) +
ggtitle("Cell Category (Tumor + TME)") +
theme_minimal(base_size = 12) +
theme(legend.position = "right")
print(p_tumor)
ggsave("umap_tumor_vs_tme.png", p_tumor, width = 10, height = 7, dpi = 300)
ggsave("umap_cell_category.png", p_category, width = 12, height = 8, dpi = 300)17. Neftel State Annotation — Tumor Cells Only
Setelah pemisahan tumor vs TME, Neftel state annotation dilakukan hanya pada tumor cells. Ini mencegah signature macrophage (MES-like palsu) atau oligodendrocyte (OPC-like palsu) mengkontaminasi interpretasi GBM cellular states.
# Subset hanya tumor cells
# subset() via @meta.data agar tidak kena masalah factor
tumor_barcodes <- rownames(gbm_10x@meta.data)[
gbm_10x@meta.data$is_tumor == TRUE & !is.na(gbm_10x@meta.data$is_tumor)
]
gbm_tumor <- subset(gbm_10x, cells = tumor_barcodes)
cat("Jumlah sel dalam tumor subset:", ncol(gbm_tumor), "\n")## Jumlah sel dalam tumor subset: 8585cat("Cluster dalam tumor subset:",
paste(sort(unique(as.character(gbm_tumor$seurat_clusters))), collapse = ", "), "\n")## Cluster dalam tumor subset: 1, 2, 3, 4, 6, 7# ── Neftel et al. 2019 — Gene Signatures ──────────────────────────────────────
# Sumber: Supplementary Table S3, Neftel et al. 2019 Cell
neftel_signatures <- list(
MES1 = c("CHI3L1", "CD44", "VIM", "ANXA2", "TGFBI", "SERPINE1",
"FN1", "MMP2", "CTGF", "ITGA5"),
MES2 = c("TIMP1", "LGALS1", "S100A4", "S100A10", "CTSB",
"LOXL2", "CCL2", "CXCL2", "MMP7", "PLAUR"),
AC = c("GFAP", "S100B", "AQP4", "ALDOC", "SLC1A3",
"GLUL", "NDRG2", "SOX9", "FABP7", "AGT"),
OPC = c("OLIG1", "OLIG2", "PDGFRA", "SOX10", "CSPG4",
"GPR17", "APOD", "CCND1", "ASCL1", "NKX2-2"),
NPC1 = c("SOX4", "SOX11", "DCX", "STMN1", "TUBA1A",
"TUBB2B", "NNAT", "MLLT11", "NOVA1", "MAP1B"),
NPC2 = c("ASCL1", "DLL3", "HES6", "EGFR", "CCND2",
"INSM1", "NEUROG1", "DLX2", "SRRM4", "PCBP4")
)
# Hitung module scores pada tumor cells
DefaultAssay(gbm_tumor) <- "RNA"
for (state in names(neftel_signatures)) {
gbm_tumor <- AddModuleScore(
object = gbm_tumor,
features = list(neftel_signatures[[state]]),
name = paste0(state, "_score"),
seed = 42
)
# AddModuleScore menambahkan suffix "1", rename agar bersih
old_name <- paste0(state, "_score1")
new_name <- paste0(state, "_score")
gbm_tumor[[new_name]] <- gbm_tumor[[old_name]]
gbm_tumor[[old_name]] <- NULL
}
cat("Module scores ditambahkan:\n")## Module scores ditambahkan:score_cols <- paste0(names(neftel_signatures), "_score")
print(summary(gbm_tumor@meta.data[, score_cols]))## MES1_score MES2_score AC_score OPC_score
## Min. :-0.6030 Min. :-0.6223 Min. :-0.94605 Min. :-0.35212
## 1st Qu.:-0.2791 1st Qu.:-0.2508 1st Qu.:-0.40037 1st Qu.:-0.14399
## Median : 0.2301 Median : 0.1161 Median :-0.14291 Median :-0.04971
## Mean : 0.2214 Mean : 0.1553 Mean :-0.08285 Mean : 0.02356
## 3rd Qu.: 0.6339 3rd Qu.: 0.4673 3rd Qu.: 0.15689 3rd Qu.: 0.08491
## Max. : 1.9514 Max. : 2.1572 Max. : 2.01690 Max. : 2.14257
## NPC1_score NPC2_score
## Min. :-0.9514 Min. :-0.495393
## 1st Qu.:-0.1370 1st Qu.:-0.197727
## Median : 0.2053 Median :-0.096281
## Mean : 0.4602 Mean :-0.003603
## 3rd Qu.: 0.9112 3rd Qu.: 0.086968
## Max. : 2.7944 Max. : 1.535921# ── Assign dominant Neftel state per cell ─────────────────────────────────────
score_matrix <- gbm_tumor@meta.data[, score_cols]
colnames(score_matrix) <- names(neftel_signatures)
# State dengan skor tertinggi = dominant state
gbm_tumor$neftel_state <- apply(score_matrix, 1, function(x) {
names(x)[which.max(x)]
})
# Hitung confidence: selisih antara skor tertinggi dan kedua tertinggi
gbm_tumor$neftel_confidence <- apply(score_matrix, 1, function(x) {
sorted <- sort(x, decreasing = TRUE)
sorted[1] - sorted[2] # margin skor tertinggi vs kedua
})
# Tandai sel dengan margin rendah sebagai "Mixed"
confidence_threshold <- 0.1
gbm_tumor$neftel_state_final <- ifelse(
gbm_tumor$neftel_confidence < confidence_threshold,
"Mixed",
gbm_tumor$neftel_state
)
cat("\nDistribusi Neftel states (raw dominant):\n")##
## Distribusi Neftel states (raw dominant):##
## AC MES1 MES2 NPC1 NPC2 OPC
## 953 2761 1284 3092 110 385##
## Distribusi Neftel states (dengan Mixed threshold = 0.1 ):##
## AC MES1 MES2 Mixed NPC1 NPC2 OPC
## 790 2163 873 1676 2779 50 254# ── Visualisasi Neftel states pada tumor cells ─────────────────────────────────
neftel_colors <- c(
"MES1" = "#D62728",
"MES2" = "#FF7F0E",
"AC" = "#2CA02C",
"OPC" = "#1F77B4",
"NPC1" = "#9467BD",
"NPC2" = "#8C564B",
"Mixed" = "#CCCCCC"
)
p_neftel <- DimPlot(
gbm_tumor,
group.by = "neftel_state_final",
reduction = "umap",
cols = neftel_colors,
label = TRUE,
repel = TRUE,
label.size = 4,
pt.size = 0.5
) +
ggtitle("Neftel et al. 2019 Cellular States\n(Tumor Cells Only)") +
theme_minimal(base_size = 12) +
theme(legend.position = "right")
print(p_neftel)
ggsave("umap_neftel_states_tumor_only.png", p_neftel,
width = 11, height = 8, dpi = 300)
# ── Score distribution per state ──────────────────────────────────────────────
score_long <- gbm_tumor@meta.data %>%
dplyr::select(neftel_state_final, dplyr::all_of(score_cols)) %>%
tidyr::pivot_longer(
cols = dplyr::all_of(score_cols),
names_to = "signature",
values_to = "score"
) %>%
dplyr::mutate(signature = gsub("_score", "", signature))
p_scores <- ggplot(score_long, aes(x = signature, y = score, fill = signature)) +
geom_violin(trim = TRUE, alpha = 0.8) +
geom_boxplot(width = 0.1, fill = "white", alpha = 0.7, outlier.size = 0.3) +
facet_wrap(~ neftel_state_final, scales = "free_y") +
scale_fill_manual(values = neftel_colors) +
labs(
title = "Neftel Signature Score Distribution per Assigned State",
x = "Signature",
y = "Module Score"
) +
theme_minimal(base_size = 11) +
theme(
legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1)
)
print(p_scores)
# ── Proporsi state per sample
sample_col <- "orig.ident"
min_cells <- 50 # eksklusi sample dengan tumor cells terlalu sedikit
if (sample_col %in% colnames(gbm_tumor@meta.data)) {
prop_df <- gbm_tumor@meta.data %>%
dplyr::group_by(.data[[sample_col]], neftel_state_final) %>%
dplyr::summarise(n = n(), .groups = "drop") %>%
dplyr::group_by(.data[[sample_col]]) %>%
dplyr::mutate(
proportion = n / sum(n),
total_cells = sum(n)
) %>%
dplyr::filter(total_cells >= min_cells)
cat("Sample yang diinklusi (≥", min_cells, "tumor cells):\n")
print(unique(prop_df[[sample_col]]))
p_prop <- ggplot(
prop_df,
aes(x = .data[[sample_col]], y = proportion, fill = neftel_state_final)
) +
geom_bar(stat = "identity") +
scale_fill_manual(values = neftel_colors) +
scale_y_continuous(labels = scales::percent) +
labs(
title = "Neftel State Proportions per Sample",
subtitle = paste("Sampel dengan <", min_cells, "tumor cells dieksklusi"),
x = "Sample",
y = "Proportion",
fill = "Neftel State"
) +
theme_minimal(base_size = 12) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
print(p_prop)
ggsave("neftel_state_proportion_per_sample.png", p_prop,
width = 12, height = 7, dpi = 300)
} else {
cat("Kolom sample '", sample_col, "' tidak ditemukan. Lewati plot proporsi.\n")
}## Sample yang diinklusi (≥ 50 tumor cells):
## [1] 102 105 105A 115 118 124 125 143
## Levels: 102 105 105A 114 115 118 124 125 126 143
##
## 102 105 105A 114 115 118 124 125 126 143
## 1236 1608 643 39 645 236 1329 649 19 2181# ── Simpan objek final ─────────────────────────────────────────────────────────
saveRDS(gbm_seurat, "gbm_seurat_corrected_annotation.rds")
saveRDS(gbm_tumor, "gbm_tumor_neftel_annotated.rds")
# Export metadata
write.csv(
gbm_seurat@meta.data,
"metadata_all_cells_corrected.csv",
row.names = TRUE
)
write.csv(
gbm_tumor@meta.data,
"metadata_tumor_neftel.csv",
row.names = TRUE
)
cat("\n=== Annotation Correction & Neftel Analysis Complete ===\n")##
## === Annotation Correction & Neftel Analysis Complete ===## Output files:## - gbm_seurat_corrected_annotation.rds : Full dataset, corrected labels## - gbm_tumor_neftel_annotated.rds : Tumor-only, Neftel states## - metadata_all_cells_corrected.csv : Metadata semua sel## - metadata_tumor_neftel.csv : Metadata tumor cells + Neftel## - umap_corrected_annotation.png## - umap_tumor_vs_tme.png## - umap_cell_category.png## - umap_neftel_states_tumor_only.png## - neftel_score_distribution.png## - neftel_state_proportion_per_sample.png18. Save Final Objects
# Save semua objek final dengan anotasi lengkap
saveRDS(gbm_10x, "gbm_10x_annotated_final.rds")
saveRDS(gbm_tumor, "gbm_tumor_neftel_annotated.rds")
# Save stability results
saveRDS(
list(
stability_summary = stability_summary,
stability_matrix = stability_matrix,
method_comparison = method_comparison,
combined_metrics = combined_metrics
),
"stability_analysis_results.rds"
)
# Export metadata
write.csv(gbm_10x@meta.data, "metadata_all_cells.csv", row.names = TRUE)
write.csv(gbm_tumor@meta.data, "metadata_tumor_neftel.csv", row.names = TRUE)
cat("\n=== Analysis Complete ===\n")##
## === Analysis Complete ===## Key output files:## - gbm_10x_annotated_final.rds : Full dataset, corrected labels## - gbm_tumor_neftel_annotated.rds : Tumor-only, Neftel states## - stability_analysis_results.rds : Bootstrap stability results## - metadata_all_cells.csv## - metadata_tumor_neftel.csv19. Session Information
## R version 4.5.2 (2025-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 10 x64 (build 19045)
##
## Matrix products: default
## LAPACK version 3.12.1
##
## locale:
## [1] LC_COLLATE=English_Indonesia.utf8 LC_CTYPE=English_Indonesia.utf8
## [3] LC_MONETARY=English_Indonesia.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_Indonesia.utf8
##
## time zone: Asia/Jakarta
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] kableExtra_1.4.0 FNN_1.1.4.1 presto_1.0.0
## [4] Rcpp_1.1.1 future.apply_1.20.1 future_1.69.0
## [7] data.table_1.18.2.1 Matrix_1.7-4 msigdbr_25.1.1
## [10] enrichplot_1.30.4 org.Hs.eg.db_3.18.0 AnnotationDbi_1.72.0
## [13] IRanges_2.44.0 S4Vectors_0.48.0 Biobase_2.70.0
## [16] BiocGenerics_0.56.0 generics_0.1.4 clusterProfiler_4.18.4
## [19] cowplot_1.2.0 pheatmap_1.0.13 RColorBrewer_1.1-3
## [22] viridis_0.6.5 viridisLite_0.4.3 ggridges_0.5.7
## [25] igraph_2.2.2 cluster_2.1.8.1 fpc_2.2-14
## [28] patchwork_1.3.2 lubridate_1.9.5 forcats_1.0.1
## [31] stringr_1.6.0 dplyr_1.2.0 purrr_1.2.1
## [34] readr_2.1.6 tidyr_1.3.2 tibble_3.3.1
## [37] ggplot2_4.0.2 tidyverse_2.0.0 Seurat_5.4.0
## [40] SeuratObject_5.3.0 sp_2.2-1
##
## loaded via a namespace (and not attached):
## [1] fs_1.6.6 matrixStats_1.5.0 spatstat.sparse_3.1-0
## [4] httr_1.4.8 prabclus_2.3-5 tools_4.5.2
## [7] sctransform_0.4.3 utf8_1.2.6 R6_2.6.1
## [10] mgcv_1.9-3 lazyeval_0.2.2 uwot_0.2.4
## [13] withr_3.0.2 gridExtra_2.3 progressr_0.18.0
## [16] textshaping_1.0.4 cli_3.6.5 spatstat.explore_3.7-0
## [19] fastDummies_1.7.5 scatterpie_0.2.6 labeling_0.4.3
## [22] sass_0.4.10 diptest_0.77-2 S7_0.2.1
## [25] robustbase_0.99-7 spatstat.data_3.1-9 pbapply_1.7-4
## [28] systemfonts_1.3.1 yulab.utils_0.2.4 gson_0.1.0
## [31] svglite_2.2.2 DOSE_4.4.0 R.utils_2.13.0
## [34] parallelly_1.46.1 rstudioapi_0.18.0 RSQLite_2.4.6
## [37] gridGraphics_0.5-1 ica_1.0-3 spatstat.random_3.4-4
## [40] GO.db_3.22.0 abind_1.4-8 R.methodsS3_1.8.2
## [43] lifecycle_1.0.5 yaml_2.3.12 qvalue_2.42.0
## [46] Rtsne_0.17 grid_4.5.2 blob_1.3.0
## [49] promises_1.5.0 crayon_1.5.3 miniUI_0.1.2
## [52] ggtangle_0.1.1 lattice_0.22-7 KEGGREST_1.50.0
## [55] pillar_1.11.1 knitr_1.51 fgsea_1.36.2
## [58] codetools_0.2-20 fastmatch_1.1-8 glue_1.8.0
## [61] ggiraph_0.9.5 ggfun_0.2.0 spatstat.univar_3.1-6
## [64] fontLiberation_0.1.0 vctrs_0.7.1 png_0.1-8
## [67] treeio_1.34.0 spam_2.11-3 gtable_0.3.6
## [70] assertthat_0.2.1 kernlab_0.9-33 cachem_1.1.0
## [73] xfun_0.56 mime_0.13 Seqinfo_1.0.0
## [76] survival_3.8-3 fitdistrplus_1.2-6 ROCR_1.0-12
## [79] nlme_3.1-168 ggtree_4.0.4 bit64_4.6.0-1
## [82] fontquiver_0.2.1 RcppAnnoy_0.0.23 bslib_0.10.0
## [85] irlba_2.3.7 KernSmooth_2.23-26 otel_0.2.0
## [88] DBI_1.2.3 nnet_7.3-20 processx_3.8.6
## [91] tidyselect_1.2.1 curl_7.0.0 bit_4.6.0
## [94] compiler_4.5.2 xml2_1.5.2 fontBitstreamVera_0.1.1
## [97] plotly_4.12.0 scales_1.4.0 DEoptimR_1.1-4
## [100] lmtest_0.9-40 callr_3.7.6 rappdirs_0.3.4
## [103] digest_0.6.39 goftest_1.2-3 spatstat.utils_3.2-1
## [106] rmarkdown_2.30 XVector_0.50.0 htmltools_0.5.9
## [109] pkgconfig_2.0.3 fastmap_1.2.0 rlang_1.1.7
## [112] htmlwidgets_1.6.4 shiny_1.12.1 farver_2.1.2
## [115] jquerylib_0.1.4 zoo_1.8-15 jsonlite_2.0.0
## [118] BiocParallel_1.44.0 mclust_6.1.2 GOSemSim_2.36.0
## [121] R.oo_1.27.1 magrittr_2.0.4 modeltools_0.2-24
## [124] ggplotify_0.1.3 dotCall64_1.2 babelgene_22.9
## [127] ape_5.8-1 ggnewscale_0.5.2 gdtools_0.5.0
## [130] reticulate_1.45.0 stringi_1.8.7 MASS_7.3-65
## [133] pkgbuild_1.4.8 plyr_1.8.9 flexmix_2.3-20
## [136] parallel_4.5.2 listenv_0.10.0 ggrepel_0.9.6
## [139] deldir_2.0-4 Biostrings_2.78.0 splines_4.5.2
## [142] tensor_1.5.1 hms_1.1.4 ps_1.9.1
## [145] spatstat.geom_3.7-0 RcppHNSW_0.6.0 reshape2_1.4.5
## [148] evaluate_1.0.5 BiocManager_1.30.27 tzdb_0.5.0
## [151] tweenr_2.0.3 httpuv_1.6.16 RANN_2.6.2
## [154] polyclip_1.10-7 scattermore_1.2 ggforce_0.5.0
## [157] xtable_1.8-4 RSpectra_0.16-2 tidytree_0.4.7
## [160] tidydr_0.0.6 later_1.4.6 ragg_1.5.0
## [163] class_7.3-23 snow_0.4-4 aplot_0.2.9
## [166] memoise_2.0.1 timechange_0.4.0 globals_0.19.0Interpretation Guide
Stability Categories
Based on Hennig (2007):
Highly Stable (>0.85): Clusters that consistently appear across bootstrap iterations with high overlap. These likely represent distinct, biologically meaningful cell populations.
Stable (0.75-0.85): Clusters that appear consistently but with some variation in composition. Still reliable for biological interpretation.
Moderately Stable (0.6-0.75): Clusters that measure some pattern but with considerable uncertainty about exact cell composition. Interpret with caution.
Unstable (<0.6): Clusters that vary substantially or dissolve frequently. May represent:
- Highly heterogeneous populations
- Transition states between other clusters
- Technical artifacts of the clustering algorithm
- Cells that don’t fit well into other clusters
Biological Interpretation
For Stable Clusters:
- High confidence in cell state assignments
- Suitable for identifying therapeutic targets
- Marker genes are reliable for validation
- Can be used to develop signatures for other datasets
For Unstable Clusters:
- May represent dynamic or plastic cell states
- Could indicate transition zones between committed states
- Requires additional validation (e.g., trajectory analysis, spatial transcriptomics)
- Consider sub-clustering or alternative analysis approaches
Clinical Implications
Therapeutic Targeting: - Prioritize stable clusters for drug target identification - Consider combination therapies for heterogeneous tumors - Unstable clusters may explain therapeutic resistance through plasticity
Biomarker Development: - Stable cluster signatures are more likely to be reproducible - Consider stability when selecting marker genes for clinical assays
Conclusion
This analysis provides a comprehensive assessment of cluster stability in glioblastoma scRNA-seq data. The results help distinguish between robust biological populations and potentially artifactual clusters, providing a foundation for reliable biological interpretation and therapeutic target identification.
The framework can be adapted for other cancer types or tissues where single-cell heterogeneity is a key feature.