# 1. Install Homebrew (macOS package manager)
/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"

# 2. Install Miniforge (Conda for bioinformatics)
brew install miniforge
conda init zsh
# ↑ Close and reopen Terminal after this!

# 3. Create dedicated RNA-seq environment
conda create -n rnaseq python=3.10 -y
conda activate rnaseq

# 4. Install all bioinformatics tools at once
conda install -c bioconda -c conda-forge \
  fastqc fastp hisat2 samtools subread sra-tools multiqc salmon -y

# 5. Verify — all should print version numbers
fastqc --version && hisat2 --version && samtools --version && salmon --version && fasterq-dump --version

install.packages("BiocManager")
BiocManager::install(c("DESeq2", "apeglm", "tximport"))
install.packages(c("ggplot2", "pheatmap", "RColorBrewer", "ggrepel", "dplyr"))

...

...

...

...

...

...

...

...

...

...

...

...

...

...

...

...

# Standard: significant genes (padj < 0.05)
sig <- res_df[res_df$padj < 0.05 & !is.na(res_df$padj), ]

# Strict: significant AND at least 2-fold change
sig_strict <- res_df[res_df$padj < 0.05 & abs(res_df$log2FoldChange) > 1 & !is.na(res_df$padj), ]

# Top 20 most significant
top20 <- head(sig[order(sig$padj), ], 20)
print(top20[, c("gene", "log2FoldChange", "padj")])

...

...

...

library(DESeq2); library(ggplot2); library(ggrepel)

# ─── Update these paths and names to match your project ───
base_path  <- "~/Desktop/RNA-Seq_yeast"   # your project folder
ctrl_name  <- "wildtype"                   # control condition label
treat_name <- "snf2_mutant"                # treatment condition label
n_ctrl     <- 3                           # number of control replicates
n_treat    <- 3                           # number of treatment replicates
n_total    <- n_ctrl + n_treat

# Load count matrix
counts_raw <- read.table(paste0(base_path, "/counts/gene_counts_all.txt"),
                          header=TRUE, skip=1, row.names=1)
counts <- counts_raw[, (ncol(counts_raw)-n_total+1):ncol(counts_raw)]

# Column names: edit to match your sample order
colnames(counts) <- c("ctrl_1","ctrl_2","ctrl_3","treat_1","treat_2","treat_3")

# ─── Define multi-factor metadata ───────────────────────────────
# EDIT batch/sex/time values to match your actual sample assignments
sample_info <- data.frame(
  condition = factor(c(ctrl_name,ctrl_name,ctrl_name, treat_name,treat_name,treat_name)),
  batch     = factor(c("batch1","batch1","batch2", "batch1","batch2","batch2")),
  # Optionally add more factors:
  # sex  = factor(c("M","F","M", "F","M","F")),
  row.names = colnames(counts)
)
print(sample_info)  # verify it looks correct

# Additive design: controls for batch, tests condition effect
dds_add <- DESeqDataSetFromMatrix(countData=counts, colData=sample_info,
                                   design=~batch + condition)
dds_add <- dds_add[rowSums(counts(dds_add)) >= 10, ]
dds_add <- DESeq(dds_add)

# Extract condition effect (batch effect is now controlled)
res_add <- results(dds_add,
                    contrast=c("condition", treat_name, ctrl_name),
                    alpha=0.05)
summary(res_add)

# Save results
res_add_df <- as.data.frame(res_add)
res_add_df$gene <- rownames(res_add_df)
res_add_df <- res_add_df[order(res_add_df$padj, na.last=TRUE), ]
write.csv(res_add_df, paste0(base_path, "/results/deseq2_additive_model.csv"), row.names=FALSE)

cat("\nAdditive model results:\n")
cat("  Significant (padj<0.05):", sum(res_add_df$padj < 0.05, na.rm=TRUE), "\n")
cat("  Up in treatment:         ", sum(res_add_df$padj < 0.05 & res_add_df$log2FoldChange > 0, na.rm=TRUE), "\n")
cat("  Down in treatment:       ", sum(res_add_df$padj < 0.05 & res_add_df$log2FoldChange < 0, na.rm=TRUE), "\n")

