JacksonFW/transcriptomics-explorer

1# =============================================================================2# edger_analysis.R3# Differential expression analysis using edgeR, with gene symbol annotation4#5# Alternative to deseq2_analysis.R — use either one, not both.6# edgeR tends to be more sensitive for experiments with fewer replicates (< 3).7# DESeq2 is generally preferred when you have 3+ replicates per group.8#9# Input:10#   counts.csv    — raw integer count matrix (genes × samples)11#   metadata.csv  — sample info with condition column12#13# Output:14#   de_results.csv   — edgeR results with gene symbols → upload to dashboard15#   ma_plot.pdf      — MA plot QC16#   volcano_plot.pdf — Volcano plot preview17#18# Usage:19#   Rscript edger_analysis.R20#   or open in RStudio and run interactively21# =============================================================================22 23# ── 0. Install packages if missing ──────────────────────────────────────────24if (!requireNamespace("BiocManager", quietly = TRUE))25  install.packages("BiocManager")26 27for (pkg in c("edgeR", "limma", "org.Hs.eg.db", "AnnotationDbi")) {28  if (!requireNamespace(pkg, quietly = TRUE))29    BiocManager::install(pkg, ask = FALSE)30}31for (pkg in c("ggplot2", "dplyr")) {32  if (!requireNamespace(pkg, quietly = TRUE))33    install.packages(pkg)34}35 36library(edgeR)37library(ggplot2)38library(dplyr)39library(org.Hs.eg.db)40library(AnnotationDbi)41 42# ── 1. USER SETTINGS — edit these ───────────────────────────────────────────43 44COUNTS_FILE   <- "counts.csv"45METADATA_FILE <- "metadata.csv"46 47# Which column in metadata contains the sample IDs48# For airway dataset: use "sample_id"49# For GEO datasets:   use "geo_accession"50SAMPLE_ID_COL <- "sample_id"51 52# Which column in metadata contains the condition/group label53# For airway dataset: use "dex"54# For GEO datasets:   check your metadata.csv column names55CONDITION_COL <- "dex"56 57# For airway dataset:58GROUP_A <- "untrt"   # reference / control59GROUP_B <- "trt"     # treatment / case60 61# Gene ID type in your count matrix row names62# "ENSEMBL"  → ENSG00000... (most common)63# "ENTREZID" → numeric IDs64# "SYMBOL"   → already gene symbols (skip annotation)65GENE_ID_TYPE <- "ENSEMBL"66 67# Annotation database (human by default)68# For mouse: install org.Mm.eg.db and set ANNOTATION_DB <- org.Mm.eg.db69ANNOTATION_DB <- org.Hs.eg.db70 71PADJ_CUTOFF <- 0.0572LFC_CUTOFF  <- 1.073 74OUTPUT_FILE <- "de_results.csv"75 76# ── 2. Load data ─────────────────────────────────────────────────────────────77cat("Loading count matrix from:", COUNTS_FILE, "\n")78counts <- read.csv(COUNTS_FILE, row.names = 1, check.names = FALSE)79counts <- round(as.matrix(counts))80counts[counts < 0] <- 081cat("Count matrix:", nrow(counts), "genes ×", ncol(counts), "samples\n")82 83cat("Loading metadata from:", METADATA_FILE, "\n")84meta <- read.csv(METADATA_FILE)85 86meta <- meta[meta[[CONDITION_COL]] %in% c(GROUP_A, GROUP_B), ]87cat("Samples:", nrow(meta),88    "(", sum(meta[[CONDITION_COL]] == GROUP_A), "×", GROUP_A,89    "/",  sum(meta[[CONDITION_COL]] == GROUP_B), "×", GROUP_B, ")\n")90 91sample_ids <- meta[[SAMPLE_ID_COL]]92counts <- counts[, colnames(counts) %in% sample_ids, drop = FALSE]93counts <- counts[, meta[[SAMPLE_ID_COL]], drop = FALSE]94 95# ── 3. Build edgeR DGEList ───────────────────────────────────────────────────96group <- factor(meta[[CONDITION_COL]], levels = c(GROUP_A, GROUP_B))97dge   <- DGEList(counts = counts, group = group)98 99# Filter lowly expressed genes100keep <- filterByExpr(dge, group = group, min.count = 10, min.total.count = 15)101dge  <- dge[keep, , keep.lib.sizes = FALSE]102cat("Genes after filtering:", nrow(dge), "\n")103 104# ── 4. Normalise and estimate dispersion ─────────────────────────────────────105dge    <- calcNormFactors(dge, method = "TMM")106design <- model.matrix(~ group)107dge    <- estimateDisp(dge, design, robust = TRUE)108cat("Common BCV:", round(sqrt(dge$common.dispersion), 3), "\n")109 110# ── 5. Fit model and test ─────────────────────────────────────────────────────111fit <- glmQLFit(dge, design, robust = TRUE)112qlf <- glmQLFTest(fit, coef = 2)113 114res    <- topTags(qlf, n = Inf, adjust.method = "BH", sort.by = "PValue")115res_df <- as.data.frame(res$table)116res_df$ensembl_id <- rownames(res_df)117rownames(res_df) <- NULL118 119# ── 6. Annotate gene IDs → symbols ───────────────────────────────────────────120if (GENE_ID_TYPE != "SYMBOL") {121  cat("\nMapping", GENE_ID_TYPE, "IDs to gene symbols...\n")122 123  symbols <- mapIds(124    ANNOTATION_DB,125    keys      = res_df$ensembl_id,126    column    = "SYMBOL",127    keytype   = GENE_ID_TYPE,128    multiVals = "first"129  )130 131  res_df$gene <- ifelse(is.na(symbols), res_df$ensembl_id, symbols)132  n_mapped <- sum(!is.na(symbols))133  cat("Mapped", n_mapped, "of", nrow(res_df), "genes to symbols\n")134} else {135  res_df$gene <- res_df$ensembl_id136}137 138# ── 7. Format and save results ────────────────────────────────────────────────139# edgeR doesn't output baseMean — use average CPM as proxy140res_df$baseMean <- rowMeans(cpm(dge, log = FALSE))[res_df$ensembl_id]141 142res_df <- res_df %>%143  dplyr::rename(log2FoldChange = logFC, pvalue = PValue, padj = FDR) %>%144  dplyr::select(gene, ensembl_id, baseMean, log2FoldChange, pvalue, padj) %>%145  dplyr::mutate(146    regulation = dplyr::case_when(147      !is.na(padj) & padj < PADJ_CUTOFF & log2FoldChange >  LFC_CUTOFF ~ "Up",148      !is.na(padj) & padj < PADJ_CUTOFF & log2FoldChange < -LFC_CUTOFF ~ "Down",149      TRUE ~ "NS"150    )151  ) %>%152  dplyr::arrange(padj, desc(abs(log2FoldChange)))153 154n_up   <- sum(res_df$regulation == "Up",   na.rm = TRUE)155n_down <- sum(res_df$regulation == "Down", na.rm = TRUE)156cat("\nSignificant DE genes:", n_up + n_down,157    "(", n_up, "up,", n_down, "down )\n")158 159write.csv(res_df, OUTPUT_FILE, row.names = FALSE)160cat("Saved", OUTPUT_FILE, "\n")161 162# ── 8. QC plots ───────────────────────────────────────────────────────────────163pdf("ma_plot.pdf", width = 7, height = 5)164plotMD(qlf, main = paste("MA plot:", GROUP_B, "vs", GROUP_A))165abline(h = c(-LFC_CUTOFF, LFC_CUTOFF), col = "blue", lty = 2)166dev.off()167cat("Saved ma_plot.pdf\n")168 169res_plot <- res_df %>% dplyr::filter(!is.na(padj))170pdf("volcano_plot.pdf", width = 7, height = 5)171ggplot(res_plot, aes(x = log2FoldChange, y = -log10(padj), colour = regulation)) +172  geom_point(alpha = 0.5, size = 0.8) +173  scale_colour_manual(values = c("Up" = "#e74c3c", "Down" = "#3498db", "NS" = "#aaaaaa")) +174  geom_vline(xintercept = c(-LFC_CUTOFF, LFC_CUTOFF), linetype = "dashed", colour = "black") +175  geom_hline(yintercept = -log10(PADJ_CUTOFF), linetype = "dashed", colour = "black") +176  labs(title = paste("Volcano:", GROUP_B, "vs", GROUP_A),177       x = "log2 Fold Change", y = "-log10(padj)") +178  theme_bw()179dev.off()180cat("Saved volcano_plot.pdf\n")181 182cat("\nAll done!\n")183cat("Upload", OUTPUT_FILE, "to the Transcriptomics Explorer dashboard.\n")184