OSCR

Integrated Multi-Tissue Transcriptomics Reveals Antagonistic Pleiotropy in Aging and Alzheimer's Disease.

Code ↔ Paper

18 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 18 matches
  1. [1] § Methods › WGCNA analysis and AP module/gene scoring ↔ 10_modules_AP.R, lines 316–342 · score 0.87 · TOMType, mergeCutHeight, minModuleSize, soft thresholding power, networkType, WGCNA
  2. [2] § Methods › Functional enrichment of AP-derived gene sets ↔ 09_go_kegg_enrichment.R, lines 67–121 · score 0.81 · clusterProfiler, p.adjust, hsa, simplify, CC, MF
  3. [3] § Methods › Cell-state enrichment analysis and cell type-adjusted differential expression ↔ 07_optional_gsva_deconv.R, lines 209–229 · score 0.81 · C8 brain cell, Brain cell state, neuron, pericyte, endothelial, excitatory
  4. [4] § Methods › Functional enrichment of AP-derived gene sets ↔ withDevelopment_section.R, lines 1318–1444 · score 0.81 · Entrez IDs, clusterProfiler, hsa, CC, MF, BP
  5. [5] § Methods › Pathway and TF analyses ↔ 05_tf_activity.R, lines 290–345 · score 0.81 · TF activities, AP overlap, AP weighted, DoRothEA, TFs, regulons
  6. [6] § Methods › WGCNA analysis and AP module/gene scoring ↔ 10_modules_AP.R, lines 316–342 · score 0.79 · merge cut height, soft thresholding power, module AP, module gene, WGCNA, network
  7. [7] § Methods › Gene classification, prioritization, and AP metrics ↔ 02_consensus_meta.R, lines 89–158 · score 0.75 · divergence score, discordant Age, shared score, log2FC, log10, consensus
  8. [8] § Methods › Pathway and TF analyses ↔ 03_pathway_antagonism.R, lines 16–39 · score 0.69 · fgseaMultilevel, fgseaSimple, permutations, NES, Pathway, antagonism
  9. [9] § Methods › WGCNA analysis and AP module/gene scoring ↔ 10_modules_AP.R, lines 807–845 · score 0.66 · ego networks, scoring genes, AP class, AP scores, hub, exported
  10. [10] § Methods › Covariate and brain-only sensitivity audits ↔ 04_ap_candidates.R, lines 665–718 · score 0.60 · AD brain meta, aging meta, AP candidates, filtered, prioritization, class
  11. [11] § Methods › Cross-condition integration and meta-analysis ↔ 01_prepare_universe.R, lines 71–122 · score 0.58 · HGNC symbols, gene symbol, log2FC, Ensembl
  12. [12] § Methods › Gene classification, prioritization, and AP metrics ↔ 02_consensus_meta.R, lines 89–158 · score 0.58 · Aging pleiotropy, aging contrasts, log2FC, consensus, prioritization, AD
  13. [13] § Methods › Cross-condition integration and meta-analysis ↔ 04_ap_candidates.R, lines 156–189 · score 0.57 · HGNC symbols, gene symbol, log2fc, Ensembl, AP
  14. [14] § Results › TF activity antagonism between healthy aging and AD ↔ 05_tf_activity.R, lines 290–345 · score 0.55 · TF activity, DoRothEA, regulons, NES, map, antagonism
  15. [15] § Results › Cell state-adjusted robustness of AD-associated AP signals ↔ 07_optional_gsva_deconv.R, lines 209–229 · score 0.53 · Cell state, endothelial, excitatory, astrocyte, inhibitory, oligodendrocyte
  16. [16] § Results › Integrated hallmark pathway across brain aging and AD ↔ utils.R, lines 109–117 · score 0.53 · SharedDown, SharedUp, AP Vulnerability, AP Resilience, quadrant, pathways
  17. [17] § Results › Integrated hallmark pathway across brain aging and AD ↔ 05_tf_activity.R, lines 147–155 · score 0.52 · SharedDown, SharedUp, AP Vulnerability, AP Resilience, quadrant, map
  18. [18] § Methods › WGCNA analysis and AP module/gene scoring ↔ 10_modules_AP.R, lines 1–60 · score 0.51 · aging samples, finite, young, GSE48350, WGCNA, Weighted

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

