############################################################
# edgeR and DESeq2 official-style DEG analysis
# for GSE121212_GeneCounts.txt
#
# Purpose:
#   Compare non-lesional vs lesional samples
#   using edgeR and DESeq2 standard workflows.
#
# Input file format:
#   Row 1: status   non-lesional / lesional
#   Row 2: patient  AD_4, AD_5, ...
#   Row 3: ID       GSM IDs
#   Row 4 onward: Gene Counts
#
# Output:
#   1. edgeR_official_all_results.txt
#   2. edgeR_official_p005_fc1.4.txt
#   3. DESeq2_official_all_results.txt
#   4. DESeq2_official_p005_fc1.4.txt
#   5. Optional FC-threshold test results
#   6. summary_gene_counts.txt
############################################################

############################
# 0. Settings
############################

INPUT_FILE <- "GSE121212_GeneCounts.txt"

GROUP_ROW <- 1
PATIENT_ROW <- 2
SAMPLE_ID_ROW <- 3
METADATA_ROWS <- 3

CONTROL_GROUP <- "non-lesional"
CASE_GROUP    <- "lesional"

FC_CUTOFF <- 1.4
P_CUTOFF  <- 0.05

OUTDIR <- "official_protocol_DEG_results"
dir.create(OUTDIR, showWarnings = FALSE)

############################
# 1. Install and load packages
############################

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

required_pkgs <- c("edgeR", "DESeq2")

for (pkg in required_pkgs) {
  if (!requireNamespace(pkg, quietly = TRUE)) {
    BiocManager::install(pkg, ask = FALSE, update = FALSE)
  }
}

library(edgeR)
library(DESeq2)

############################
# 2. Read count table
############################

raw <- read.delim(
  INPUT_FILE,
  header = FALSE,
  check.names = FALSE,
  stringsAsFactors = FALSE,
  na.strings = c("NaN", "NA", "")
)

# Metadata rows
group_original <- as.character(unlist(raw[GROUP_ROW, -1]))
patient <- as.character(unlist(raw[PATIENT_ROW, -1]))
sample_names <- as.character(unlist(raw[SAMPLE_ID_ROW, -1]))

# Convert group names to R-safe names
group <- make.names(group_original)

CONTROL_GROUP <- make.names(CONTROL_GROUP)
CASE_GROUP    <- make.names(CASE_GROUP)

if (!all(c(CONTROL_GROUP, CASE_GROUP) %in% group)) {
  stop(
    "CONTROL_GROUP or CASE_GROUP was not found in the group row.\n\n",
    "CONTROL_GROUP after make.names(): ", CONTROL_GROUP, "\n",
    "CASE_GROUP after make.names(): ", CASE_GROUP, "\n",
    "Groups found after make.names(): ", paste(unique(group), collapse = ", "), "\n\n",
    "Original groups found: ", paste(unique(group_original), collapse = ", ")
  )
}

# Count table starts after metadata rows
count_data <- raw[(METADATA_ROWS + 1):nrow(raw), ]

gene_id <- as.character(count_data[[1]])

counts <- count_data[, -1, drop = FALSE]
counts <- as.data.frame(lapply(counts, function(x) as.numeric(as.character(x))))

rownames(counts) <- gene_id
colnames(counts) <- sample_names

# Remove genes with missing or empty Gene IDs
valid_gene_id <- !is.na(rownames(counts)) & rownames(counts) != ""
counts <- counts[valid_gene_id, , drop = FALSE]

# Remove genes containing NaN / NA values
# Do not replace NaN with 0.
valid_counts <- rowSums(is.na(counts)) == 0
counts <- counts[valid_counts, , drop = FALSE]

# If duplicated Gene IDs exist, sum their counts
if (any(duplicated(rownames(counts)))) {
  counts <- rowsum(as.matrix(counts), group = rownames(counts))
} else {
  counts <- as.matrix(counts)
}

# Counts must be integer-like
counts <- round(counts)

sample_info <- data.frame(
  sample = colnames(counts),
  group = factor(group, levels = c(CONTROL_GROUP, CASE_GROUP)),
  patient = patient,
  stringsAsFactors = FALSE
)

rownames(sample_info) <- sample_info$sample

# Keep only the two groups used in this comparison
keep_samples <- sample_info$group %in% c(CONTROL_GROUP, CASE_GROUP)
counts <- counts[, keep_samples, drop = FALSE]
sample_info <- sample_info[keep_samples, , drop = FALSE]
sample_info$group <- droplevels(sample_info$group)

cat("Samples used:\n")
print(table(sample_info$group))

