library(edgeR)
library(limma)
library(readxl)


# ---- load from Excel ----
infile <- "SupplementaryData_02_transcriptomics.xlsx"
sheet  <- "mapped_read_counts"

all_data.df <- read_excel(infile, sheet = sheet)
all_data.df <- as.data.frame(all_data.df, check.names = FALSE)

# first column = gene IDs (same as your CSV)
counts <- all_data.df[, -1, drop = FALSE]
rownames(counts) <- all_data.df[[1]]

# coerce to numeric defensively (Excel can be sneaky)
counts <- as.data.frame(lapply(counts, function(x) as.numeric(x)))
counts[is.na(counts)] <- 0

group <- c(
  "HLpDay","HLpDay","HLpDay",
  "HLmDay","HLmDay","HLmDay",
  "HLd00","HLd00","HLd00",
  "HLr00","HLr00","HLr00",
  "LLpDay","LLpDay","LLpDay",
  "LLmDay","LLmDay","LLmDay",
  "LLd00","LLd00","LLd00",
  "LLr00","LLr00","LLr00",
  "HLpNight","HLpNight","HLpNight",
  "HLmNight","HLmNight","HLmNight",
  "HLdNight","HLdNight","HLdNight",
  "HLrNight","HLrNight","HLrNight",
  "LLpNight","LLpNight","LLpNight",
  "LLmNight","LLmNight","LLmNight",
  "LLdNight","LLdNight","LLdNight",
  "LLrNight","LLrNight","LLrNight",
  "HLpNight","HLpNight","HLpNight",
  "HLmNight","HLmNight","HLmNight",
  "HLdNight","HLdNight","HLdNight",
  "HLrNight","HLrNight","HLrNight",
  "LLpNight","LLpNight","LLpNight",
  "LLmNight","LLmNight","LLmNight",
  "LLdNight","LLdNight","LLdNight",
  "LLrNight","LLrNight","LLrNight",
  "HLpDay","HLpDay","HLpDay",
  "HLmDay","HLmDay","HLmDay",
  "HLdDay","HLdDay","HLdDay",
  "HLrDay","HLrDay","HLrDay",
  "LLpDay","LLpDay","LLpDay",
  "LLmDay","LLmDay","LLmDay",
  "LLdDay","LLdDay","LLdDay",
  "LLrDay","LLrDay","LLrDay",
  "HLpDay","HLpDay","HLpDay",
  "HLmDay","HLmDay","HLmDay",
  "HLd03","HLd03","HLd03",
  "HLr03","HLr03","HLr03",
  "LLpDay","LLpDay","LLpDay",
  "LLmDay","LLmDay","LLmDay",
  "LLd03","LLd03","LLd03",
  "LLr03","LLr03","LLr03"
)

cds <- DGEList(counts, group = group)
cds <- calcNormFactors(cds)
design1 <- model.matrix(~0 + group)
cds <- estimateGLMCommonDisp(cds, design1)
cds <- estimateGLMTrendedDisp(cds, design1)
cds <- estimateGLMTagwiseDisp(cds, design1)
fit <- glmFit(cds, design1)

my.contrasts <- makeContrasts(
  HLpDayvsLLpDay     = groupHLpDay     - groupLLpDay,
  HLpNightvsLLpNight = groupHLpNight   - groupLLpNight,
  HLmDayvsLLmDay     = groupHLmDay     - groupLLmDay,
  HLmNightvsLLmNight = groupHLmNight   - groupLLmNight,
  HLpDayvsHLmDay     = groupHLpDay     - groupHLmDay,
  HLpNightvsHLmNight = groupHLpNight   - groupHLmNight,
  LLpDayvsLLmDay     = groupLLpDay     - groupLLmDay,
  LLpNightvsLLmNight = groupLLpNight   - groupLLmNight,
  HLmDayvsHLmNight   = groupHLmDay     - groupHLmNight,
  HLpDayvsHLpNight   = groupHLpDay     - groupHLpNight,
  LLmDayvsLLmNight   = groupLLmDay     - groupLLmNight,
  LLpDayvsLLpNight   = groupLLpDay     - groupLLpNight,
  levels = design1
)