R · 1,242 lines · 59 KB · no license · 4 matches

  1. #!/usr/bin/env Rscript
  2. # =============================================================================
  3. # 10_modules_AP_consensus.R — UPDATED FULL SCRIPT (robust covariates)
  4. # Uses ONLY:
  5. # - results_aging_brain/GSE48350_samples_ALL_tissues_Young_vs_Old.csv
  6. # (all brain healthy aging samples; includes young + old; 150 total)
  7. # - results_AD/samples_ALL_tissues_ageGT59.csv (take only rows with group == "AD")
  8. # for sample annotation.
  9. # =============================================================================
  11. suppressPackageStartupMessages({
  12. library(dplyr); library(readr); library(stringr); library(tidyr)
  13. library(purrr); library(tibble); library(ggplot2); library(scales)
  14. library(WGCNA)
  15. suppressWarnings({
  16. ok_igraph <- requireNamespace("igraph", quietly=TRUE)
  17. ok_ggraph <- requireNamespace("ggraph", quietly=TRUE) && requireNamespace("tidygraph", quietly=TRUE)
  18. })
  19. })
  21. set.seed(2371)
  22. allowWGCNAThreads()
  23. options(stringsAsFactors = FALSE)
  25. # ------------------ Config / IO ----------------------------------------------
  26. source("utils.R")
  27. load(".config_env.RData") # expects: AGE_BRAIN_EXPR_RDS, AD_BRAIN_EXPR_RDS, outdir, AP_COLORS, PNG_DPI, REGION_P_FDR
  29. dir.create(file.path(outdir, "wgcna"), FALSE, TRUE)
  30. .dirnets <- file.path(outdir, "wgcna", "networks"); dir.create(.dirnets, FALSE, TRUE)
  32. # CSVs (override via env)
  33. AD_SAMPLES_CSV <- Sys.getenv("AD_SAMPLES_CSV",
  34. unset = file.path("results_AD","samples_ALL_tissues_ageGT59.csv"))
  35. HEALTHY_SAMPLES_CSV <- Sys.getenv("HEALTHY_SAMPLES_CSV",
  36. unset = file.path("results_aging_brain","GSE48350_samples_ALL_tissues_Young_vs_Old.csv"))
  38. # Tunables
  39. CONSENSUS_TOPK_GENES <- as.integer(Sys.getenv("CONSENSUS_TOPK_GENES", unset="15000"))
  40. MIN_FINITE_PROP <- as.numeric(Sys.getenv("MIN_FINITE_PROP", unset="0.80"))
  41. MIN_FINITE_MIN <- as.integer(Sys.getenv("MIN_FINITE_MIN", unset="12"))
  43. NETWORK_MAX_GENES <- as.integer(Sys.getenv("NETWORK_MAX_GENES", unset="400"))
  44. EDGE_TOPK_PER_GENE <- as.integer(Sys.getenv("EDGE_TOPK_PER_GENE", unset="5"))
  45. EDGE_MIN_TOM <- as.numeric(Sys.getenv("EDGE_MIN_TOM", unset="0.05"))
  46. EDGE_MIN_ADJ <- as.numeric(Sys.getenv("EDGE_MIN_ADJ", unset="0.10"))
  47. TOP_HUBS_PER_MODULE <- as.integer(Sys.getenv("TOP_HUBS_PER_MODULE", unset="10"))
  48. MIN_MODULE_SIZE <- as.integer(Sys.getenv("MIN_MODULE_SIZE", unset="120"))
  49. MIN_EDGES_TARGET <- as.integer(Sys.getenv("MIN_EDGES_TARGET", unset="40"))
  50. MIN_EDGE_MIN_TOM <- as.numeric(Sys.getenv("MIN_EDGE_MIN_TOM", unset="0.01"))
  51. MIN_EDGE_MIN_ADJ <- as.numeric(Sys.getenv("MIN_EDGE_MIN_ADJ", unset="0.02"))
  52. FORCE_EXPORT_IF_EMPTY<- as.logical(Sys.getenv("FORCE_EXPORT_IF_EMPTY",unset="TRUE"))
  53. TOP_EXPORT <- as.integer(Sys.getenv("AP_TOP_MODULE_EXPORT", unset="6"))
  54. TOP_EGO_PER_CLASS <- as.integer(Sys.getenv("TOP_EGO_PER_CLASS", unset="15"))
  56. # Optional weightings
  57. AP_SCORE_USE_KME <- as.logical(Sys.getenv("AP_SCORE_USE_KME", unset="FALSE"))
  58. AP_MODULE_USE_KME_WEIGHTS <- as.logical(Sys.getenv("AP_MODULE_USE_KME_WEIGHTS", unset="TRUE"))
  60. `%||%` <- function(a, b) if (!is.null(a) && length(a)>0 && !all(is.na(a))) a else b
  61. log_msg <- function(...) cat(format(Sys.time(), "[%Y-%m-%d %H:%M:%S]"), sprintf(...), "\n")
  62. psafe <- function(p) pmax(pmin(as.numeric(p), 1), .Machine$double.xmin)
  64. stopifnot(file.exists(AGE_BRAIN_EXPR_RDS), file.exists(AD_BRAIN_EXPR_RDS))
  65. stopifnot(file.exists(AD_SAMPLES_CSV), file.exists(HEALTHY_SAMPLES_CSV))
  67. # ------------------ Helpers: IO & metadata -----------------------------------
  68. read_samples_table <- function(path){
  69. tryCatch(readr::read_csv(path, show_col_types=FALSE),
  70. error=function(e) tryCatch(readr::read_tsv(path, show_col_types=FALSE),
  71. error=function(e2) readr::read_delim(path, delim="\t", show_col_types=FALSE)))
  72. }
  74. load_expr_meta <- function(path) {
  75. obj <- readRDS(path)
  76. if (is.list(obj) && all(c("expr","meta") %in% names(obj))) {
  77. expr <- as.matrix(obj$expr); meta <- as.data.frame(obj$meta)
  78. } else if (is.matrix(obj) || is.data.frame(obj)) {
  79. expr <- as.matrix(obj); meta <- data.frame(sample_id=colnames(expr))
  80. } else stop("Unsupported RDS format at ", path)
  81. if (is.null(rownames(expr))) stop("Expression must have rownames=gene identifiers.")
  82. if (!"sample_id" %in% names(meta)) meta$sample_id <- rownames(meta) %||% colnames(expr)
  83. meta <- meta %>% distinct(sample_id, .keep_all=TRUE)
  84. keep <- intersect(colnames(expr), meta$sample_id)
  85. list(expr=expr[, keep, drop=FALSE], meta=meta[match(keep, meta$sample_id), , drop=FALSE])
  86. }
  88. coerce_meta <- function(meta, context=c("healthy","ad")) {
  89. context <- match.arg(context)
  90. nm <- tolower(names(meta))
  91. pick <- function(cands) { i <- which(nm %in% cands); if (length(i)) names(meta)[i[1]] else NA }
  92. age_col <- pick(c("age","age_years","ages","age_at_death","age_at_sampling"))
  93. dz_col <- pick(c("disease","diagnosis","dx","ad_status","group","phenotype"))
  94. sex_col <- pick(c("sex","gender"))
  95. bat_col <- pick(c("batch","plate","run","study_batch"))
  96. reg_col <- pick(c("region","brain_region","tissue","area","region_raw"))
  97. meta$age <- if (!is.na(age_col)) suppressWarnings(as.numeric(meta[[age_col]])) else NA_real_
  98. if (!is.na(dz_col)) {
  99. v <- as.character(meta[[dz_col]])
  100. is_ad <- grepl("\\b(ad|alz|alzheimer|case|patient|disease)\\b", v, TRUE)
  101. is_ctrl <- grepl("\\b(ctrl|control|cn|healthy|normal)\\b", v, TRUE)
  102. meta$disease <- ifelse(is_ad, 1L, ifelse(is_ctrl, 0L, NA_integer_))
  103. } else meta$disease <- if (context=="ad") NA_integer_ else 0L
  104. meta$sex <- if (!is.na(sex_col)) as.character(meta[[sex_col]]) else NA_character_
  105. meta$batch <- if (!is.na(bat_col)) as.character(meta[[bat_col]]) else NA_character_
  106. meta$region <- if (!is.na(reg_col)) as.character(meta[[reg_col]]) else NA_character_
  107. meta$context <- context
  108. meta
  109. }
  111. merge_external_info <- function(meta, csv_path, dataset_label=""){
  112. if (!file.exists(csv_path)) return(meta)
  113. ext <- read_samples_table(csv_path)
  114. nm <- tolower(names(ext))
  115. pick <- function(...) { cands <- c(...); i <- which(nm %in% tolower(cands)); if (length(i)) names(ext)[i[1]] else NA }
  116. idc <- pick("sample_id","gsm","id","sample","geo_accession")
  117. if (is.na(idc)) { log_msg("[%s] No sample_id-like column in %s; skipping merge.", dataset_label, csv_path); return(meta) }
  118. ext <- ext %>% mutate(sample_id = as.character(.data[[idc]]))
  119. agec <- pick("age","age_years","Age","AgeYears")
  120. sexc <- pick("sex","gender")
  121. regc <- pick("region","brain_region","tissue","area","region_raw")
  122. tisc <- pick("tissue","Tissue")
  123. aggc <- pick("age_group","Age_group","AGE_GROUP")
  125. if (!is.na(agec)) ext$age <- suppressWarnings(as.numeric(ext[[agec]]))
  126. if (!is.na(sexc)) ext$sex <- as.character(ext[[sexc]])
  127. if (!is.na(regc)) ext$region <- as.character(ext[[regc]])
  128. if (!is.na(tisc)) ext$tissue <- as.character(ext[[tisc]])
  129. if (!is.na(aggc)) ext$age_group <- as.character(ext[[aggc]])
  131. meta2 <- meta %>%
  132. left_join(ext %>% select(sample_id, any_of(c("age","sex","region","tissue","age_group"))), by="sample_id")
  133. for (col in c("age","sex","region","tissue","age_group")) {
  134. new <- paste0(col, ".y"); old <- paste0(col, ".x")
  135. if (new %in% names(meta2)) { meta2[[col]] <- meta2[[new]] %||% meta2[[old]]; meta2[[new]] <- NULL; meta2[[old]] <- NULL }
  136. }
  137. meta2
  138. }
  140. # ------------------ Cleaning & math ------------------------------------------
  141. harmonize_gene_names <- function(E) {
  142. rn <- rownames(E)
  143. rn <- gsub("\\s+", "", rn)
  144. rn <- gsub("\\.\\d+$", "", rn)
  145. rownames(E) <- toupper(rn)
  146. E
  147. }
  148. dedup_by_maxvar <- function(E) {
  149. rn <- rownames(E)
  150. if (anyDuplicated(rn)) {
  151. idxs <- split(seq_along(rn), rn)
  152. pick <- vapply(idxs, function(ix){ vv <- apply(E[ix,,drop=FALSE], 1, stats::var, na.rm=TRUE); ix[which.max(vv)][1] }, integer(1))
  153. E <- E[pick,,drop=FALSE]; rownames(E) <- names(pick)
  154. }
  155. E
  156. }
  157. median_impute_rows <- function(E) {
  158. for (i in seq_len(nrow(E))) {
  159. v <- E[i,]; ok <- is.finite(v)
  160. if (!all(ok)) { if (any(ok)) { med <- stats::median(v[ok]); v[!ok] <- med; E[i,] <- v } else { E[i,] <- 0 } }
  161. }
  162. E
  163. }
  164. clean_expression <- function(E, name="") {
  165. log_msg("[%s] raw: %s genes x %s samples", name, nrow(E), ncol(E))
  166. E <- harmonize_gene_names(E) |> dedup_by_maxvar()
  167. finite_counts <- rowSums(is.finite(E))
  168. min_needed <- max(MIN_FINITE_MIN, ceiling(MIN_FINITE_PROP * ncol(E)))
  169. keep <- finite_counts >= min_needed
  170. if (!any(keep)) stop(sprintf("[%s] No genes pass finite filter.", name))
  171. E <- E[keep,,drop=FALSE] |> median_impute_rows()
  172. v <- apply(E, 1, stats::var)
  173. E <- E[v > 0,, drop=FALSE]
  174. log_msg("[%s] cleaned: %s genes x %s samples (min finite per gene %s)", name, nrow(E), ncol(E), min_needed)
  175. E
  176. }
  177. sel_top_var_consensus <- function(E1, E2, k=15000, tag="CONS") {
  178. v1 <- apply(E1, 1, stats::var); v2 <- apply(E2, 1, stats::var)
  179. r1 <- rank(-v1, ties.method="average"); r2 <- rank(-v2, ties.method="average")
  180. rs <- r1 + r2; k <- min(k, length(rs))
  181. keep <- names(sort(rs))[seq_len(k)]
  182. log_msg("[%s] choosing top %s / %s shared genes by consensus variance", tag, k, length(rs))
  183. list(E1=E1[keep,,drop=FALSE], E2=E2[keep,,drop=FALSE], genes=keep)
  184. }
  185. std_me_names <- function(x) paste0("ME", sprintf("%02d", as.integer(gsub("^ME","", as.character(x)))))
  187. # ===== NEW: robust covariate utilities =======================================
  189. sanitize_covars <- function(covs, sample_order=NULL) {
  190. if (is.null(covs) || ncol(covs)==0) return(NULL)
  191. covs <- as.data.frame(covs, stringsAsFactors = FALSE)
  192. if (!is.null(sample_order)) {
  193. covs <- covs[match(sample_order, rownames(covs)), , drop=FALSE]
  194. }
  195. for (j in names(covs)) {
  196. v <- covs[[j]]
  197. if (is.character(v) || is.logical(v)) covs[[j]] <- droplevels(factor(v))
  198. }
  199. keep <- vapply(covs, function(col) any(is.finite(col) | !is.na(col)), logical(1))
  200. covs <- covs[, keep, drop=FALSE]
  201. if (ncol(covs)==0) return(NULL)
  202. drop_j <- logical(ncol(covs))
  203. for (j in seq_along(covs)) {
  204. v <- covs[[j]]
  205. if (is.factor(v)) {
  206. v <- droplevels(v)
  207. if (nlevels(v) < 2) drop_j[j] <- TRUE else covs[[j]] <- v
  208. } else if (is.numeric(v)) {
  209. if (all(!is.finite(v)) || stats::var(v[is.finite(v)])==0) drop_j[j] <- TRUE
  210. }
  211. }
  212. covs <- covs[, !drop_j, drop=FALSE]
  213. if (ncol(covs)==0) return(NULL)
  214. covs
  215. }
  217. drop_collinear <- function(y, covs, thr=0.98) {
  218. if (is.null(covs) || ncol(covs)==0) return(NULL)
  219. mm <- tryCatch(model.matrix(~ . , data=cbind.data.frame(covs)), error=function(e) NULL)
  220. if (is.null(mm)) return(NULL)
  221. mm <- mm[, colnames(mm)!="(Intercept)", drop=FALSE]
  222. if (ncol(mm) <= 1) return(as.data.frame(covs))
  223. C <- suppressWarnings(cor(mm, use="pairwise.complete.obs"))
  224. diag(C) <- 0
  225. to_drop <- c()
  226. while (TRUE) {
  227. mx <- suppressWarnings(max(abs(C), na.rm=TRUE))
  228. if (!is.finite(mx) || mx < thr) break
  229. ij <- which(abs(C) == mx, arr.ind=TRUE)[1,]
  230. jdrop <- colnames(C)[ij[2]]
  231. to_drop <- c(to_drop, jdrop)
  232. keep <- setdiff(colnames(C), to_drop)
  233. if (length(keep) < 2) break
  234. C <- C[keep, keep, drop=FALSE]
  235. }
  236. if (length(to_drop)) {
  237. keep <- setdiff(colnames(mm), to_drop)
  238. return(as.data.frame(covs)) # keep sanitized original; mm used only for checks
  239. }
  240. as.data.frame(covs)
  241. }
  243. # ===== Associations ===========================================================
  244. assoc_linear <- function(y, x, covars=NULL, min_n=4) {
  245. df <- data.frame(y=y, x=x)
  246. if (!is.null(covars)) df <- cbind(df, covars)
  247. df <- df[complete.cases(df), ]
  248. if (nrow(df) < min_n || sd(df$x, na.rm=TRUE) == 0) return(c(beta=NA, se=NA, p=NA, n=nrow(df)))
  249. fit <- tryCatch(lm(y ~ x + ., df), error=function(e) NULL)
  250. if (is.null(fit)) return(c(beta=NA, se=NA, p=NA, n=nrow(df)))
  251. s <- summary(fit)$coefficients
  252. if (!"x" %in% rownames(s)) return(c(beta=NA, se=NA, p=NA, n=nrow(df)))
  253. c(beta=unname(s["x","Estimate"]), se=unname(s["x","Std. Error"]), p=unname(s["x","Pr(>|t|)"]), n=nrow(df))
  254. }
  255. assoc_binary <- function(y, x, covars=NULL, min_n=6) {
  256. df <- data.frame(y=y, x=x)
  257. if (!is.null(covars)) df <- cbind(df, covars)
  258. df <- df[complete.cases(df), ]
  259. ux <- unique(df$x); ux <- ux[is.finite(ux)]
  260. if (length(ux) < 2 || nrow(df) < min_n) return(c(beta=NA, se=NA, p=NA, n=nrow(df)))
  261. fit <- tryCatch(lm(y ~ x + ., df), error=function(e) NULL)
  262. if (is.null(fit)) return(c(beta=NA, se=NA, p=NA, n=nrow(df)))
  263. s <- summary(fit)$coefficients
  264. if (!"x" %in% rownames(s)) return(c(beta=NA, se=NA, p=NA, n=nrow(df)))
  265. c(beta=unname(s["x","Estimate"]), se=unname(s["x","Std. Error"]), p=unname(s["x","Pr(>|t|)"]), n=unname(nrow(df)))
  266. }
  267. fisher_combine_p <- function(pvals){ pvals <- psafe(pvals); 1 - pchisq(-2*sum(log(pvals)), df=2*length(pvals)) }
  269. # ------------------ Load datasets --------------------------------------------
  270. age_ds <- load_expr_meta(AGE_BRAIN_EXPR_RDS); age_ds$meta <- coerce_meta(age_ds$meta, "healthy")
  271. ad_ds <- load_expr_meta(AD_BRAIN_EXPR_RDS); ad_ds$meta <- coerce_meta(ad_ds$meta, "ad")
  273. # Merge external info
  274. # Healthy aging annotations from GSE48350_samples_ALL_tissues_Young_vs_Old.csv (brain)
  275. age_ds$meta <- merge_external_info(age_ds$meta, HEALTHY_SAMPLES_CSV, "GSE48350_healthy")
  276. # AD annotations from samples_ALL_tissues_ageGT59.csv (filter to AD rows only)
  277. ad_ds$meta <- merge_external_info(ad_ds$meta, AD_SAMPLES_CSV, "AD>59")
  279. # AD ids from CSV — take only group == "AD" if 'group' column exists; else assume all rows are AD
  280. ad_csv_raw <- read_samples_table(AD_SAMPLES_CSV)
  281. stopifnot("sample_id" %in% names(ad_csv_raw))
  282. grp_col <- names(ad_csv_raw)[tolower(names(ad_csv_raw))=="group"] %||% NA
  283. ad_csv <- if (is.na(grp_col)) {
  284. log_msg("Note: 'group' not found in AD CSV; assuming all rows are AD.")
  285. ad_csv_raw
  286. } else {
  287. dplyr::filter(ad_csv_raw, grepl("^\\s*AD\\s*$", .data[[grp_col]], ignore.case = TRUE))
  288. }
  289. stopifnot(nrow(ad_csv)>0)
  290. ad_ids <- intersect(colnames(ad_ds$expr), as.character(ad_csv$sample_id))
  291. stopifnot(length(ad_ids)>0)
  293. # ------------------ Healthy-only for CONSENSUS -------------------------------
  294. # Healthy aging cohort from GSE48350 CSV (all brain healthy samples; young + old)
  295. expr_age_healthy <- age_ds$expr[, intersect(colnames(age_ds$expr),
  296. age_ds$meta$sample_id[(age_ds$meta$disease %||% 0) == 0]), drop=FALSE]
  298. # Controls from AD dataset (explicitly exclude AD ids)
  299. expr_ctrl_healthy <- ad_ds$expr[, intersect(colnames(ad_ds$expr),
  300. ad_ds$meta$sample_id[(ad_ds$meta$disease %||% 0) == 0 &
  301. !(ad_ds$meta$sample_id %in% ad_ids)]), drop=FALSE]
  303. EA0 <- clean_expression(expr_age_healthy, "GSE48350_healthy")
  304. EC0 <- clean_expression(expr_ctrl_healthy, "GSE48350_controls")
  306. # Shared gene universe and top-K consensus
  307. genes_shared <- intersect(rownames(EA0), rownames(EC0))
  308. EA_shared <- EA0[genes_shared,,drop=FALSE]
  309. EC_shared <- EC0[genes_shared,,drop=FALSE]
  310. cons <- sel_top_var_consensus(EA_shared, EC_shared, k=CONSENSUS_TOPK_GENES, tag="CONS")
  311. EA <- cons$E1; EC <- cons$E2; genes_cons <- cons$genes
  313. # WGCNA consensus on healthy brains
  314. DatA <- as.data.frame(t(EA)); DatC <- as.data.frame(t(EC))
  315. pick_power <- function(X){ sft <- tryCatch(pickSoftThreshold(X, networkType="signed", verbose=0), error=function(e) NULL); pe <- if (!is.null(sft)) sft$powerEstimate else NA_real_; if (!is.finite(pe)) 6 else pe }
  316. pA <- pick_power(DatA); pC <- pick_power(DatC)
  317. cons_power <- floor(median(c(pA,pC), na.rm=TRUE)); if (!is.finite(cons_power)) cons_power <- 6
  318. log_msg("Consensus soft-threshold power: %s (A=%s, C=%s) with %s genes", cons_power, pA, pC, length(genes_cons))
  320. multiExpr <- list(HealthyA= list(data=DatA), HealthyC= list(data=DatC))
  321. net_cons <- tryCatch(
  322. blockwiseConsensusModules(
  323. multiExpr,
  324. power=cons_power, networkType="signed", TOMType="signed",
  325. minModuleSize=MIN_MODULE_SIZE, reassignThreshold=0,
  326. mergeCutHeight=0.25, numericLabels=TRUE, pamRespectsDendro=FALSE, verbose=2
  327. ),
  328. error=function(e){
  329. log_msg("Consensus failed: %s. Falling back to single-dataset (HealthyA).", conditionMessage(e))
  330. single <- blockwiseModules(DatA, power=cons_power, networkType="signed", TOMType="signed",
  331. minModuleSize=MIN_MODULE_SIZE, reassignThreshold=0,
  332. mergeCutHeight=0.25, numericLabels=TRUE, pamRespectsDendro=FALSE, verbose=2)
  333. list(colors = single$colors, dendrograms = single$dendrograms, blockGenes = single$blockGenes)
  334. }
  335. )
  336. cons_colors <- net_cons$colors; names(cons_colors) <- colnames(DatA)
  337. module_df <- tibble(
  338. gene = names(cons_colors),
  339. module_label = as.integer(cons_colors),
  340. module = std_me_names(as.integer(cons_colors))
  341. )
  342. write_csv(module_df, file.path(outdir, "wgcna", "Consensus_modules_genes.csv"))
  344. # ------------------ Eigengenes helper ----------------------------------------
  345. compute_MEs <- function(EgxS, colors_vec, prefix="") {
  346. if (is.null(EgxS) || ncol(EgxS) < 2) { log_msg("%s compute_MEs: not enough samples (%s).", prefix, ncol(EgxS)); return(NULL) }
  347. genes <- intersect(rownames(EgxS), names(colors_vec))
  348. cols <- colors_vec[genes]; cols <- cols[!is.na(cols) & cols != 0]
  349. genes <- names(cols)
  350. if (length(genes) < 3) { log_msg("%s compute_MEs: <3 non-grey genes.", prefix); return(NULL) }
  351. X <- as.data.frame(t(EgxS[genes, , drop=FALSE]))
  352. MEs <- tryCatch(orderMEs(moduleEigengenes(X, colors=cols, excludeGrey=TRUE)$eigengenes),
  353. error=function(e){ log_msg("%s compute_MEs: %s", prefix, conditionMessage(e)); NULL })
  354. if (!is.null(MEs)) { rownames(MEs) <- rownames(X); colnames(MEs) <- std_me_names(colnames(MEs)) }
  355. MEs
  356. }
  358. # ------------------ Cohorts (age metadata, continuous) -----------------------
  359. age_meta1 <- age_ds$meta %>%
  360. mutate(disease = 0L) %>%
  361. mutate(age = suppressWarnings(as.numeric(age))) %>% filter(is.finite(age))
  362. age_meta2 <- ad_ds$meta %>%
  363. mutate(disease = ifelse(sample_id %in% ad_ids, 1L, disease)) %>%
  364. mutate(disease = ifelse(is.na(disease), 0L, disease)) %>% filter(disease==0) %>%
  365. mutate(age = suppressWarnings(as.numeric(age))) %>% filter(is.finite(age))
  367. E_age1 <- EA[, intersect(colnames(EA), age_meta1$sample_id), drop=FALSE]
  368. E_age2 <- EC[, intersect(colnames(EC), age_meta2$sample_id), drop=FALSE]
  370. # 3) AD vs Control (>59) in GSE48350
  371. ctrl_59_meta <- ad_ds$meta %>%
  372. mutate(disease = ifelse(sample_id %in% ad_ids, 1L, disease),
  373. disease = ifelse(is.na(disease), 0L, disease),
  374. age = suppressWarnings(as.numeric(age))) %>%
  375. filter(disease==0, is.finite(age), age > 59)
  376. ad_59_ids <- intersect(colnames(ad_ds$expr), ad_ids)
  378. region_pick <- function(df){
  379. if ("region_raw" %in% names(df)) df$region_raw else if ("region" %in% names(df)) df$region else if ("tissue" %in% names(df)) df$tissue else rep(NA_character_, nrow(df))
  380. }
  381. ad_meta <- tibble(sample_id = ad_59_ids, disease=1L) %>%
  382. left_join(ad_csv %>% transmute(sample_id=as.character(sample_id),
  383. age = suppressWarnings(as.numeric(age)),
  384. sex = as.character(sex),
  385. region = region_pick(ad_csv),
  386. tissue = as.character(tissue)), by="sample_id")
  387. ctrl_meta <- ctrl_59_meta %>% select(sample_id, age, sex, region, tissue, disease)
  388. ad_ctrl_meta <- bind_rows(ad_meta %>% mutate(context="ad"),
  389. ctrl_meta %>% mutate(context="control")) %>%
  390. mutate(age = suppressWarnings(as.numeric(age)))
  392. AD_full0 <- ad_ds$expr %>% harmonize_gene_names() %>% { .[intersect(rownames(.), genes_cons), , drop=FALSE] }
  393. E_ad_ctrl59 <- AD_full0[, intersect(colnames(AD_full0), ad_ctrl_meta$sample_id), drop=FALSE] |> median_impute_rows()
  394. log_msg("[AD>59] genes=%s samples=%s (AD=%s, CTRL=%s)",
  395. nrow(E_ad_ctrl59), ncol(E_ad_ctrl59),
  396. sum(ad_ctrl_meta$disease==1, na.rm=TRUE),
  397. sum(ad_ctrl_meta$disease==0, na.rm=TRUE))
  399. # ------------------ Module-trait heatmaps (eigengene view) -------------------
  400. mod_trait_heatmap <- function(MEs, meta, mode=c("age","ad"), prefix=""){
  401. if (is.null(MEs) || ncol(MEs)==0) { log_msg("%s heatmap: no eigengenes; skipping.", prefix); return(invisible(NULL)) }
  402. mode <- match.arg(mode)
  403. traits <- if (mode=="age") data.frame(age = suppressWarnings(as.numeric(meta$age)))
  404. else data.frame(disease = suppressWarnings(as.numeric(meta$disease)))
  405. rownames(traits) <- meta$sample_id
  406. keep <- intersect(rownames(MEs), rownames(traits))
  407. if (length(keep) < 3) { log_msg("%s heatmap: <3 overlapping samples. Skipping.", prefix); return(invisible(NULL)) }
  408. MEs <- MEs[keep,,drop=FALSE]; traits <- traits[keep,,drop=FALSE]
  409. safe_cor_p <- function(x, y){
  410. ok <- is.finite(x) & is.finite(y); n <- sum(ok); if (n < 3) return(c(cor=NA_real_, p=NA_real_))
  411. x <- x[ok]; y <- y[ok]; if (sd(x)==0 || sd(y)==0) return(c(cor=NA_real_, p=NA_real_))
  412. r <- suppressWarnings(cor(x, y, method="pearson")); if (!is.finite(r)) return(c(cor=NA_real_, p=NA_real_))
  413. tstat <- r * sqrt((n-2)/(1 - r^2)); p <- 2 * stats::pt(-abs(tstat), df=n-2); c(cor=r, p=p)
  414. }
  415. C <- matrix(NA_real_, nrow=ncol(MEs), ncol=ncol(traits), dimnames=list(colnames(MEs), colnames(traits))); P <- C
  416. for (i in seq_len(ncol(MEs))) for (j in seq_len(ncol(traits))) { res <- safe_cor_p(MEs[,i], traits[,j]); C[i,j] <- res[1]; P[i,j] <- res[2] }
  417. if (all(!is.finite(C))) { log_msg("%s heatmap: all correlations NA. Skipping.", prefix); return(invisible(NULL)) }
  418. lab <- paste0(sprintf("%.2f", C), "\n(", ifelse(is.finite(P), sprintf("%.1e", P), "NA"), ")")
  419. df <- as.data.frame(as.table(C)); names(df) <- c("ME","Trait","Cor"); df$P <- as.vector(P); df$lab <- as.vector(lab)
  420. p <- ggplot(df, aes(Trait, ME, fill=Cor, label=lab)) +
  421. geom_tile() + geom_text(size=3) +
  422. scale_fill_gradient2(low="#2166AC", mid="white", high="#B2182B", midpoint=0, na.value="grey85") +
  423. labs(title=paste0(prefix, " module–trait relationships (", mode, ")"), x=NULL, y=NULL) +
  424. theme_minimal(base_size=12)
  425. safe_ggsave(file.path(outdir, "wgcna", paste0(prefix, "_module_trait_heatmap_", mode)), p, w=6.8, h=0.5 + 0.25*max(1, nrow(C)), dpi=PNG_DPI)
  426. }
  428. # ------------------ Module associations (eigengene) --------------------------
  429. MEs_age1 <- compute_MEs(E_age1, cons_colors, prefix="GSE48350_AgingCont")
  430. MEs_age2 <- compute_MEs(E_age2, cons_colors, prefix="GSE48350Ctrl_AgingCont")
  431. MEs_ad59 <- compute_MEs(E_ad_ctrl59, cons_colors, prefix="GSE48350_AD>59")
  433. mod_trait_heatmap(MEs_age1, age_meta1, mode="age", prefix="GSE48350_AgingCont")
  434. mod_trait_heatmap(MEs_age2, age_meta2, mode="age", prefix="GSE48350Ctrl_AgingCont")
  435. mod_trait_heatmap(MEs_ad59, ad_ctrl_meta, mode="ad", prefix="GSE48350_AD")
  437. # ======== KEY FIX: robust covariate building for association tests ===========
  438. build_covars <- function(meta, include=c("sex","batch","region","tissue","age")) {
  439. use_cols <- intersect(include, names(meta))
  440. if (!length(use_cols)) return(NULL)
  441. covs <- meta[, use_cols, drop=FALSE]
  442. rownames(covs) <- meta$sample_id
  443. sanitize_covars(covs)
  444. }
  446. test_binary_trait <- function(MEs, meta, trait, include_covars=c("sex","batch","region","tissue","age")){
  447. if (is.null(MEs) || ncol(MEs)==0) return(tibble(module=character(), beta=numeric(), p=numeric(), FDR=numeric()))
  448. stopifnot(trait %in% names(meta))
  449. keep <- intersect(rownames(MEs), meta$sample_id)
  450. if (length(keep) < 6) return(tibble(module=character(), beta=numeric(), p=numeric(), FDR=numeric()))
  451. MEs <- MEs[keep,,drop=FALSE]
  452. mm <- meta[match(keep, meta$sample_id), , drop=FALSE]
  454. covs_raw <- build_covars(mm, include=include_covars)
  455. covs <- sanitize_covars(covs_raw, sample_order=rownames(MEs))
  456. if (!is.null(covs)) {
  457. covs <- drop_collinear(as.numeric(mm[[trait]]), covs, thr=0.98)
  458. }
  460. out <- bind_rows(lapply(colnames(MEs), function(m){
  461. y <- as.numeric(MEs[,m])
  462. a <- assoc_binary(y, as.numeric(mm[[trait]]), covs, min_n=6)
  463. tibble(module=m, beta=as.numeric(a["beta"]), p=as.numeric(a["p"]))
  464. }))
  465. out$p <- psafe(out$p); out$FDR <- p.adjust(out$p, "BH"); out
  466. }
  468. res_age_mod1 <- test_binary_trait(MEs_age1, age_meta1 %>% mutate(age_group = ifelse(age > 59,1L,0L)), "age_group", include_covars=c("sex","batch","region","tissue"))
  469. res_age_mod2 <- test_binary_trait(MEs_age2, age_meta2 %>% mutate(age_group = ifelse(age > 59,1L,0L)), "age_group", include_covars=c("sex","batch","region","tissue"))
  470. res_ad_mod <- test_binary_trait(MEs_ad59, ad_ctrl_meta, "disease", include_covars=c("sex","batch","region","tissue","age"))
  472. res_age_mod1$module <- std_me_names(res_age_mod1$module)
  473. res_age_mod2$module <- std_me_names(res_age_mod2$module)
  474. res_ad_mod$module <- std_me_names(res_ad_mod$module)
  476. combine_age_modules <- function(a1, a2){
  477. if (nrow(a1)>0 && nrow(a2)>0) {
  478. a1 %>% rename(beta1=beta, p1=p, FDR1=FDR) %>%
  479. full_join(a2 %>% rename(beta2=beta, p2=p, FDR2=FDR), by="module") %>%
  480. mutate(beta = rowMeans(cbind(beta1, beta2), na.rm=TRUE),
  481. p = fisher_combine_p(na.omit(c(psafe(p1), psafe(p2)))),
  482. FDR = p.adjust(psafe(p), "BH")) %>%
  483. select(module, beta, p, FDR, beta1, beta2, p1, p2, FDR1, FDR2)
  484. } else {
  485. a <- if (nrow(a1)>0) a1 else a2
  486. a %>% rename(beta=beta, p=p, FDR=FDR)
  487. }
  488. }
  489. age_combined_mod <- combine_age_modules(res_age_mod1, res_age_mod2)
  491. ap_direction <- function(a, b) {
  492. res <- rep(NA_character_, length(a))
  493. res[a > 0 & b > 0] <- "AgeUp-ADDown (AP_Resilience)"
  494. res[a < 0 & b < 0] <- "AgeDown-ADUp (AP_Vulnerability)"
  495. res[is.na(res)] <- "Other"; res
  496. }
  497. call_ap_modules <- function(df_age, df_ad) {
  498. df <- full_join(
  499. df_age %>% select(module, beta_age=beta, p_age=p, FDR_age=FDR),
  500. df_ad %>% select(module, beta_ad=beta, p_ad=p, FDR_ad=FDR),
  501. by="module"
  502. ) %>%
  503. mutate(p_age = psafe(p_age), p_ad = psafe(p_ad),
  504. direction = ap_direction(beta_age, -beta_ad),
  505. AP_score = (sign(beta_age) * -sign(beta_ad)) * ( -log10(p_age) + -log10(p_ad) ),
  506. age_sig = FDR_age < REGION_P_FDR,
  507. ad_sig = FDR_ad < REGION_P_FDR)
  508. df
  509. }
  510. AP_modules <- call_ap_modules(age_combined_mod, res_ad_mod)
  511. sizes <- module_df %>% count(module, name="module_size")
  512. AP_modules <- AP_modules %>%
  513. left_join(sizes, by="module") %>%
  514. mutate(AP_class = case_when(
  515. direction == "AgeUp-ADDown (AP_Resilience)" ~ "AP_Resilience",
  516. direction == "AgeDown-ADUp (AP_Vulnerability)" ~ "AP_Vulnerability",
  517. TRUE ~ "Other")) %>%
  518. arrange(desc(AP_score))
  519. write_csv(AP_modules, file.path(outdir, "AP_modules.csv"))
  521. # ------------------ Hubs (kME) -----------------------------------------------
  522. hub_scores <- function(module_df, EgxS) {
  523. split(module_df$gene, module_df$module) %>%
  524. imap_dfr(function(genes, m) {
  525. g <- intersect(genes, rownames(EgxS))
  526. if (length(g) < 3 || ncol(EgxS) < 2) return(tibble(module=m, gene=NA_character_, kME=NA_real_))
  527. datExpr <- t(EgxS[g, , drop=FALSE])
  528. ME <- WGCNA::moduleEigengenes(datExpr, rep(1, ncol(datExpr)))$eigengenes[,1]
  529. kME <- suppressWarnings(cor(datExpr, ME, use="pairwise.complete.obs"))
  530. tibble(module=m, gene=colnames(datExpr), kME=as.numeric(kME)) %>% arrange(desc(kME)) %>% head(TOP_HUBS_PER_MODULE)
  531. })
  532. }
  533. age_hubs <- hub_scores(module_df, EA)
  534. ctrl_hubs <- hub_scores(module_df, EC)
  535. write_csv(age_hubs, file.path(outdir, "AP_module_hubs_GSE48350_healthy.csv"))
  536. write_csv(ctrl_hubs, file.path(outdir, "AP_module_hubs_GSE48350_controls.csv"))
  538. # -------------------- export summary -----------------------------------------
  539. write_csv(AP_modules, file.path(outdir, "AP_modules_summary_from_genes.csv"))
  540. cat("Quadrant plot and module summary saved in:", outdir, "\n")
  542. # ============================================================================#
  543. # ================= AP gene-level scoring (AGE = CONTINUOUS) =================#
  544. # ============================================================================#
  546. log_msg("Scoring genes for AP importance (age=continuous per cohort, meta-analyzed)...")
  548. # Covariates
  549. cov_age1 <- build_covars(age_meta1, include=c("sex","batch","region","tissue"))
  550. cov_age2 <- build_covars(age_meta2, include=c("sex","batch","region","tissue"))
  551. cov_ad <- build_covars(ad_ctrl_meta, include=c("sex","batch","region","tissue","age"))
  552. cov_ad <- drop_collinear(as.numeric(setNames(ad_ctrl_meta$disease, ad_ctrl_meta$sample_id)[colnames(E_ad_ctrl59)]), cov_ad, thr=0.98)
  554. # Align matrices & traits
  555. align_cols <- function(E, ids) E[, intersect(colnames(E), ids), drop=FALSE]
  556. Eg_age1 <- align_cols(E_age1, age_meta1$sample_id)
  557. Eg_age2 <- align_cols(E_age2, age_meta2$sample_id)
  558. Eg_ad <- align_cols(E_ad_ctrl59, ad_ctrl_meta$sample_id)
  560. x_age1 <- setNames(scale(age_meta1$age)[,1], age_meta1$sample_id) # scaled
  561. x_age2 <- setNames(scale(age_meta2$age)[,1], age_meta2$sample_id)
  562. x_ad <- setNames(ad_ctrl_meta$disease, ad_ctrl_meta$sample_id)
  564. # Association matrices
  565. do_assoc_matrix_cont <- function(EgxS, x, covars, min_n=4) {
  566. if (is.null(EgxS) || ncol(EgxS) < min_n) {
  567. log_msg("Assoc matrix (continuous): too few samples (%s).", ncol(EgxS));
  568. return(tibble(gene=character(), beta=numeric(), se=numeric(), p=numeric()))
  569. }
  570. if (!is.null(covars)) {
  571. covars <- sanitize_covars(covars, sample_order=colnames(EgxS))
  572. }
  573. genes <- rownames(EgxS)
  574. purrr::map_dfr(genes, function(g){
  575. y <- EgxS[g, ]
  576. a <- assoc_linear(y, x[colnames(EgxS)], covars, min_n=min_n)
  577. tibble(gene=g, beta=as.numeric(a["beta"]), se=as.numeric(a["se"]), p=psafe(a["p"]))
  578. })
  579. }
  580. do_assoc_matrix_bin <- function(EgxS, x, covars, min_n=6) {
  581. if (is.null(EgxS) || ncol(EgxS) < min_n) {
  582. log_msg("Assoc matrix (binary): too few samples (%s).", ncol(EgxS));
  583. return(tibble(gene=character(), beta=numeric(), se=numeric(), p=numeric()))
  584. }
  585. if (!is.null(covars)) {
  586. covars <- sanitize_covars(covars, sample_order=colnames(EgxS))
  587. }
  588. genes <- rownames(EgxS)
  589. purrr::map_dfr(genes, function(g){
  590. y <- EgxS[g, ]
  591. a <- assoc_binary(y, x[colnames(EgxS)], covars, min_n=min_n)
  592. tibble(gene=g, beta=as.numeric(a["beta"]), se=as.numeric(a["se"]), p=psafe(a["p"]))
  593. })
  594. }
  596. ga1 <- do_assoc_matrix_cont(Eg_age1, x_age1, cov_age1, min_n=4) %>% rename(beta_age1=beta, se_age1=se, p_age1=p)
  597. ga2 <- do_assoc_matrix_cont(Eg_age2, x_age2, cov_age2, min_n=4) %>% rename(beta_age2=beta, se_age2=se, p_age2=p)
  598. gad <- do_assoc_matrix_bin (Eg_ad, x_ad, cov_ad, min_n=6) %>% rename(beta_ad=beta, se_ad=se, p_ad=p)
  600. # Combine aging cohorts (inverse-variance for beta, Fisher for p)
  601. combine_age_genes <- function(ga1, ga2){
  602. full_join(ga1, ga2, by="gene") %>%
  603. mutate(
  604. beta_age = case_when(
  605. is.finite(beta_age1) & is.finite(beta_age2) & is.finite(se_age1) & is.finite(se_age2) ~ {
  606. w1 <- 1/(se_age1^2); w2 <- 1/(se_age2^2); (beta_age1*w1 + beta_age2*w2)/(w1+w2)
  607. },
  608. is.finite(beta_age1) ~ beta_age1,
  609. is.finite(beta_age2) ~ beta_age2,
  610. TRUE ~ NA_real_
  611. ),
  612. p_age = {
  613. pvec <- c(psafe(p_age1), psafe(p_age2)); pvec <- pvec[is.finite(pvec)]
  614. if (length(pvec)==0) NA_real_ else fisher_combine_p(pvec)
  615. }
  616. ) %>% select(gene, beta_age, p_age)
  617. }
  618. age_comb <- combine_age_genes(ga1, ga2)
  620. # Build gene table
  621. gene_tbl <- tibble(gene=genes_cons) %>%
  622. left_join(age_comb, by="gene") %>%
  623. left_join(gad %>% select(gene, beta_ad, p_ad), by="gene")
  625. # FDR (p.adjust keeps NA if p is NA)
  626. gene_tbl$p_age <- as.numeric(gene_tbl$p_age)
  627. gene_tbl$p_ad <- as.numeric(gene_tbl$p_ad)
  628. gene_tbl$FDR_age <- p.adjust(gene_tbl$p_age, "BH")
  629. gene_tbl$FDR_ad <- p.adjust(gene_tbl$p_ad, "BH")
  631. # Attach module + kME from healthy matrices
  632. all_kME <- function(E, module_df) {
  633. split(module_df$gene, module_df$module) %>%
  634. imap_dfr(function(gset, m){
  635. g <- intersect(gset, rownames(E))
  636. if (length(g) < 3 || ncol(E) < 2) return(tibble(gene=character(), module=m, kME=numeric()))
  637. datExpr <- t(E[g, , drop=FALSE]); ME <- WGCNA::moduleEigengenes(datExpr, rep(1, ncol(datExpr)))$eigengenes[,1]
  638. k <- suppressWarnings(cor(datExpr, ME, use="pairwise.complete.obs"))
  639. tibble(gene=colnames(datExpr), module=m, kME=as.numeric(k))
  640. })
  641. }
  642. kME_A <- all_kME(EA, module_df)
  643. kME_C <- all_kME(EC, module_df)
  644. kME_tbl <- full_join(kME_A %>% rename(kME_A=kME), kME_C %>% rename(kME_C=kME), by=c("gene","module")) %>%
  645. mutate(kME_mean = rowMeans(cbind(kME_A, kME_C), na.rm=TRUE))
  647. gene_tbl <- gene_tbl %>%
  648. left_join(module_df %>% select(gene, module), by="gene") %>%
  649. left_join(kME_tbl %>% select(gene, kME_A, kME_C, kME_mean), by="gene")
  651. # ------------------ AP gene metrics (normalized) -----------------------------
  652. p_or_one <- function(p) ifelse(is.na(p), 1, psafe(p))
  654. AP_antagonism <- sign(gene_tbl$beta_age) * -sign(gene_tbl$beta_ad)
  655. AP_antagonism[is.na(AP_antagonism)] <- 0
  657. AP_sig_sum <- (-log10(p_or_one(gene_tbl$p_age))) + (-log10(p_or_one(gene_tbl$p_ad)))
  658. effect_sum <- abs(gene_tbl$beta_age) + abs(gene_tbl$beta_ad)
  660. kme_w <- if (isTRUE(AP_SCORE_USE_KME)) pmax(gene_tbl$kME_mean, 0) else 1
  661. AP_score_raw <- AP_antagonism * effect_sum * AP_sig_sum * kme_w
  663. norm_minmax01 <- function(x){
  664. xf <- x[is.finite(x)]
  665. if (!length(xf) || length(unique(xf)) == 1) return(rep(0.5, length(x)))
  666. rng <- range(xf)
  667. (x - rng[1]) / (rng[2] - rng[1])
  668. }
  669. robust_z <- function(x){
  670. m <- stats::median(x[is.finite(x)], na.rm=TRUE)
  671. s <- stats::mad(x[is.finite(x)], constant=1.4826, na.rm=TRUE)
  672. if (!is.finite(s) || s==0) return(rep(0, length(x)))
  673. (x - m) / s
  674. }
  675. percentile_0_100 <- function(x){
  676. r <- rank(x, na.last="keep", ties.method="average")
  677. 100 * (r - 1) / (sum(is.finite(x)) - 1)
  678. }
  680. AP_score_norm01 <- norm_minmax01(AP_score_raw)
  681. AP_score_z <- robust_z(AP_score_raw)
  682. AP_score_percentile <- percentile_0_100(AP_score_raw)
  684. # Per-gene AP Z (for module aggregation; sign-aware)
  685. z_from_beta_p <- function(beta, p) {
  686. pp <- p_or_one(p)
  687. zz <- sign(beta) * qnorm(1 - pp/2)
  688. zz[!is.finite(zz)] <- 0
  689. zz
  690. }
  691. Z_age <- z_from_beta_p(gene_tbl$beta_age, gene_tbl$p_age)
  692. Z_ad <- z_from_beta_p(-gene_tbl$beta_ad, gene_tbl$p_ad) # minus sign encodes AP (AD down = resilience)
  694. Z_AP_gene <- (Z_age + Z_ad) / sqrt(2) # Stouffer combine (equal weights)
  696. # Store gene outputs
  697. gene_tbl$AP_antagonism <- AP_antagonism
  698. gene_tbl$AP_sig_sum <- AP_sig_sum
  699. gene_tbl$AP_score_gene_raw <- AP_score_raw
  700. gene_tbl$AP_score_gene_norm01 <- AP_score_norm01
  701. gene_tbl$AP_score_gene_z <- AP_score_z
  702. gene_tbl$AP_score_gene_percentile <- AP_score_percentile
  703. gene_tbl$Z_age <- Z_age
  704. gene_tbl$Z_ad <- Z_ad
  705. gene_tbl$Z_AP_gene <- Z_AP_gene
  706. gene_tbl$AP_score_gene <- AP_score_norm01 # alias used by downstream code
  708. ap_class_gene <- function(b_age, b_ad) {
  709. aa <- sign(b_age) * -sign(b_ad)
  710. ifelse(aa > 0 & sign(b_age) > 0, "AP_Resilience",
  711. ifelse(aa > 0 & sign(b_age) < 0, "AP_Vulnerability", "Other"))
  712. }
  713. gene_tbl$AP_class_gene <- ap_class_gene(gene_tbl$beta_age, gene_tbl$beta_ad)
  715. AP_genes_ranked <- gene_tbl %>%
  716. arrange(desc(AP_score_gene_norm01)) %>%
  717. select(gene, module, kME_A, kME_C, kME_mean,
  718. beta_age, p_age, FDR_age,
  719. beta_ad, p_ad, FDR_ad,
  720. AP_antagonism, AP_sig_sum, AP_class_gene,
  721. AP_score_gene_raw, AP_score_gene_norm01, AP_score_gene_z, AP_score_gene_percentile,
  722. Z_age, Z_ad, Z_AP_gene,
  723. AP_score_gene)
  725. readr::write_csv(AP_genes_ranked, file.path(outdir, "wgcna", "AP_genes_ranked.csv"))
  727. # ------------------ Gene-aggregated AP per module ----------------------------
  728. log_msg("Aggregating AP evidence across genes within modules (Stouffer Z)...")
  730. w_gene <- if (isTRUE(AP_MODULE_USE_KME_WEIGHTS)) pmax(gene_tbl$kME_mean, 0) else rep(1, nrow(gene_tbl))
  731. w_gene[!is.finite(w_gene)] <- 0
  733. agg_stouffer <- function(z, w){
  734. ok <- is.finite(z) & is.finite(w) & w > 0
  735. if (!any(ok)) return(NA_real_)
  736. sum(w[ok] * z[ok]) / sqrt(sum((w[ok])^2))
  737. }
  739. AP_modules_geneAgg <- gene_tbl %>%
  740. group_by(module) %>%
  741. summarise(
  742. n_genes = n(),
  743. n_antagonism = sum(AP_antagonism > 0, na.rm=TRUE),
  744. frac_antagonism= n_antagonism / n_genes,
  745. mean_AP_norm01 = mean(AP_score_gene_norm01, na.rm=TRUE),
  746. top10_AP_norm01= {
  747. k <- ceiling(0.10 * n_genes)
  748. if (k < 1) k <- 1
  749. mean(head(sort(AP_score_gene_norm01, decreasing=TRUE), k), na.rm=TRUE)
  750. },
  751. Z_from_genes = agg_stouffer(Z_AP_gene, w_gene[match(cur_data_all()$gene, gene_tbl$gene)]),
  752. p_from_genes = 2*pnorm(-abs(Z_from_genes))
  753. ) %>%
  754. ungroup() %>%
  755. mutate(FDR_from_genes = p.adjust(p_from_genes, "BH")) %>%
  756. arrange(p_from_genes)
  758. readr::write_csv(AP_modules_geneAgg, file.path(outdir, "wgcna", "AP_modules_geneAggregated.csv"))
  760. AP_modules_combined <- AP_modules %>%
  761. select(module, beta_age, p_age, FDR_age, beta_ad, p_ad, FDR_ad, AP_score, AP_class, module_size) %>%
  762. rename(EG_beta_age=beta_age, EG_p_age=p_age, EG_FDR_age=FDR_age,
  763. EG_beta_ad=beta_ad, EG_p_ad=p_ad, EG_FDR_ad=FDR_ad, EG_AP_score=AP_score, EG_AP_class=AP_class) %>%
  764. full_join(AP_modules_geneAgg, by="module")
  766. readr::write_csv(AP_modules_combined, file.path(outdir, "wgcna", "AP_modules_combined.csv"))
  768. # ---------- Ego networks for top genes ---------------------------------------
  769. .sim_cache <- new.env(parent=emptyenv())
  770. get_similarity_for_module <- function(mod, wobj) {
  771. key <- paste0(mod, "_", wobj$name)
  772. if (!is.null(.sim_cache[[key]])) return(.sim_cache[[key]])
  773. genes_m <- (wobj$module_df %>% filter(module==mod))$gene
  774. genes_m <- intersect(genes_m, rownames(wobj$expr))
  775. if (length(genes_m) < 5) return(NULL)
  776. datExpr <- t(wobj$expr[genes_m, , drop=FALSE])
  777. S <- NULL; ok <- TRUE
  778. tryCatch({ TOM <- WGCNA::TOMsimilarityFromExpr(datExpr, power=wobj$power, networkType=wobj$networkType); colnames(TOM) <- rownames(TOM) <- colnames(datExpr); S <<- TOM }, error=function(e){ ok <<- FALSE })
  779. if (!ok || is.null(S)) { Adj <- WGCNA::adjacency(datExpr, power=wobj$power, type=wobj$networkType); colnames(Adj) <- rownames(Adj) <- colnames(datExpr); S <- Adj }
  780. .sim_cache[[key]] <- S; S
  781. }
  782. export_gene_ego <- function(gene, mod, wobj, topk=15) {
  783. S <- get_similarity_for_module(mod, wobj)
  784. if (is.null(S) || !(gene %in% colnames(S))) return(invisible(NULL))
  785. v <- S[, gene]; v <- v[names(v) != gene]; v <- sort(v, decreasing=TRUE)
  786. nbrs <- head(names(v[v > 0]), min(topk, sum(v > 0)))
  787. if (!length(nbrs)) return(invisible(NULL))
  788. nodes <- unique(c(gene, nbrs))
  789. sub <- S[nodes, nodes, drop=FALSE]
  790. edges <- as.data.frame(as.table(sub), stringsAsFactors=FALSE)
  791. names(edges) <- c("from","to","weight")
  792. edges <- edges %>% filter(from < to, weight > 0)
  793. base <- file.path(outdir, "wgcna", "networks", paste0(wobj$name, "_", mod, "_ego_", make.names(gene)))
  794. write_csv(edges, paste0(base, "_edges.csv"))
  795. readr::write_csv(tibble(gene=nodes), paste0(base, "_nodes.csv"))
  796. if (ok_igraph && ok_ggraph && nrow(edges) > 0) {
  797. g <- igraph::graph_from_data_frame(edges, directed=FALSE, vertices=tibble(name=nodes))
  798. p <- ggraph::ggraph(g, layout="fr") + ggraph::geom_edge_link(alpha=0.3) +
  799. ggraph::geom_node_point(size=3, alpha=0.9) +
  800. ggraph::geom_node_text(aes(label=name), repel=TRUE, size=3) +
  801. ggplot2::labs(title=paste0("Ego network: ", gene, " (", mod, ")")) +
  802. ggplot2::theme_void(base_size=12)
  803. safe_ggsave(paste0(base, "_plot"), p, w=6.5, h=5.5, dpi=PNG_DPI)
  804. }
  805. }
  807. topR <- AP_genes_ranked %>% filter(AP_class_gene=="AP_Resilience", module!="ME00") %>% head(200)
  808. topV <- AP_genes_ranked %>% filter(AP_class_gene=="AP_Vulnerability", module!="ME00") %>% head(200)
  809. if (nrow(topR)==0) topR <- AP_genes_ranked %>% filter(module!="ME00") %>% arrange(desc(AP_sig_sum)) %>% head(200)
  810. if (nrow(topV)==0) topV <- AP_genes_ranked %>% filter(module!="ME00") %>% arrange(desc(AP_sig_sum)) %>% head(200)
  812. topR_pick <- topR %>% arrange(desc(AP_score_gene_norm01)) %>% head(TOP_EGO_PER_CLASS)
  813. topV_pick <- topV %>% arrange(desc(AP_score_gene_norm01)) %>% head(TOP_EGO_PER_CLASS)
  814. if (nrow(topR_pick)+nrow(topV_pick) == 0) {
  815. log_msg("No top AP genes to export ego networks for.")
  816. } else {
  817. wobj_export <- list(name="ConsensusOn_GSE48350", module_df=module_df, expr=EC, power=cons_power, networkType="signed")
  818. for (row in seq_len(nrow(topR_pick))) export_gene_ego(topR_pick$gene[row], topR_pick$module[row], wobj_export)
  819. for (row in seq_len(nrow(topV_pick))) export_gene_ego(topV_pick$gene[row], topV_pick$module[row], wobj_export)
  820. log_msg("Exported ego networks for top AP genes (per class).")
  821. }
  823. # ---- Save combined objects for reproducibility ----
  824. saveRDS(list(
  825. modules = module_df,
  826. AP_modules = AP_modules,
  827. AP_modules_geneAgg = AP_modules_geneAgg,
  828. AP_modules_combined = AP_modules_combined,
  829. AP_genes = AP_genes_ranked,
  830. hubs_age = age_hubs,
  831. hubs_ctrl = ctrl_hubs,
  832. power = cons_power,
  833. colors = cons_colors,
  834. genes_cons = genes_cons,
  835. AP_SCORE_USE_KME = AP_SCORE_USE_KME,
  836. AP_MODULE_USE_KME_WEIGHTS = AP_MODULE_USE_KME_WEIGHTS
  837. ), file = file.path(outdir, "wgcna", "AP_wgcna_objects.rds"))
  839. # ------------------ Quadrant plot (eigengene view) ---------------------------
  842. suppressPackageStartupMessages({ library(ggrepel); library(scales); library(readr); library(dplyr); library(ggplot2) })
  844. # Load the combined module table produced above
  845. mods <- readr::read_csv(file.path(outdir, "wgcna", "AP_modules_combined.csv"), show_col_types = FALSE)
  847. # If EG_AP_class already exists, use it; otherwise infer from EG betas
  848. infer_class <- function(b_age, b_ad) {
  849. if (is.finite(b_age) && is.finite(b_ad)) {
  850. if (b_age > 0 && b_ad < 0) "AP_Resilience"
  851. else if (b_age < 0 && b_ad > 0) "AP_Vulnerability"
  852. else "Other"
  853. } else "Other"
  854. }
  856. AP_COLORS <- c("AP_Resilience"="#009E73", "AP_Vulnerability"="#D55E00", "Other"="grey70")
  858. AP_modules_plot <- mods %>%
  859. mutate(
  860. AP_class = if (!is.null(EG_AP_class)) EG_AP_class else mapply(infer_class, EG_beta_age, EG_beta_ad),
  861. x = EG_beta_age,
  862. y = -EG_beta_ad,
  863. size_raw = pmax(0,
  864. -log10(pmax(EG_p_age, .Machine$double.xmin)) +
  865. -log10(pmax(EG_p_ad, .Machine$double.xmin))
  866. ),
  867. size_raw = ifelse(is.finite(size_raw), size_raw, 0),
  868. size_s = scales::rescale(
  869. pmin(size_raw, stats::quantile(size_raw[size_raw>0], 0.99, na.rm=TRUE)),
  870. to = c(2.5, 9),
  871. from = range(size_raw[size_raw>0], na.rm=TRUE)
  872. )
  873. )
  875. plt_quad <- ggplot(AP_modules_plot, aes(x = x, y = y)) +
  876. geom_hline(yintercept = 0, linetype = 2, linewidth = 0.3) +
  877. geom_vline(xintercept = 0, linetype = 2, linewidth = 0.3) +
  878. geom_point(aes(color = AP_class, size = size_s), alpha = 0.9) +
  879. ggrepel::geom_text_repel(aes(label = module),
  880. min.segment.length = 0, seed = 42,
  881. box.padding = 0.25, point.padding = 0.2,
  882. size = 3.2, color = "grey20", max.overlaps = Inf) +
  883. scale_color_manual(values = AP_COLORS, name = "AP Class", drop = FALSE) +
  884. scale_size_identity(guide = "legend", name = "-log10(p_age)+-log10(p_ad)") +
  885. labs(
  886. x = "Module effect (Healthy Aging, β_age; eigengene)",
  887. y = "Module effect in AD (−β_ad; eigengene)",
  888. title = "AP Module Quadrant Plot (Aging vs AD)"
  889. ) +
  890. coord_equal(expand = TRUE) +
  891. theme_minimal(base_size = 12) +
  892. theme(legend.position = "right")
  894. dir.create(file.path(outdir, "wgcna"), showWarnings = FALSE, recursive = TRUE)
  895. safe_ggsave(file.path(outdir, "wgcna", "fig_modules_quadrant_AP_modules_from_EG"), plt_quad, w = 7.5, h = 6.5, dpi = PNG_DPI)
  897. # ==================== High-AP modules: plots + gene composition ==============
  898. suppressPackageStartupMessages({ library(dplyr); library(readr); library(ggplot2); library(ggrepel); library(scales); })
  900. mods <- read_csv(file.path(outdir, "wgcna", "AP_modules_combined.csv"), show_col_types = FALSE)
  901. genes <- read_csv(file.path(outdir, "wgcna", "AP_genes_ranked.csv"), show_col_types = FALSE)
  903. # --- Selection knobs (override via env) --------------------------------------
  904. AP_HIGHSCORE_MIN <- as.numeric(Sys.getenv("AP_HIGHSCORE_MIN", unset = NA)) # e.g. 4
  905. AP_TOP_N <- as.integer(Sys.getenv("AP_TOP_N", unset = 6)) # fallback if none pass threshold
  907. # Helper: pick high-AP modules by eigengene-level score (antagonism-weighted)
  908. pick_high_ap <- function(df, min_score = NA, top_n = 6){
  909. d <- df %>% mutate(EG_AP_abs = abs(EG_AP_score))
  910. if (is.finite(min_score)) {
  911. keep <- d %>% filter(EG_AP_abs >= min_score)
  912. if (nrow(keep) > 0) return(keep)
  913. }
  914. d %>% slice_max(order_by = EG_AP_abs, n = top_n, with_ties = FALSE)
  915. }
  917. mods_sel <- pick_high_ap(mods, AP_HIGHSCORE_MIN, AP_TOP_N)
  919. if (nrow(mods_sel) == 0) {
  920. warning("No modules selected for plotting; check thresholds.")
  921. } else {
  922. # -------------- Plot 1: Quadrant for selected modules (eigengene stats) ---
  923. infer_class <- function(b_age, b_ad) {
  924. if (is.finite(b_age) && is.finite(b_ad)) {
  925. if (b_age > 0 && b_ad < 0) "AP_Resilience"
  926. else if (b_age < 0 && b_ad > 0) "AP_Vulnerability"
  927. else "Other"
  928. } else "Other"
  929. }
  931. AP_COLORS <- c("AP_Resilience" = "#009E73",
  932. "AP_Vulnerability" = "#D55E00",
  933. "Other" = "grey70")
  935. plot_df <- mods_sel %>%
  936. mutate(
  937. AP_class = if ("EG_AP_class" %in% names(.)) EG_AP_class else mapply(infer_class, EG_beta_age, EG_beta_ad),
  938. x = EG_beta_age,
  939. y = -EG_beta_ad,
  940. size_raw = pmax(0,
  941. -log10(pmax(EG_p_age, .Machine$double.xmin)) +
  942. -log10(pmax(EG_p_ad, .Machine$double.xmin))
  943. ),
  944. size_raw = ifelse(is.finite(size_raw), size_raw, 0),
  945. size_s = scales::rescale(
  946. pmin(size_raw, quantile(size_raw[size_raw>0], 0.99, na.rm=TRUE)),
  947. to = c(2.5, 9),
  948. from = range(size_raw[size_raw>0], na.rm=TRUE)
  949. )
  950. )
  952. plt_quad_hi <- ggplot(plot_df, aes(x = x, y = y)) +
  953. geom_hline(yintercept = 0, linetype = 2, linewidth = 0.3) +
  954. geom_vline(xintercept = 0, linetype = 2, linewidth = 0.3) +
  955. geom_point(aes(color = AP_class, size = size_s), alpha = 0.9) +
  956. ggrepel::geom_text_repel(aes(label = module),
  957. seed = 42, box.padding = 0.25, point.padding = 0.2,
  958. size = 3.2, color = "grey20", min.segment.length = 0) +
  959. scale_color_manual(values = AP_COLORS, name = "AP Class", drop = FALSE) +
  960. scale_size_identity(guide = "legend", name = "-log10(p_age)+-log10(p_ad)") +
  961. labs(
  962. title = "High-AP Modules (eigengene view)",
  963. x = "β_age (Healthy aging, eigengene)",
  964. y = "−β_ad (AD vs Control, eigengene)",
  965. subtitle = sprintf("Selected by |AP score|%s",
  966. if (is.finite(AP_HIGHSCORE_MIN)) paste0(" ≥ ", AP_HIGHSCORE_MIN) else
  967. paste0(" — top ", AP_TOP_N))
  968. ) +
  969. coord_equal(expand = TRUE) +
  970. theme_minimal(base_size = 12) +
  971. theme(legend.position = "right")
  973. safe_ggsave(file.path(outdir, "wgcna", "fig_quadrant_highAP_modules"), plt_quad_hi, w = 7.5, h = 6.5, dpi = PNG_DPI)
  975. # --------- Composition: #genes per module by AP_class_gene (gene-level) ----
  976. comp_long <- genes %>%
  977. semi_join(mods_sel %>% select(module), by = "module") %>%
  978. mutate(AP_class_gene = factor(AP_class_gene, levels = c("AP_Resilience","AP_Vulnerability","Other"))) %>%
  979. count(module, AP_class_gene, name = "n_genes") %>%
  980. group_by(module) %>%
  981. mutate(total_genes = sum(n_genes), frac = n_genes / total_genes) %>%
  982. ungroup()
  984. # Save a tidy composition table (answers your “are you determining the numbers?” → yes)
  985. write_csv(comp_long, file.path(outdir, "wgcna", "highAP_modules_gene_composition.csv"))
  987. # -------------- Plot 2: Stacked bars by gene AP class (counts & fraction) --
  988. plt_counts <- ggplot(comp_long, aes(x = module, y = n_genes, fill = AP_class_gene)) +
  989. geom_col(width = 0.8, color = "white") +
  990. scale_fill_manual(values = c("AP_Resilience"="#009E73","AP_Vulnerability"="#D55E00","Other"="grey70"),
  991. name = "Gene AP class") +
  992. labs(title = "Gene composition of high-AP modules",
  993. subtitle = "Counts by AP_Resilience / AP_Vulnerability / Other",
  994. x = "Module", y = "# genes") +
  995. theme_minimal(base_size = 12) +
  996. theme(legend.position = "right")
  997. safe_ggsave(file.path(outdir, "wgcna", "fig_highAP_modules_gene_counts"), plt_counts, w = 8, h = 5, dpi = PNG_DPI)
  999. plt_frac <- ggplot(comp_long, aes(x = module, y = frac, fill = AP_class_gene)) +
  1000. geom_col(width = 0.8, color = "white") +
  1001. scale_y_continuous(labels = percent_format()) +
  1002. scale_fill_manual(values = c("AP_Resilience"="#009E73","AP_Vulnerability"="#D55E00","Other"="grey70"),
  1003. name = "Gene AP class") +
  1004. labs(title = "Gene composition of high-AP modules",
  1005. subtitle = "Fraction within each module",
  1006. x = "Module", y = "Fraction") +
  1007. theme_minimal(base_size = 12) +
  1008. theme(legend.position = "right")
  1009. safe_ggsave(file.path(outdir, "wgcna", "fig_highAP_modules_gene_fractions"), plt_frac, w = 8, h = 5, dpi = PNG_DPI)
  1011. # ================== SINGLE FIGURE: Lollipop (top) + Donuts (bottom) ==================
  1012. suppressPackageStartupMessages({
  1013. library(dplyr); library(readr); library(ggplot2); library(scales)
  1014. })
  1015. has_patchwork <- requireNamespace("patchwork", quietly = TRUE)
  1017. mods <- readr::read_csv(file.path(outdir, "wgcna", "AP_modules_combined.csv"), show_col_types = FALSE)
  1018. genes <- readr::read_csv(file.path(outdir, "wgcna", "AP_genes_ranked.csv"), show_col_types = FALSE)
  1020. # -------- selection knobs (change via env if you like) --------
  1021. AP_HIGHSCORE_MIN <- as.numeric(Sys.getenv("AP_HIGHSCORE_MIN", unset = NA)) # e.g., 4
  1022. AP_TOP_N <- as.integer(Sys.getenv("AP_TOP_N", unset = 6)) # top-N if no cutoff
  1024. pick_high_ap <- function(df, min_score = NA, top_n = 6){
  1025. d <- df %>% mutate(EG_AP_abs = abs(EG_AP_score))
  1026. if (is.finite(min_score)) {
  1027. keep <- d %>% filter(EG_AP_abs >= min_score)
  1028. if (nrow(keep) > 0) return(keep)
  1029. }
  1030. d %>% slice_max(order_by = EG_AP_abs, n = top_n, with_ties = FALSE)
  1031. }
  1033. mods_sel <- pick_high_ap(mods, AP_HIGHSCORE_MIN, AP_TOP_N) %>%
  1034. arrange(desc(abs(EG_AP_score))) %>%
  1035. mutate(module = factor(module, levels = unique(module)))
  1037. if (nrow(mods_sel) == 0) {
  1038. warning("No modules selected — lower AP_HIGHSCORE_MIN or increase AP_TOP_N.");
  1039. } else {
  1040. AP_COLORS <- c("AP_Resilience"="#009E73","AP_Vulnerability"="#D55E00","Other"="grey70")
  1042. # -------------------- TOP: Lollipop rank by |AP score| ---------------------
  1043. p_top <- ggplot(mods_sel, aes(y = module, x = abs(EG_AP_score))) +
  1044. geom_segment(aes(yend = module, x = 0, xend = abs(EG_AP_score)), linewidth = 1, alpha = 0.45) +
  1045. geom_point(aes(color = EG_AP_class,
  1046. size = -log10(pmax(EG_p_age, .Machine$double.xmin)) +
  1047. -log10(pmax(EG_p_ad, .Machine$double.xmin))),
  1048. alpha = 0.95) +
  1049. scale_color_manual(values = AP_COLORS, drop = FALSE, name = "AP class") +
  1050. scale_size_continuous(range = c(3, 9), name = expression(-log[10](p[age])+-log[10](p[AD]))) +
  1051. labs(title = "High-AP modules", subtitle = "Ranked by |AP score| (eigengene)",
  1052. x = "|AP score|", y = NULL) +
  1053. theme_minimal(base_size = 12) +
  1054. theme(legend.position = "right", panel.grid.major.y = element_blank())
  1056. # --------------- BOTTOM: Donut row of gene-class composition ---------------
  1057. comp <- genes %>%
  1058. semi_join(mods_sel %>% select(module), by = "module") %>%
  1059. mutate(AP_class_gene = factor(AP_class_gene, levels = c("AP_Resilience","AP_Vulnerability","Other")),
  1060. module = factor(module, levels = levels(mods_sel$module))) %>%
  1061. count(module, AP_class_gene, name = "n") %>%
  1062. group_by(module) %>% mutate(frac = n/sum(n)) %>% ungroup()
  1064. p_bottom <- ggplot(comp, aes(x = 2, y = frac, fill = AP_class_gene)) +
  1065. geom_col(width = 0.6, color = "white") +
  1066. coord_polar(theta = "y") +
  1067. facet_wrap(~ module, nrow = 1) +
  1068. scale_fill_manual(values = AP_COLORS, name = "Gene AP class") +
  1069. xlim(0.5, 2.6) +
  1070. labs(title = NULL, subtitle = NULL) +
  1071. theme_void(base_size = 12) +
  1072. theme(legend.position = "right",
  1073. strip.text = element_text(face = "bold"))
  1075. # ---------------------- Stack (donut at the bottom) ------------------------
  1076. dir.create(file.path(outdir, "wgcna"), recursive = TRUE, showWarnings = FALSE)
  1077. if (has_patchwork) {
  1078. fig <- p_top / p_bottom + patchwork::plot_layout(heights = c(2, 1), guides = "collect")
  1079. safe_ggsave(file.path(outdir, "wgcna", "fig_highAP_lollipop_plus_donuts"), fig, w = 11, h = 7.5, dpi = PNG_DPI)
  1080. } else {
  1081. # Save separately if patchwork isn’t installed
  1082. safe_ggsave(file.path(outdir, "wgcna", "fig_highAP_lollipop"), p_top, w = 11, h = 5.0, dpi = PNG_DPI)
  1083. safe_ggsave(file.path(outdir, "wgcna", "fig_highAP_donuts"), p_bottom, w = 11, h = 3.0, dpi = PNG_DPI)
  1084. }
  1086. # Also write the composition table (counts + fractions) for the figure
  1087. readr::write_csv(comp, file.path(outdir, "wgcna", "highAP_modules_gene_composition_for_donuts.csv"))
  1088. }
  1090. # ---------------- Top genes per module (10 highest AP_score_gene) ----------
  1091. TOP_GENES_PER_MODULE <- as.integer(Sys.getenv("TOP_GENES_PER_MODULE", unset = "10"))
  1093. top_genes <- genes %>%
  1094. semi_join(mods_sel %>% dplyr::select(module), by = "module") %>%
  1095. dplyr::group_by(module) %>%
  1096. dplyr::slice_max(order_by = AP_score_gene, n = TOP_GENES_PER_MODULE, with_ties = FALSE) %>%
  1097. dplyr::ungroup() %>%
  1098. # build a per-module ordering key so each facet sorts independently
  1099. dplyr::group_by(module) %>%
  1100. dplyr::arrange(AP_score_gene, .by_group = TRUE) %>%
  1101. dplyr::mutate(key = paste(module, gene, sep = "::")) %>%
  1102. dplyr::ungroup()
  1104. # lock factor levels using the per-module order
  1105. key_levels <- top_genes %>%
  1106. dplyr::group_by(module) %>%
  1107. dplyr::arrange(AP_score_gene, .by_group = TRUE) %>%
  1108. dplyr::pull(key)
  1109. top_genes$key <- factor(top_genes$key, levels = unique(key_levels))
  1111. p_topgenes <- ggplot(top_genes,
  1112. aes(x = AP_score_gene, y = key, color = AP_class_gene)) +
  1113. geom_segment(aes(x = 0, xend = AP_score_gene, y = key, yend = key),
  1114. alpha = 0.45, linewidth = 0.9) +
  1115. geom_point(size = 2.6, alpha = 0.95) +
  1116. scale_color_manual(values = AP_COLORS,
  1117. name = "Gene AP class",
  1118. breaks = c("AP_Resilience","AP_Vulnerability","Other"),
  1119. drop = FALSE) +
  1120. scale_y_discrete(labels = function(x) sub(".*::", "", x)) +
  1121. facet_wrap(~ module, scales = "free_y") +
  1122. labs(
  1123. title = sprintf("Top %d genes by AP score (per module)", TOP_GENES_PER_MODULE),
  1124. x = "AP_score_gene (normalized 0–1)", y = NULL
  1125. ) +
  1126. coord_cartesian(clip = "off") +
  1127. theme_minimal(base_size = 12) +
  1128. theme(
  1129. legend.position = "right",
  1130. strip.text = element_text(face = "bold"),
  1131. panel.grid.major.y = element_blank()
  1132. )
  1134. safe_ggsave(file.path(outdir, "wgcna", "fig_highAP_top_genes_per_module"),
  1135. p_topgenes, w = 11, h = 6.5, dpi = PNG_DPI)
  1137. # export the exact list that appears in the plot
  1138. readr::write_csv(top_genes %>%
  1139. dplyr::transmute(module, gene,
  1140. AP_score_gene,
  1141. AP_class_gene,
  1142. beta_age, p_age, FDR_age,
  1143. beta_ad, p_ad, FDR_ad,
  1144. kME_mean),
  1145. file.path(outdir, "wgcna", "top_genes_per_module_for_figure.csv"))
  1148. # Also print a quick text summary to the log
  1149. log_msg("High-AP modules selected: %s", paste(sort(unique(mods_sel$module)), collapse=", "))
  1150. log_msg("Composition table written to: %s",
  1151. file.path(outdir, "wgcna", "highAP_modules_gene_composition.csv"))
  1154. # ======================= ALL MODULES: where are the strong genes? =======================
  1155. suppressPackageStartupMessages({ library(dplyr); library(readr); library(ggplot2); library(ggrepel); library(scales) })
  1157. # knobs (override via env)
  1158. AP_GENE_MIN_NORM <- as.numeric(Sys.getenv("AP_GENE_MIN_NORM", unset = "0.70")) # threshold on AP_score_gene_norm01
  1159. TOP_GENES_PER_MOD <- as.integer(Sys.getenv("TOP_GENES_PER_MOD", unset = "3")) # labels per module
  1161. # Ensure factor order: rank modules globally by |EG_AP_score| when available
  1162. mods_ranked <- mods %>%
  1163. mutate(EG_AP_abs = abs(EG_AP_score)) %>%
  1164. arrange(desc(EG_AP_abs)) %>%
  1165. mutate(module = factor(module, levels = unique(module)))
  1167. genes2 <- genes %>%
  1168. mutate(AP_class_gene = factor(AP_class_gene, levels = c("AP_Resilience","AP_Vulnerability","Other"))) %>%
  1169. left_join(mods_ranked %>% select(module, EG_AP_abs), by = "module") %>%
  1170. mutate(module = factor(module, levels = levels(mods_ranked$module)))
  1172. # ---------- (1) Heatmap: counts of high-AP genes per module × gene AP class ----------
  1173. hi_counts <- genes2 %>%
  1174. filter(is.finite(AP_score_gene_norm01)) %>%
  1175. mutate(is_high = AP_score_gene_norm01 >= AP_GENE_MIN_NORM) %>%
  1176. group_by(module, AP_class_gene) %>%
  1177. summarise(n_high = sum(is_high, na.rm = TRUE), .groups = "drop")
  1179. # Fill missing combos with zero for a clean heatmap
  1180. all_mods <- levels(genes2$module)
  1181. all_class <- levels(genes2$AP_class_gene)
  1182. hi_counts <- tidyr::complete(hi_counts, module = all_mods, AP_class_gene = all_class, fill = list(n_high = 0))
  1184. p_heat <- ggplot(hi_counts, aes(x = AP_class_gene, y = module, fill = n_high)) +
  1185. geom_tile(color = "white", linewidth = 0.3) +
  1186. scale_fill_gradient(low = "grey95", high = "black", name = "# high-AP genes") +
  1187. scale_x_discrete(position = "top") +
  1188. labs(
  1189. title = sprintf("High-AP genes across ALL modules (threshold: ≥ %.2f)", AP_GENE_MIN_NORM),
  1190. x = "Gene AP class", y = "Module (ranked by |EG AP|)"
  1191. ) +
  1192. theme_minimal(base_size = 12) +
  1193. theme(panel.grid = element_blank(), axis.text.x = element_text(face = "bold"))
  1194. safe_ggsave(file.path(outdir, "wgcna", "fig_allModules_highAP_heatmap"), p_heat, w = 6.5, h = max(6, 0.25*length(all_mods)), dpi = PNG_DPI)
  1196. # ---------- (2) Small multiples: top K genes per module (labels) ---------------
  1197. topK_per_module <- genes2 %>%
  1198. group_by(module) %>%
  1199. arrange(desc(AP_score_gene_norm01), .by_group = TRUE) %>%
  1200. slice_head(n = TOP_GENES_PER_MOD) %>%
  1201. ungroup() %>%
  1202. mutate(label = paste0(gene, " (", scales::number(AP_score_gene_norm01, accuracy = 0.01), ")"))
  1204. # If some modules are tiny, ensure at least one row
  1205. topK_per_module <- topK_per_module %>% filter(!is.na(module))
  1207. p_smallmult <- ggplot(topK_per_module,
  1208. aes(x = AP_score_gene_norm01, y = reorder(label, AP_score_gene_norm01))) +
  1209. geom_col(aes(fill = AP_class_gene), width = 0.8, color = "white") +
  1210. geom_vline(xintercept = AP_GENE_MIN_NORM, linetype = 2, linewidth = 0.3) +
  1211. facet_wrap(~ module, ncol = 3, scales = "free_y") +
  1212. scale_fill_manual(values = AP_COLORS, name = "Gene AP class") +
  1213. scale_x_continuous(limits = c(0, 1), labels = label_number(accuracy = 0.1)) +
  1214. labs(
  1215. title = sprintf("Top %d genes per module by AP score (normalized 0–1)", TOP_GENES_PER_MOD),
  1216. subtitle = "ALL modules; dotted line = high-AP threshold",
  1217. x = "AP_score_gene (norm.)", y = NULL
  1218. ) +
  1219. theme_minimal(base_size = 12) +
  1220. theme(legend.position = "right",
  1221. strip.text = element_text(face = "bold"),
  1222. panel.grid.major.y = element_blank())
  1223. safe_ggsave(file.path(outdir, "wgcna", "fig_allModules_topK_genes"), p_smallmult, w = 11, h = 10, dpi = PNG_DPI)
  1225. # Export tidy tables used in these figures
  1226. readr::write_csv(hi_counts, file.path(outdir, "wgcna", "allModules_highAP_geneClass_counts.csv"))
  1227. readr::write_csv(topK_per_module, file.path(outdir, "wgcna", "allModules_topK_genes.csv"))
  1229. # Log a quick heads-up
  1230. log_msg("All-modules figures saved: high-AP heatmap + top-%d genes per module. Threshold=%.2f",
  1231. TOP_GENES_PER_MOD, AP_GENE_MIN_NORM)
  1233. }
  1236. log_msg("Done. Wrote:\n - %s\n - %s\n - %s\n - %s\n - %s\n - %s",
  1237. file.path(outdir, "AP_modules.csv"),
  1238. file.path(outdir, "wgcna", "AP_modules_geneAggregated.csv"),
  1239. file.path(outdir, "wgcna", "AP_modules_combined.csv"),
  1240. file.path(outdir, "wgcna"),
  1241. file.path(outdir, "wgcna", "networks"),
  1242. file.path(outdir, "wgcna", "AP_genes_ranked.csv"))