cat("\nGenes after removing NA/NaN rows:\n")
print(nrow(counts))

############################
# 3. Helper functions
############################

add_fc_columns <- function(df, logfc_col = "logFC") {
  df$FoldChange <- ifelse(
    df[[logfc_col]] >= 0,
    2 ^ df[[logfc_col]],
    1 / (2 ^ df[[logfc_col]])
  )
  
  df$Direction <- ifelse(df[[logfc_col]] >= 0, "Up", "Down")
  
  return(df)
}

write_result <- function(df, filename) {
  write.table(
    df,
    file = file.path(OUTDIR, filename),
    sep = "\t",
    quote = FALSE,
    row.names = FALSE
  )
}

############################
# 4. edgeR official-style workflow
############################
#
# Workflow:
#   DGEList
#   filterByExpr
#   TMM normalization
#   estimateDisp
#   glmQLFit
#   glmQLFTest
#
# This is a non-paired two-group comparison.
############################

group_edger <- sample_info$group

design_edger <- model.matrix(~ group_edger)

y <- DGEList(
  counts = counts,
  group = group_edger
)

# Low-expression gene filtering
keep_edger <- filterByExpr(
  y,
  design = design_edger
)

y_filtered <- y[keep_edger, , keep.lib.sizes = FALSE]

cat("\nedgeR genes after filterByExpr:\n")
print(nrow(y_filtered))

# TMM normalization
y_filtered <- calcNormFactors(
  y_filtered,
  method = "TMM"
)

# Dispersion estimation
y_filtered <- estimateDisp(
  y_filtered,
  design_edger
)

# Quasi-likelihood GLM
fit_edger <- glmQLFit(
  y_filtered,
  design_edger,
  robust = TRUE
)

# Coefficient 2 corresponds to lesional vs non-lesional
qlf_edger <- glmQLFTest(
  fit_edger,
  coef = 2
)

edger_all <- topTags(
  qlf_edger,
  n = Inf,
  sort.by = "none"
)$table

edger_all$GeneID <- rownames(edger_all)

edger_all <- edger_all[, c(
  "GeneID",
  "logFC",
  "logCPM",
  "F",
  "PValue",
  "FDR"
)]

edger_all <- add_fc_columns(
  edger_all,
  logfc_col = "logFC"
)

edger_all <- edger_all[order(edger_all$PValue), ]

edger_p005_fc14 <- subset(
  edger_all,
  PValue < P_CUTOFF & FoldChange >= FC_CUTOFF
)

write_result(
  edger_all,
  "edgeR_official_all_results.txt"
)

write_result(
  edger_p005_fc14,
  "edgeR_official_p005_fc1.4.txt"
)

############################
# 5. DESeq2 official-style workflow
############################
#
# Workflow:
#   DESeqDataSetFromMatrix
#   pre-filtering based on the DESeq2 vignette example
#   DESeq
#   results with independentFiltering = TRUE
#
# This is a non-paired two-group comparison.
############################

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData = sample_info,
  design = ~ group
)

# Pre-filtering based on the DESeq2 vignette example
#
# Keep genes that have at least 10 counts
# in at least as many samples as the smallest group size.
#
# In this dataset:
#   non-lesional: 20 samples
#   lesional:     20 samples
#
# Therefore, genes are kept if they have >= 10 counts
# in at least 20 samples among the 40 samples.

smallestGroupSize <- min(table(sample_info$group))

keep_deseq2 <- rowSums(counts(dds) >= 10) >= smallestGroupSize
dds <- dds[keep_deseq2, ]

cat("\nDESeq2 smallest group size:\n")
print(smallestGroupSize)

cat("\nDESeq2 genes after vignette-style pre-filtering:\n")
print(nrow(dds))

dds <- DESeq(dds)

res_deseq2 <- results(
  dds,
  contrast = c("group", CASE_GROUP, CONTROL_GROUP),
  alpha = P_CUTOFF,
  independentFiltering = TRUE
)

deseq2_all <- as.data.frame(res_deseq2)
deseq2_all$GeneID <- rownames(deseq2_all)

deseq2_all <- deseq2_all[, c(
  "GeneID",
  "baseMean",
  "log2FoldChange",
  "lfcSE",
  "stat",
  "pvalue",
  "padj"
)]

names(deseq2_all) <- c(
  "GeneID",
  "baseMean",
  "logFC",
  "lfcSE",
  "stat",
  "PValue",
  "FDR"
)

deseq2_all <- add_fc_columns(
  deseq2_all,
  logfc_col = "logFC"
)

# Exclude genes with NA p-values
deseq2_all_valid <- subset(
  deseq2_all,
  !is.na(PValue)
)

