###############################################################################
# concordance_analysis.R — MetaOmixTools method concordance
#
# Single script: runs all analyses, saves concordance_results.csv,
# and generates publication-quality TABLE figures (PDF + PNG).
#
# INPUT (place in same folder):
# data1_1.tsv data2_2.tsv data3_3.tsv data4_4.tsv <- RankProd
# gene1.txt gene2.txt gene3.txt gene4.txt <- RRA / MetaEnrich
#
# OUTPUT:
# concordance_results.csv <- full data table
# concordance_tableA.pdf/png <- MetaRank symmetric table
# concordance_tableB.pdf/png <- MetaEnrich symmetric table
# concordance_combined.pdf/png <- both panels side by side
###############################################################################
suppressPackageStartupMessages({
library(RankProd)
library(RobustRankAggreg)
library(clusterProfiler)
library(org.Hs.eg.db)
library(metap)
library(ggplot2)
library(dplyr)
library(patchwork)
})
set.seed(123)
OUTPUT_DIR <- "."
# ─────────────────────────────────────────────────────────────────────────────
# HELPERS
# ─────────────────────────────────────────────────────────────────────────────
jaccard <- function(a, b) {
a <- a[!is.na(a)]; b <- b[!is.na(b)]
u <- length(union(a, b))
if (u == 0) return(NA_real_)
round(length(intersect(a, b)) / u, 3)
}
build_matrix <- function(files) {
dfs <- lapply(files, function(f)
read.table(f, header = FALSE, sep = "\t",
col.names = c("Gene","Score"),
stringsAsFactors = FALSE, fill = TRUE))
genes <- unique(unlist(lapply(dfs, `[[`, "Gene")))
mat <- sapply(dfs, function(df) {
v <- setNames(as.numeric(df$Score), df$Gene); v[genes]
})
mat <- as.data.frame(mat); rownames(mat) <- genes
for (j in seq_len(ncol(mat)))
mat[[j]][is.na(mat[[j]])] <- mean(mat[[j]], na.rm = TRUE)
mat
}
run_rp <- function(mat) {
r <- RankProducts(mat, cl = rep(1, ncol(mat)), na.rm = TRUE, logged = TRUE,
gene.names = rownames(mat), plot = FALSE,
MinNumOfValidPairs = NA)
d <- data.frame(Gene = rownames(mat), Rank = r$RPrank[,1])
d[order(d$Rank), ]
}
run_rp_adv <- function(mat, origin) {
r <- RP.advance(mat, cl = rep(1, ncol(mat)), origin = origin,
na.rm = TRUE, logged = TRUE,
gene.names = rownames(mat), plot = FALSE)
d <- data.frame(Gene = rownames(mat), Rank = r$RPrank[,1])
d[order(d$Rank), ]
}
run_rra <- function(lists, method) {
r <- as.data.frame(RobustRankAggreg::aggregateRanks(
glist = lists, method = method, full = TRUE, exact = FALSE))
colnames(r)[1:2] <- c("Gene","Score")
r$Rank <- seq_len(nrow(r))
r[order(r$Rank), ]
}
comb_pval <- function(pv, method) {
pv <- pv[!is.na(pv) & pv > 0]
if (length(pv) < 2) return(NA_real_)
tryCatch(switch(method,
Fisher = metap::sumlog(pv)$p,
Stouffer = metap::sumz(pv)$p,
Tippett = metap::minimump(pv)$p,
Wilkinson = metap::wilkinsonp(pv)$p),
error = function(e) NA_real_)
}
# ─────────────────────────────────────────────────────────────────────────────
# SECTION 1 — MetaRank
# ─────────────────────────────────────────────────────────────────────────────
mat <- build_matrix(c("data1_1.tsv","data2_2.tsv","data3_3.tsv","data4_4.tsv"))
cat(sprintf(" Matrix: %d genes × %d datasets\n", nrow(mat), ncol(mat)))
cat(" Running RP ...\n"); res_rp <- run_rp(mat)
cat(" Running RP.advance ...\n"); res_rpa <- run_rp_adv(mat, c(1,2,1,2))
rra_lists <- lapply(c("gene1.txt","gene2.txt","gene3.txt","gene4.txt"),
function(f) { g <- readLines(f); trimws(g[nzchar(trimws(g))]) })
rra_methods <- c("RRA","min","median","mean","geom.mean","stuart")
cat(" Running 6 RRA variants ...\n")
rra_res <- setNames(lapply(rra_methods, function(m) run_rra(rra_lists, m)),
rra_methods)
top_ns <- c(50, 100, 200)
rows <- list()
# RP vs RP.advance
for (n in top_ns)
rows[[length(rows)+1]] <- data.frame(
Module="MetaRank", Group="RankProd",
Method_A="RP", Method_B="RP.advance", TopN=n,
Jaccard=jaccard(head(res_rp$Gene,n), head(res_rpa$Gene,n)))
# RP vs every RRA (cross-library)
for (m in rra_methods)
for (n in top_ns)
rows[[length(rows)+1]] <- data.frame(
Module="MetaRank", Group="RP vs RRA",
Method_A="RP", Method_B=m, TopN=n,
Jaccard=jaccard(head(res_rp$Gene,n), head(rra_res[[m]]$Gene,n)))
# ALL pairwise RRA (all 15 combinations)
cat(sprintf(" Computing all %d RRA-internal pairs ...\n",
choose(length(rra_methods), 2)))
for (pair in combn(rra_methods, 2, simplify=FALSE))
for (n in top_ns)
rows[[length(rows)+1]] <- data.frame(
Module="MetaRank", Group="RRA internal",
Method_A=pair[1], Method_B=pair[2], TopN=n,
Jaccard=jaccard(head(rra_res[[pair[1]]]$Gene,n),
head(rra_res[[pair[2]]]$Gene,n)))
mr_df <- do.call(rbind, rows)
# ─────────────────────────────────────────────────────────────────────────────
# SECTION 2 — MetaEnrich
# ─────────────────────────────────────────────────────────────────────────────
ora_full <- lapply(seq_along(rra_lists), function(i) {
genes <- head(rra_lists[[i]], 500)
cat(sprintf(" Dataset %d (%d genes) ... ", i, length(genes)))
res <- tryCatch(
suppressMessages(enrichGO(
gene=genes, OrgDb=org.Hs.eg.db, keyType="SYMBOL", ont="BP",
pvalueCutoff=1, pAdjustMethod="BH", readable=FALSE)),
error=function(e) NULL)
if (is.null(res)) { cat("no result\n"); return(NULL) }
df <- as.data.frame(res@result); df <- df[df$Count>0, ]
cat(sprintf("%d terms\n", nrow(df))); df
})
all_ids <- unique(unlist(lapply(ora_full, function(d)
if (!is.null(d)) d$ID else character(0))))
cat(sprintf(" Unique GO:BP terms: %d\n", length(all_ids)))
pval_mat <- sapply(ora_full, function(d) {
pv <- rep(1, length(all_ids)); names(pv) <- all_ids
if (!is.null(d)) pv[d$ID] <- d$pvalue; pv
}); rownames(pval_mat) <- all_ids
cat(" Combining p-values ...\n")
comb_methods <- c("Fisher","Stouffer","Tippett","Wilkinson")
comb_pvals <- setNames(lapply(comb_methods, function(m)
apply(pval_mat, 1, comb_pval, method=m)), comb_methods)
sig_terms <- lapply(comb_pvals, function(pv) {
adj <- p.adjust(pv,"BH"); names(adj)[!is.na(adj) & adj<0.05] })
for (m in comb_methods)
cat(sprintf(" %-12s : %d significant terms\n", m, length(sig_terms[[m]])))
me_rows <- list()
for (pair in combn(comb_methods, 2, simplify=FALSE))
me_rows[[length(me_rows)+1]] <- data.frame(
Module="MetaEnrich", Group="p-value combination",
Method_A=pair[1], Method_B=pair[2], TopN=NA_integer_,
Jaccard=jaccard(sig_terms[[pair[1]]], sig_terms[[pair[2]]]))
me_df <- do.call(rbind, me_rows)
# ─────────────────────────────────────────────────────────────────────────────
# SAVE CSV
# ─────────────────────────────────────────────────────────────────────────────
all_res <- rbind(mr_df, me_df)
write.csv(all_res, file.path(OUTPUT_DIR,"concordance_results.csv"), row.names=FALSE)
# ─────────────────────────────────────────────────────────────────────────────
# FIGURE HELPERS — symmetric coloured table via ggplot
# ─────────────────────────────────────────────────────────────────────────────
BLUE_LOW <- "#E6F1FB"
BLUE_MID <- "#378ADD"
BLUE_HIGH <- "#0C447C"
# Map Jaccard value → fill colour
jcol <- function(j) {
if (is.na(j)) return("grey93")
scales::col_numeric(c(BLUE_LOW, BLUE_MID, BLUE_HIGH), c(0,1))(j)
}
# Map Jaccard value → legible text colour
tcol <- function(j) if (!is.na(j) && j >= 0.55) "white" else "grey15"
# Build a symmetric matrix from a long data frame, then draw as a ggplot table
draw_table <- function(methods,
long_df,
title, subtitle,
note = NULL) {
n <- length(methods)
# Fill symmetric matrix
mat <- matrix(NA_real_, n, n, dimnames = list(methods, methods))
diag(mat) <- 1
for (i in seq_len(nrow(long_df))) {
a <- long_df$A[i]; b <- long_df$B[i]; j <- long_df$Jaccard[i]
if (a %in% methods && b %in% methods) {
mat[a, b] <- j; mat[b, a] <- j
}
}
# Long format for ggplot
cells <- expand.grid(row = seq_len(n), col = seq_len(n),
stringsAsFactors = FALSE)
cells$val <- mapply(function(r,c) mat[r,c], cells$row, cells$col)
cells$label <- ifelse(is.na(cells$val), "\u2014",
ifelse(cells$val==1 & cells$row==cells$col,
"1.00", sprintf("%.2f", cells$val)))
cells$fill <- sapply(cells$val, jcol)
cells$tcol <- sapply(cells$val, tcol)
col_labels <- methods
row_labels <- methods
p <- ggplot() +
# tiles
geom_tile(data = cells, aes(x = col, y = -row),
fill = cells$fill, color = "white", linewidth = 0.6,
width = 0.94, height = 0.84) +
# cell text
geom_text(data = cells, aes(x = col, y = -row, label = label),
color = cells$tcol, size = 3.1, family = "sans") +
# column headers
annotate("text", x = seq_len(n), y = 0.65,
label = col_labels, fontface = "bold",
size = 3.1, color = "grey15", hjust = 0.5) +
# row labels
annotate("text", x = 0.22, y = -seq_len(n),
label = row_labels, size = 3.1,
color = "grey15", hjust = 1) +
# axis titles (italic)
annotate("text", x = mean(seq_len(n)), y = 1.55,
label = "Method B", size = 3.0, fontface = "italic",
color = "grey45", hjust = 0.5) +
annotate("text", x = -0.65, y = -mean(seq_len(n)),
label = "Method A", size = 3.0, fontface = "italic",
color = "grey45", angle = 90, hjust = 0.5) +
scale_x_continuous(limits = c(-1.1, n + 0.7), expand = c(0,0)) +
scale_y_continuous(limits = c(-n - 0.65, 2.1), expand = c(0,0)) +
labs(title = title, subtitle = subtitle,
caption = if (!is.null(note)) note else NULL) +
theme_void(base_size = 11) +
theme(
plot.title = element_text(face = "bold", size = 12, hjust = 0,
margin = margin(b = 3)),
plot.subtitle = element_text(size = 9, color = "grey45", hjust = 0,
margin = margin(b = 10)),
plot.caption = element_text(size = 7.5, color = "grey50",
hjust = 0, margin = margin(t = 8)),
plot.margin = margin(8, 12, 8, 28)
)
p
}
# Colour scale legend as a thin ggplot strip
legend_strip <- function() {
d <- data.frame(x = seq(0, 1, length.out = 200), y = 1)
ggplot(d, aes(x = x, y = y, fill = x)) +
geom_tile(height = 0.8, width = 0.006, show.legend = FALSE) +
scale_fill_gradientn(colors = c(BLUE_LOW, BLUE_MID, BLUE_HIGH),
limits = c(0,1)) +
annotate("text", x = c(0, 0.25, 0.5, 0.75, 1),
y = 0.05, label = c("0","0.25","0.50","0.75","1.00"),
size = 2.7, color = "grey30", vjust = 1) +
annotate("text", x = 0.5, y = 1.9,
label = "Jaccard similarity",
size = 2.9, color = "grey30", hjust = 0.5) +
scale_x_continuous(limits = c(-0.05, 1.05), expand = c(0,0)) +
scale_y_continuous(limits = c(-0.3, 2.4), expand = c(0,0)) +
theme_void() +
theme(plot.margin = margin(2, 40, 2, 40))
}
# ─────────────────────────────────────────────────────────────────────────────
# TABLE A — MetaRank (top-100)
# ─────────────────────────────────────────────────────────────────────────────
mr100 <- all_res %>%
filter(Module == "MetaRank", TopN == 100) %>%
select(A = Method_A, B = Method_B, Jaccard)
# For the RP-vs-RRA cross-library rows the Method_B was stored as e.g. "RRA"
# (the run_rra method name), so it already matches the rra_methods labels.
# RP.advance has Method_B = "RP.advance" — keep as is.
all_mr_methods <- c("RP","RP.advance","RRA","min","median","mean","geom.mean","stuart")
short <- c("RP"="RP","RP.advance"="RP.adv","RRA"="RRA","min"="min",
"median"="median","mean"="mean","geom.mean"="geo.m","stuart"="Stuart")
mr100$A <- ifelse(mr100$A %in% names(short), short[mr100$A], mr100$A)
mr100$B <- ifelse(mr100$B %in% names(short), short[mr100$B], mr100$B)
display_methods <- unname(short[all_mr_methods])
fig_A <- draw_table(
methods = display_methods,
long_df = mr100,
title = "A MetaRank \u2014 pairwise ranking agreement",
subtitle = paste0("Jaccard similarity, top-100 genes | ",
"Grey cells = pair not applicable or self-comparison"),
note = paste("RP.advance uses origin vector (1,2,1,2) on four lung cancer",
"studies (GSE10082, GSE19188, GSE63459, GSE75037).")
)
# ─────────────────────────────────────────────────────────────────────────────
# TABLE B — MetaEnrich
# ─────────────────────────────────────────────────────────────────────────────
me_long <- me_df %>% select(A = Method_A, B = Method_B, Jaccard)
fig_B <- draw_table(
methods = comb_methods,
long_df = me_long,
title = "B MetaEnrich \u2014 term-set agreement",
subtitle = "Jaccard similarity on significant GO:BP terms (BH < 0.05)",
note = "ORA input: top-500 genes per ranked list; org.Hs.eg.db 3.23.1."
)
# ─────────────────────────────────────────────────────────────────────────────
# COMBINE
# ─────────────────────────────────────────────────────────────────────────────
leg <- legend_strip()
dataset_caption <- paste(
"Datasets: four lung cancer studies (GSE10082, GSE19188, GSE63459, GSE75037).",
"MetaRank: RankProd input 22k\u201355k genes/study; RRA input 13k\u201321k genes/list.",
"MetaEnrich: top-500 genes/list; ORA against GO:BP (org.Hs.eg.db 3.23.1).",
"All analyses: impute NA, BH correction, set.seed(123).",
sep = " "
)
combined <- (fig_A | fig_B) /
(leg + plot_spacer() + plot_layout(widths = c(1, 2.5))) +
plot_layout(heights = c(9, 1)) +
plot_annotation(
caption = dataset_caption,
theme = theme(
plot.caption = element_text(size = 7.5, color = "grey45",
lineheight = 1.35, hjust = 0,
margin = margin(t = 6)),
plot.margin = margin(6, 10, 6, 6)
)
)
# ─────────────────────────────────────────────────────────────────────────────
# SAVE
# ─────────────────────────────────────────────────────────────────────────────
save_fig <- function(plot, base, w, h) {
ggsave(file.path(OUTPUT_DIR, paste0(base,".pdf")),
plot, width=w, height=h, device=cairo_pdf)
ggsave(file.path(OUTPUT_DIR, paste0(base,".png")),
plot, width=w, height=h, dpi=300)
cat(sprintf(" %s.pdf / .png (%g x %g in)\n", base, w, h))
}
save_fig(fig_A, "concordance_tableA", 8.0, 5.8)
save_fig(fig_B, "concordance_tableB", 4.8, 3.8)
save_fig(combined, "concordance_combined", 13.0, 7.5)
