needed <- c("DESeq2", "SummarizedExperiment", "clusterProfiler",
            "AnnotationDbi", "org.Hs.eg.db")
available <- vapply(needed, requireNamespace, logical(1), quietly = TRUE)
if (!all(available)) stop("Prepare packages first: ",
                         paste(needed[!available], collapse = ", "))
suppressPackageStartupMessages({
  library(DESeq2)
  library(SummarizedExperiment)
  library(clusterProfiler)
  library(org.Hs.eg.db)
})
args <- commandArgs(trailingOnly = TRUE)
if (length(args) != 1L || !file.exists(args[1])) {
  stop("Usage: Rscript airway_qc_enrichment_learning.R AIRWAY_BUNDLE_PATH")
}
bundle <- readRDS(args[1])
stopifnot(identical(bundle$input, "Public airway teaching dataset"))
dds <- bundle$dds
res <- bundle$res
stopifnot(inherits(dds, "DESeqDataSet"), inherits(res, "DESeqResults"),
          identical(rownames(dds), rownames(res)),
          identical(colnames(dds), rownames(colData(dds))))
out <- tempfile("perry-airway-qc-", tmpdir = getwd())
stopifnot(dir.create(out))

qc <- data.frame(
  sample_id = colnames(dds),
  condition = colData(dds)$condition,
  subject = colData(dds)$subject,
  modeled_gene_count_sum = colSums(counts(dds)),
  modeled_genes_with_nonzero_count = colSums(counts(dds) > 0),
  size_factor = sizeFactors(dds)
)
write.table(qc, file.path(out, "sample_count_qc.tsv"),
             sep = "\t", row.names = FALSE, quote = FALSE)

vsd <- varianceStabilizingTransformation(dds, blind = FALSE)
png(file.path(out, "pca.png"), width = 1200, height = 850, res = 120)
print(plotPCA(vsd, intgroup = c("condition", "subject")))
dev.off()
vmat <- assay(vsd)
distances <- as.matrix(dist(t(vmat)))
png(file.path(out, "sample_distance.png"), width = 1100, height = 1100)
heatmap(distances, symm = TRUE, scale = "none",
         main = "Airway teaching data: sample distance")
dev.off()

variances <- apply(vmat, 1, var)
top <- head(order(variances, decreasing = TRUE), 30L)
centered <- sweep(vmat[top, , drop = FALSE], 1,
                  rowMeans(vmat[top, , drop = FALSE]), "-")
png(file.path(out, "high_variance_genes.png"), width = 1200, height = 1200)
heatmap(centered, scale = "none", Colv = NA,
         main = "High-variance genes: row-centered VST")
dev.off()
capture.output(sessionInfo(), file = file.path(out, "sessionInfo.txt"))

ids <- rownames(res)
stopifnot(!anyDuplicated(ids), all(grepl("^ENSG[0-9]+$", ids)))
orgdb <- org.Hs.eg.db
stopifnot("ENSEMBL" %in% AnnotationDbi::keytypes(orgdb))
known <- intersect(ids, AnnotationDbi::keys(orgdb, keytype = "ENSEMBL"))
stopifnot(length(known) > 0L)
annotation <- AnnotationDbi::select(
  orgdb, keys = known, keytype = "ENSEMBL",
  columns = c("GOALL", "ONTOLOGYALL")
)
bp_rows <- !is.na(annotation$GOALL) & !is.na(annotation$ONTOLOGYALL) &
           annotation$ONTOLOGYALL == "BP"
bp_ids <- unique(annotation$ENSEMBL[bp_rows])
fdr_eligible <- is.finite(res$pvalue) & is.finite(res$padj)
universe <- intersect(ids[fdr_eligible], bp_ids)
candidate_up <- is.finite(res$padj) & res$padj < 0.05 &
                is.finite(res$log2FoldChange) & res$log2FoldChange > 0
foreground <- intersect(ids[candidate_up], universe)
stopifnot(length(universe) > 0L, all(foreground %in% universe))
writeLines(universe, file.path(out, "ora_background_ensembl.txt"))
writeLines(foreground, file.path(out, "ora_positive_foreground_ensembl.txt"))
writeLines(setdiff(ids, known), file.path(out, "not_in_annotation_keyspace.txt"))
capture.output(orgdb, file = file.path(out, "annotation_database_info.txt"))
audit <- data.frame(
  definition = c("modeled_genes", "finite_p_and_padj", "BP_ORA_background",
                 "positive_FDR_candidates_in_background"),
  n = c(length(ids), sum(fdr_eligible), length(universe), length(foreground))
)
write.table(audit, file.path(out, "enrichment_input_audit.tsv"),
             sep = "\t", quote = FALSE, row.names = FALSE)

if (length(foreground) > 0L) {
  ora <- enrichGO(gene = foreground, universe = universe, OrgDb = orgdb,
                   keyType = "ENSEMBL", ont = "BP", pAdjustMethod = "BH",
                   pvalueCutoff = 1, qvalueCutoff = 1,
                   minGSSize = 10, maxGSSize = 500, readable = FALSE)
  ora_table <- as.data.frame(ora)
  write.table(ora_table, file.path(out, "ora_positive_BP_returned_terms.tsv"),
               sep = "\t", quote = FALSE, row.names = FALSE)
} else {
  message("No foreground genes: ORA skipped, not an enrichment result")
}

# A separate eligibility rule for preranked GSEA, not the ORA foreground.
# Keep finite test statistics even if independent filtering made padj NA.
rank_ok <- is.finite(res$stat) & is.finite(res$pvalue) & ids %in% bp_ids
ranking <- setNames(res$stat[rank_ok], ids[rank_ok])
ranking <- sort(ranking, decreasing = TRUE)
stopifnot(length(ranking) >= 10L, !anyDuplicated(names(ranking)))
write.table(data.frame(gene_id = names(ranking), Wald_stat = unname(ranking)),
             file.path(out, "gsea_ranked_genes.tsv"),
             sep = "\t", quote = FALSE, row.names = FALSE)
# Use gseGO's fixed default seed; preserve the actual package version.
gsea <- gseGO(geneList = ranking, OrgDb = orgdb, keyType = "ENSEMBL",
               ont = "BP", minGSSize = 10, maxGSSize = 500,
               pvalueCutoff = 1, pAdjustMethod = "BH",
               verbose = FALSE, seed = TRUE)
write.table(as.data.frame(gsea), file.path(out, "gsea_BP_returned_terms.tsv"),
             sep = "\t", quote = FALSE, row.names = FALSE)
saveRDS(list(background = universe, foreground = foreground,
             ranking = ranking, gsea = gsea,
             gsea_seed_setting = TRUE,
             contrast = bundle$contrast, dataset = bundle$input),
         file.path(out, "enrichment_bundle.rds"))
message("Teaching QC/enrichment outputs: ", normalizePath(out))