From ddb6a201ef8ae2e05b68e2ca86fb374daa0826c5 Mon Sep 17 00:00:00 2001 From: zhihankshi Date: Tue, 7 Jul 2026 18:05:07 -0400 Subject: [PATCH] updated the syn2580853 mitochondria folder --- Mitochondria/syn2580853/Proteomics/README.md | 76 ++ .../syn2580853/Proteomics/TauProteomics.R | 838 ++++++++++++ Mitochondria/syn2580853/README.md | 23 +- Mitochondria/syn2580853/RNA_Seq/README.md | 74 ++ Mitochondria/syn2580853/RNA_Seq/TauAnalysis.R | 1164 +++++++++++++++++ 5 files changed, 2174 insertions(+), 1 deletion(-) create mode 100644 Mitochondria/syn2580853/Proteomics/README.md create mode 100644 Mitochondria/syn2580853/Proteomics/TauProteomics.R create mode 100644 Mitochondria/syn2580853/RNA_Seq/README.md create mode 100644 Mitochondria/syn2580853/RNA_Seq/TauAnalysis.R diff --git a/Mitochondria/syn2580853/Proteomics/README.md b/Mitochondria/syn2580853/Proteomics/README.md new file mode 100644 index 0000000..4f1b7b8 --- /dev/null +++ b/Mitochondria/syn2580853/Proteomics/README.md @@ -0,0 +1,76 @@ +# FOLDER PATHS +``` +syn2580853/ +└── Proteomics/ + ├── TMT_all/ + │ └── baylor_proteomics_fly_proteinoutput.txt + ├── RESULTS_syn2580853/ + │ ├── tables/ + │ └── figs/ + ├── TauProteomics.R + ├── EmoryDrosophilaTau_biospecimen_metadata.csv + └── README.md +``` + +# Data Collection + +The following files from Synapse project **syn2580853** (EmoryDrosophilaTau) are used for analysis: + +1. **Baylor proteomics output** — `baylor_proteomics_fly_proteinoutput.txt` (MaxQuant LFQ intensities, 27 samples) +2. **Biospecimen metadata** — `EmoryDrosophilaTau_biospecimen_metadata.csv` (genotype, age, assay) + +Place the Baylor file in a local folder called `TMT_all/`. Place the biospecimen metadata in the `Proteomics/` folder alongside `TauProteomics.R`. + +# Genotypes & Timepoints + +| Genotype | Construct | Samples | +|---|---|---| +| Control | Elav-GAL4/+ | Days 1, 10, 20 (n=3 each) | +| TauWT | Elav-GAL4/+; UAS-TauWT/+ | Days 1, 10, 20 (n=3 each) | +| TauR406W | Elav-GAL4/+; UAS-TauR406W/+ | Days 1, 10, 20 (n=3 each) | + +# QC & Downstream Analysis + +The full proteomics pathway is in [TauProteomics.R](TauProteomics.R). + +The following processes are run on the data: + +1. Quality Control (PCA) +2. limma differential expression (9 contrasts) +3. Visualization (volcano plots, heatmaps) +4. Mitochondrial protein subset +5. enrichR GO ORA +6. GSEA (fgsea, supplementary contrasts) + +## Contrasts + +| Contrast | Description | +|---|---| +| `DE_TauWT_vs_Control_Day*` | Wild-type Tau vs Control | +| `DE_TauR406W_vs_Control_Day*` | R406W mutant vs Control | +| `DE_TauR406W_vs_TauWT_Day*` | Mutant vs wild-type Tau (supplementary) | + +## Key Output Tables + +| File | Description | +|---|---| +| `proteomics_metadata.csv` | Sample genotype and day | +| `DE_TauWT_vs_Control_Day*.csv` | TauWT DE results | +| `DE_TauR406W_vs_Control_Day*.csv` | TauR406W DE results | +| `DE_mitochondrial_proteins_*.csv` | Mitochondrial protein subsets | +| `GO_proteomics_*.csv` | enrichR GO results | + +## How to Run + +1. Download data files into `TMT_all/` and `Proteomics/` +2. Update `setwd()` in `TauProteomics.R` to your local `Proteomics/` path +3. Run the script in RStudio + +```r +setwd("Your_Working_Directory/Proteomics") +source("TauProteomics.R") +``` + +## Proteomics + +The RNA-seq processing for this dataset can be found [here](../RNA_Seq). diff --git a/Mitochondria/syn2580853/Proteomics/TauProteomics.R b/Mitochondria/syn2580853/Proteomics/TauProteomics.R new file mode 100644 index 0000000..b6b96a6 --- /dev/null +++ b/Mitochondria/syn2580853/Proteomics/TauProteomics.R @@ -0,0 +1,838 @@ +# ============================================================ +# Emory Drosophila Tau Model - PROTEOMICS Analysis Pipeline +# Adapted from OSD-514 Proteomics Procedure +# ============================================================ + +# -------------------------------------------------------- +# STEP 1: Load Packages +# -------------------------------------------------------- +library(dplyr) +library(tidyr) +library(readr) +library(stringr) +library(ggplot2) +library(ggrepel) +library(limma) +library(pheatmap) +library(matrixStats) +library(tibble) + +message("✓ Packages loaded") + +# -------------------------------------------------------- +# STEP 2: Set Up Directories +# -------------------------------------------------------- +# Point this to your local syn2580853/Proteomics folder after downloading +# the Baylor TMT file and biospecimen metadata from Synapse. +setwd("Your_Working_Directory/Proteomics") + +OUT_DIR <- "RESULTS_syn2580853" +EXTRAS_DIR <- "extras" +GSEA_DIR <- file.path(EXTRAS_DIR, "GSEA") +dir.create(file.path(OUT_DIR, "figs"), recursive = TRUE, showWarnings = FALSE) +dir.create(file.path(OUT_DIR, "tables"), recursive = TRUE, showWarnings = FALSE) +dir.create(file.path(EXTRAS_DIR, "tables"), recursive = TRUE, showWarnings = FALSE) +dir.create(file.path(GSEA_DIR, "tables"), recursive = TRUE, showWarnings = FALSE) +dir.create(file.path(GSEA_DIR, "figs"), recursive = TRUE, showWarnings = FALSE) + +message("✓ Output directory: ", OUT_DIR) + +# -------------------------------------------------------- +# STEP 3: Read Proteomics Data +# -------------------------------------------------------- +message("Reading proteomics data...") + +prot_raw <- read.table("TMT_all/baylor_proteomics_fly_proteinoutput.txt", + header = TRUE, + sep = "\t", + quote = "", + stringsAsFactors = FALSE) + +message("Raw data: ", nrow(prot_raw), " proteins x ", ncol(prot_raw), " columns") + +# -------------------------------------------------------- +# STEP 4: Clean Protein IDs and Extract LFQ Intensities +# -------------------------------------------------------- + +lfq_cols <- grep("^LFQ", colnames(prot_raw), ignore.case = TRUE, value = TRUE) +message("Found ", length(lfq_cols), " LFQ intensity columns") + +id_col <- if ("Protein.IDs" %in% colnames(prot_raw)) { + "Protein.IDs" +} else { + "Protein IDs" +} + +# Extract expression matrix +expr_raw <- prot_raw[, lfq_cols] +rownames(expr_raw) <- prot_raw[[id_col]] + +# Clean column names to sample numbers (01, 02, ...) +colnames(expr_raw) <- sub("^LFQ(?:[.]|\\s+)?intensity(?:[.]|\\s+)?", "", colnames(expr_raw), ignore.case = TRUE, perl = TRUE) +colnames(expr_raw) <- sprintf("%02d", as.integer(colnames(expr_raw))) + +# Convert to numeric matrix +expr_raw <- as.matrix(expr_raw) +mode(expr_raw) <- "numeric" + +# Replace 0 with NA (MaxQuant uses 0 for missing values) +expr_raw[expr_raw == 0] <- NA + +message("Expression matrix: ", nrow(expr_raw), " proteins x ", ncol(expr_raw), " samples") + +# -------------------------------------------------------- +# STEP 5: Filter Contaminants and Reverse Hits +# -------------------------------------------------------- +message("Filtering contaminants...") + +# Remove contaminants (CON__) and reverse hits (REV__) +is_contaminant <- grepl("^CON__|^REV__", rownames(expr_raw)) +expr_filt <- expr_raw[!is_contaminant, ] + +message("Removed ", sum(is_contaminant), " contaminant/reverse proteins") +message("Remaining: ", nrow(expr_filt), " proteins") + +# -------------------------------------------------------- +# STEP 6: Build Metadata (CORRECTED) +# -------------------------------------------------------- +message("Building metadata...") + +# Read biospecimen metadata +bio_meta <- read.csv("EmoryDrosophilaTau_biospecimen_metadata.csv") + +# Filter to proteomics samples (LC-MSMS) +prot_meta <- bio_meta %>% + filter(assay == "LC-MSMS") %>% + dplyr::select(specimenID, genotype, ageDays) + +# Extract sample number from specimenID ("LFQ intensity 01" -> "01") +prot_meta$sample <- str_extract(prot_meta$specimenID, "\\d+$") + +# Verify the mapping +message("Expression columns: ", paste(colnames(expr_filt), collapse = ", ")) +message("Metadata samples: ", paste(prot_meta$sample, collapse = ", ")) + +# Simplify genotype names +prot_meta <- prot_meta %>% + mutate( + genotype_simple = case_when( + genotype == "Elav-GAL4/+" ~ "Control", + genotype == "Elav-GAL4/+;UAS-TauWT/+" ~ "TauWT", + genotype == "Elav-GAL4/+;UAS-TauR406W/+" ~ "TauR406W" + ), + day = paste0("Day", sprintf("%02d", ageDays)), + group = paste(genotype_simple, day, sep = "_") + ) + +# Set factor levels +prot_meta$genotype_simple <- factor(prot_meta$genotype_simple, + levels = c("Control", "TauWT", "TauR406W")) +prot_meta$day <- factor(prot_meta$day, + levels = c("Day01", "Day10", "Day20")) +prot_meta$group <- factor(prot_meta$group) + +# --- ALIGN: Keep only common samples --- +common_samples <- intersect(colnames(expr_filt), prot_meta$sample) +message("\nCommon samples found: ", length(common_samples)) + +if (length(common_samples) == 0) { + stop("ERROR: No matching samples!") +} + +# Subset and order +expr_filt <- expr_filt[, common_samples] +prot_meta <- prot_meta %>% + filter(sample %in% common_samples) %>% + arrange(match(sample, colnames(expr_filt))) + +rownames(prot_meta) <- prot_meta$sample + +# --- VERIFY --- +stopifnot(all(colnames(expr_filt) == prot_meta$sample)) +message("✓ Alignment verified! ", length(common_samples), " samples matched.") + +message("✓ Metadata created") +print(table(prot_meta$genotype_simple, prot_meta$day)) + +# -------------------------------------------------------- +# STEP 7: Log2 Transform and Filter +# -------------------------------------------------------- +message("Log2 transforming and filtering...") + +# Log2 transform +expr_log <- log2(expr_filt) + +# Filter: keep proteins with values in at least 70% of samples +min_samples <- ceiling(0.7 * ncol(expr_log)) +keep <- rowSums(!is.na(expr_log)) >= min_samples +expr_filt2 <- expr_log[keep, ] + +# Also require some variance +expr_filt2 <- expr_filt2[rowVars(expr_filt2, na.rm = TRUE) > 0, ] + +message("Proteins before filtering: ", nrow(expr_log)) +message("Proteins after filtering: ", nrow(expr_filt2)) + +# -------------------------------------------------------- +# STEP 8: Impute Missing Values (minimal imputation) +# -------------------------------------------------------- +message("Imputing missing values...") + +# For limma, impute missing values with row minimum - 1 (downshifted minimum) +expr_imputed <- expr_filt2 +for (i in 1:nrow(expr_imputed)) { + row_min <- min(expr_imputed[i, ], na.rm = TRUE) + expr_imputed[i, is.na(expr_imputed[i, ])] <- row_min - 1 +} + +# -------------------------------------------------------- +# STEP 9: Quality Control - PCA +# -------------------------------------------------------- +message("Running PCA...") + +pca <- prcomp(t(expr_imputed), center = TRUE, scale = FALSE) +pca_df <- as.data.frame(pca$x) %>% + rownames_to_column("sample") %>% + left_join(prot_meta, by = "sample") + +percentVar <- round(100 * summary(pca)$importance[2, 1:2]) + +png(file.path(OUT_DIR, "figs", "01_PCA_proteomics.png"), + width = 1000, height = 800, res = 150) +ggplot(pca_df, aes(x = PC1, y = PC2, color = genotype_simple, shape = day)) + + geom_point(size = 4) + + scale_color_manual(values = c("Control" = "steelblue", + "TauWT" = "salmon", + "TauR406W" = "darkred")) + + theme_minimal(base_size = 14) + + labs( + x = paste0("PC1: ", percentVar[1], "% variance"), + y = paste0("PC2: ", percentVar[2], "% variance"), + title = "PCA: Proteomics Samples", + color = "Genotype", + shape = "Day" + ) +dev.off() + +message("✓ PCA plot saved") + +# -------------------------------------------------------- +# STEP 10: Differential Analysis with limma +# -------------------------------------------------------- +message("Running limma differential analysis...") + +# Create design matrix +design <- model.matrix(~ 0 + group, data = prot_meta) +colnames(design) <- gsub("group", "", colnames(design)) + +# Fit linear model +fit <- lmFit(expr_imputed, design) + +# Define contrasts +contrast_matrix <- makeContrasts( + # TauWT vs Control (overall across all days) + TauWT_vs_Control_Day01 = TauWT_Day01 - Control_Day01, + TauWT_vs_Control_Day10 = TauWT_Day10 - Control_Day10, + TauWT_vs_Control_Day20 = TauWT_Day20 - Control_Day20, + + # TauR406W vs Control + TauR406W_vs_Control_Day01 = TauR406W_Day01 - Control_Day01, + TauR406W_vs_Control_Day10 = TauR406W_Day10 - Control_Day10, + TauR406W_vs_Control_Day20 = TauR406W_Day20 - Control_Day20, + + # TauR406W vs TauWT (mutant vs wild-type tau) + TauR406W_vs_TauWT_Day01 = TauR406W_Day01 - TauWT_Day01, + TauR406W_vs_TauWT_Day10 = TauR406W_Day10 - TauWT_Day10, + TauR406W_vs_TauWT_Day20 = TauR406W_Day20 - TauWT_Day20, + + levels = design +) + +# Fit contrasts +fit2 <- contrasts.fit(fit, contrast_matrix) +fit2 <- eBayes(fit2) + +message("✓ Limma analysis complete") + +# -------------------------------------------------------- +# STEP 11: Extract and Save Results +# -------------------------------------------------------- +message("Extracting results...") + +# Function to extract and save results +save_limma_results <- function(fit, coef_name) { + tt <- topTable(fit, coef = coef_name, number = Inf, adjust.method = "BH") %>% + rownames_to_column("protein_id") %>% + arrange(adj.P.Val) + + table_dir <- if (grepl("^TauR406W_vs_TauWT_", coef_name)) { + file.path(EXTRAS_DIR, "tables") + } else { + file.path(OUT_DIR, "tables") + } + + write.csv(tt, + file.path(table_dir, paste0("DE_", coef_name, ".csv")), + row.names = FALSE) + + n_sig <- sum(tt$adj.P.Val < 0.05, na.rm = TRUE) + n_up <- sum(tt$adj.P.Val < 0.05 & tt$logFC > 0.5, na.rm = TRUE) + n_down <- sum(tt$adj.P.Val < 0.05 & tt$logFC < -0.5, na.rm = TRUE) + + message(sprintf(" %s: %d sig (↑%d, ↓%d)", coef_name, n_sig, n_up, n_down)) + + return(tt) +} + +# Extract all contrasts +contrasts_list <- colnames(contrast_matrix) +results_list <- list() + +for (contrast in contrasts_list) { + results_list[[contrast]] <- save_limma_results(fit2, contrast) +} + +# -------------------------------------------------------- +# STEP 12: Volcano Plots +# -------------------------------------------------------- +message("Creating volcano plots...") + +save_volcano <- function(tt, contrast_name, top_n = 15) { + tt <- tt %>% + mutate( + sig = case_when( + adj.P.Val < 0.05 & logFC > 0.5 ~ "Up", + adj.P.Val < 0.05 & logFC < -0.5 ~ "Down", + TRUE ~ "NS" + ) + ) + + # Top genes to label + top_up <- tt %>% filter(sig == "Up") %>% arrange(adj.P.Val) %>% head(top_n/2) + top_down <- tt %>% filter(sig == "Down") %>% arrange(adj.P.Val) %>% head(top_n/2) + top_labels <- bind_rows(top_up, top_down) + + p <- ggplot(tt, aes(x = logFC, y = -log10(P.Value), color = sig)) + + geom_point(alpha = 0.6, size = 1.5) + + geom_text_repel(data = top_labels, + aes(label = protein_id), + size = 2.5, + max.overlaps = 20) + + scale_color_manual(values = c("Down" = "blue", "Up" = "red", "NS" = "grey70")) + + geom_vline(xintercept = c(-0.5, 0.5), linetype = "dashed", color = "gray40") + + geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray40") + + theme_minimal(base_size = 12) + + labs( + x = "log2 Fold Change", + y = "-log10(p-value)", + title = paste("Volcano:", gsub("_", " ", contrast_name)) + ) + + ggsave(file.path(OUT_DIR, "figs", paste0("volcano_", contrast_name, ".png")), + p, width = 10, height = 8, dpi = 200) +} + +# Create volcano for key comparisons +save_volcano(results_list[["TauWT_vs_Control_Day01"]], "TauWT_vs_Control_Day01") +save_volcano(results_list[["TauWT_vs_Control_Day10"]], "TauWT_vs_Control_Day10") +save_volcano(results_list[["TauWT_vs_Control_Day20"]], "TauWT_vs_Control_Day20") +save_volcano(results_list[["TauR406W_vs_Control_Day20"]], "TauR406W_vs_Control_Day20") + +message("✓ Volcano plots saved") + +# -------------------------------------------------------- +# STEP 13: Heatmaps for All Key Contrasts +# -------------------------------------------------------- +message("Creating heatmaps for all contrasts...") + +create_heatmap <- function(contrast_name, results_df, expr_mat, prot_meta, + top_n = 50, padj_cut = 0.1) { + + # Get top DE proteins + top_proteins <- results_df %>% + dplyr::filter(adj.P.Val < padj_cut) %>% + arrange(adj.P.Val) %>% + head(top_n) %>% + pull(protein_id) + + if (length(top_proteins) < 5) { + # If not enough significant, take top by p-value + top_proteins <- results_df %>% + arrange(P.Value) %>% + head(top_n) %>% + pull(protein_id) + message(" ", contrast_name, ": Using top ", length(top_proteins), " by p-value (few significant)") + } else { + message(" ", contrast_name, ": ", length(top_proteins), " significant proteins") + } + + # Subset expression matrix + heat_mat <- expr_mat[top_proteins, ] + + # Z-score normalize + heat_mat_z <- t(scale(t(heat_mat))) + heat_mat_z[heat_mat_z > 2.5] <- 2.5 + heat_mat_z[heat_mat_z < -2.5] <- -2.5 + + # Create simplified column labels (Genotype_Day_Rep) + # Create replicate numbers within each group + prot_meta_ordered <- prot_meta %>% + group_by(genotype_simple, day) %>% + mutate(rep = row_number()) %>% + ungroup() + + col_labels <- paste0( + substr(prot_meta_ordered$genotype_simple, 1, 4), # Ctrl, TauW, TauR + "_D", prot_meta_ordered$ageDays, + "_", prot_meta_ordered$rep + ) + + # Annotation + ann_col <- data.frame( + Genotype = prot_meta$genotype_simple, + Day = prot_meta$day, + row.names = prot_meta$sample + ) + + ann_colors <- list( + Genotype = c(Control = "steelblue", TauWT = "salmon", TauR406W = "darkred"), + Day = c(Day01 = "#fee8c8", Day10 = "#fdbb84", Day20 = "#e34a33") + ) + + # Order columns by genotype then day + col_order <- order(prot_meta$genotype_simple, prot_meta$day) + heat_mat_z <- heat_mat_z[, col_order] + ann_col <- ann_col[col_order, ] + col_labels <- col_labels[col_order] + + # Simplify row labels (protein names) + row_labels <- sapply(rownames(heat_mat_z), function(x) { + # Extract gene name from protein ID (e.g., sp|P02515|HSP22_DROME -> HSP22) + parts <- strsplit(x, "\\|")[[1]] + if (length(parts) >= 3) { + gene <- gsub("_DROME.*", "", parts[3]) + return(gene) + } + # If format different, truncate + if (nchar(x) > 20) return(paste0(substr(x, 1, 17), "...")) + return(x) + }) + + # Save heatmap + png(file.path(OUT_DIR, "figs", paste0("heatmap_", contrast_name, ".png")), + width = 1400, height = 1200, res = 150) + + pheatmap(heat_mat_z, + annotation_col = ann_col, + annotation_colors = ann_colors, + labels_col = col_labels, + labels_row = row_labels, + cluster_cols = FALSE, + cluster_rows = TRUE, + show_rownames = TRUE, + show_colnames = TRUE, + fontsize_row = 7, + fontsize_col = 8, + border_color = NA, + main = paste0("Top DE Proteins: ", gsub("_", " ", contrast_name))) + + dev.off() + message(" ✓ Saved: heatmap_", contrast_name, ".png") +} + +# Create heatmaps for all key contrasts +key_contrasts <- c( + "TauWT_vs_Control_Day01", + "TauWT_vs_Control_Day10", + "TauWT_vs_Control_Day20", + "TauR406W_vs_Control_Day01", + "TauR406W_vs_Control_Day10", + "TauR406W_vs_Control_Day20", + "TauR406W_vs_TauWT_Day20" +) + +for (contrast in key_contrasts) { + if (contrast %in% names(results_list)) { + create_heatmap( + contrast_name = contrast, + results_df = results_list[[contrast]], + expr_mat = expr_imputed, + prot_meta = prot_meta + ) + } +} + +message("✓ All heatmaps created!") + +# -------------------------------------------------------- +# STEP 14: Summary +# -------------------------------------------------------- +message("\n========================================") +message("PROTEOMICS ANALYSIS COMPLETE!") +message("========================================") +message("\nOutput in: ", OUT_DIR) +message("\nDE Summary:") +for (contrast in contrasts_list) { + n_sig <- sum(results_list[[contrast]]$adj.P.Val < 0.05, na.rm = TRUE) + message(sprintf(" %s: %d significant proteins", contrast, n_sig)) +} +message("\nFiles generated:") +message(" - tables/proteomics_metadata.csv") +message(" - tables/DE_*.csv (all contrasts)") +message(" - figs/01_PCA_proteomics.png") +message(" - figs/volcano_*.png") +message(" - figs/heatmap_top50_DE_proteins.png") + +# ============================================================ +# STEP 15: GSEA Analysis for Proteomics +# ============================================================ + +message("\n========================================") +message("Starting GSEA Analysis...") +message("========================================") + +# --- Load additional packages --- +if (!requireNamespace("fgsea", quietly = TRUE)) BiocManager::install("fgsea") +if (!requireNamespace("org.Dm.eg.db", quietly = TRUE)) BiocManager::install("org.Dm.eg.db") +if (!requireNamespace("GO.db", quietly = TRUE)) BiocManager::install("GO.db") + +library(fgsea) +library(org.Dm.eg.db) +library(AnnotationDbi) +library(GO.db) +library(data.table) + +# -------------------------------------------------------- +# -------------------------------------------------------- +# STEP 16: Map Protein IDs to Gene Symbols (FIXED) +# -------------------------------------------------------- +message("Mapping protein IDs to FlyBase gene symbols...") + +# Extract UniProt accession from protein ID +# Format: sp|P02515|HSP22_DROME -> P02515 +# Format: tr|Q9VZR2|Q9VZR2_DROME -> Q9VZR2 + +extract_uniprot_id <- function(protein_id) { + parts <- strsplit(protein_id, "\\|")[[1]] + if (length(parts) >= 2) { + return(parts[2]) + } + return(NA) +} + +# Create mapping with UniProt IDs +protein_ids <- rownames(expr_imputed) + +uniprot_map <- data.frame( + protein_id = protein_ids, + uniprot_id = sapply(protein_ids, extract_uniprot_id), + stringsAsFactors = FALSE +) + +message(" Extracted ", sum(!is.na(uniprot_map$uniprot_id)), " UniProt accessions") + +# Get unique UniProt IDs +uniprot_ids <- unique(na.omit(uniprot_map$uniprot_id)) + +# Map to SYMBOL using org.Dm.eg.db +uniprot_to_symbol <- AnnotationDbi::select( + org.Dm.eg.db, + keys = uniprot_ids, + keytype = "UNIPROT", + columns = c("UNIPROT", "SYMBOL") +) + +message(" Successfully mapped: ", sum(!is.na(uniprot_to_symbol$SYMBOL)), " proteins to FlyBase symbols") + +# Merge back with protein IDs +gene_map <- uniprot_map %>% + left_join(uniprot_to_symbol, by = c("uniprot_id" = "UNIPROT")) %>% + dplyr::rename(gene_symbol = SYMBOL) %>% + dplyr::filter(!is.na(gene_symbol) & gene_symbol != "") %>% + distinct(protein_id, .keep_all = TRUE) + +message(" Final gene map: ", nrow(gene_map), " proteins with FlyBase symbols") + +# Check overlap with GO pathways +genes_in_pathways <- unique(unlist(pathways)) +overlap <- sum(gene_map$gene_symbol %in% genes_in_pathways) +message(" Genes overlapping with GO pathways: ", overlap) + +# -------------------------------------------------------- +# STEP 17: Build Ranked Gene Lists for Each Contrast +# -------------------------------------------------------- +message("Building ranked gene lists...") + +build_ranks <- function(de_results, gene_map) { + # Check if protein_id is already a column + if ("protein_id" %in% colnames(de_results)) { + de_with_genes <- de_results %>% + left_join(gene_map, by = "protein_id") + } else { + de_with_genes <- de_results %>% + rownames_to_column("protein_id") %>% + left_join(gene_map, by = "protein_id") + } + + # Filter valid entries + de_with_genes <- de_with_genes %>% + dplyr::filter(!is.na(gene_symbol), gene_symbol != "") %>% + dplyr::filter(!is.na(logFC), !is.na(P.Value)) + + # For duplicate genes, keep the one with lowest p-value + de_unique <- de_with_genes %>% + group_by(gene_symbol) %>% + slice_min(P.Value, n = 1) %>% + ungroup() + + # Create rank score: sign(logFC) * -log10(pvalue) + de_unique <- de_unique %>% + mutate(rank_score = sign(logFC) * -log10(P.Value)) + + # Named vector for fgsea + ranks <- setNames(de_unique$rank_score, de_unique$gene_symbol) + ranks <- sort(ranks, decreasing = TRUE) + ranks <- ranks[is.finite(ranks)] + + return(ranks) +} + +# -------------------------------------------------------- +# STEP 18: Build GO Gene Sets +# -------------------------------------------------------- +message("Building GO Biological Process gene sets...") + +# Get all GO BP terms for Drosophila +go_bp <- AnnotationDbi::select( + org.Dm.eg.db, + keys = keys(org.Dm.eg.db, keytype = "SYMBOL"), + columns = c("SYMBOL", "GOALL", "ONTOLOGYALL"), + keytype = "SYMBOL" +) %>% + filter(ONTOLOGYALL == "BP") %>% + dplyr::select(SYMBOL, GOALL) %>% + distinct() + +# Convert to list format for fgsea +pathways <- split(go_bp$SYMBOL, go_bp$GOALL) + +# Filter by size (10-500 genes per pathway) +pathways <- pathways[sapply(pathways, length) >= 10 & sapply(pathways, length) <= 500] +message("GO BP pathways loaded: ", length(pathways)) + +# Get GO term names +go_terms <- AnnotationDbi::select(GO.db, keys = names(pathways), + columns = "TERM", keytype = "GOID") +term_map <- setNames(go_terms$TERM, go_terms$GOID) + +# -------------------------------------------------------- +# STEP 19: Run fGSEA for Each Contrast +# -------------------------------------------------------- +run_gsea_for_contrast <- function(contrast_name, de_results, gene_map, pathways, term_map) { + + message("\n--- Running GSEA for: ", contrast_name, " ---") + + # Build ranks using the updated function + ranks <- build_ranks(de_results, gene_map) + message(" Genes in ranked list: ", length(ranks)) + + if (length(ranks) < 50) { + message(" WARNING: Too few genes for GSEA") + return(NULL) + } + + # Run fgsea + set.seed(42) + fgsea_res <- fgsea( + pathways = pathways, + stats = ranks, + minSize = 10, + maxSize = 500, + nPermSimple = 10000 + ) + + # Add GO term names + fgsea_res <- fgsea_res %>% + mutate(term = term_map[pathway]) %>% + arrange(padj) + + n_sig <- sum(fgsea_res$padj < 0.05, na.rm = TRUE) + message(" Significant pathways (padj < 0.05): ", n_sig) + + # Save results table + fgsea_out <- fgsea_res %>% + mutate(leadingEdge = sapply(leadingEdge, paste, collapse = ";")) %>% + dplyr::select(pathway, term, size, NES, pval, padj, leadingEdge) + + write.csv(fgsea_out, + file.path(GSEA_DIR, "tables", paste0("fgsea_GO_BP_", contrast_name, ".csv")), + row.names = FALSE) + + # --- Barplot of top pathways --- + top_paths <- fgsea_res %>% + dplyr::filter(is.finite(padj)) %>% + head(20) %>% + mutate(label = ifelse(is.na(term) | term == "", pathway, term)) + + if (nrow(top_paths) > 0) { + p <- ggplot(top_paths, aes(x = reorder(label, NES), y = NES, fill = -log10(padj))) + + geom_col() + + coord_flip() + + scale_fill_gradient(low = "steelblue", high = "darkred") + + labs( + title = paste0("GSEA GO Biological Process\n", gsub("_", " ", contrast_name)), + x = NULL, + y = "Normalized Enrichment Score (NES)", + fill = "-log10(FDR)" + ) + + theme_minimal(base_size = 11) + + theme(plot.title = element_text(hjust = 0.5)) + + ggsave(file.path(GSEA_DIR, "figs", paste0("fgsea_GO_BP_", contrast_name, ".png")), + p, width = 10, height = 8, dpi = 300) + } + + # --- Mito-specific pathways --- + mito_pat <- paste( + "mitochond", "oxidative phosph", "electron transport", + "respiratory chain", "ATP synthase", "tricarboxylic", + "TCA", "beta-oxid", "mitophagy", "fission", "fusion", + sep = "|" + ) + + mito_paths <- fgsea_res %>% + dplyr::filter(grepl(mito_pat, term, ignore.case = TRUE)) %>% + arrange(padj) + + if (nrow(mito_paths) > 0) { + write.csv( + mito_paths %>% mutate(leadingEdge = sapply(leadingEdge, paste, collapse = ";")), + file.path(GSEA_DIR, "tables", paste0("fgsea_mito_", contrast_name, ".csv")), + row.names = FALSE + ) + + message(" Mito-related pathways: ", nrow(mito_paths)) + + # Mito barplot + mito_plot <- mito_paths %>% + head(15) %>% + mutate(label = ifelse(is.na(term) | term == "", pathway, term)) + + if (nrow(mito_plot) > 0) { + p_mito <- ggplot(mito_plot, aes(x = reorder(label, NES), y = NES, fill = -log10(padj))) + + geom_col() + + coord_flip() + + scale_fill_gradient(low = "steelblue", high = "darkred") + + labs( + title = paste0("Mitochondrial Pathways\n", gsub("_", " ", contrast_name)), + x = NULL, + y = "NES", + fill = "-log10(FDR)" + ) + + theme_minimal(base_size = 11) + + theme(plot.title = element_text(hjust = 0.5)) + + ggsave(file.path(GSEA_DIR, "figs", paste0("fgsea_mito_", contrast_name, ".png")), + p_mito, width = 10, height = 6, dpi = 300) + } + } else { + message(" No mito-related pathways found") + } + + return(fgsea_res) +} + +# -------------------------------------------------------- +# STEP 20: Run GSEA for All Contrasts +# -------------------------------------------------------- +message("\nRunning GSEA for all contrasts...") + +gsea_results <- list() + +for (contrast in names(results_list)) { + gsea_results[[contrast]] <- run_gsea_for_contrast( + contrast_name = contrast, + de_results = results_list[[contrast]], + gene_map = gene_map, + pathways = pathways, + term_map = term_map + ) +} + +# -------------------------------------------------------- +# STEP 21: Summary +# -------------------------------------------------------- +message("\n========================================") +message("GSEA ANALYSIS COMPLETE!") +message("========================================") +message("\nOutput in: ", normalizePath(GSEA_DIR)) +message("\nGSEA Summary:") +for (contrast in names(gsea_results)) { + if (!is.null(gsea_results[[contrast]])) { + n_sig <- sum(gsea_results[[contrast]]$padj < 0.05, na.rm = TRUE) + n_mito <- sum(grepl("mitochond|oxidative|electron|respiratory|ATP|TCA", + gsea_results[[contrast]]$term, ignore.case = TRUE) & + gsea_results[[contrast]]$padj < 0.05, na.rm = TRUE) + message(sprintf(" %s: %d significant pathways (%d mito-related)", + contrast, n_sig, n_mito)) + } +} +message("\nFiles generated:") +message(" - tables/fgsea_GO_BP_*.csv (all pathways)") +message(" - tables/fgsea_mito_*.csv (mito pathways)") +message(" - figs/fgsea_GO_BP_*.png (top pathway barplots)") +message(" - figs/fgsea_mito_*.png (mito pathway barplots)") + +# --- Create GSEA trend plots --- +message("\nCreating GSEA trend plots...") + +for (contrast in names(gsea_results)) { + if (!is.null(gsea_results[[contrast]]) && nrow(gsea_results[[contrast]]) > 0) { + + # Convert to data.frame explicitly + fg <- as.data.frame(gsea_results[[contrast]]) + + # Get top 10 positive and top 10 negative NES + top_up <- fg[fg$NES > 0, ] + top_up <- top_up[order(top_up$pval), ] + top_up <- head(top_up, 10) + + top_down <- fg[fg$NES < 0, ] + top_down <- top_down[order(top_down$pval), ] + top_down <- head(top_down, 10) + + top_paths <- rbind(top_up, top_down) + + if (nrow(top_paths) > 0) { + # Create short term labels + top_paths$term_short <- ifelse( + is.na(top_paths$term) | top_paths$term == "", + top_paths$pathway, + ifelse(nchar(top_paths$term) > 45, + paste0(substr(top_paths$term, 1, 42), "..."), + top_paths$term) + ) + + p <- ggplot(top_paths, + aes(x = reorder(term_short, NES), y = NES, fill = -log10(padj))) + + geom_col() + + coord_flip() + + scale_fill_gradient(low = "steelblue", high = "darkred", + name = "-log10(FDR)") + + theme_minimal(base_size = 10) + + labs( + x = "", + y = "Normalized Enrichment Score (NES)", + title = paste0("GSEA: ", gsub("_", " ", contrast)), + subtitle = "Top pathways | Positive NES = Upregulated" + ) + + geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") + + theme(axis.text.y = element_text(size = 8)) + + out_file <- file.path(GSEA_DIR, "figs", paste0("fgsea_trends_", contrast, ".png")) + ggsave(out_file, p, width = 11, height = 8, dpi = 200) + } + } +} + +message("✓ GSEA trend plots created") \ No newline at end of file diff --git a/Mitochondria/syn2580853/README.md b/Mitochondria/syn2580853/README.md index cde982c..74f22ec 100644 --- a/Mitochondria/syn2580853/README.md +++ b/Mitochondria/syn2580853/README.md @@ -1 +1,22 @@ -Coming Soon! +# syn2580853 — Emory Drosophila Tau Model + +Drosophila melanogaster model of Tau-driven neurodegeneration (EmoryTau / AMP-AD). This directory holds analysis code and documentation for the Synapse project **syn2580853**. + +## Study Overview + +Human Tau (wild-type and R406W FTLD mutant) is expressed in neurons via Elav-GAL4. RNA-seq and LC-MS/MS proteomics were collected at Day 1, Day 10, and Day 20. + + +| Genotype | Construct | +| -------- | --------------------------- | +| Control | Elav-GAL4/+ | +| TauWT | Elav-GAL4/+; UAS-TauWT/+ | +| TauR406W | Elav-GAL4/+; UAS-TauR406W/+ | + + +## Subdirectories + +- [RNA_Seq](RNA_Seq) — transcriptomics pipeline (`TauAnalysis.R`) +- [Proteomics](Proteomics) — proteomics pipeline (`TauProteomics.R`) + +Raw data and analysis outputs are stored on Synapse and generated locally when the scripts are run; they are not committed to this repository. \ No newline at end of file diff --git a/Mitochondria/syn2580853/RNA_Seq/README.md b/Mitochondria/syn2580853/RNA_Seq/README.md new file mode 100644 index 0000000..e6c59c3 --- /dev/null +++ b/Mitochondria/syn2580853/RNA_Seq/README.md @@ -0,0 +1,74 @@ +# FOLDER PATHS +``` +syn2580853/ +└── RNA_Seq/ + ├── Collapsed_Counts/ + │ ├── Trim_*_counts.txt + │ └── metadata_from_filenames.csv + ├── RESULTS_syn2580853/ + │ ├── tables/ + │ └── figs/ + ├── GSEA/ + │ ├── tables/ + │ └── figs/ + ├── TauAnalysis.R + └── README.md +``` + +# Data Collection + +The following files from Synapse project **syn2580853** (EmoryDrosophilaTau) are used for analysis: + +1. **18 RNA-seq count files** — `Trim_*_counts.txt` (e.g. `Trim_DE01D01_counts.txt`) +2. **RNA-seq metadata** — `EmoryDrosophilaTau_RNAseq_metadata.csv` (optional reference) + +Download the 18 `Trim_*_counts.txt` files into a local folder called `Collapsed_Counts/` before running the pipeline. + +# Sample Naming + +| Code | Meaning | +|---|---| +| `DE` | Control (Elav-GAL4 driver only) | +| `DT` | TauWT | +| `01` / `10` / `20` | Day 1, 10, or 20 | +| `D01`–`D03` | Biological replicate | + +Example: `Trim_DT20D02_counts.txt` = TauWT, Day 20, replicate 2. + +`TauAnalysis.R` parses sample names from filenames and writes `metadata_from_filenames.csv` into `Collapsed_Counts/`. + +# QC & Downstream Analysis + +The full RNA-seq pathway is in [TauAnalysis.R](TauAnalysis.R). + +The following processes are run on the data: + +1. Quality Control (library sizes, PCA) +2. DESeq2 (`~ day + genotype`, Tau vs Control) +3. Visualization (volcano plots, heatmaps) +4. Timepoint-specific contrasts (Day01, Day10, Day20) +5. Mitochondrial gene subset +6. enrichR GO ORA +7. GSEA (fgsea, GO Biological Process) + +## Key Output Tables + +| File | Description | +|---|---| +| `DE_Tau_vs_Control.csv` | Overall Tau vs Control | +| `DE_Tau_vs_Control_shrunk.csv` | apeglm-shrunk LFC | +| `DE_Tau_vs_Control_Day*.csv` | Per-timepoint contrasts | +| `DE_mitochondrial_genes.csv` | Curated mitochondrial gene subset | +| `sample_metadata.csv` | Sample annotations | +| `fgsea_GO_BP_Tau_vs_Control.csv` | GSEA pathway results | + +## How to Run + +1. Download count files into `Collapsed_Counts/` +2. Update `setwd()` in `TauAnalysis.R` to your local `RNA_Seq/` path +3. Run the script in RStudio + +```r +setwd("Your_Working_Directory/RNA_Seq") +source("TauAnalysis.R") +``` diff --git a/Mitochondria/syn2580853/RNA_Seq/TauAnalysis.R b/Mitochondria/syn2580853/RNA_Seq/TauAnalysis.R new file mode 100644 index 0000000..fa2e065 --- /dev/null +++ b/Mitochondria/syn2580853/RNA_Seq/TauAnalysis.R @@ -0,0 +1,1164 @@ +# ============================================================ +# Emory Drosophila Tau Model - RNA-seq Analysis Pipeline +# ============================================================ + +# STEP 1: Load packages +# -------------------------------------------------------- +# These are like "apps" that give R extra abilities + +library(DESeq2) # For differential expression +library(ggplot2) # For plotting +library(pheatmap) # For heatmaps +library(dplyr) # For data manipulation +library(tibble) # For data frames +library(stringr) # For text manipulation +library(dplyr) +select <- dplyr::select # Force dplyr's select to be default + +message("✓ Packages loaded successfully!") + +# STEP 2: Set working directory +# -------------------------------------------------------- +# Point this to your local syn2580853/RNA_Seq folder after downloading +# count files from Synapse into Collapsed_Counts/. +setwd("Your_Working_Directory/RNA_Seq") + +message("Working directory: ", getwd()) + +# STEP 3: Create output folders +# -------------------------------------------------------- +OUT_DIR <- "RESULTS_syn2580853" +GSEA_DIR <- "GSEA" +EXTRAS_DIR <- "extras" +dir.create(file.path(OUT_DIR, "figs"), recursive = TRUE, showWarnings = FALSE) +dir.create(file.path(OUT_DIR, "tables"), recursive = TRUE, showWarnings = FALSE) +dir.create(file.path(GSEA_DIR, "figs"), recursive = TRUE, showWarnings = FALSE) +dir.create(file.path(GSEA_DIR, "tables"), recursive = TRUE, showWarnings = FALSE) +dir.create(file.path(EXTRAS_DIR, "tables"), recursive = TRUE, showWarnings = FALSE) +message("✓ Output folders created") + +# STEP 4: Find your counts files +# -------------------------------------------------------- +count_files <- list.files("Collapsed_Counts", pattern = "_counts\\.txt$", full.names = TRUE) +message("Found ", length(count_files), " count files") + +# Check: Did it find 18 files? +if (length(count_files) == 0) { + stop("ERROR: No count files found! Check that Collapsed_Counts/ contains the _counts.txt files") +} + +# Show the files found +print(basename(count_files)) + +# STEP 5: Read and merge all counts into one table +# -------------------------------------------------------- +message("Reading count files...") + +# Function to read one counts file +read_counts <- function(file) { + df <- read.table(file, header = FALSE, sep = "\t", + col.names = c("gene", "count"), + stringsAsFactors = FALSE) + return(df) +} + +# Get sample names from filenames (e.g., "DE01D01" from "Trim_DE01D01_counts.txt") +sample_names <- basename(count_files) %>% + str_remove("^Trim_") %>% + str_remove("_counts\\.txt$") + +message("Sample names: ", paste(sample_names, collapse = ", ")) + +# Read all files into a list +counts_list <- lapply(count_files, read_counts) +names(counts_list) <- sample_names + +# Merge all into one data frame +counts_merged <- counts_list[[1]] %>% rename(!!sample_names[1] := count) + +for (i in 2:length(counts_list)) { + temp <- counts_list[[i]] %>% rename(!!sample_names[i] := count) + counts_merged <- full_join(counts_merged, temp, by = "gene") +} + +# Convert to matrix format +rownames(counts_merged) <- counts_merged$gene +counts_mat <- as.matrix(counts_merged[, -1]) +mode(counts_mat) <- "integer" +counts_mat[is.na(counts_mat)] <- 0 + +message("✓ Count matrix created: ", nrow(counts_mat), " genes x ", ncol(counts_mat), " samples") + +# STEP 6: Create metadata from sample names +# -------------------------------------------------------- +# Sample naming: DE01D01 +# DE = ELAV Control (E), DT = Tau (T) +# 01/10/20 = Day 1, 10, 20 +# D01/D02/D03 = replicate number + +meta <- data.frame( + sample = sample_names, + genotype = ifelse(str_detect(sample_names, "^DE"), "Control", "Tau_WT"), + day = case_when( + str_detect(sample_names, "^..01") ~ "Day01", + str_detect(sample_names, "^..10") ~ "Day10", + str_detect(sample_names, "^..20") ~ "Day20" + ), + replicate = str_extract(sample_names, "D0[0-9]$"), + row.names = sample_names +) + +# Set factor levels (Control is reference) +meta$genotype <- factor(meta$genotype, levels = c("Control", "Tau_WT")) +meta$day <- factor(meta$day, levels = c("Day01", "Day10", "Day20")) + +# Display metadata +message("✓ Metadata created:") +print(meta) + +# Save metadata +write.csv(meta, file.path(OUT_DIR, "tables", "sample_metadata.csv"), row.names = FALSE) +write.csv(meta, file.path("Collapsed_Counts", "metadata_from_filenames.csv"), row.names = FALSE) + +# STEP 7: Quality Control - Library Sizes +# -------------------------------------------------------- +message("Running QC...") + +lib_sizes <- colSums(counts_mat) + +# Create bar plot +png(file.path(OUT_DIR, "figs", "01_library_sizes.png"), width = 1000, height = 600, res = 150) +par(mar = c(8, 4, 3, 1)) +barplot(lib_sizes / 1e6, + las = 2, + ylab = "Library Size (millions)", + main = "Library Sizes per Sample", + col = ifelse(meta$genotype == "Control", "steelblue", "salmon"), + cex.names = 0.8) +legend("topright", legend = c("Control", "Tau_WT"), fill = c("steelblue", "salmon")) +dev.off() + +message("✓ Library size plot saved") + +# STEP 8: Filter low-count genes +# -------------------------------------------------------- +# Keep genes with at least 10 counts in at least 20% of samples +min_samples <- ceiling(ncol(counts_mat) * 0.2) +keep <- rowSums(counts_mat >= 10) >= min_samples +counts_filt <- counts_mat[keep, ] + +message("Genes before filtering: ", nrow(counts_mat)) +message("Genes after filtering: ", nrow(counts_filt)) +message("Removed ", nrow(counts_mat) - nrow(counts_filt), " low-count genes") + +# STEP 9: Create DESeq2 object +# -------------------------------------------------------- +message("Creating DESeq2 dataset...") + +dds <- DESeqDataSetFromMatrix( + countData = counts_filt, + colData = meta, + design = ~ day + genotype # Test genotype effect, accounting for day +) + +# STEP 10: Run differential expression analysis +# -------------------------------------------------------- +message("Running DESeq2 (this may take a minute)...") + +dds <- DESeq(dds) + +message("✓ DESeq2 complete!") + +# STEP 11: Transform data for visualization +# -------------------------------------------------------- +vsd <- vst(dds, blind = TRUE) + +# STEP 12: PCA Plot +# -------------------------------------------------------- +pca_data <- plotPCA(vsd, intgroup = c("genotype", "day"), returnData = TRUE) +percentVar <- round(100 * attr(pca_data, "percentVar")) + +png(file.path(OUT_DIR, "figs", "02_PCA_plot.png"), width = 1000, height = 800, res = 150) +ggplot(pca_data, aes(x = PC1, y = PC2, color = genotype, shape = day)) + + geom_point(size = 5) + + xlab(paste0("PC1: ", percentVar[1], "% variance")) + + ylab(paste0("PC2: ", percentVar[2], "% variance")) + + scale_color_manual(values = c("Control" = "steelblue", "Tau_WT" = "salmon")) + + theme_minimal(base_size = 14) + + ggtitle("PCA: Samples by Genotype and Day") + + theme(legend.position = "right") +dev.off() + +message("✓ PCA plot saved") + +# STEP 13: Get differential expression results +# -------------------------------------------------------- +message("Extracting results: Tau vs Control...") + +res <- results(dds, + contrast = c("genotype", "Tau_WT", "Control"), + alpha = 0.05) + +# View summary +summary(res) + +# Convert to data frame and sort by adjusted p-value +res_df <- as.data.frame(res) %>% + rownames_to_column("gene") %>% + arrange(padj) + +# Save full results +write.csv(res_df, + file.path(OUT_DIR, "tables", "DE_Tau_vs_Control.csv"), + row.names = FALSE) + +# Count significant genes +n_sig <- sum(res_df$padj < 0.05, na.rm = TRUE) +n_up <- sum(res_df$padj < 0.05 & res_df$log2FoldChange > 0, na.rm = TRUE) +n_down <- sum(res_df$padj < 0.05 & res_df$log2FoldChange < 0, na.rm = TRUE) + +message("✓ Significant DE genes (padj < 0.05): ", n_sig) +message(" - Upregulated in Tau: ", n_up) +message(" - Downregulated in Tau: ", n_down) + +# STEP 14: Volcano Plot +# -------------------------------------------------------- +res_plot <- res_df %>% + mutate( + significance = case_when( + is.na(padj) ~ "NS", + padj < 0.05 & log2FoldChange > 1 ~ "Up", + padj < 0.05 & log2FoldChange < -1 ~ "Down", + TRUE ~ "NS" + ) + ) + +png(file.path(OUT_DIR, "figs", "03_volcano_plot.png"), width = 1000, height = 800, res = 150) +ggplot(res_plot, aes(x = log2FoldChange, y = -log10(pvalue), color = significance)) + + geom_point(alpha = 0.5, size = 1.5) + + scale_color_manual(values = c("Up" = "red", "Down" = "blue", "NS" = "gray70")) + + theme_minimal(base_size = 14) + + geom_vline(xintercept = c(-1, 1), linetype = "dashed", color = "gray40") + + geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "gray40") + + ggtitle("Volcano Plot: Tau vs Control") + + xlab("log2 Fold Change") + + ylab("-log10(p-value)") +dev.off() + +message("✓ Volcano plot saved") + +# STEP 15: Top DE genes heatmap +# -------------------------------------------------------- +top_genes <- res_df %>% + filter(padj < 0.05) %>% + arrange(padj) %>% + head(50) %>% + pull(gene) + +if (length(top_genes) >= 5) { + mat <- assay(vsd)[top_genes, ] + mat_z <- t(scale(t(mat))) # Z-score normalize + + anno_col <- data.frame( + Genotype = meta$genotype, + Day = meta$day, + row.names = rownames(meta) + ) + + png(file.path(OUT_DIR, "figs", "04_heatmap_top50.png"), width = 1000, height = 1200, res = 150) + pheatmap(mat_z, + annotation_col = anno_col, + cluster_cols = TRUE, + cluster_rows = TRUE, + show_rownames = TRUE, + fontsize_row = 8, + main = "Top 50 DE Genes") + dev.off() + + message("✓ Heatmap saved") +} else { + message("Not enough significant genes for heatmap") +} + +# ============================================================ +# DONE! +# ============================================================ +message("\n========================================") +message("ANALYSIS COMPLETE!") +message("========================================") +message("Results saved to: ", OUT_DIR) +message("\nOutput files:") +message(" - tables/sample_metadata.csv") +message(" - tables/DE_Tau_vs_Control.csv") +message(" - figs/01_library_sizes.png") +message(" - figs/02_PCA_plot.png") +message(" - figs/03_volcano_plot.png") +message(" - figs/04_heatmap_top50.png") +message("\nTop 10 DE genes:") +print(head(res_df, 10)) + +# ============================================================ +# PART 2: TIMEPOINT-SPECIFIC ANALYSIS & ENHANCED VISUALIZATION +# Following OSD-514 Procedure +# ============================================================ + +# Install additional packages if needed +if (!requireNamespace("apeglm", quietly = TRUE)) { + BiocManager::install("apeglm", ask = FALSE) +} +if (!requireNamespace("ggrepel", quietly = TRUE)) { + install.packages("ggrepel") +} +if (!requireNamespace("enrichR", quietly = TRUE)) { + install.packages("enrichR") +} + +library(apeglm) +library(ggrepel) + +# -------------------------------------------------------- +# STEP 16: LFC Shrinkage (apeglm) for better fold change estimates +# -------------------------------------------------------- +message("Applying LFC shrinkage...") + +# Get the coefficient name for Tau vs Control +resultsNames(dds) + +res_shrunk <- lfcShrink(dds, + coef = "genotype_Tau_WT_vs_Control", + type = "apeglm") + +# Save shrunk results +res_shrunk_df <- as.data.frame(res_shrunk) %>% + rownames_to_column("gene") %>% + arrange(padj) + +write.csv(res_shrunk_df, + file.path(OUT_DIR, "tables", "DE_Tau_vs_Control_shrunk.csv"), + row.names = FALSE) + +message("✓ LFC shrinkage complete") + +# -------------------------------------------------------- +# STEP 17: Timepoint-Specific Comparisons +# -------------------------------------------------------- +message("Running timepoint-specific analyses...") + +# Create a combined factor for genotype:day interaction +meta$group <- factor(paste(meta$genotype, meta$day, sep = "_")) + +# Rebuild DESeq2 with interaction design +dds_time <- DESeqDataSetFromMatrix( + countData = counts_filt, + colData = meta, + design = ~ group +) +dds_time <- DESeq(dds_time) + +# Function to run contrast and save results +run_contrast <- function(dds, contrast_name, num, denom) { + res <- results(dds, contrast = c("group", num, denom), alpha = 0.05) + res_df <- as.data.frame(res) %>% + rownames_to_column("gene") %>% + arrange(padj) + + # Save full results + write.csv(res_df, + file.path(OUT_DIR, "tables", paste0("DE_", contrast_name, ".csv")), + row.names = FALSE) + + # Count significant + n_sig <- sum(res_df$padj < 0.05, na.rm = TRUE) + n_up <- sum(res_df$padj < 0.05 & res_df$log2FoldChange > 1, na.rm = TRUE) + n_down <- sum(res_df$padj < 0.05 & res_df$log2FoldChange < -1, na.rm = TRUE) + + message(sprintf("%s: %d significant (↑%d, ↓%d)", + contrast_name, n_sig, n_up, n_down)) + + return(list(res = res, res_df = res_df)) +} + +# Run timepoint-specific contrasts +res_Day01 <- run_contrast(dds_time, "Tau_vs_Control_Day01", + "Tau_WT_Day01", "Control_Day01") +res_Day10 <- run_contrast(dds_time, "Tau_vs_Control_Day10", + "Tau_WT_Day10", "Control_Day10") +res_Day20 <- run_contrast(dds_time, "Tau_vs_Control_Day20", + "Tau_WT_Day20", "Control_Day20") + +message("✓ Timepoint-specific analyses complete") + +# -------------------------------------------------------- +# STEP 18: Enhanced Volcano Plots (with gene labels) +# -------------------------------------------------------- +message("Creating enhanced volcano plots...") + +save_volcano_enhanced <- function(res_df, tag, top_n_labels = 15, + thr_p = 0.05, thr_fc = 1) { + + vdf <- res_df %>% + mutate( + group = case_when( + is.na(padj) ~ "NS", + padj < thr_p & log2FoldChange >= thr_fc ~ "Up (sig)", + padj < thr_p & log2FoldChange <= -thr_fc ~ "Down (sig)", + TRUE ~ "NS" + ), + group = factor(group, levels = c("Down (sig)", "NS", "Up (sig)")) + ) + + # Select top genes to label + top_genes <- vdf %>% + filter(padj < thr_p) %>% + arrange(padj) %>% + head(top_n_labels) %>% + pull(gene) + + p <- ggplot(vdf, aes(x = log2FoldChange, y = -log10(pvalue), color = group)) + + geom_point(alpha = 0.6, size = 1.5) + + geom_vline(xintercept = c(-thr_fc, thr_fc), linetype = "dashed", color = "gray40") + + geom_hline(yintercept = -log10(thr_p), linetype = "dashed", color = "gray40") + + geom_text_repel( + data = subset(vdf, gene %in% top_genes), + aes(label = gene), + size = 3, + max.overlaps = 20, + box.padding = 0.5 + ) + + scale_color_manual(values = c("Down (sig)" = "#2C7BB6", + "NS" = "grey70", + "Up (sig)" = "#D7191C")) + + theme_classic(base_size = 12) + + theme(legend.position = "top") + + labs( + x = "log2 Fold Change", + y = "-log10(p-value)", + title = paste("Volcano:", gsub("_", " ", tag)), + subtitle = "Dashed lines: |LFC| = 1, p = 0.05" + ) + + ggsave(file.path(OUT_DIR, "figs", paste0("volcano_", tag, ".png")), + p, width = 10, height = 8, dpi = 200) + + message(" Saved: volcano_", tag, ".png") +} + +# Create volcano plots for each comparison +save_volcano_enhanced(res_df, "Tau_vs_Control_Overall") +save_volcano_enhanced(res_Day01$res_df, "Tau_vs_Control_Day01") +save_volcano_enhanced(res_Day10$res_df, "Tau_vs_Control_Day10") +save_volcano_enhanced(res_Day20$res_df, "Tau_vs_Control_Day20") + +message("✓ Volcano plots saved") + +# -------------------------------------------------------- +# STEP 19: Enhanced Heatmaps (with sample annotations) +# -------------------------------------------------------- +message("Creating enhanced heatmaps...") + +save_heatmap_enhanced <- function(res_df, tag, topN = 40, cap_z = 2.5) { + + # Get top genes + top_genes <- res_df %>% + filter(padj < 0.05) %>% + arrange(padj, -abs(log2FoldChange)) %>% + head(topN) %>% + pull(gene) + + if (length(top_genes) < 5) { + message(" Not enough significant genes for heatmap: ", tag) + return(invisible(NULL)) + } + + # Get expression matrix + mat <- assay(vsd)[top_genes, ] + + # Z-score normalize and cap + mat_z <- t(scale(t(mat))) + mat_z[mat_z > cap_z] <- cap_z + mat_z[mat_z < -cap_z] <- -cap_z + + # Annotation + ann <- data.frame( + Genotype = meta$genotype, + Day = meta$day, + row.names = rownames(meta) + ) + + # Order columns by genotype then day + col_order <- order(meta$genotype, meta$day) + mat_z <- mat_z[, col_order] + ann <- ann[col_order, , drop = FALSE] + + # Create short labels for columns + short_labels <- paste0( + ifelse(meta$genotype == "Control", "C", "T"), + "_", + gsub("Day", "D", meta$day), + "_", + meta$replicate + ) + names(short_labels) <- rownames(meta) + + # Color palette + ann_colors <- list( + Genotype = c(Control = "steelblue", Tau_WT = "salmon"), + Day = c(Day01 = "#fee8c8", Day10 = "#fdbb84", Day20 = "#e34a33") + ) + + png(file.path(OUT_DIR, "figs", paste0("heatmap_", tag, ".png")), + width = 1400, height = 1200, res = 150) + + pheatmap(mat_z, + annotation_col = ann, + annotation_colors = ann_colors, + labels_col = short_labels[colnames(mat_z)], + cluster_cols = FALSE, + cluster_rows = TRUE, + show_rownames = TRUE, + fontsize_row = 8, + fontsize_col = 9, + border_color = NA, + main = paste0("Top ", length(top_genes), " DE Genes: ", + gsub("_", " ", tag), " (Z-score, cap ±", cap_z, ")")) + + dev.off() + + message(" Saved: heatmap_", tag, ".png") +} + +# Create heatmaps +save_heatmap_enhanced(res_df, "Tau_vs_Control_Overall", topN = 50) +save_heatmap_enhanced(res_Day01$res_df, "Tau_vs_Control_Day01", topN = 40) +save_heatmap_enhanced(res_Day10$res_df, "Tau_vs_Control_Day10", topN = 40) +save_heatmap_enhanced(res_Day20$res_df, "Tau_vs_Control_Day20", topN = 40) + +message("✓ Heatmaps saved") + +# -------------------------------------------------------- +# STEP 20: Mitochondrial Gene Analysis +# -------------------------------------------------------- +message("Running mitochondrial gene analysis...") + +# Install org.Dm.eg.db if needed +if (!requireNamespace("org.Dm.eg.db", quietly = TRUE)) { + BiocManager::install("org.Dm.eg.db", ask = FALSE) +} +library(org.Dm.eg.db) +library(AnnotationDbi) +library(GO.db) + +# Curated mitochondrial GO terms (from OSD-514 procedure) +mito_go_bp <- c( + "GO:0006119", # oxidative phosphorylation + "GO:0022900", # electron transport chain + "GO:0006099", # tricarboxylic acid (TCA) cycle + "GO:0006635", # fatty acid beta-oxidation + "GO:0000422" # mitophagy +) + +mito_go_cc <- c( + "GO:0005739", # mitochondrion + "GO:0005743", # mitochondrial inner membrane + "GO:0005747", # respiratory chain complex I + "GO:0005753", # ATP synthase complex + "GO:0005759" # mitochondrial matrix +) + +# Get all genes annotated with mito GO terms +# Note: Your data uses gene symbols, so we'll query by symbol +all_mito_go <- c(mito_go_bp, mito_go_cc) + +# Get genes for each GO term +mito_genes_list <- lapply(all_mito_go, function(go_id) { + tryCatch({ + genes <- AnnotationDbi::select(org.Dm.eg.db, + keys = go_id, + keytype = "GO", + columns = "SYMBOL") + unique(na.omit(genes$SYMBOL)) + }, error = function(e) character(0)) +}) + +# Combine all mito genes +mito_genes_curated <- unique(unlist(mito_genes_list)) +message("Curated mitochondrial gene set: ", length(mito_genes_curated), " genes") + +# Find overlap with our DE genes +mito_in_data <- intersect(mito_genes_curated, rownames(counts_filt)) +message("Mito genes in our dataset: ", length(mito_in_data)) + +# Check which are significant +mito_sig <- res_df %>% + filter(gene %in% mito_in_data & padj < 0.05) + +message("Significantly DE mito genes: ", nrow(mito_sig)) + +# Save mito gene results +mito_results <- res_df %>% + filter(gene %in% mito_in_data) %>% + arrange(padj) + +write.csv(mito_results, + file.path(OUT_DIR, "tables", "DE_mitochondrial_genes.csv"), + row.names = FALSE) + +# -------------------------------------------------------- +# STEP 21: Mitochondrial Gene Volcano Plot +# -------------------------------------------------------- +save_volcano_mito <- function(res_df, mito_genes, tag, thr_p = 0.05, thr_fc = 1) { + + vdf <- res_df %>% + mutate( + is_mito = gene %in% mito_genes, + group = case_when( + is.na(padj) ~ "NS", + padj < thr_p & log2FoldChange >= thr_fc ~ "Up (sig)", + padj < thr_p & log2FoldChange <= -thr_fc ~ "Down (sig)", + TRUE ~ "NS" + ), + group = factor(group, levels = c("Down (sig)", "NS", "Up (sig)")) + ) + + # Get mito genes to label (significant ones) + mito_to_label <- vdf %>% + filter(is_mito & padj < 0.1) %>% + arrange(padj) %>% + head(20) %>% + pull(gene) + + p <- ggplot(vdf, aes(x = log2FoldChange, y = -log10(pvalue), color = group)) + + geom_point(alpha = 0.5, size = 1.5) + + # Highlight mito genes with black outline + geom_point(data = subset(vdf, is_mito), + shape = 21, stroke = 0.8, size = 2.5, + aes(fill = group), color = "black") + + geom_vline(xintercept = c(-thr_fc, thr_fc), linetype = "dashed", color = "gray40") + + geom_hline(yintercept = -log10(thr_p), linetype = "dashed", color = "gray40") + + geom_text_repel( + data = subset(vdf, gene %in% mito_to_label), + aes(label = gene), + size = 3, + max.overlaps = 30, + box.padding = 0.4 + ) + + scale_color_manual(values = c("Down (sig)" = "#2C7BB6", + "NS" = "grey70", + "Up (sig)" = "#D7191C")) + + scale_fill_manual(values = c("Down (sig)" = "#2C7BB6", + "NS" = "grey70", + "Up (sig)" = "#D7191C")) + + theme_classic(base_size = 12) + + theme(legend.position = "top") + + labs( + x = "log2 Fold Change", + y = "-log10(p-value)", + title = paste("Volcano (Mito Focus):", gsub("_", " ", tag)), + subtitle = "Black-outlined points = mitochondrial genes" + ) + + ggsave(file.path(OUT_DIR, "figs", paste0("volcano_mito_", tag, ".png")), + p, width = 10, height = 8, dpi = 200) + + message(" Saved: volcano_mito_", tag, ".png") +} + +save_volcano_mito(res_df, mito_in_data, "Tau_vs_Control") + +# -------------------------------------------------------- +# STEP 22: Mitochondrial Gene Heatmap +# -------------------------------------------------------- +if (length(mito_in_data) >= 10) { + # Get DE mito genes (relax threshold slightly) + mito_de <- res_df %>% + filter(gene %in% mito_in_data & padj < 0.1) %>% + arrange(padj) %>% + head(50) %>% + pull(gene) + + if (length(mito_de) >= 5) { + mat <- assay(vsd)[mito_de, ] + mat_z <- t(scale(t(mat))) + mat_z[mat_z > 2.5] <- 2.5 + mat_z[mat_z < -2.5] <- -2.5 + + ann <- data.frame( + Genotype = meta$genotype, + Day = meta$day, + row.names = rownames(meta) + ) + + col_order <- order(meta$genotype, meta$day) + mat_z <- mat_z[, col_order] + ann <- ann[col_order, , drop = FALSE] + + ann_colors <- list( + Genotype = c(Control = "steelblue", Tau_WT = "salmon"), + Day = c(Day01 = "#fee8c8", Day10 = "#fdbb84", Day20 = "#e34a33") + ) + + png(file.path(OUT_DIR, "figs", "heatmap_mitochondrial_genes.png"), + width = 1400, height = 1200, res = 150) + + pheatmap(mat_z, + annotation_col = ann, + annotation_colors = ann_colors, + cluster_cols = FALSE, + cluster_rows = TRUE, + show_rownames = TRUE, + fontsize_row = 8, + border_color = NA, + main = paste0("Mitochondrial DE Genes (", length(mito_de), " genes)")) + + dev.off() + + message("✓ Mitochondrial heatmap saved") + } +} + +message("\n✓ Mitochondrial analysis complete") + +# -------------------------------------------------------- +# STEP 23: GO Enrichment Analysis (using enrichR) +# -------------------------------------------------------- +message("Running GO enrichment analysis...") + +library(enrichR) + +# Set up enrichR databases +dbs <- c("GO_Biological_Process_2021", + "GO_Cellular_Component_2021", + "GO_Molecular_Function_2021", + "KEGG_2021_Drosophila") + +# Get significant gene lists +sig_up <- res_df %>% + filter(padj < 0.05 & log2FoldChange > 1) %>% + pull(gene) + +sig_down <- res_df %>% + filter(padj < 0.05 & log2FoldChange < -1) %>% + pull(gene) + +sig_all <- res_df %>% + filter(padj < 0.05 & abs(log2FoldChange) > 1) %>% + pull(gene) + +message("Upregulated genes: ", length(sig_up)) +message("Downregulated genes: ", length(sig_down)) + +# Run enrichment +if (length(sig_all) >= 10) { + + enrich_all <- enrichr(sig_all, dbs) + + # Save results + for (db_name in names(enrich_all)) { + write.csv(enrich_all[[db_name]], + file.path(EXTRAS_DIR, "tables", + paste0("GO_", gsub(" ", "_", db_name), "_all_DE.csv")), + row.names = FALSE) + } + + # Plot top GO Biological Process terms + if (nrow(enrich_all[["GO_Biological_Process_2021"]]) > 0) { + top_bp <- enrich_all[["GO_Biological_Process_2021"]] %>% + filter(Adjusted.P.value < 0.05) %>% + arrange(Adjusted.P.value) %>% + head(20) %>% + mutate(Term = substr(Term, 1, 50)) # Truncate long names + + if (nrow(top_bp) > 0) { + p <- ggplot(top_bp, aes(x = -log10(Adjusted.P.value), + y = reorder(Term, -log10(Adjusted.P.value)))) + + geom_bar(stat = "identity", fill = "steelblue") + + theme_minimal(base_size = 11) + + labs(x = "-log10(FDR)", + y = "", + title = "GO Biological Process Enrichment", + subtitle = "Tau vs Control DE Genes") + + theme(axis.text.y = element_text(size = 9)) + + ggsave(file.path(OUT_DIR, "figs", "GO_BP_enrichment.png"), + p, width = 10, height = 8, dpi = 200) + + message("✓ GO enrichment plot saved") + } + } + + # Identify mito-related enriched terms + mito_terms <- enrich_all[["GO_Biological_Process_2021"]] %>% + filter(grepl("mitochond|respiratory|oxidative|electron|ATP", + Term, ignore.case = TRUE)) + + if (nrow(mito_terms) > 0) { + write.csv(mito_terms, + file.path(EXTRAS_DIR, "tables", "GO_mitochondrial_terms.csv"), + row.names = FALSE) + message("Found ", nrow(mito_terms), " mitochondria-related GO terms") + } + +} else { + message("Not enough DE genes for enrichment analysis") +} + +# ============================================================ +# ANALYSIS COMPLETE +# ============================================================ +message("\n========================================") +message("FULL ANALYSIS COMPLETE!") +message("========================================") +message("\nOutput files in: ", OUT_DIR) +message("\nTables generated:") +message(" - DE_Tau_vs_Control.csv (overall)") +message(" - DE_Tau_vs_Control_shrunk.csv (with LFC shrinkage)") +message(" - DE_Tau_vs_Control_Day01/10/20.csv (timepoint-specific)") +message(" - DE_mitochondrial_genes.csv") +message(" - GO enrichment tables") +message("\nFigures generated:") +message(" - Volcano plots (overall + timepoint-specific + mito focus)") +message(" - Heatmaps (overall + timepoint-specific + mito)") +message(" - GO enrichment barplot") + +# Print summary +message("\n===== SUMMARY =====") +message("Total DE genes (padj < 0.05): ", sum(res_df$padj < 0.05, na.rm = TRUE)) +message(" Day01: ", sum(res_Day01$res_df$padj < 0.05, na.rm = TRUE)) +message(" Day10: ", sum(res_Day10$res_df$padj < 0.05, na.rm = TRUE)) +message(" Day20: ", sum(res_Day20$res_df$padj < 0.05, na.rm = TRUE)) +message("Mito genes in dataset: ", length(mito_in_data)) +message("Significant mito genes: ", nrow(mito_sig)) + +# ============================================================ +# PART 3: GSEA (Gene Set Enrichment Analysis) with fgsea +# Following OSD-514 Procedure +# ============================================================ + +message("Starting GSEA analysis...") + +# Install fgsea if needed +if (!requireNamespace("fgsea", quietly = TRUE)) { + BiocManager::install("fgsea", ask = FALSE) +} +if (!requireNamespace("data.table", quietly = TRUE)) { + install.packages("data.table") +} + +library(fgsea) +library(data.table) + +# -------------------------------------------------------- +# STEP 24: Build Gene Ranks +# -------------------------------------------------------- +# GSEA uses ALL genes ranked by a metric (not just significant ones) +# We'll use: sign(log2FC) * -log10(pvalue) as the ranking metric + +build_ranks <- function(res_df) { + # Remove NAs + df <- res_df %>% + filter(!is.na(pvalue) & !is.na(log2FoldChange) & pvalue > 0) + + # Create ranking metric: signed -log10(p-value) + ranks <- sign(df$log2FoldChange) * -log10(df$pvalue) + names(ranks) <- df$gene + + # Sort by rank + ranks <- sort(ranks, decreasing = TRUE) + + return(ranks) +} + +# Build ranks for overall comparison +ranks_overall <- build_ranks(res_df) +message("Ranked ", length(ranks_overall), " genes for GSEA") + +# -------------------------------------------------------- +# STEP 25: Build GO Biological Process Gene Sets +# -------------------------------------------------------- +message("Building GO BP gene sets...") + +# Get all GO BP terms for genes in our dataset +gene_list <- names(ranks_overall) + +# Query GO annotations for our genes +go_annotations <- AnnotationDbi::select( + org.Dm.eg.db, + keys = gene_list, + keytype = "SYMBOL", + columns = c("SYMBOL", "GO", "ONTOLOGY") +) + +# Filter to Biological Process only - WITH dplyr:: prefix +go_bp <- go_annotations %>% + dplyr::filter(ONTOLOGY == "BP") %>% + dplyr::filter(!is.na(GO)) %>% + dplyr::select(GO, SYMBOL) %>% + dplyr::distinct() + +# Create pathway list (GO term -> gene vector) +pathways_bp <- split(go_bp$SYMBOL, go_bp$GO) + +# Filter pathways by size (10-500 genes) +pathway_sizes <- sapply(pathways_bp, length) +pathways_bp <- pathways_bp[pathway_sizes >= 10 & pathway_sizes <= 500] + +message("Built ", length(pathways_bp), " GO BP gene sets (size 10-500)") + +# Get GO term names for labeling +go_term_names <- AnnotationDbi::select( + GO.db, + keys = names(pathways_bp), + keytype = "GOID", + columns = "TERM" +) +term_map <- setNames(go_term_names$TERM, go_term_names$GOID) + +# -------------------------------------------------------- +# STEP 26: Run fgsea +# -------------------------------------------------------- +message("Running fgsea (this may take a minute)...") + +set.seed(42) +fgsea_results <- fgsea( + pathways = pathways_bp, + stats = ranks_overall, + minSize = 10, + maxSize = 500, + nPermSimple = 10000 +) + +# Add term names +fgsea_results$term <- term_map[fgsea_results$pathway] + +# Sort by adjusted p-value +fgsea_results <- fgsea_results[order(padj)] + +message("GSEA complete: ", sum(fgsea_results$padj < 0.05, na.rm = TRUE), + " significant pathways (FDR < 0.05)") + +# -------------------------------------------------------- +# STEP 27: Save GSEA Results +# -------------------------------------------------------- + +# Save full results (collapse leadingEdge to string) +fgsea_export <- fgsea_results %>% + as.data.frame() %>% + mutate(leadingEdge = sapply(leadingEdge, paste, collapse = ";")) + +write.csv(fgsea_export, + file.path(GSEA_DIR, "tables", "fgsea_GO_BP_Tau_vs_Control.csv"), + row.names = FALSE) + +# Save significant results +fgsea_sig <- fgsea_export %>% filter(padj < 0.05) +write.csv(fgsea_sig, + file.path(GSEA_DIR, "tables", "fgsea_GO_BP_significant.csv"), + row.names = FALSE) + +message("Saved GSEA tables to: ", file.path(GSEA_DIR, "tables")) + +# -------------------------------------------------------- +# STEP 28: GSEA Barplot (Top Pathways by NES) +# -------------------------------------------------------- +message("Creating GSEA plots...") + +# Get top positive and negative NES pathways +top_up <- fgsea_results %>% + filter(padj < 0.1 & NES > 0) %>% + arrange(desc(NES)) %>% + head(10) + +top_down <- fgsea_results %>% + filter(padj < 0.1 & NES < 0) %>% + arrange(NES) %>% + head(10) + +top_pathways <- rbind(top_up, top_down) %>% + mutate(term_short = ifelse(nchar(term) > 50, + paste0(substr(term, 1, 47), "..."), + term)) + +if (nrow(top_pathways) > 0) { + p <- ggplot(top_pathways, aes(x = reorder(term_short, NES), y = NES, + fill = -log10(padj))) + + geom_col() + + coord_flip() + + scale_fill_gradient(low = "blue", high = "red") + + theme_minimal(base_size = 11) + + labs( + x = "", + y = "Normalized Enrichment Score (NES)", + fill = "-log10(FDR)", + title = "GSEA: GO Biological Process", + subtitle = "Tau vs Control - Top Enriched Pathways" + ) + + theme(axis.text.y = element_text(size = 9)) + + ggsave(file.path(GSEA_DIR, "figs", "fgsea_barplot_top_pathways.png"), + p, width = 12, height = 8, dpi = 200) + + message("✓ Saved: fgsea_barplot_top_pathways.png") +} + +# -------------------------------------------------------- +# STEP 29: Identify Mito-Related GSEA Pathways +# -------------------------------------------------------- +message("Identifying mitochondrial pathways in GSEA results...") + +# Pattern to match mito-related terms +mito_pattern <- paste( + "mitochond", "oxidative phosph", "electron transport", + "respiratory chain", "ATP synthase", "tricarboxylic", "TCA", + "beta-oxid", "mitophagy", "cristae", + sep = "|" +) + +fgsea_mito <- fgsea_results %>% + filter(grepl(mito_pattern, term, ignore.case = TRUE)) + +message("Found ", nrow(fgsea_mito), " mitochondria-related pathways") +message(" Significant (FDR < 0.05): ", sum(fgsea_mito$padj < 0.05, na.rm = TRUE)) + +# Save mito pathways +if (nrow(fgsea_mito) > 0) { + fgsea_mito_export <- fgsea_mito %>% + as.data.frame() %>% + mutate(leadingEdge = sapply(leadingEdge, paste, collapse = ";")) + + write.csv(fgsea_mito_export, + file.path(GSEA_DIR, "tables", "fgsea_mitochondrial_pathways.csv"), + row.names = FALSE) + + # Print top mito pathways + message("\nTop Mitochondrial Pathways:") + print(fgsea_mito %>% + as.data.frame() %>% + dplyr::select(term, NES, padj) %>% + head(10)) +} + +# -------------------------------------------------------- +# STEP 30: Enrichment Curves for Top Mito Pathways +# -------------------------------------------------------- + +# Plot enrichment curves for top mito-related pathways +if (nrow(fgsea_mito) > 0) { + + # Get top 4 mito pathways (by significance) + top_mito <- fgsea_mito %>% + arrange(padj) %>% + head(4) %>% + pull(pathway) + + for (pw in top_mito) { + pw_name <- term_map[pw] + if (is.na(pw_name)) pw_name <- pw + + # Create enrichment plot + p <- plotEnrichment(pathways_bp[[pw]], ranks_overall) + + ggtitle(paste0(pw_name, "\nNES = ", + round(fgsea_results[pathway == pw, NES], 2), + ", FDR = ", + signif(fgsea_results[pathway == pw, padj], 2))) + + # Clean filename + fname <- gsub("[^A-Za-z0-9]", "_", substr(pw_name, 1, 40)) + + ggsave(file.path(GSEA_DIR, "figs", paste0("enrichment_curve_", fname, ".png")), + p, width = 8, height = 5, dpi = 200) + } + + message("✓ Saved enrichment curves for top mito pathways") +} + +# -------------------------------------------------------- +# STEP 31: GSEA Table Plot (fgsea style) +# -------------------------------------------------------- + +# Create table plot for top pathways +if (nrow(top_pathways) > 0) { + top_pw_ids <- top_pathways$pathway + + png(file.path(GSEA_DIR, "figs", "fgsea_table_top_pathways.png"), + width = 1600, height = 900, res = 150) + + plotGseaTable( + pathways_bp[top_pw_ids], + ranks_overall, + fgsea_results[pathway %in% top_pw_ids], + gseaParam = 0.5 + ) + + dev.off() + + message("✓ Saved: fgsea_table_top_pathways.png") +} + +# -------------------------------------------------------- +# STEP 32: Run GSEA for Each Timepoint +# -------------------------------------------------------- +message("\nRunning GSEA for each timepoint...") + +run_gsea_timepoint <- function(res_df, tag) { + ranks <- build_ranks(res_df) + + if (length(ranks) < 100) { + message(" ", tag, ": Too few ranked genes, skipping") + return(NULL) + } + + set.seed(42) + fg <- fgsea( + pathways = pathways_bp, + stats = ranks, + minSize = 10, + maxSize = 500, + nPermSimple = 10000 + ) + + fg$term <- term_map[fg$pathway] + fg <- fg[order(padj)] + + n_sig <- sum(fg$padj < 0.05, na.rm = TRUE) + message(" ", tag, ": ", n_sig, " significant pathways") + + # Save results + fg_export <- fg %>% + as.data.frame() %>% + mutate(leadingEdge = sapply(leadingEdge, paste, collapse = ";")) + + write.csv(fg_export, + file.path(GSEA_DIR, "tables", paste0("fgsea_GO_BP_", tag, ".csv")), + row.names = FALSE) + + return(fg) +} + +fgsea_Day01 <- run_gsea_timepoint(res_Day01$res_df, "Day01") +fgsea_Day10 <- run_gsea_timepoint(res_Day10$res_df, "Day10") +fgsea_Day20 <- run_gsea_timepoint(res_Day20$res_df, "Day20") + +# ============================================================ +# GSEA COMPLETE +# ============================================================ +message("\n========================================") +message("GSEA ANALYSIS COMPLETE!") +message("========================================") +message("\nOutput in: ", GSEA_DIR) +message("\nTables:") +message(" - fgsea_GO_BP_Tau_vs_Control.csv (all pathways)") +message(" - fgsea_GO_BP_significant.csv (FDR < 0.05)") +message(" - fgsea_mitochondrial_pathways.csv") +message(" - fgsea_GO_BP_Day01/10/20.csv (timepoint-specific)") +message("\nFigures:") +message(" - fgsea_barplot_top_pathways.png") +message(" - fgsea_table_top_pathways.png") +message(" - enrichment_curve_*.png (mito pathways)") + +message("\n===== GSEA SUMMARY =====") +message("Overall Tau vs Control:") +message(" Significant pathways (FDR < 0.05): ", + sum(fgsea_results$padj < 0.05, na.rm = TRUE)) +message(" Mito-related pathways found: ", nrow(fgsea_mito)) +message(" Mito-related significant: ", + sum(fgsea_mito$padj < 0.05, na.rm = TRUE)) \ No newline at end of file