# Interaction model: tests if treatment effect DIFFERS by batch
dds_int <- DESeqDataSetFromMatrix(countData=counts, colData=sample_info,
                                   design=~batch * condition)
dds_int <- dds_int[rowSums(counts(dds_int)) >= 10, ]
dds_int <- DESeq(dds_int)

# List all available result names (choose the interaction term)
resultsNames(dds_int)
# Look for a name like: "batchbatch2.conditiontreat_name"

# Extract the interaction term (genes where batch modifies treatment effect)
interaction_term <- grep("batch.*condition", resultsNames(dds_int), value=TRUE)[1]
res_int <- results(dds_int, name=interaction_term, alpha=0.05)
summary(res_int)

res_int_df <- as.data.frame(res_int)
res_int_df$gene <- rownames(res_int_df)
write.csv(res_int_df, paste0(base_path, "/results/deseq2_interaction.csv"), row.names=FALSE)
cat("Genes with significant interaction:", sum(res_int_df$padj < 0.05, na.rm=TRUE), "\n")

# LRT: finds genes where REMOVING condition significantly worsens the model fit
dds_lrt <- DESeqDataSetFromMatrix(countData=counts, colData=sample_info,
                                   design=~batch + condition)
dds_lrt <- dds_lrt[rowSums(counts(dds_lrt)) >= 10, ]
dds_lrt <- DESeq(dds_lrt, test="LRT", reduced=~batch)  # reduced removes condition

res_lrt <- results(dds_lrt, alpha=0.05)
summary(res_lrt)

res_lrt_df <- as.data.frame(res_lrt)
res_lrt_df$gene <- rownames(res_lrt_df)
res_lrt_df <- res_lrt_df[order(res_lrt_df$padj, na.last=TRUE), ]
write.csv(res_lrt_df, paste0(base_path, "/results/deseq2_LRT.csv"), row.names=FALSE)
cat("Genes significant by LRT (padj<0.05):", sum(res_lrt_df$padj < 0.05, na.rm=TRUE), "\n")

# Paired design: each control is matched to a treatment from the same source
sample_info_paired <- data.frame(
  condition = factor(c(ctrl_name,ctrl_name,ctrl_name, treat_name,treat_name,treat_name)),
  patient   = factor(c("P1","P2","P3", "P1","P2","P3")),  # same patient IDs
  row.names = colnames(counts)
)

# Patient goes first — this blocks for individual variation
dds_paired <- DESeqDataSetFromMatrix(countData=counts, colData=sample_info_paired,
                                      design=~patient + condition)
dds_paired <- dds_paired[rowSums(counts(dds_paired)) >= 10, ]
dds_paired <- DESeq(dds_paired)
res_paired <- results(dds_paired, contrast=c("condition", treat_name, ctrl_name), alpha=0.05)
summary(res_paired)
cat("Significant (paired, padj<0.05):", sum(res_paired$padj < 0.05, na.rm=TRUE), "\n")

library(DESeq2); library(ggplot2); library(ggrepel)

# Run this AFTER loading counts and sample_info from the DESeq2 step
# (sample_info must have both 'condition' and 'batch' columns)
base_path <- "~/Desktop/RNA-Seq_yeast"

dds_raw <- DESeqDataSetFromMatrix(countData=counts, colData=sample_info, design=~condition)
dds_raw <- dds_raw[rowSums(counts(dds_raw)) >= 10, ]
dds_raw <- estimateSizeFactors(dds_raw)
vst_raw  <- vst(dds_raw, blind=TRUE)

# PCA: color by condition, shape by batch — batch effect = shapes separate
pca_data <- plotPCA(vst_raw, intgroup=c("condition","batch"), returnData=TRUE)
pct_var  <- round(100 * attr(pca_data, "percentVar"))

p_check <- ggplot(pca_data, aes(PC1, PC2, color=condition, shape=batch, label=name)) +
  geom_point(size=4) +
  geom_text_repel(size=3.5, show.legend=FALSE) +
  labs(title="PCA: check for batch effects",
       subtitle="If shapes separate on PC1/PC2 — batch effect present",
       x=paste0("PC1: ", pct_var[1], "% variance"),
       y=paste0("PC2: ", pct_var[2], "% variance")) +
  theme_minimal(base_size=13)
