diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..93e92c0 --- /dev/null +++ b/.gitignore @@ -0,0 +1,13 @@ +# R +.Rhistory +.RData +.Rproj.user + +# macOS +.DS_Store + +# Large / regenerable outputs +results/**/*.csv.gz +results/**/*.pdf +diagnostic_output/** +priority_audit_output/** diff --git a/config.R.example b/config.R.example new file mode 100644 index 0000000..b4a55ec --- /dev/null +++ b/config.R.example @@ -0,0 +1,46 @@ +# ============================================================================ +# COMPASS configuration +# Copy this file to config.R and fill in the actual paths / values before +# running 00_preflight.R or 01_run_final_benchmark.R. +# +# cp config.R.example config.R +# +# All analysis parameters are defined here so that every script shares the +# same values without duplication. +# ============================================================================ + +# ---- Input paths ---- +COMPASS_CSV <- "data/human_normal_all.csv" +RNA_RDA <- "data/1.5-rna_final_score.rda" +CHIP_RDA <- "data/4-chip_final_score.rda" +RNA_FREQ_RDA <- "data/1.3-rna_freq.rda" + +GMT <- c( + PRC2 = "data/BENPORATH_PRC2_TARGETS.v2026.1.Hs.gmt", + SUZ12 = "data/BENPORATH_SUZ12_TARGETS.v2026.1.Hs.gmt", + EED = "data/BENPORATH_EED_TARGETS.v2026.1.Hs.gmt" +) + +# ---- Output directories ---- +OUT_DIR <- "results/benchmark_output" +DIAG_OUT_DIR <- "diagnostic_output" +AUDIT_OUT_DIR <- "priority_audit_output" + +# ---- Analysis parameters ---- +GENE_BIOTYPE <- "protein_coding" +EXCLUDE_TISSUE <- "esc" +TOP_K <- c(100L, 250L, 500L, 1000L) +THESIS_ROC_BOOT_N <- 1000L +JOURNAL_ROC_BOOT_N <- 2000L +N_BOOT <- 2000L +PRIMARY_REFERENCE <- "PRC2" +PRIMARY_K <- 500L +SEED <- 20260727L +WRITE_FULL_CONTEXT_TABLE <- FALSE +RUN_MODE <- "final_journal" + +# ---- Sanity-check thresholds (reproduced from thesis) ---- +EXPECTED_N <- 18642L +EXPECTED_POSITIVES <- c(PRC2 = 626L, SUZ12 = 998L, EED = 1017L) +EXPECTED_AUC <- c(PRC2 = 0.882, SUZ12 = 0.862, EED = 0.813) +AUC_TOLERANCE <- 0.0015 diff --git a/scripts/00_preflight.R b/scripts/00_preflight.R index 2af5408..a484d34 100644 --- a/scripts/00_preflight.R +++ b/scripts/00_preflight.R @@ -1,6 +1,12 @@ #!/usr/bin/env Rscript # Reproducibility preflight for final journal run. +if (!file.exists("config.R")) { + stop( + "Missing config.R. Copy config.R.example to config.R and fill in your ", + "local paths / parameters.\n cp config.R.example config.R" + ) +} source("config.R") required <- c( diff --git a/scripts/01_run_final_benchmark.R b/scripts/01_run_final_benchmark.R index 74d6e93..3f684bc 100644 --- a/scripts/01_run_final_benchmark.R +++ b/scripts/01_run_final_benchmark.R @@ -15,7 +15,14 @@ # check fails materially, baseline comparisons should NOT be interpreted yet. # ============================================================================ +if (!file.exists("config.R")) { + stop( + "Missing config.R. Copy config.R.example to config.R and fill in your ", + "local paths / parameters.\n cp config.R.example config.R" + ) +} source("config.R") +source("scripts/utils.R") required_pkgs <- c("data.table", "dplyr", "tidyr", "readr", "pROC", "ggplot2") missing_pkgs <- required_pkgs[ @@ -53,65 +60,9 @@ die_if_missing <- function(path, label = basename(path)) { if (!file.exists(path)) stop("Missing required file: ", label, "\nExpected: ", path) } -normalize_gene <- function(x) toupper(trimws(as.character(x))) - -normalize_species <- function(x) { - z <- tolower(trimws(as.character(x))) - z[z %in% c("hs", "hg38", "homo sapiens", "homo_sapiens")] <- "human" - z[z %in% c("mm", "mm10", "mm39", "mus musculus", "mus_musculus")] <- "mouse" - z -} - -normalize_status <- function(x) tolower(trimws(as.character(x))) - -normalize_tissue <- function(x) { - z <- tolower(trimws(as.character(x))) - z <- gsub("_", "-", z, fixed = TRUE) - z <- gsub("[[:space:]]+", "-", z) - z <- gsub("-+", "-", z) - z[z == "testis"] <- "testicle" - z[z == "head-and neck"] <- "head-and-neck" - z -} - -load_rda_list <- function(path) { - e <- new.env(parent = emptyenv()) - nm <- load(path, envir = e) - mget(nm, envir = e, inherits = FALSE) -} - -require_cols <- function(df, cols, label) { - miss <- setdiff(cols, names(df)) - if (length(miss) > 0) { - stop(label, " missing columns: ", paste(miss, collapse = ", ")) - } -} - -read_gmt_genes <- function(path) { - ln <- readLines(path, warn = FALSE) - ln <- ln[nchar(ln) > 0][1] - normalize_gene(strsplit(ln, "\t")[[1]][-(1:2)]) -} - -rna_conf_score <- function(x) { - dplyr::case_when( - tolower(as.character(x)) == "high" ~ 1.0, - tolower(as.character(x)) == "medium" ~ 0.4, - tolower(as.character(x)) == "low" ~ 0.2, - TRUE ~ 0 - ) -} - -chip_bind_score <- function(x) { - dplyr::case_when( - as.character(x) == "High confidence binding" ~ 1.0, - as.character(x) == "Medium-high confidence" ~ 0.8, - as.character(x) == "Medium confidence" ~ 0.4, - as.character(x) == "Low-medium confidence" ~ 0.3, - as.character(x) == "Low confidence" ~ 0.15, - TRUE ~ 0 - ) -} +# ============================================================================ +# Script-specific helpers +# ============================================================================ auc_ci_thesis <- function(labels, scores, n = THESIS_ROC_BOOT_N) { r <- pROC::roc( @@ -137,46 +88,41 @@ auc_ci_thesis <- function(labels, scores, n = THESIS_ROC_BOOT_N) { ) } -average_precision_grouped <- function(y, score) { - ok <- is.finite(score) & !is.na(y) - y <- as.integer(y[ok]) - score <- score[ok] - npos <- sum(y == 1) - if (npos == 0) return(NA_real_) - - o <- order(score, decreasing = TRUE) - y <- y[o] - score <- score[o] - - grp <- cumsum(c(TRUE, diff(score) != 0)) - pos_by <- as.numeric(tapply(y, grp, sum)) - n_by <- as.numeric(tapply(y, grp, length)) - - cum_pos <- cumsum(pos_by) - cum_n <- cumsum(n_by) - precision <- cum_pos / cum_n - recall <- cum_pos / npos - delta_recall <- c(recall[1], diff(recall)) - - sum(delta_recall * precision) +auc_ci_journal <- function(labels, scores, n = JOURNAL_ROC_BOOT_N) { + r <- pROC::roc( + labels, scores, + levels = c(0, 1), + direction = "<", + quiet = TRUE + ) + ci <- as.numeric( + pROC::ci.auc( + r, + method = "bootstrap", + boot.n = n, + boot.stratified = FALSE, + progress = "none" + ) + ) + list( + auc = as.numeric(pROC::auc(r)), + lo = ci[1], + hi = ci[3] + ) } topk_metrics <- function(gene, y, score, k) { kk <- min(as.integer(k), length(score)) - # deterministic tie-breaking only for Top-K extraction o <- order(-score, gene) idx <- o[seq_len(kk)] tp <- sum(y[idx] == 1) prevalence <- mean(y == 1) - precision <- tp / kk - recall <- tp / sum(y == 1) - enrichment <- precision / prevalence tibble( k = kk, true_positives = tp, - precision = precision, - recall = recall, - enrichment = enrichment + precision = tp / kk, + recall = tp / sum(y == 1), + enrichment = (tp / kk) / prevalence ) } @@ -614,7 +560,7 @@ context <- full_join( # Scale-robust rank baselines within tissue. MeanRank_missing0 = 0.50 * rna_rank + 0.50 * chip_rank, - RankProduct_missing0 = sqrt(pmax(rna_rank, 0) * pmax(chip_rank, 0)), + RankGeometricMean_missing0 = sqrt(pmax(rna_rank, 0) * pmax(chip_rank, 0)), # Very permissive single-best-evidence baseline. MaxEvidence = pmax(rna_zero, chip_zero), @@ -636,7 +582,7 @@ baseline_cols <- c( "EqualMean_missing0", "AvailableMean", "MeanRank_missing0", - "RankProduct_missing0", + "RankGeometricMean_missing0", "MaxEvidence", "StrictConcordance_proxy" ) diff --git a/scripts/02_run_diagnostic_2000.R b/scripts/02_run_diagnostic_2000.R index 92cb017..6a6c242 100644 --- a/scripts/02_run_diagnostic_2000.R +++ b/scripts/02_run_diagnostic_2000.R @@ -24,25 +24,19 @@ suppressPackageStartupMessages({ library(readr) }) -DATA_DIR <- "data" -BENCH_DIR <- "benchmark_output" -OUT_DIR <- "diagnostic_output" - -GMT <- c( - PRC2 = file.path(DATA_DIR, "BENPORATH_PRC2_TARGETS.v2026.1.Hs.gmt"), - SUZ12 = file.path(DATA_DIR, "BENPORATH_SUZ12_TARGETS.v2026.1.Hs.gmt"), - EED = file.path(DATA_DIR, "BENPORATH_EED_TARGETS.v2026.1.Hs.gmt") -) - -COMPASS_CSV <- file.path(DATA_DIR, "human_normal_all.csv") -CONTEXT_CSV <- file.path(BENCH_DIR, "04_context_level_baseline_inputs.csv.gz") -GENE_TSV <- file.path(BENCH_DIR, "03_gene_level_scores_COMPASS_and_baselines.tsv") +if (!file.exists("config.R")) { + stop( + "Missing config.R. Copy config.R.example to config.R and fill in your ", + "local paths / parameters.\n cp config.R.example config.R" + ) +} +source("config.R") +source("scripts/utils.R") -TOP_K <- c(100L, 250L, 500L, 1000L) -N_BOOT <- 2000L -SEED <- 20260727L +dir.create(DIAG_OUT_DIR, recursive = TRUE, showWarnings = FALSE) -dir.create(OUT_DIR, recursive = TRUE, showWarnings = FALSE) +CONTEXT_CSV <- file.path(OUT_DIR, "04_context_level_baseline_inputs.csv.gz") +GENE_TSV <- file.path(OUT_DIR, "03_gene_level_scores_COMPASS_and_baselines.tsv") must_exist <- c(COMPASS_CSV, CONTEXT_CSV, GENE_TSV, unname(GMT)) missing <- must_exist[!file.exists(must_exist)] @@ -50,24 +44,6 @@ if (length(missing) > 0) { stop("Missing required files:\n", paste(missing, collapse = "\n")) } -normalize_gene <- function(x) toupper(trimws(as.character(x))) - -normalize_tissue <- function(x) { - z <- tolower(trimws(as.character(x))) - z <- gsub("_", "-", z, fixed = TRUE) - z <- gsub("[[:space:]]+", "-", z) - z <- gsub("-+", "-", z) - z[z == "testis"] <- "testicle" - z[z == "head-and neck"] <- "head-and-neck" - z -} - -read_gmt_genes <- function(path) { - ln <- readLines(path, warn = FALSE) - ln <- ln[nchar(ln) > 0][1] - normalize_gene(strsplit(ln, "\t")[[1]][-(1:2)]) -} - auc_rank <- function(y, score) { ok <- is.finite(score) & !is.na(y) y <- as.integer(y[ok]) @@ -79,28 +55,6 @@ auc_rank <- function(y, score) { (sum(r[y == 1]) - n1 * (n1 + 1) / 2) / (n1 * n0) } -average_precision_grouped <- function(y, score) { - ok <- is.finite(score) & !is.na(y) - y <- as.integer(y[ok]) - score <- score[ok] - npos <- sum(y == 1) - if (npos == 0) return(NA_real_) - - o <- order(score, decreasing = TRUE) - y <- y[o] - score <- score[o] - - grp <- cumsum(c(TRUE, diff(score) != 0)) - pos_by <- as.numeric(tapply(y, grp, sum)) - n_by <- as.numeric(tapply(y, grp, length)) - cp <- cumsum(pos_by) - cn <- cumsum(n_by) - precision <- cp / cn - recall <- cp / npos - delta_recall <- c(recall[1], diff(recall)) - sum(delta_recall * precision) -} - calc_metrics <- function(df, methods, sets, analysis_name) { out <- list() top_out <- list() @@ -261,7 +215,7 @@ full_diag <- gene_scores %>% EqualMean_missing0, AvailableMean, MeanRank_missing0, - RankProduct_missing0, + RankGeometricMean_missing0, MaxEvidence, StrictConcordance_proxy ) %>% @@ -285,8 +239,8 @@ full_methods <- c( ) full_res <- calc_metrics(full_diag, full_methods, sets, "full_canonical_universe") -write_tsv(full_res$metrics, file.path(OUT_DIR, "01_full_metrics_with_single_modality.tsv")) -write_tsv(full_res$topk, file.path(OUT_DIR, "02_full_topK_with_single_modality.tsv")) +write_tsv(full_res$metrics, file.path(DIAG_OUT_DIR, "01_full_metrics_with_single_modality.tsv")) +write_tsv(full_res$topk, file.path(DIAG_OUT_DIR, "02_full_topK_with_single_modality.tsv")) # --------------------------------------------------------------------------- # B. Provenance of the winning tissue for each baseline @@ -340,7 +294,7 @@ for (m in prov_methods) { } provenance <- bind_rows(prov_rows) -write_tsv(provenance, file.path(OUT_DIR, "03_topK_winning_context_provenance.tsv")) +write_tsv(provenance, file.path(DIAG_OUT_DIR, "03_topK_winning_context_provenance.tsv")) # --------------------------------------------------------------------------- # C. Matched both-modality context benchmark @@ -382,8 +336,8 @@ matched_res <- calc_metrics( "matched_both_modalities" ) -write_tsv(matched_res$metrics, file.path(OUT_DIR, "04_matched_both_metrics.tsv")) -write_tsv(matched_res$topk, file.path(OUT_DIR, "05_matched_both_topK.tsv")) +write_tsv(matched_res$metrics, file.path(DIAG_OUT_DIR, "04_matched_both_metrics.tsv")) +write_tsv(matched_res$topk, file.path(DIAG_OUT_DIR, "05_matched_both_topK.tsv")) matched_summary <- tibble( item = c( @@ -401,7 +355,7 @@ matched_summary <- tibble( sum(matched_gene$gene %in% sets$EED) ) ) -write_tsv(matched_summary, file.path(OUT_DIR, "06_matched_both_summary.tsv")) +write_tsv(matched_summary, file.path(DIAG_OUT_DIR, "06_matched_both_summary.tsv")) # --------------------------------------------------------------------------- # D. Paired bootstrap for all three reference sets @@ -426,7 +380,7 @@ boot_full <- bootstrap_deltas( B = N_BOOT ) -write_tsv(boot_full, file.path(OUT_DIR, "07_bootstrap_all_sets_full.tsv")) +write_tsv(boot_full, file.path(DIAG_OUT_DIR, "07_bootstrap_all_sets_full.tsv")) matched_boot_methods <- c( "COMPASS", @@ -444,7 +398,7 @@ boot_matched <- bootstrap_deltas( B = N_BOOT ) -write_tsv(boot_matched, file.path(OUT_DIR, "08_bootstrap_all_sets_matched_both.tsv")) +write_tsv(boot_matched, file.path(DIAG_OUT_DIR, "08_bootstrap_all_sets_matched_both.tsv")) # --------------------------------------------------------------------------- # E. Compact interpretation table @@ -463,7 +417,7 @@ key <- full_res$metrics %>% select(reference_set, method, ROC_AUC, AP, AP_over_prevalence) %>% arrange(reference_set, desc(ROC_AUC)) -write_tsv(key, file.path(OUT_DIR, "09_key_diagnostic_comparison.tsv")) +write_tsv(key, file.path(DIAG_OUT_DIR, "09_key_diagnostic_comparison.tsv")) cat("\nDiagnostic v3 completed.\n") cat("Please return the entire diagnostic_output folder.\n") diff --git a/scripts/03_run_priority_audit.R b/scripts/03_run_priority_audit.R index 37a4b1b..136217b 100644 --- a/scripts/03_run_priority_audit.R +++ b/scripts/03_run_priority_audit.R @@ -21,30 +21,22 @@ suppressPackageStartupMessages({ library(readr) }) -DATA_DIR <- "data" -BENCH_DIR <- "benchmark_output" -OUT_DIR <- "priority_audit_output" +if (!file.exists("config.R")) { + stop( + "Missing config.R. Copy config.R.example to config.R and fill in your ", + "local paths / parameters.\n cp config.R.example config.R" + ) +} +source("config.R") +source("scripts/utils.R") -COMPASS_CSV <- file.path(DATA_DIR, "human_normal_all.csv") -CONTEXT_CSV <- file.path(BENCH_DIR, "04_context_level_baseline_inputs.csv.gz") -GENE_TSV <- file.path(BENCH_DIR, "03_gene_level_scores_COMPASS_and_baselines.tsv") -METRICS_TSV <- file.path(BENCH_DIR, "06_reference_metrics.tsv") +dir.create(AUDIT_OUT_DIR, recursive = TRUE, showWarnings = FALSE) -TOP_K <- c(100L, 250L, 500L, 1000L) EPS <- 1e-10 -dir.create(OUT_DIR, recursive = TRUE, showWarnings = FALSE) - -normalize_gene <- function(x) toupper(trimws(as.character(x))) -normalize_tissue <- function(x) { - z <- tolower(trimws(as.character(x))) - z <- gsub("_", "-", z, fixed = TRUE) - z <- gsub("[[:space:]]+", "-", z) - z <- gsub("-+", "-", z) - z[z == "testis"] <- "testicle" - z[z == "head-and neck"] <- "head-and-neck" - z -} +CONTEXT_CSV <- file.path(OUT_DIR, "04_context_level_baseline_inputs.csv.gz") +GENE_TSV <- file.path(OUT_DIR, "03_gene_level_scores_COMPASS_and_baselines.tsv") +METRICS_TSV <- file.path(OUT_DIR, "06_reference_metrics.tsv") required <- c(COMPASS_CSV, CONTEXT_CSV, GENE_TSV, METRICS_TSV) if (any(!file.exists(required))) { @@ -116,7 +108,7 @@ methods <- c( "EqualMean_missing0", "AvailableMean", "MeanRank_missing0", - "RankProduct_missing0", + "RankGeometricMean_missing0", "MaxEvidence", "StrictConcordance_proxy" ) @@ -211,7 +203,7 @@ for (m in methods) { k=k, top_genes=nrow(pg), winner_both_fraction=mean(pg$winner_any_both), - winner_single_only_fraction=mean(!pg$winner_any_both), + winner_no_both_fraction=mean(!pg$winner_any_both), winner_convergent_fraction=mean(pg$winner_any_convergent), winner_convergent_highRNA_fraction=mean(pg$winner_any_convergent_highRNA), metadata_join_failure_fraction=mean(pg$winning_context_n==0) @@ -227,8 +219,8 @@ for (m in methods) { priority_summary <- bind_rows(summary_rows) priority_detail <- bind_rows(detail_rows) -write_tsv(priority_summary, file.path(OUT_DIR,"01_topK_winning_context_quality.tsv")) -write_tsv(priority_detail, file.path(OUT_DIR,"02_top500_winning_context_detail.tsv")) +write_tsv(priority_summary, file.path(AUDIT_OUT_DIR,"01_topK_winning_context_quality.tsv")) +write_tsv(priority_detail, file.path(AUDIT_OUT_DIR,"02_top500_winning_context_detail.tsv")) # ------------------------ # "Any concordant context" (less strict than winning-context concordance) @@ -262,7 +254,7 @@ for (m in methods) { ) } } -write_tsv(bind_rows(any_rows), file.path(OUT_DIR,"03_topK_any_concordant_context.tsv")) +write_tsv(bind_rows(any_rows), file.path(AUDIT_OUT_DIR,"03_topK_any_concordant_context.tsv")) # ------------------------ # Top-K Jaccard overlap with COMPASS @@ -293,7 +285,7 @@ for (k in TOP_K) { ) } } -write_tsv(bind_rows(jrows), file.path(OUT_DIR,"04_topK_jaccard_vs_COMPASS.tsv")) +write_tsv(bind_rows(jrows), file.path(AUDIT_OUT_DIR,"04_topK_jaccard_vs_COMPASS.tsv")) # ------------------------ # Performance / evidence-quality trade-off table @@ -303,7 +295,7 @@ q500 <- priority_summary %>% select( method, winner_both_fraction, - winner_single_only_fraction, + winner_no_both_fraction, winner_convergent_fraction, winner_convergent_highRNA_fraction ) @@ -312,7 +304,7 @@ tradeoff <- metrics %>% inner_join(q500, by="method") %>% arrange(reference_set, desc(ROC_AUC)) -write_tsv(tradeoff, file.path(OUT_DIR,"05_performance_vs_evidence_quality.tsv")) +write_tsv(tradeoff, file.path(AUDIT_OUT_DIR,"05_performance_vs_evidence_quality.tsv")) cat("\nPriority audit completed.\n") cat("Return the entire priority_audit_output folder.\n") diff --git a/scripts/utils.R b/scripts/utils.R new file mode 100644 index 0000000..19cd46f --- /dev/null +++ b/scripts/utils.R @@ -0,0 +1,79 @@ +# ============================================================================ +# COMPASS shared utility functions +# Sourced by all analysis scripts to avoid duplication. +# ============================================================================ + +normalize_gene <- function(x) toupper(trimws(as.character(x))) + +normalize_tissue <- function(x) { + z <- tolower(trimws(as.character(x))) + z <- gsub("_", "-", z, fixed = TRUE) + z <- gsub("[[:space:]]+", "-", z) + z <- gsub("-+", "-", z) + z[z == "testis"] <- "testicle" + z[z == "head-and neck"] <- "head-and-neck" + z +} + +read_gmt_genes <- function(path) { + ln <- readLines(path, warn = FALSE) + ln <- ln[nchar(ln) > 0][1] + normalize_gene(strsplit(ln, "\t")[[1]][-(1:2)]) +} + +load_rda_list <- function(path) { + e <- new.env(parent = emptyenv()) + nm <- load(path, envir = e) + mget(nm, envir = e, inherits = FALSE) +} + +average_precision_grouped <- function(y, score) { + ok <- is.finite(score) & !is.na(y) + y <- as.integer(y[ok]) + score <- score[ok] + npos <- sum(y == 1) + if (npos == 0) return(NA_real_) + + o <- order(score, decreasing = TRUE) + y <- y[o] + score <- score[o] + + grp <- cumsum(c(TRUE, diff(score) != 0)) + pos_by <- as.numeric(tapply(y, grp, sum)) + n_by <- as.numeric(tapply(y, grp, length)) + + cum_pos <- cumsum(pos_by) + cum_n <- cumsum(n_by) + precision <- cum_pos / cum_n + recall <- cum_pos / npos + delta_recall <- c(recall[1], diff(recall)) + + sum(delta_recall * precision) +} + +require_cols <- function(df, cols, label) { + miss <- setdiff(cols, names(df)) + if (length(miss) > 0) { + stop(label, " missing columns: ", paste(miss, collapse = ", ")) + } +} + +rna_conf_score <- function(x) { + dplyr::case_when( + tolower(as.character(x)) == "high" ~ 1.0, + tolower(as.character(x)) == "medium" ~ 0.4, + tolower(as.character(x)) == "low" ~ 0.2, + TRUE ~ 0 + ) +} + +chip_bind_score <- function(x) { + dplyr::case_when( + as.character(x) == "High confidence binding" ~ 1.0, + as.character(x) == "Medium-high confidence" ~ 0.8, + as.character(x) == "Medium confidence" ~ 0.4, + as.character(x) == "Low-medium confidence" ~ 0.3, + as.character(x) == "Low confidence" ~ 0.15, + TRUE ~ 0 + ) +}