SINGLE-CELL RNA-SEQ, DE ANALYSIS + PATHWAY ENRICHMENT
This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.
1.0 PACKAGE INSTALLATIONS
################################################################################
# scRNA-seq Differential Expression Pipeline
# Dataset : GSE161277 — Human Colorectal Carcinogenesis (Zheng et al., 2022)
# 4 patients | matched carcinoma, normal, adenoma, para-cancer, blood
# 10x Genomics | sparse MTX format (~338 MB download)
# Goal : Abnormal vs. Normal DE within epithelial cells
# Approach: (1) Seurat QC + normalization + clustering
# (2) Pseudobulk DE with DESeq2 <- statistically rigorous
# (3) Single-cell DE with MAST <- accounts for dropout
# Citation: Zheng X et al. (2022) Adv Sci. PMID: 35221332
################################################################################
# ==============================================================================
# 0. INSTALL / LOAD PACKAGES
# ==============================================================================
##install.packages(c("Seurat", "ggplot2", "dplyr", "pheatmap", "ggrepel",
## "patchwork", "Matrix", 'glmGamPoi'))
# -----------------------------
# Run ONCE on a new computer
# -----------------------------
print("Resource: Zheng X, Song J, Yu C, Zhou Z et al. Single-cell transcriptomic profiling unravels the adenoma-initiation role of protein tyrosine kinases during colorectal tumorigenesis. Signal Transduct Target Ther 2022 Feb 28;7(1):60. PMID: 35221332")## [1] "Resource: Zheng X, Song J, Yu C, Zhou Z et al. Single-cell transcriptomic profiling unravels the adenoma-initiation role of protein tyrosine kinases during colorectal tumorigenesis. Signal Transduct Target Ther 2022 Feb 28;7(1):60. PMID: 35221332"if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install(
c(
"GEOquery", "Seurat", "ggplot2", "dplyr", "pheatmap",
"ggrepel", "patchwork", "Matrix", "glmGamPoi",
"DESeq2", "MAST", "fgsea", "org.Hs.eg.db",
"clusterProfiler"
),
update = TRUE,
ask = FALSE
)
# -----------------------------
# Run EVERY session
# -----------------------------
suppressPackageStartupMessages({
library(Seurat)
library(GEOquery)
library(DESeq2)
library(MAST)
library(tibble)
library(ggplot2)
library(ggrepel)
library(dplyr)
library(pheatmap)
library(patchwork)
library(Matrix)
library(glmGamPoi)
library(clusterProfiler)
library(org.Hs.eg.db)
})
set.seed(42)2.0 DATASET UPLOAD AND METADATA SETUP
# ==============================================================================
# 1. DATA DOWNLOAD FROM GEO
# ==============================================================================
# GSE161277 provides per-sample 10x MEX files (matrix.mtx.gz, barcodes.tsv.gz,
# features.tsv.gz) bundled as GSE161277_RAW.tar (~338 MB).
# Sparse format: loads directly into memory without large temp files.
#
# Sample map (13 samples across 4 patients):
# GSM4904234 Patient0 carcinoma
# GSM4904235 Patient1 adenoma
# GSM4904236 Patient1 carcinoma
# GSM4904237 Patient1 normal
# GSM4904238 Patient2 adenoma
# GSM4904239 Patient2 carcinoma
# GSM4904240 Patient2 normal
# GSM4904241 Patient2 para-cancer
# GSM4904242 Patient3 adenoma_1
# GSM4904243 Patient3 adenoma_2
# GSM4904244 Patient3 blood
# GSM4904245 Patient3 carcinoma
# GSM4904246 Patient3 normal
#
# For Tumor vs. Normal DE we use the 3 matched pairs:
# Patient1: GSM4904236 (carcinoma) + GSM4904237 (normal)
# Patient2: GSM4904239 (carcinoma) + GSM4904240 (normal)
# Patient3: GSM4904245 (carcinoma) + GSM4904246 (normal)
# ==============================================================================
# 1. DATA DOWNLOAD FROM GEO
# ==============================================================================
# Expand ~ to the full system path (/Users/tinacole/GSE161277_data)
outdir <- path.expand("~/GSE161277_data")
tar_file <- file.path(outdir, "GSE161277_RAW.tar")
dir.create(outdir, showWarnings = FALSE, recursive = TRUE)
# Download only if the archive doesn't already exist locally
if (!file.exists(tar_file)) {
message("Downloading GSE161277_RAW.tar (~338 MB)...")
options(timeout = 900)
download.file(
url = "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE161277&format=file",
destfile = tar_file,
mode = "wb"
)
}
message("Extracting TAR...")
untar(tar_file, exdir = outdir)
# ==============================================================================
# 2. BUILD SAMPLE METADATA TABLE
# ==============================================================================
sample_meta <- data.frame(
gsm = c("GSM4904234", "GSM4904235", "GSM4904236", "GSM4904237",
"GSM4904238", "GSM4904239", "GSM4904240", "GSM4904241",
"GSM4904242", "GSM4904243", "GSM4904244", "GSM4904245",
"GSM4904246"),
patient = c("Patient0", "Patient1", "Patient1", "Patient1",
"Patient2", "Patient2", "Patient2", "Patient2",
"Patient3", "Patient3", "Patient3", "Patient3", "Patient3"),
condition = c("Tumor", "Adenoma", "Tumor", "Normal",
"Adenoma", "Tumor", "Normal", "Para-cancer",
"Adenoma", "Adenoma", "Blood", "Tumor", "Normal"),
stringsAsFactors = FALSE
)
de_samples <- sample_meta %>%
mutate(group = if_else(condition == "Normal", "Normal", "Abnormal"))3.0 LOAD & MERGE COUNT MATRICES
message("Loading per-sample 10x matrices...")
seurat_list <- lapply(seq_len(nrow(de_samples)), function(i) {
meta <- de_samples[i, ]
path <- file.path(outdir, meta$gsm)
mat <- Read10X(data.dir = path)
seu <- CreateSeuratObject(
counts = mat,
project = paste0(meta$patient, "_", meta$condition),
min.cells = 3,
min.features = 200
)
seu$patient <- meta$patient
seu$condition <- meta$condition
seu$group <- meta$group
seu$gsm <- meta$gsm
seu$sample_id <- paste0(meta$patient, "_", meta$condition)
seu <- RenameCells(seu, add.cell.id = meta$gsm)
seu
})
seu <- merge(seurat_list[[1]], y = seurat_list[-1], project = "GSE161277_CRC")
message(sprintf("Merged object: %d genes x %d cells", nrow(seu), ncol(seu)))
seu <- subset(seu, cells = sample(Cells(seu), 5000)) #subset after merger4.0: QC & NORMALIZATION
# ==============================================================================
# 4. QUALITY CONTROL
# ==============================================================================
seu[["percent.mt"]] <- PercentageFeatureSet(seu, pattern = "^MT-")
p_qc <- VlnPlot(seu,
features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
ncol = 3, pt.size = 0, group.by = "sample_id")
p_qc
ggsave("QC_violin.pdf", p_qc, width = 14, height = 5)
seu <- subset(seu,
subset = nFeature_RNA > 200 &
nFeature_RNA < 6000 &
percent.mt < 20)
message(sprintf("After QC: %d cells remain", ncol(seu)))
# ==============================================================================
# 5. NORMALIZATION & FEATURE SELECTION
# ==============================================================================
message("Running SCTransform...")
seu <- SCTransform(seu, vars.to.regress = "percent.mt",
variable.features.n = 3000,
ncells = 5000,
verbose = FALSE)6.0 DIMENSIONALITY REDUCTION & CLUSTERING
# ==============================================================================
# 6. DIMENSIONALITY REDUCTION, INTEGRATION & CLUSTERING
# ==============================================================================
# PCA
seu <- RunPCA(seu, npcs = 50, verbose = FALSE)
N_DIMS <- 30
# Harmony batch correction by patient
library(harmony)
seu <- RunHarmony(seu, group.by.vars = "sample_id")## Initializing centroids# UMAP, neighbors, clusters on Harmony embeddings
seu <- RunUMAP(seu, reduction = "harmony", dims = 1:N_DIMS, verbose = FALSE)
seu <- FindNeighbors(seu, reduction = "harmony", dims = 1:N_DIMS, verbose = FALSE)
seu <- FindClusters(seu, resolution = 0.5, verbose = FALSE)
# Visualize
DimPlot(seu, reduction = "umap", label = TRUE, pt.size = 0.5)
DimPlot(seu, reduction = "umap", group.by = "condition", pt.size = 0.5)
DimPlot(seu, reduction = "umap", group.by = "sample_id", pt.size = 0.5)
# Find markers on corrected clusters
DefaultAssay(seu) <- "SCT"
seu <- PrepSCTFindMarkers(seu)
markers <- FindAllMarkers(seu, only.pos = TRUE, min.pct = 0.1, logfc.threshold = 0.1)
markers %>% group_by(cluster) %>% top_n(5, avg_log2FC)7.0 ANNOTATION
# ==============================================================================
# 7. CELL-TYPE ANNOTATION (AUTOMATED FALLBACK INCLUDED)
# ==============================================================================
# Scoring canonical epithelial expression for robust automatic subsetting
epi_genes <- c("EPCAM", "KRT18", "KRT19", "CDH1")
valid_epi <- intersect(epi_genes, rownames(seu))
seu <- AddModuleScore(seu, features = list(valid_epi), name = "Epithelial_Score")
# Define cell clusters as Epithelial if mean score > 0
epi_clusters <- AverageExpression(seu,
features = valid_epi,
group.by = "seurat_clusters",
assay = "SCT",
layer = "data")$SCT
epi_cluster_ids <- colnames(epi_clusters)[colMeans(epi_clusters) > 0] # was 0.5
seu$Cell_type <- if_else(Idents(seu) %in% epi_cluster_ids, "Epithelial", "Other")
message(sprintf("Identified %d epithelial cells", sum(seu$Cell_type == "Epithelial")))
####Checking counts. ####
# Check to confirm
epi_cluster_ids # you'll see "g0", "g1" etc.## [1] "g0" "g1" "g2" "g3" "g4" "g5" "g6" "g7" "g8" "g9" "g10" "g11"# Fix
epi_cluster_ids_fixed <- sub("^g", "", epi_cluster_ids)
seu$Cell_type <- if_else(as.character(Idents(seu)) %in% epi_cluster_ids_fixed,
"Epithelial", "Other")
message(sprintf("Identified %d epithelial cells", sum(seu$Cell_type == "Epithelial")))
# Check the actual values
colMeans(epi_clusters)## g0 g1 g2 g3 g4 g5 g6
## 0.07916325 3.40343750 0.07942708 1.65858209 0.10218254 0.14975248 2.76977401
## g7 g8 g9 g10 g11
## 0.02192982 0.03437500 1.92549669 0.20220588 0.05128205# Check module score distribution
hist(seu$Epithelial_Score1, breaks = 50)
summary(seu$Epithelial_Score1)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -0.6512 -0.5319 -0.4157 -0.1373 0.2345 1.7927epi_cluster_ids_fixed <- sub("^g", "", colnames(epi_clusters)[colMeans(epi_clusters) > 1])
seu$Cell_type <- if_else(as.character(Idents(seu)) %in% epi_cluster_ids_fixed,
"Epithelial", "Other")
message(sprintf("Identified %d epithelial cells", sum(seu$Cell_type == "Epithelial")))8.1 PSEUDOBULK DE WITH DESeq2 (8.1)
# ==============================================================================
# 8A. PSEUDOBULK DE WITH DESeq2
# ==============================================================================
message("Running pseudobulk DE (Abnormal vs. Normal, epithelial cells)...")
seu_epi <- subset(seu, subset = Cell_type == "Epithelial")
DefaultAssay(seu_epi) <- "RNA"
# Ensure raw count extraction for pseudo-bulk aggregation
pb_counts <- AggregateExpression(
seu_epi,
assays = "RNA",
return.seurat = FALSE,
group.by = c("patient", "group"),
#layer = "counts" #replaced slot = "counts"
)$RNA
col_meta <- data.frame(
sample = colnames(pb_counts),
patient = factor(sub("_.*", "", colnames(pb_counts))),
group = factor(sub(".*_", "", colnames(pb_counts)), levels = c("Normal", "Abnormal")),
stringsAsFactors = FALSE
)
dds <- DESeqDataSetFromMatrix(
countData = round(pb_counts),
colData = col_meta,
design = ~ patient + group
)
keep_genes <- rowSums(counts(dds) >= 10) >= 2
dds <- dds[keep_genes, ]
dds <- DESeq(dds)
# Get specific coefficient name dynamically
res_coef <- resultsNames(dds)[grep("group_Abnormal_vs_Normal", resultsNames(dds))]
res_deseq_shrunk <- lfcShrink(dds, coef = res_coef, type = "apeglm")
sig_deseq <- as.data.frame(res_deseq_shrunk) %>%
tibble::rownames_to_column("gene") %>%
filter(!is.na(padj), padj < 0.05, abs(log2FoldChange) > 1) %>%
arrange(padj)
write.csv(sig_deseq, "DESeq2_results_Tumor_vs_Normal.csv", row.names = FALSE)8.2 DE: SINGLE-CELL WITH MAST
# ==============================================================================
# 8B. SINGLE-CELL DE WITH MAST
# ==============================================================================
message("Running MAST DE (Abnormal vs. Normal, epithelial cells)...")
DefaultAssay(seu_epi) <- "RNA"
seu_epi <- NormalizeData(seu_epi, verbose = FALSE)
# Seurat's internal MAST wrapper for stability across versions
Idents(seu_epi) <- "group"
mast_de <- FindMarkers(
seu_epi,
ident.1 = "Abnormal",
ident.2 = "Normal",
test.use = "MAST",
latent.vars = "patient",
logfc.threshold = 0.25,
min.pct = 0.1
)
sig_mast <- mast_de %>%
tibble::rownames_to_column("gene") %>%
filter(p_val_adj < 0.05, abs(avg_log2FC) > 0.5)
write.csv(sig_mast, "MAST_results_Tumor_vs_Normal.csv", row.names = FALSE)9.0 VISUALIZATIONS
# ==============================================================================
# 9. VISUALIZATIONS
# ==============================================================================
volcano_data <- as.data.frame(res_deseq_shrunk) %>%
tibble::rownames_to_column("gene") %>%
filter(!is.na(padj)) %>%
mutate(
sig = case_when(
padj < 0.05 & log2FoldChange > 1 ~ "Up in Abnormal",
padj < 0.05 & log2FoldChange < -1 ~ "Down in Abnormal",
TRUE ~ "NS"
),
label = ifelse(sig != "NS" & abs(log2FoldChange) > 2 & -log10(padj) > 5, gene, NA)
)
p_volcano <- ggplot(volcano_data, aes(x = log2FoldChange, y = -log10(padj), colour = sig, label = label)) +
geom_point(alpha = 0.6, size = 1.2) +
geom_text_repel(size = 2.5, max.overlaps = 15, colour = "black") +
scale_colour_manual(values = c("Up in Abnormal" = "#d62728", "Down in Abnormal" = "#1f77b4", "NS" = "grey60")) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed", colour = "grey40") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed", colour = "grey40") +
theme_bw(base_size = 12)
p_volcano
ggsave("Volcano_DESeq2.pdf", p_volcano, width = 7, height = 6)BONUS: PATHWAY ENRICHMENT AND CROSS METHOD COMPARISON
# ==============================================================================
# 10. PATHWAY ENRICHMENT (GO & KEGG)
# ==============================================================================
up_genes <- sig_deseq %>% filter(log2FoldChange > 1) %>% pull(gene)
down_genes <- sig_deseq %>% filter(log2FoldChange < -1) %>% pull(gene)
background <- rownames(dds)
convert_to_entrez <- function(genes) {
res <- bitr(genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
return(res$ENTREZID)
}
up_entrez <- convert_to_entrez(up_genes)
dn_entrez <- convert_to_entrez(down_genes)
bg_entrez <- convert_to_entrez(background)
# GO - Upregulated
if (length(up_entrez) > 0) {
go_up <- enrichGO(gene = up_entrez, universe = bg_entrez, OrgDb = org.Hs.eg.db,
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05)
write.csv(as.data.frame(go_up), "GO_BP_upregulated.csv", row.names = FALSE)
dotplot(go_up, showCategory = 15, title = "GO BP - Up in Abnormal")
}
# GO - Downregulated
if (length(dn_entrez) > 0) {
go_dn <- enrichGO(gene = dn_entrez, universe = bg_entrez, OrgDb = org.Hs.eg.db,
ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05)
write.csv(as.data.frame(go_dn), "GO_BP_downregulated.csv", row.names = FALSE)
dotplot(go_dn, showCategory = 15, title = "GO BP - Down in Abnormal")
}
# KEGG - Upregulated
if (length(up_entrez) > 0) {
kegg_up <- enrichKEGG(gene = up_entrez, universe = bg_entrez,
organism = "hsa", pvalueCutoff = 0.05)
write.csv(as.data.frame(kegg_up), "KEGG_upregulated.csv", row.names = FALSE)
dotplot(kegg_up, showCategory = 15, title = "KEGG - Up in Abnormal")
}
# KEGG - Downregulated
if (length(dn_entrez) > 0) {
kegg_dn <- enrichKEGG(gene = dn_entrez, universe = bg_entrez,
organism = "hsa", pvalueCutoff = 0.05)
write.csv(as.data.frame(kegg_dn), "KEGG_downregulated.csv", row.names = FALSE)
dotplot(kegg_dn, showCategory = 15, title = "KEGG - Down in Abnormal")
}
# GSEA
# Use full DESeq2 results (all genes ranked, not just significant)
ranked_genes <- as.data.frame(res_deseq_shrunk) %>%
tibble::rownames_to_column("gene") %>%
filter(!is.na(log2FoldChange)) %>%
arrange(desc(log2FoldChange)) %>%
dplyr::select(gene, log2FoldChange) %>%
deframe()
gsea_go <- gseGO(geneList = ranked_genes,
OrgDb = org.Hs.eg.db,
keyType = "SYMBOL",
ont = "BP",
pAdjustMethod = "BH",
pvalueCutoff = 0.05)
dotplot(gsea_go, showCategory = 15, split = ".sign") +
facet_grid(. ~ .sign) +
theme(
axis.text.y = element_text(size = 7),
axis.text.x = element_text(size = 8),
strip.text = element_text(size = 10),
plot.margin = margin(10, 10, 10, 120)
)
write.csv(as.data.frame(gsea_go), "GSEA_GO_results.csv", row.names = FALSE)
# Use full DESeq2 results, not sig_deseq
all_results <- as.data.frame(res_deseq_shrunk) %>%
tibble::rownames_to_column("gene") %>%
filter(!is.na(log2FoldChange))
ranked_genes <- all_results %>%
arrange(desc(log2FoldChange)) %>%
dplyr::select(gene, log2FoldChange) %>%
deframe()
# ==============================================================================
# 11. SAVE SESSION
# ==============================================================================
saveRDS(seu, "seurat_full.rds")
saveRDS(seu_epi, "seurat_epithelial.rds")
saveRDS(dds, "deseq2_dds.rds")
message("Pipeline completed successfully.")