print(p_check)
ggsave(paste0(base_path, "/results/pca_batch_check.png"), p_check, width=9, height=7, dpi=300)

# Install if needed: BiocManager::install("sva")
library(sva); library(DESeq2); library(ggplot2); library(ggrepel)

# batch and group vectors must be in the same order as your count columns
# 1 = batch1, 2 = batch2  |  1 = control, 2 = treatment
batch_vec <- as.integer(sample_info$batch)   # e.g. c(1,1,2, 1,2,2)
group_vec <- as.integer(sample_info$condition) # e.g. c(1,1,1, 2,2,2)

# Apply ComBat-seq: corrects raw counts, preserves count nature
counts_mat       <- as.matrix(counts)
counts_corrected <- ComBat_seq(counts_mat, batch=batch_vec, group=group_vec)

cat("ComBat-seq done. Matrix dimensions:", dim(counts_corrected), "\n")

# Save corrected count matrix
write.csv(as.data.frame(counts_corrected),
          paste0(base_path, "/counts/gene_counts_batch_corrected.csv"))

# Re-run DESeq2 with corrected counts — simple design, no batch term needed
sample_info_simple <- data.frame(
  condition = sample_info$condition,
  row.names = colnames(counts_corrected)
)
dds_bc <- DESeqDataSetFromMatrix(countData=counts_corrected, colData=sample_info_simple,
                                  design=~condition)
dds_bc <- dds_bc[rowSums(counts(dds_bc)) >= 10, ]
dds_bc <- DESeq(dds_bc)
res_bc <- results(dds_bc, contrast=c("condition", treat_name, ctrl_name), alpha=0.05)
summary(res_bc)
res_bc_df <- as.data.frame(res_bc); res_bc_df$gene <- rownames(res_bc_df)
write.csv(res_bc_df, paste0(base_path, "/results/deseq2_combatseq_corrected.csv"), row.names=FALSE)
cat("Significant after ComBat-seq:", sum(res_bc_df$padj < 0.05, na.rm=TRUE), "\n")

# Install if needed: BiocManager::install("limma")
library(limma); library(DESeq2)

# Get VST-normalized data from your DESeq2 object
vst_data <- vst(dds_bc, blind=FALSE)
mat       <- assay(vst_data)

# Remove batch effect from the VST matrix (visualization only)
mat_corrected <- removeBatchEffect(mat, batch=batch_vec)

# Replace the assay data with corrected values for plotting
assay(vst_data) <- mat_corrected

# PCA after batch removal
pca_after <- plotPCA(vst_data, intgroup="condition") +
  ggtitle("PCA AFTER limma batch removal") +
  theme_minimal(base_size=13)
print(pca_after)
ggsave(paste0(base_path, "/results/pca_after_batch_removal.png"), pca_after, width=8, height=6, dpi=300)
cat("PCA plot saved! Use the corrected matrix for heatmaps and other visualizations.\n")

library(sva); library(DESeq2)

# Normalize raw counts for SVA input
dds_sva <- DESeqDataSetFromMatrix(countData=counts, colData=sample_info, design=~condition)
dds_sva <- dds_sva[rowSums(counts(dds_sva)) >= 10, ]
dds_sva <- estimateSizeFactors(dds_sva)
log_norm <- log2(counts(dds_sva, normalized=TRUE) + 1)

# Define full model (includes condition) and null model (no condition)
mod_full <- model.matrix(~condition, data=colData(dds_sva))
mod_null <- model.matrix(~1,          data=colData(dds_sva))

# Estimate number of surrogate variables
n_sv <- num.sv(log_norm, mod_full, method="leek")
cat("Estimated surrogate variables:", n_sv, "\n")

svobj <- sva(log_norm, mod_full, mod_null, n.sv=n_sv)

# Add surrogate variables to DESeq2 metadata
for (i in 1:n_sv) {
  colData(dds_sva)[[paste0("SV", i)]] <- svobj$sv[, i]
}
sv_terms    <- paste(paste0("SV", 1:n_sv), collapse=" + ")
design_sva  <- as.formula(paste("~", sv_terms, "+ condition"))
design(dds_sva) <- design_sva
dds_sva <- DESeq(dds_sva)
res_sva <- results(dds_sva, contrast=c("condition", treat_name, ctrl_name), alpha=0.05)
summary(res_sva)
cat("Significant after SVA correction:", sum(res_sva$padj < 0.05, na.rm=TRUE), "\n")