res_list <- list()
for (cn in colnames(my.contrasts)) {
  lrt <- glmLRT(fit, contrast = my.contrasts[, cn])
  tt <- topTags(lrt, n = Inf)$table
  sub <- tt[, c("logFC","logCPM","LR","PValue","FDR")]
  names(sub) <- paste(c("logFC","logCPM","LR","PValue","FDR"), cn, sep = "_")
  res_list[[cn]] <- sub
}

final <- do.call(cbind, res_list)
write.csv(final, "all_contrasts_with_FDR.csv")
## ===================== Fe × Light interaction (edgeR) =====================
suppressPackageStartupMessages({library(dplyr); library(stringr)})

# Vectorized parsers
parse_light <- function(x) {
  out <- rep(NA_character_, length(x))
  out[startsWith(x, "HL")] <- "HL"
  out[startsWith(x, "LL")] <- "LL"
  out
}

parse_fe <- function(x) {
  # 3rd character in labels like "HLpDay", "LLmNight", "HLd00", "HLr03"
  substr(x, 3, 3)
}

parse_time <- function(x) {
  out <- rep(NA_character_, length(x))
  out[grepl("Night", x, ignore.case = FALSE)] <- "Night"
  out[grepl("Day",   x, ignore.case = FALSE)] <- "Day"
  # Treat coded times (00, 23, 03) as Day if still NA
  coded_day <- grepl("(00|23|03)", x)
  out[is.na(out) & coded_day] <- "Day"
  out
}

# Rebuild meta (then continue as before)
meta <- data.frame(sample = colnames(counts), group = group, stringsAsFactors = FALSE)
meta$Light     <- factor(parse_light(meta$group), levels = c("LL","HL"))
meta$Fe        <- factor(parse_fe(meta$group),    levels = c("p","m","d","r"))
meta$TimeOfDay <- factor(parse_time(meta$group),  levels = c("Day","Night"))


# Keep only +Fe (p) vs −Fe (m); drop resupply (r) and DFOB (d)
keep_idx <- meta$Fe %in% c("p","m") & !is.na(meta$Light) & !is.na(meta$TimeOfDay)
counts_int <- counts[, keep_idx, drop = FALSE]
meta_int   <- droplevels(meta[keep_idx, ])

# DGE object + filtering + TMM
dge <- DGEList(counts = counts_int)
keep_genes <- filterByExpr(dge, group = interaction(meta_int$Fe, meta_int$Light, meta_int$TimeOfDay, drop = TRUE))
dge <- dge[keep_genes,, keep.lib.sizes = FALSE]
dge <- calcNormFactors(dge, method = "TMM")

# Design with interaction and TimeOfDay covariate
design_int <- model.matrix(~ Fe * Light + TimeOfDay, data = meta_int)

# QL pipeline
dge <- estimateDisp(dge, design_int)
fit <- glmQLFit(dge, design_int)

# Find the Fe:Light interaction coefficient
int_col <- grep("Fe.*:Light|Fe:Light", colnames(design_int))
if (length(int_col) != 1) stop("Could not uniquely identify Fe:Light term.")
qlf_int <- glmQLFTest(fit, coef = int_col)

# Full table + FDR
tt_int <- topTags(qlf_int, n = Inf)$table
tt_int$gene_id <- rownames(tt_int)
tt_int$FDR <- p.adjust(tt_int$PValue, "BH")
tt_int <- tt_int[, c("gene_id","logFC","logCPM","F","PValue","FDR")]

dir.create("interaction_rna", showWarnings = FALSE)
write.csv(tt_int, file.path("interaction_rna","rna_FeXLight_interaction_edgeR.csv"), row.names = FALSE)

summary_rna <- data.frame(
  domain   = "RNA-seq",
  n_tested = nrow(tt_int),
  n_FDR_5  = sum(tt_int$FDR < 0.05, na.rm = TRUE),
  n_FDR_10 = sum(tt_int$FDR < 0.10, na.rm = TRUE)
)
write.csv(summary_rna, file.path("interaction_rna","rna_FeXLight_interaction_summary.csv"), row.names = FALSE)