fix: harden pipeline checks

Author: Ekin

scripts/00_get_data.R

@@ -60,6 +60,19 @@ if (!force_download && file.exists(counts_out) && file.exists(metadata_out)) {
   }
 
   sample_ids <- sub(".*\\b(POS_\\d+|NEG_\\d+)\\b.*", "\\1", pheno$title)
+  sample_ids[!grepl("^(POS|NEG)_\\d+$", sample_ids)] <- NA_character_
+
+  # Fall back to GEO accessions if title parsing fails for a subset of records.
+  if ("geo_accession" %in% colnames(pheno)) {
+    sample_ids[is.na(sample_ids)] <- pheno$geo_accession[is.na(sample_ids)]
+  }
+
+  if (any(is.na(sample_ids))) {
+    stop("Could not derive sample IDs for all records from GEO metadata.")
+  }
+  if (any(duplicated(sample_ids))) {
+    stop("Duplicate sample IDs detected in GEO metadata.")
+  }
 
   metadata <- data.frame(
     sample_id = sample_ids,
@@ -71,9 +84,28 @@ if (!force_download && file.exists(counts_out) && file.exists(metadata_out)) {
   metadata$condition <- factor(metadata$condition, levels = c("negative", "positive"))
 
   common <- intersect(colnames(raw_counts), rownames(metadata))
+  if (length(common) == 0) {
+    stop("No overlapping sample IDs between count matrix and metadata.")
+  }
   raw_counts <- raw_counts[, common, drop = FALSE]
   metadata <- metadata[common, , drop = FALSE]
 
+  missing_condition <- is.na(metadata$condition)
+  if (any(missing_condition)) {
+    warning(
+      sprintf(
+        "Dropping %d samples with unknown condition labels after alignment.",
+        sum(missing_condition)
+      )
+    )
+    metadata <- metadata[!missing_condition, , drop = FALSE]
+    raw_counts <- raw_counts[, rownames(metadata), drop = FALSE]
+  }
+
+  if (ncol(raw_counts) == 0) {
+    stop("No samples left after metadata alignment/condition filtering.")
+  }
+
   stopifnot(is.data.frame(metadata))
   stopifnot("condition" %in% colnames(metadata))
 

scripts/01_qc.R

@@ -17,6 +17,12 @@ counts <- readRDS("data/counts_raw.rds")
 metadata <- readRDS("data/metadata.rds")
 
 metadata <- metadata[colnames(counts), , drop = FALSE]
+if (!all(colnames(counts) == rownames(metadata))) {
+  stop("Counts and metadata are not aligned by sample ID.")
+}
+if (any(is.na(metadata$condition))) {
+  stop("Metadata contains missing condition labels. Re-run scripts/00_get_data.R.")
+}
 
 message("Starting QC: ", nrow(counts), " genes, ", ncol(counts), " samples")
 
@@ -47,9 +53,18 @@ set.seed(123)
 n <- 30
 pos_samples <- rownames(metadata)[metadata$condition == "positive"]
 neg_samples <- rownames(metadata)[metadata$condition == "negative"]
+target_n <- min(n, length(pos_samples), length(neg_samples))
+
+if (target_n == 0) {
+  stop("At least one condition has zero samples after QC filtering.")
+}
+if (target_n < n) {
+  message("Using ", target_n, " samples per condition (limited by available data).")
+}
+
 balanced <- c(
-  sample(pos_samples, min(n, length(pos_samples))),
-  sample(neg_samples, min(n, length(neg_samples)))
+  sample(pos_samples, target_n),
+  sample(neg_samples, target_n)
 )
 counts <- counts[, balanced, drop = FALSE]
 metadata <- metadata[balanced, , drop = FALSE]
@@ -76,4 +91,4 @@ dds <- DESeqDataSetFromMatrix(counts, metadata, design = ~ condition)
 vsd <- varianceStabilizingTransformation(dds, blind = TRUE)
 saveRDS(vsd, "data/vst_data.rds")
 
-message("QC complete")
+message("QC complete")

scripts/06_enrichment.R

@@ -21,11 +21,15 @@ message("Running enrichment analysis...")
 sig_genes <- res_df %>%
   filter(padj < 0.05, abs(log2FoldChange) > 1) %>%
   pull(gene)
+sig_genes <- unique(sig_genes)
 
 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(length(sig_genes), " DE genes (", n_up, " up, ", n_down, " down)")
+if (length(sig_genes) == 0) {
+  stop("No DE genes passed padj/log2FC thresholds; enrichment cannot proceed.")
+}
 
 gene_map <- tryCatch({
   bitr(sig_genes, fromType = "ENSEMBL",
@@ -55,8 +59,8 @@ go_bp <- enrichGO(
 )
 
 go_df <- as.data.frame(go_bp)
+write.csv(go_df, "results/tables/go_biological_process.csv", row.names = FALSE)
 if (nrow(go_df) > 0) {
-  write.csv(go_df, "results/tables/go_biological_process.csv", row.names = FALSE)
   message("  ", nrow(go_df), " GO terms enriched")
 
   p1 <- dotplot(go_bp, showCategory = 20, font.size = 9) +
@@ -65,26 +69,35 @@ if (nrow(go_df) > 0) {
   ggsave("results/figures/go_dotplot.png", p1, width = 10, height = 9, dpi = 300)
 }
 
-kegg <- enrichKEGG(
-  gene = gene_map$ENTREZID,
-  organism = "hsa",
-  pAdjustMethod = "BH",
-  pvalueCutoff = 0.05,
-  qvalueCutoff = 0.1
+kegg <- tryCatch(
+  enrichKEGG(
+    gene = gene_map$ENTREZID,
+    organism = "hsa",
+    pAdjustMethod = "BH",
+    pvalueCutoff = 0.05,
+    qvalueCutoff = 0.1
+  ),
+  error = function(e) {
+    warning("KEGG enrichment unavailable (network/service issue): ", conditionMessage(e))
+    NULL
+  }
 )
 
-kegg_df <- as.data.frame(kegg)
-if (nrow(kegg_df) > 0) {
-  kegg_read <- setReadable(kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
-  kegg_df <- as.data.frame(kegg_read)
-  write.csv(kegg_df, "results/tables/kegg_pathways.csv", row.names = FALSE)
-  message("  ", nrow(kegg_df), " KEGG pathways enriched")
-  
-  p2 <- dotplot(kegg_read, showCategory = 20, font.size = 9) +
-    ggtitle("KEGG Pathway Enrichment") +
-    theme(plot.title = element_text(face = "bold", hjust = 0.5))
-  ggsave("results/figures/kegg_dotplot.png", p2, width = 10, height = 9, dpi = 300)
+kegg_df <- data.frame()
+if (!is.null(kegg)) {
+  kegg_df <- as.data.frame(kegg)
+  if (nrow(kegg_df) > 0) {
+    kegg_read <- setReadable(kegg, OrgDb = org.Hs.eg.db, keyType = "ENTREZID")
+    kegg_df <- as.data.frame(kegg_read)
+    message("  ", nrow(kegg_df), " KEGG pathways enriched")
+
+    p2 <- dotplot(kegg_read, showCategory = 20, font.size = 9) +
+      ggtitle("KEGG Pathway Enrichment") +
+      theme(plot.title = element_text(face = "bold", hjust = 0.5))
+    ggsave("results/figures/kegg_dotplot.png", p2, width = 10, height = 9, dpi = 300)
+  }
 }
+write.csv(kegg_df, "results/tables/kegg_pathways.csv", row.names = FALSE)
 
 if (nrow(go_df) > 0) {
   message("\nTop GO terms:")
@@ -96,4 +109,4 @@ if (nrow(kegg_df) > 0) {
   print(head(kegg_df[, c("Description", "p.adjust", "Count")], 5), row.names = FALSE)
 }
 
-message("\nEnrichment analysis complete")
+message("\nEnrichment analysis complete")

tests/testthat/test-smoke.R

@@ -27,6 +27,10 @@ test_that("DESeq2 results CSV exists and is well-formed", {
   expect_gt(nrow(res), 0)
   expect_true(any(!is.na(res$padj)))
   expect_true(all(nchar(res$gene) > 0))
+
+  pvals <- res$pvalue[!is.na(res$pvalue)]
+  expect_true(length(pvals) > 0)
+  expect_true(all(pvals >= 0 & pvals <= 1))
 })
 
 test_that("key figures exist", {
@@ -36,3 +40,13 @@ test_that("key figures exist", {
   expect_true(file.exists(volcano))
   expect_true(file.exists(pca))
 })
+
+test_that("key derived tables exist", {
+  top_genes <- file.path(root, "results/tables/top_genes.csv")
+  go_terms <- file.path(root, "results/tables/go_biological_process.csv")
+  kegg_terms <- file.path(root, "results/tables/kegg_pathways.csv")
+
+  expect_true(file.exists(top_genes))
+  expect_true(file.exists(go_terms))
+  expect_true(file.exists(kegg_terms))
+})