# Install if needed: install.packages("patchwork")
library(patchwork)

# Before: VST from uncorrected counts
dds_before <- DESeqDataSetFromMatrix(countData=counts, colData=sample_info, design=~condition)
dds_before <- dds_before[rowSums(counts(dds_before)) >= 10, ]
dds_before <- estimateSizeFactors(dds_before)
vst_before <- vst(dds_before, blind=TRUE)

# After: VST from ComBat-seq corrected counts
dds_after <- DESeqDataSetFromMatrix(countData=counts_corrected, colData=sample_info_simple, design=~condition)
dds_after <- dds_after[rowSums(counts(dds_after)) >= 10, ]
dds_after <- estimateSizeFactors(dds_after)
vst_after <- vst(dds_after, blind=TRUE)

p_before <- plotPCA(vst_before, intgroup="condition") +
  ggtitle("BEFORE batch correction") + theme_minimal(base_size=12)
p_after  <- plotPCA(vst_after,  intgroup="condition") +
  ggtitle("AFTER batch correction")  + theme_minimal(base_size=12)

ggsave(paste0(base_path, "/results/pca_before_after_batch.png"),
       p_before + p_after, width=14, height=6, dpi=300)
cat("Before/after PCA comparison saved!\n")

if (!requireNamespace("BiocManager", quietly=TRUE)) install.packages("BiocManager")

# Core enrichment packages
BiocManager::install(c("clusterProfiler", "enrichplot", "ReactomePA", "fgsea"), ask=FALSE)

# Organism annotation databases (install the ones relevant to your organism)
BiocManager::install(c(
  "org.Hs.eg.db",   # Human
  "org.Mm.eg.db",   # Mouse
  "org.Rn.eg.db",   # Rat
  "org.Sc.sgd.db",  # Yeast (S. cerevisiae)
  "org.At.tair.db", # Arabidopsis
  "org.Dm.eg.db",   # Fruit fly
  "org.Ce.eg.db",   # C. elegans
  "org.Dr.eg.db"    # Zebrafish
), ask=FALSE)

install.packages(c("ggplot2", "dplyr"))  # if not already installed

library(clusterProfiler); library(enrichplot); library(ggplot2); library(dplyr)
library(org.Sc.sgd.db)  # ← CHANGE to your organism's OrgDb package

# ─── Edit these ───────────────────────────────────────────────────────────
base_path  <- "~/Desktop/RNA-Seq_yeast"
OrgDb_pkg  <- org.Sc.sgd.db  # your organism OrgDb
key_type   <- "GENENAME"      # "SYMBOL" for most; "TAIR" for arabidopsis
kegg_org   <- "sce"           # KEGG organism code
react_org  <- "yeast"         # Reactome organism name
# ──────────────────────────────────────────────────────────────────────────

# Load DESeq2 results
res_df <- read.csv(paste0(base_path, "/results/deseq2_results.csv"))

# 1. Significant gene list for ORA (padj<0.05, |LFC|>1)
sig_genes <- res_df$gene[!is.na(res_df$padj) & res_df$padj < 0.05 &
                          abs(res_df$log2FoldChange) > 1]
cat("Significant genes (ORA input):", length(sig_genes), "\n")

# 2. Ranked gene list for GSEA (all genes, ranked by signed log10 p-value)
res_ranked <- res_df[!is.na(res_df$padj) & !is.na(res_df$log2FoldChange), ]
res_ranked$score <- sign(res_ranked$log2FoldChange) * (-log10(res_ranked$padj + 1e-300))
res_ranked <- res_ranked[order(res_ranked$score, decreasing=TRUE), ]
gene_list <- setNames(res_ranked$score, res_ranked$gene)
cat("Ranked gene list (GSEA input):", length(gene_list), "genes\n")