10_modules_AP.R at commit 377f6ae, no license · at the source

Overview

Authors: Rana Salihoglu1, Şehnaz Can2, Thomas Dandekar1, Elena Bencurova1
ORCID iDs: Rana Salihoglu, Şehnaz Can, Thomas Dandekar, Elena Bencurova
  1. Department of Bioinformatics, Julius-Maximilians-Universität Würzburg, Würzburg, Germany
  2. Department of Bioengineering, Marmara University, Istanbul, Turkey
Institutions: University of Würzburg (Germany); Marmara University (Türkiye)
Journal: Computational and structural biotechnology journal, volume 35, issue 1, article 0134
Dates: received 30 March 2026; accepted 21 May 2026; published online 8 June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.34133/csbj.0134 · PMID 42267139 · PMCID PMC13243799 · OpenAlex W7162071788
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), Alzheimer's / dementia (population), cellular / molecular (subfield)
Methods: Statistics, Smoothing, state filtering, decompositions, Preprocessing, Connectivity
Topic: Single-cell and spatial transcriptomics (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Citations: not cited yet (Europe PMC); 95 references in the paper

Abstract

Aging is the strongest risk factor for Alzheimer’s disease (AD); however, some individuals age without major cognitive decline, suggesting that resilience and vulnerability may be associated with distinct molecular trajectories. To investigate these trajectories, we performed an integrated transcriptomic analysis of human dermal fibroblasts (GSE113957) and multi-region brain profiles (GSE48350), extending previous dataset-specific studies that focused primarily on age prediction, regional variation, or synaptic/immune signatures. Healthy aging and AD were compared within a novel antagonistic pleiotropy (AP) framework. This approach prioritized genes and candidate transcriptional regulators with opposing age and disease-associated expression patterns. Across tissues, healthy aging was associated with relative preservation of metabolic, mitochondrial, and lipid-homeostatic programs, whereas AD was associated with suppression of these programs alongside greater inflammatory and immune pathway activity. AP-Vulnerability genes (Age↓/AD↑), including TAC1, FREM3, and SLC25A46, declined with age but were induced in AD. Conversely, AP-Resilience genes (Age↑/AD↓), including PTH2, PPDPF, and NEFH, increased during healthy aging but were reduced in AD. Pathway analyses suggested an association between metabolic programs and resilience, and between immune activation and vulnerability. Transcription-factor inference prioritized PPARG, NFE2L2, and TEAD4 as candidate resilience-associated regulators, showing directionally opposite patterns relative to immune- and developmental-related regulators in AD.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 18 matches between paragraphs and lines of code.

salihoglu/Alzheimer_AP

License: none: the authors keep all their rights
State: the link answers, verified on 27 September 2026
Evidence: files inventoried
Commit: 377f6aed20f1984d011e84b99e5ed5cd2db0c811, 30 March 2026
Languages: R (14)
Size: 15 files, 14 scripts
Software Heritage: not archived
Found in: “Data Availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (10 files), ggplot2 (9 files), clusterProfiler (2 files), igraph (1 file), patchwork (1 file), pheatmap (1 file), SingleCellExperiment (1 file), WGCNA (1 file)
Availability: 1 check, the latest on 27 September 2026: the link answers
  • 27 September 2026: the link answers
15 files
At the source: github.com/salihoglu/Alzheimer_AP

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 14 scripts, each with its path and the digest of its content;
  • 18 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Data Availability

https://github.com/salihoglu/Alzheimer_AP

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 4 authors, 88 references.

Cite

This paper

Salihoglu, R., Can, Ş., Dandekar, T., & Bencurova, E. (2026). Integrated Multi-Tissue Transcriptomics Reveals Antagonistic Pleiotropy in Aging and Alzheimer's Disease. Computational and structural biotechnology journal, 35(1), 0134. https://doi.org/10.34133/csbj.0134

BibTeX

@article{salihoglu2026integrated,
author = {Salihoglu, Rana and Can, Şehnaz and Dandekar, Thomas and Bencurova, Elena},
title = {{Integrated Multi-Tissue Transcriptomics Reveals Antagonistic Pleiotropy in Aging and Alzheimer's Disease}},
journal = {Computational and structural biotechnology journal},
year = {2026},
month = jun,
volume = {35},
number = {1},
pages = {0134},
publisher = {AAAS Science Partner Journal Program},
issn = {2001-0370},
doi = {10.34133/csbj.0134},
url = {https://doi.org/10.34133/csbj.0134},
pmid = {42267139},
pmcid = {PMC13243799}
}

RIS

TY - JOUR
AU - Salihoglu, Rana
AU - Can, Şehnaz
AU - Dandekar, Thomas
AU - Bencurova, Elena
TI - Integrated Multi-Tissue Transcriptomics Reveals Antagonistic Pleiotropy in Aging and Alzheimer's Disease
T2 - Computational and structural biotechnology journal
J2 - Comput Struct Biotechnol J
PY - 2026
DA - 2026/06/08
VL - 35
IS - 1
SP - 0134
SN - 2001-0370
PB - AAAS Science Partner Journal Program
DO - 10.34133/csbj.0134
UR - https://doi.org/10.34133/csbj.0134
LA - en
ER -

CSL-JSON

{
"id": "10.34133/csbj.0134",
"type": "article-journal",
"title": "Integrated Multi-Tissue Transcriptomics Reveals Antagonistic Pleiotropy in Aging and Alzheimer's Disease",
"container-title": "Computational and structural biotechnology journal",
"author": [
{
"family": "Salihoglu",
"given": "Rana"
},
{
"family": "Can",
"given": "Şehnaz"
},
{
"family": "Dandekar",
"given": "Thomas"
},
{
"family": "Bencurova",
"given": "Elena"
}
],
"container-title-short": "Comput Struct Biotechnol J",
"volume": "35",
"issue": "1",
"page": "0134",
"DOI": "10.34133/csbj.0134",
"PMID": "42267139",
"PMCID": "PMC13243799",
"ISSN": "2001-0370",
"publisher": "AAAS Science Partner Journal Program",
"URL": "https://doi.org/10.34133/csbj.0134",
"language": "en",
"issued": {
"date-parts": [
[
2026,
6,
8
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s41593-026-02367-0 [code]
A reproducible three-dimensional model of human brain tissue to investigate physiological and disease-associated microglia phenotypes.
Journal: Nature neuroscience
In common: WGCNA, SingleCellExperiment, igraph, 5 other tools, Alzheimer's / dementia, cellular / molecular, 3 references
[2] doi:10.1038/s41380-026-03629-w [code]
Maternal fasting during early gestation induces epigenetic alterations and schizophrenia-related phenotypes.
Journal: Molecular psychiatry
In common: WGCNA, igraph, clusterProfiler, 4 other tools, genetics / omics, cellular / molecular, 4 references
[3] doi:10.1016/j.xcrm.2026.102766 [code]
A longitudinal single-cell and spatial multiomic atlas of pediatric high-grade glioma.
Journal: Cell reports. Medicine
In common: WGCNA, SingleCellExperiment, igraph, 5 other tools, genetics / omics, cellular / molecular, 1 reference
[4] doi:10.34133/csbj.0108 [code]
Latent Factor Modeling Reveals Unexpected Spatial Heterogeneity in Human Alzheimer's Disease Brain Transcriptomes.
Journal: Computational and structural biotechnology journal
In common: clusterProfiler, pheatmap, ggplot2, 1 other tool, Alzheimer's / dementia, genetics / omics, cellular / molecular, 6 references
[5] doi:10.1038/s41593-026-02384-z [code]
cGAS-mediated type I IFN signaling contributes to disease progression in drug-refractory epilepsy.
Journal: Nature neuroscience
In common: WGCNA, SingleCellExperiment, igraph, 5 other tools, cellular / molecular, 1 reference
[6] doi:10.1038/s44318-026-00806-z [code]
Interspecific diversity in the neuronal composition of the mammalian cortex arises from heterochrony in neurogenesis.
Journal: The EMBO journal
In common: WGCNA, SingleCellExperiment, igraph, 5 other tools, 1 reference
[7] doi:10.1002/alz.71804
A meta-refined human Alzheimer's disease-associated gene subset shows partial mouse-model pathway correspondence and limited cross-cohort machine-learning transportability.
Journal: Alzheimer's & dementia : the journal of the Alzheimer's Association
In common: Alzheimer's / dementia, genetics / omics, cellular / molecular, 10 references
[8] doi:10.1038/s41467-026-72598-z [code]
Functional impact of genetic background on variable expressivity in neurodevelopmental disorders.
Journal: Nature communications
In common: WGCNA, clusterProfiler, pheatmap, 3 other tools, 4 references
[9] doi:10.1038/s41398-026-04200-5 [code]
Postmortem brain single-nucleus and bulk gene expression analyses identify shared and distinct abnormalities in bipolar disorder and major depressive disorder.
Journal: Translational psychiatry
In common: WGCNA, SingleCellExperiment, igraph, 4 other tools, genetics / omics, cellular / molecular, 1 reference
[10] doi:10.1038/s41467-026-73305-8 [code]
Comparative analysis of the cellular landscape in mammalian striatum.
Journal: Nature communications
In common: WGCNA, SingleCellExperiment, clusterProfiler, 4 other tools, genetics / omics, cellular / molecular, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.