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"

install.packages("harmony", repos = "https://cloud.r-project.org")

## 
## The downloaded binary packages are in
##  /var/folders/rg/x_7b05fn3sj3v_jq8q367xzm0000gn/T//RtmpSh1Rcb/downloaded_packages

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("org.Hs.eg.db", force = TRUE)

BiocManager::install(
  c(
    "GEOquery", "Seurat", "ggplot2", "dplyr", "pheatmap",
    "ggrepel", "patchwork", "Matrix", "glmGamPoi",
    "DESeq2", "MAST", "fgsea", "org.Hs.eg.db",
    "clusterProfiler", "ReactomePA"
  ),
  update = TRUE,
  ask = FALSE, force = TRUE
)

## 
## The downloaded binary packages are in
##  /var/folders/rg/x_7b05fn3sj3v_jq8q367xzm0000gn/T//RtmpSh1Rcb/downloaded_packages

install.packages("harmony")  # harmony is on CRAN not Bioconductor

## 
## The downloaded binary packages are in
##  /var/folders/rg/x_7b05fn3sj3v_jq8q367xzm0000gn/T//RtmpSh1Rcb/downloaded_packages

suppressPackageStartupMessages({
  library(Seurat)
  library(GEOquery)
  library(DESeq2)
  library(MAST)
  library(tibble)
  library(ggplot2)
  library(ggrepel)
  library(dplyr)
  library(pheatmap)
  library(ReactomePA)
  library(patchwork)
  library(Matrix)
  library(glmGamPoi)
  library(clusterProfiler)
  library(org.Hs.eg.db)
  library(harmony)  # add this
})
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 merger

4.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

# QC distributions with cutoff lines
par(mfrow = c(1,3))
hist(seu$nFeature_RNA, breaks = 100, main = "Genes per Cell", xlab = "nFeature_RNA")
abline(v = c(200, 6000), col = "red", lty = 2)

hist(seu$nCount_RNA, breaks = 100, main = "UMIs per Cell", xlab = "nCount_RNA")

hist(seu$percent.mt, breaks = 100, main = "% Mitochondrial", xlab = "percent.mt")
abline(v = 20, col = "red", lty = 2)

par(mfrow = c(1,1))

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
# ==============================================================================

# Cell cycle scoring (before PCA so scores inform embeddings)
s.genes <- cc.genes$s.genes
g2m.genes <- cc.genes$g2m.genes
seu <- CellCycleScoring(seu, s.features = s.genes, g2m.features = g2m.genes, set.ident = FALSE)
seu$CC.Difference <- seu$S.Score - seu$G2M.Score

# PCA + Harmony (correct for patient, not condition)
seu <- RunPCA(seu, npcs = 50, verbose = FALSE)
N_DIMS <- 30
seu <- RunHarmony(seu, group.by.vars = "patient", verbose = FALSE)

# UMAP, neighbors, clusters
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)

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)

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

epi_genes <- c("EPCAM", "KRT18", "KRT19", "CDH1")
valid_epi <- intersect(epi_genes, rownames(seu))
seu <- AddModuleScore(seu, features = list(valid_epi), name = "Epithelial_Score")

epi_clusters <- AverageExpression(seu,
                                  features = valid_epi,
                                  group.by = "seurat_clusters",
                                  assay = "SCT",
                                  layer = "data")$SCT

# Threshold > 1; strip leading "g" from column names if present
epi_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)


# Heatmap of top DE genes
top_genes <- sig_deseq %>% 
  arrange(padj) %>% 
  head(30) %>% 
  pull(gene)

DoHeatmap(seu_epi, 
          features = top_genes, 
          group.by = "group",
          assay = "SCT",
          size = 3) +
  theme(axis.text.y = element_text(size = 7))

### Looks like we only have 8 genes, but we need all 30 ### 

#library(Seurat)
#library(ggplot2)

# 1. Clean your top 30 genes list to ensure they exist in the SCT assay
top30_genes <- intersect(top_genes, rownames(seu_epi[["SCT"]]))

# 2. Scale all 30 genes explicitly so DoHeatmap can render them
seu_epi <- ScaleData(seu_epi, features = top30_genes, assay = "SCT")

# 3. Generate the single-cell heatmap
DoHeatmap(
  seu_epi,
  features = top30_genes,
  group.by = "group",
  assay = "SCT",
  slot = "scale.data",
  size = 4,              # Gene name text size
  angle = 45             # Group label angle at the top
) + 
  scale_fill_gradientn(colors = c("magenta", "black", "yellow")) # Pink-Black-Yellow palette shown in reference