# GO Biological Process ORA
go_bp <- enrichGO(gene         = sig_genes,
                   OrgDb        = OrgDb_pkg,
                   keyType      = key_type,
                   ont          = "BP",        # "BP", "MF", or "CC"
                   pAdjustMethod = "BH",
                   pvalueCutoff = 0.05,
                   qvalueCutoff = 0.20,
                   readable     = TRUE)

go_df <- as.data.frame(go_bp)
cat("GO BP terms enriched:", nrow(go_df), "\n")
write.csv(go_df, paste0(base_path, "/results/GO_BP_ORA.csv"), row.names=FALSE)

# Dot plot (top 20 terms, sized by gene count, coloured by p-value)
p_dot <- dotplot(go_bp, showCategory=20, title="GO Biological Process — ORA") +
  theme_minimal(base_size=12)
ggsave(paste0(base_path, "/results/GO_BP_dotplot.png"), p_dot, width=10, height=10, dpi=300)
cat("GO dotplot saved!\n")

# Also run Molecular Function and Cellular Component
go_mf <- enrichGO(gene=sig_genes, OrgDb=OrgDb_pkg, keyType=key_type, ont="MF",
                   pAdjustMethod="BH", pvalueCutoff=0.05, readable=TRUE)
go_cc <- enrichGO(gene=sig_genes, OrgDb=OrgDb_pkg, keyType=key_type, ont="CC",
                   pAdjustMethod="BH", pvalueCutoff=0.05, readable=TRUE)
write.csv(as.data.frame(go_mf), paste0(base_path, "/results/GO_MF_ORA.csv"), row.names=FALSE)
write.csv(as.data.frame(go_cc), paste0(base_path, "/results/GO_CC_ORA.csv"), row.names=FALSE)
cat("GO MF terms:", nrow(as.data.frame(go_mf)), " | CC terms:", nrow(as.data.frame(go_cc)), "\n")

# GO GSEA (uses all ranked genes — no threshold cutoff needed)
gsea_bp <- gseGO(geneList      = gene_list,
                  OrgDb         = OrgDb_pkg,
                  keyType       = key_type,
                  ont           = "BP",
                  minGSSize     = 10,
                  maxGSSize     = 500,
                  pvalueCutoff  = 0.05,
                  pAdjustMethod = "BH",
                  verbose       = FALSE,
                  seed          = 42)

gsea_df <- as.data.frame(gsea_bp)
cat("GSEA GO BP terms significant:", nrow(gsea_df), "\n")
cat("  Activated (NES>0):",  sum(gsea_df$NES > 0), "\n")
cat("  Suppressed (NES<0):", sum(gsea_df$NES < 0), "\n")
write.csv(gsea_df, paste0(base_path, "/results/GSEA_GO_BP.csv"), row.names=FALSE)

# Ridge plot — shows distribution of ranked genes within each term
p_ridge <- ridgeplot(gsea_bp, showCategory=15) +
  theme_minimal(base_size=11) +
  labs(title="GSEA GO Biological Process — Ridge Plot")
ggsave(paste0(base_path, "/results/GSEA_ridgeplot.png"), p_ridge, width=11, height=10, dpi=300)

# Enrichment score plot for top terms
gseaplot2(gsea_bp, geneSetID=1:min(3,nrow(gsea_df)),
          title="Top GSEA enrichment scores")
cat("GSEA plots saved!\n")

# ─── KEGG ORA ────────────────────────────────────────────────────────────
# For most organisms, clusterProfiler can use gene symbols directly.
# For human/mouse, convert SYMBOL → ENTREZID first (see note below).
kegg_ora <- enrichKEGG(gene          = sig_genes,
                        organism      = kegg_org,
                        pAdjustMethod = "BH",
                        pvalueCutoff  = 0.05)

kegg_df <- as.data.frame(kegg_ora)
cat("KEGG pathways enriched:", nrow(kegg_df), "\n")
write.csv(kegg_df, paste0(base_path, "/results/KEGG_ORA.csv"), row.names=FALSE)

# KEGG bar plot
p_kegg <- barplot(kegg_ora, showCategory=20, title="KEGG Pathway Enrichment") +
  theme_minimal(base_size=12)
ggsave(paste0(base_path, "/results/KEGG_barplot.png"), p_kegg, width=10, height=8, dpi=300)