deseq2_all_valid <- deseq2_all_valid[order(deseq2_all_valid$PValue), ]

deseq2_p005_fc14 <- subset(
  deseq2_all_valid,
  PValue < P_CUTOFF & FoldChange >= FC_CUTOFF
)

write_result(
  deseq2_all_valid,
  "DESeq2_official_all_results.txt"
)

write_result(
  deseq2_p005_fc14,
  "DESeq2_official_p005_fc1.4.txt"
)

############################
# 6. Optional: FC-threshold tests
############################
#
# These are stricter tests.
#
# Main article-style condition:
#   p < 0.05 AND observed fold change >= 1.4
#
# FC-threshold test:
#   tests whether the expression change is significantly
#   greater than 1.4-fold.
#
# Therefore, these results should be treated as supplemental.
############################

LOG2_FC_THRESHOLD <- log2(FC_CUTOFF)

############################
# 6-1. edgeR glmTreat
############################

treat_edger <- glmTreat(
  fit_edger,
  coef = 2,
  lfc = LOG2_FC_THRESHOLD
)

edger_treat_all <- topTags(
  treat_edger,
  n = Inf,
  sort.by = "none"
)$table

edger_treat_all$GeneID <- rownames(edger_treat_all)

edger_treat_all <- edger_treat_all[, c(
  "GeneID",
  "logFC",
  "logCPM",
  "PValue",
  "FDR"
)]

edger_treat_all <- add_fc_columns(
  edger_treat_all,
  logfc_col = "logFC"
)

edger_treat_all <- edger_treat_all[order(edger_treat_all$PValue), ]

edger_treat_p005 <- subset(
  edger_treat_all,
  PValue < P_CUTOFF
)

write_result(
  edger_treat_all,
  "edgeR_glmTreat_fc1.4_all_results.txt"
)

write_result(
  edger_treat_p005,
  "edgeR_glmTreat_fc1.4_p005.txt"
)

############################
# 6-2. DESeq2 lfcThreshold test
############################

res_deseq2_lfc <- results(
  dds,
  contrast = c("group", CASE_GROUP, CONTROL_GROUP),
  alpha = P_CUTOFF,
  independentFiltering = TRUE,
  lfcThreshold = LOG2_FC_THRESHOLD,
  altHypothesis = "greaterAbs"
)

deseq2_lfc_all <- as.data.frame(res_deseq2_lfc)
deseq2_lfc_all$GeneID <- rownames(deseq2_lfc_all)

deseq2_lfc_all <- deseq2_lfc_all[, c(
  "GeneID",
  "baseMean",
  "log2FoldChange",
  "lfcSE",
  "stat",
  "pvalue",
  "padj"
)]

names(deseq2_lfc_all) <- c(
  "GeneID",
  "baseMean",
  "logFC",
  "lfcSE",
  "stat",
  "PValue",
  "FDR"
)

deseq2_lfc_all <- add_fc_columns(
  deseq2_lfc_all,
  logfc_col = "logFC"
)

deseq2_lfc_all <- subset(
  deseq2_lfc_all,
  !is.na(PValue)
)

deseq2_lfc_all <- deseq2_lfc_all[order(deseq2_lfc_all$PValue), ]

deseq2_lfc_p005 <- subset(
  deseq2_lfc_all,
  PValue < P_CUTOFF
)

write_result(
  deseq2_lfc_all,
  "DESeq2_lfcThreshold_fc1.4_all_results.txt"
)

write_result(
  deseq2_lfc_p005,
  "DESeq2_lfcThreshold_fc1.4_p005.txt"
)

############################
# 7. Summary
############################

summary_table <- data.frame(
  Method = c(
    "edgeR official-style: p < 0.05 and FC >= 1.4",
    "DESeq2 official-style: p < 0.05 and FC >= 1.4",
    "edgeR glmTreat: FC-threshold test",
    "DESeq2 lfcThreshold: FC-threshold test"
  ),
  TestedGenes = c(
    nrow(edger_all),
    nrow(deseq2_all_valid),
    nrow(edger_treat_all),
    nrow(deseq2_lfc_all)
  ),
  SelectedGenes = c(
    nrow(edger_p005_fc14),
    nrow(deseq2_p005_fc14),
    nrow(edger_treat_p005),
    nrow(deseq2_lfc_p005)
  )
)

write_result(
  summary_table,
  "summary_gene_counts.txt"
)

cat("\nSummary:\n")
print(summary_table)

cat("\nFinished.\n")
cat("Results were written to: ", OUTDIR, "\n")