BONUS: PATHWAY ENRICHMENT AND CROSS METHOD COMPARISON

options(clusterProfiler.download.method = "auto")


# ==============================================================================
# 10. PATHWAY ENRICHMENT (GO, KEGG, REACTOME, GSEA)
# ==============================================================================
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)

# Full DESeq2 results for GSEA — defined FIRST so all methods can use it
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()

ranked_entrez <- all_results %>%
  mutate(entrez = mapIds(org.Hs.eg.db, keys = gene, column = "ENTREZID",
                         keytype = "SYMBOL", multiVals = "first")) %>%
  filter(!is.na(entrez), !is.na(log2FoldChange)) %>%
  arrange(desc(log2FoldChange)) %>%
  dplyr::select(entrez, log2FoldChange) %>%
  deframe()

# 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)
  print(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)
  print(dotplot(go_dn, showCategory = 15, title = "GO BP - Down in Abnormal"))
}

# KEGG - Upregulated
if (length(up_entrez) > 0) {
  tryCatch({
    kegg_up <- enrichKEGG(gene = up_entrez, universe = bg_entrez,
                          organism = "hsa", pvalueCutoff = 0.05,
                          use_internal_data = TRUE)
    write.csv(as.data.frame(kegg_up), "KEGG_upregulated.csv", row.names = FALSE)
    print(dotplot(kegg_up, showCategory = 15, title = "KEGG - Up in Abnormal"))
  }, error = function(e) message("KEGG up unavailable: ", e$message))
}

# KEGG - Downregulated
if (length(dn_entrez) > 0) {
  tryCatch({
    kegg_dn <- enrichKEGG(gene = dn_entrez, universe = bg_entrez,
                          organism = "hsa", pvalueCutoff = 0.05,
                          use_internal_data = TRUE)
    write.csv(as.data.frame(kegg_dn), "KEGG_downregulated.csv", row.names = FALSE)
    print(dotplot(kegg_dn, showCategory = 15, title = "KEGG - Down in Abnormal"))
  }, error = function(e) message("KEGG down unavailable: ", e$message))
}

# REACTOME - Upregulated
if (length(up_entrez) > 0) {
  tryCatch({
    reactome_up <- enrichPathway(gene = up_entrez, universe = bg_entrez,
                                 organism = "human", pAdjustMethod = "BH",
                                 pvalueCutoff = 0.05, readable = TRUE)
    write.csv(as.data.frame(reactome_up), "Reactome_upregulated.csv", row.names = FALSE)
    print(dotplot(reactome_up, showCategory = 15, title = "Reactome - Up in Abnormal"))
  }, error = function(e) message("Reactome up unavailable: ", e$message))
}

# REACTOME - Downregulated
if (length(dn_entrez) > 0) {
  tryCatch({
    reactome_dn <- enrichPathway(gene = dn_entrez, universe = bg_entrez,
                                 organism = "human", pAdjustMethod = "BH",
                                 pvalueCutoff = 0.05, readable = TRUE)
    write.csv(as.data.frame(reactome_dn), "Reactome_downregulated.csv", row.names = FALSE)
    print(dotplot(reactome_dn, showCategory = 15, title = "Reactome - Down in Abnormal"))
  }, error = function(e) message("Reactome down unavailable: ", e$message))
}

# GSEA - GO
gsea_go <- gseGO(geneList      = ranked_genes,
                 OrgDb         = org.Hs.eg.db,
                 keyType       = "SYMBOL",
                 ont           = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff  = 0.05)
print(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)

# GSEA - Reactome
tryCatch({
  gsea_reactome <- gsePathway(geneList = ranked_entrez, organism = "human",
                               pAdjustMethod = "BH", pvalueCutoff = 0.05)
  print(dotplot(gsea_reactome, showCategory = 15, split = ".sign") +
    facet_grid(. ~ .sign) +
    theme(axis.text.y = element_text(size = 7)))
  write.csv(as.data.frame(gsea_reactome), "GSEA_Reactome_results.csv", row.names = FALSE)
}, error = function(e) message("GSEA Reactome unavailable: ", e$message))

# ==============================================================================
# 11. SAVE SESSION
# ==============================================================================
saveRDS(seu,     "seurat_full.rds")
saveRDS(seu_epi, "seurat_epithelial.rds")
saveRDS(dds,     "deseq2_dds.rds")
message("Pipeline completed successfully.")