# ─── KEGG GSEA ───────────────────────────────────────────────────────────
gsea_kegg <- gseKEGG(geneList      = gene_list,
                      organism      = kegg_org,
                      minGSSize     = 10,
                      pvalueCutoff  = 0.05,
                      pAdjustMethod = "BH",
                      verbose       = FALSE)
write.csv(as.data.frame(gsea_kegg), paste0(base_path, "/results/KEGG_GSEA.csv"), row.names=FALSE)
cat("KEGG GSEA terms:", nrow(as.data.frame(gsea_kegg)), "\n")

# Note for human/mouse — must convert gene IDs to ENTREZ first:
# entrez <- bitr(sig_genes, fromType="SYMBOL", toType="ENTREZID", OrgDb=org.Hs.eg.db)$ENTREZID
# Then use entrez as the gene= argument and gene_list with ENTREZ names for GSEA

library(ReactomePA)

# Convert gene names to Entrez IDs (required by ReactomePA)
id_map <- bitr(sig_genes, fromType=key_type, toType="ENTREZID", OrgDb=OrgDb_pkg)
entrez_sig <- id_map$ENTREZID
cat("Mapped to Entrez:", length(entrez_sig), "/", length(sig_genes), "genes\n")

# Reactome ORA
react_ora <- enrichPathway(gene          = entrez_sig,
                             organism      = react_org,
                             pAdjustMethod = "BH",
                             pvalueCutoff  = 0.05,
                             readable      = TRUE)

react_df <- as.data.frame(react_ora)
cat("Reactome pathways enriched:", nrow(react_df), "\n")
write.csv(react_df, paste0(base_path, "/results/Reactome_ORA.csv"), row.names=FALSE)

p_react <- dotplot(react_ora, showCategory=20, title="Reactome Pathway ORA") +
  theme_minimal(base_size=12)
ggsave(paste0(base_path, "/results/Reactome_dotplot.png"), p_react, width=10, height=10, dpi=300)

# Reactome GSEA — convert ranked list to Entrez IDs
id_map_all <- bitr(res_ranked$gene, fromType=key_type, toType="ENTREZID", OrgDb=OrgDb_pkg)
ranked_entrez <- setNames(
  res_ranked$score[match(id_map_all[[key_type]], res_ranked$gene)],
  id_map_all$ENTREZID
)
ranked_entrez <- sort(ranked_entrez, decreasing=TRUE)
ranked_entrez <- ranked_entrez[!duplicated(names(ranked_entrez))]

gsea_react <- gsePathway(geneList=ranked_entrez, organism=react_org,
                           pvalueCutoff=0.05, pAdjustMethod="BH", verbose=FALSE)
write.csv(as.data.frame(gsea_react), paste0(base_path, "/results/Reactome_GSEA.csv"), row.names=FALSE)
cat("Reactome GSEA terms:", nrow(as.data.frame(gsea_react)), "\n")

library(enrichplot)

# ─── Enrichment map: shows overlap between GO terms ───
go_bp_sim <- pairwise_termsim(go_bp)  # compute term similarity
p_emap <- emapplot(go_bp_sim, showCategory=30, layout="kk") +
  ggtitle("GO BP Enrichment Map")
ggsave(paste0(base_path, "/results/GO_enrichment_map.png"), p_emap,
       width=13, height=11, dpi=300)

# ─── Concept network: genes connected to their GO terms ───
# Shows which specific genes drive each enriched term
p_cnet <- cnetplot(go_bp, showCategory=8, colorEdge=TRUE,
                    node_label="gene",
                    foldChange=setNames(res_df$log2FoldChange, res_df$gene)) +
  ggtitle("Gene-Concept Network (top 8 GO terms)") +
  scale_color_gradient2(name="log2FC", low="#58a6ff", mid="white", high="#f0883e")
ggsave(paste0(base_path, "/results/GO_cnetplot.png"), p_cnet,
       width=13, height=11, dpi=300)

# ─── Upset plot: genes shared across multiple GO terms ───
p_upset <- upsetplot(go_bp, n=10)
ggsave(paste0(base_path, "/results/GO_upsetplot.png"), p_upset,
       width=12, height=6, dpi=300)
cat("Network and upset plots saved!\n")

...

...