OSCR

The aging epigenome: integrative analyses reveal intersection with Alzheimer's disease.

Code ↔ Paper

6 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 6 matches
  1. [1] § Methods › Differentially methylated regions analysis ↔ code/utility/coMethDMR_aux.R, lines 151–199 · score 0.63 · linear regression, co methylated, medians, errors, models, CpGs
  2. [2] § Methods › Pathway analysis ↔ code/utility/pathway.R, lines 69–184 · score 0.61 · methylRRA, methylGSA, pathways, UCSC, enriched, enrichment
  3. [3] § Methods › Inflation assessment and correction ↔ code/utility/annotation_and_bacon.R, lines 483–626 · score 0.57 · inflation factors, bacon corrected, bias, genomic
  4. [4] § Methods › Inflation assessment and correction ↔ code/utility/annotation_and_bacon.R, lines 483–626 · score 0.57 · Genomic inflation factors, lambda, bias, bacon
  5. [5] § Methods › Association of DNA methylation at individual CpGs with chronological age ↔ code/dmr/coMethDMR.Rmd, lines 212–240 · score 0.54 · CD4T, Gran, Mono, NK, sex, beta
  6. [6] § Methods › Pre-processing of DNA methylation data ↔ code/dmr/coMethDMR.Rmd, lines 187–210 · score 0.50 · bisulfite conversion, status, EPIC, FHS, DNAm

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 · 629 lines · 24 KB · MIT · 2 matches

  1. #######################################################################################################
  2. # =================================================================================================== #
  3. # Function for annotation and inflation adjusted
  4. # =================================================================================================== #
  5. #######################################################################################################
  6. library(SummarizedExperiment)
  7. library(tidyverse)
  8. library(bacon)
  9. library(GWASTools)
  10. library(minfi)
  11. library(rGREAT)
  12. # ===================================================================================================
  13. # Annotation Single CpGs
  14. # ===================================================================================================
  15. # Function: get_anno_gr
  16. # ------------------------------------------------------------------------------
  17. # Description:
  18. # Retrieves annotation information as a GRanges object for the specified methylation array and genome.
  19. #
  20. # Parameters:
  21. # array : Character. An input array (e.g., methylation array) for which genomic annotations
  22. # are required. Supported options include "HM450", "EPICv1", "EPIC", and "EPICv2".
  23. # genome : Character. A string specifying the genome build used for annotation, such as "hg19"
  24. # or "hg38".
  25. # dir.data.aux : (Optional) Character. A directory path containing auxiliary data used by the
  26. # annotation function, particularly for reading external manifest files when using
  27. # genome "hg38".
  28. #
  29. # Returns:
  30. # A GRanges object containing the genomic annotation data.
  32. get_anno_gr <- function(array = "HM450",
  33. genome = "hg19",
  34. dir.data.aux = NULL) {
  36. if(array == "HM450"){
  37. library(IlluminaHumanMethylation450kanno.ilmn12.hg19)
  38. anno <- minfi::getAnnotation(IlluminaHumanMethylation450kanno.ilmn12.hg19)
  39. }
  40. if(array %in% c("EPICv1", "EPIC")){
  41. if(genome == "hg19") {
  42. library(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
  43. anno <- minfi::getAnnotation(IlluminaHumanMethylationEPICanno.ilm10b4.hg19)
  44. }
  45. if(genome == "hg38") {
  46. anno <- read_csv(
  47. file.path(dir.data.aux, "infinium-methylationepic-v-1-0-b5-manifest-file.csv"),
  48. show_col_types = F,
  49. skip = 7
  50. )
  51. }
  52. }
  53. if(array == "EPICv2") {
  54. library(IlluminaHumanMethylationEPICv2anno.20a1.hg38)
  55. anno <- minfi::getAnnotation(IlluminaHumanMethylationEPICv2anno.20a1.hg38)
  56. }
  58. if(array %in% c("EPICv1", "EPIC") & genome == "hg38") {
  59. anno.gr <- anno %>% makeGRangesFromDataFrame(
  60. seqnames.field = "CHR_hg38",
  61. start.field = "Start_hg38", end.field = "End_hg38",
  62. strand.field = "Strand_hg38",
  63. keep.extra.columns = T,
  64. na.rm = T
  65. )
  66. } else {
  67. anno.gr <- anno %>% makeGRangesFromDataFrame(
  68. start.field = "pos", end.field = "pos", keep.extra.columns = T
  69. )
  70. }
  71. anno.gr
  72. }
  73. # ------------------------------------------------------------------------------
  74. # Function: annotate_results
  75. # ------------------------------------------------------------------------------
  76. # Description:
  77. # Annotates a CpG results data frame by merging genomic annotation (from get_anno_gr) with
  78. # additional GREAT analysis results.
  79. #
  80. # Parameters:
  81. # result : Data frame. A table containing CpG IDs and associated statistics to annotate.
  82. # array : Character. The methylation array type (e.g., "HM450", "EPICv1", "EPIC").
  83. # genome : Character. The genome version ("hg19" or "hg38").
  84. # dir.data.aux : Character. Directory path for auxiliary data files containing precomputed annotations.
  85. # save : Logical. Whether to save the annotated result as a CSV file.
  86. # dir.save : Character. Directory path where the annotated result file will be saved.
  87. # prefix : Character. A prefix used for naming the output file.
  88. #
  89. # Returns:
  90. # A data frame with additional columns for genomic annotation and GREAT analysis.
  92. annotate_results <- function(result,
  93. array = "HM450",
  94. genome = "hg19",
  95. dir.data.aux = NULL,
  96. save = T,
  97. dir.save = NULL,
  98. prefix = "Framingham"){
  100. # load(file.path(dir.data.aux,"E073_15_coreMarks_segments.rda"))
  101. # load(file.path(dir.data.aux,"meta_analysis_cpgs.rda"))
  103. anno.gr <- get_anno_gr(array = array, genome = genome, dir.data.aux = dir.data.aux)
  104. anno_df <- as.data.frame(anno.gr)
  106. if(genome == "hg19") {
  108. if(array == "HM450"){
  109. load(file.path(dir.data.aux,"great_HM450_array_annotation.rda"))
  110. }
  111. if(array %in% c("EPICv1", "EPIC")){
  112. load(file.path(dir.data.aux,"great_EPIC_array_annotation.rda"))
  113. }
  114. result <- cbind(
  115. result,
  116. anno_df[result$cpg,c("seqnames", "start", "end", "width", "Relation_to_Island", "UCSC_RefGene_Name", "UCSC_RefGene_Group")]
  117. )
  118. result <- dplyr::left_join(result, great, by = c("seqnames","start","end","cpg"))
  119. }
  121. # Add annotation
  123. if(genome == "hg38"){
  125. if(array %in% c("EPICv1", "EPIC")) {
  126. load(file.path(dir.data.aux,"great_EPIC_array_annotation.rda"))
  127. anno_df <- anno_df %>%
  128. mutate(cpg = Name) %>%
  129. dplyr::select(cpg, seqnames, start, end, width,
  130. UCSC_RefGene_Group, UCSC_RefGene_Name, Relation_to_UCSC_CpG_Island) %>%
  131. unique()
  132. } else {
  133. load(file.path(dir.data.aux,"great_EPICv2_array_annotation.rda"))
  134. anno_df <- anno_df %>%
  135. mutate(cpg = gsub("_.*", "",anno_df$Name)) %>%
  136. dplyr::select(cpg, seqnames, start, end, width,
  137. UCSC_RefGene_Group, UCSC_RefGene_Name, Relation_to_Island) %>%
  138. unique()
  139. }
  141. result <- left_join(
  142. result,
  143. anno_df
  144. )
  146. result <- dplyr::left_join(result, great %>% ungroup %>% dplyr::select(cpg, GREAT_annotation))
  147. }
  149. if(save){
  150. write_csv(
  151. result,
  152. file.path(dir.save, paste0(prefix, "_annotated_results.csv"))
  153. )
  154. }
  156. return(result)
  157. }
  158. # ===================================================================================================
  159. # Annotation DMR
  160. # ===================================================================================================
  161. # Annotate the regions with E073 15-core marks segmentation states
  162. # ------------------------------------------------------------------------------
  163. # Description:
  164. # Annotates regions with ChromHMM segmentation states based on E073 15-core marks segmentation data.
  165. #
  166. # Parameters:
  167. # result : Data frame. Regions to be annotated.
  168. # ChmmModels.gr: GRanges object. ChromHMM segmentation data containing state information.
  169. # region_var : Character. The column name in 'result' representing region identifiers.
  170. #
  171. # Returns:
  172. # The original data frame with an added column "E073_15_coreMarks_segments_state" for segmentation states.
  174. annotate_coreMarks_segments <- function(result, ChmmModels.gr, region_var) {
  175. # Create a GRanges object from the result data frame
  176. result.gr <- result %>%
  177. makeGRangesFromDataFrame(start.field = "start",
  178. end.field = "end",
  179. seqnames.field = "seqnames")
  181. # Find overlaps between the result regions and the ChromHMM models
  182. hits <- findOverlaps(result.gr, ChmmModels.gr) %>% as.data.frame()
  183. hits$state <- ChmmModels.gr$state[hits$subjectHits]
  184. hits$region <- result[[region_var]][hits$queryHits]
  186. # Match the state information to each region based on the 'region' string
  187. result$E073_15_coreMarks_segments_state <- hits$state[match(result[[region_var]], hits$region)]
  189. return(result)
  190. }
  191. # ------------------------------------------------------------------------------
  192. # Function: annotate_great
  193. # ------------------------------------------------------------------------------
  194. # Description:
  195. # Annotates regions with gene associations using GREAT analysis.
  196. #
  197. # Parameters:
  198. # result : Data frame. Regions to be annotated.
  199. # genome : Character. The genome version ("hg19" or "hg38") used in GREAT analysis.
  200. #
  201. # Returns:
  202. # A data frame with an additional column "GREAT_annotation" containing GREAT-derived gene information.
  203. # ------------------------------------------------------------------------------
  204. annotate_great <- function(result, genome) {
  205. # Create a GRanges object from the result data frame
  206. result.gr <- result %>%
  207. makeGRangesFromDataFrame(start.field = "start",
  208. end.field = "end",
  209. seqnames.field = "seqnames")
  211. # Submit the GREAT job and retrieve gene associations
  212. job <- submitGreatJob(result.gr, species = genome)
  213. regionsToGenes_gr <- rGREAT::getRegionGeneAssociations(job)
  214. regionsToGenes <- as.data.frame(regionsToGenes_gr)
  216. # Create annotation strings for each region
  217. GREAT_annotation <- lapply(seq_len(length(regionsToGenes$annotated_genes)), function(i) {
  218. g <- ifelse(regionsToGenes$dist_to_TSS[[i]] > 0,
  219. paste0(regionsToGenes$annotated_genes[[i]], " (+", regionsToGenes$dist_to_TSS[[i]], ")"),
  220. paste0(regionsToGenes$annotated_genes[[i]], " (", regionsToGenes$dist_to_TSS[[i]], ")"))
  221. paste0(g, collapse = ";")
  222. })
  224. # Select key columns from GREAT output and combine with the annotation strings
  225. great <- dplyr::select(regionsToGenes, seqnames, start, end, width)
  226. great <- data.frame(great, GREAT_annotation = unlist(GREAT_annotation))
  228. # Merge the GREAT annotation with the original result data frame
  229. result <- dplyr::left_join(result, great, by = c("seqnames", "start", "end"))
  231. return(result)
  232. }
  233. # ------------------------------------------------------------------------------
  234. # Function: annotate_region
  235. # ------------------------------------------------------------------------------
  236. # Description:
  237. # Annotates regions with CpG UCSC information by summarizing overlapping CpG details.
  238. #
  239. # Parameters:
  240. # result : Data frame. Regions to be annotated.
  241. # array : Character. The methylation array type (e.g., "HM450", "EPIC").
  242. # genome : Character. The genome version ("hg19" or "hg38").
  243. # dir.data.aux: Character. Directory path for auxiliary annotation data.
  244. # cpg : Vector. List of CpG identifiers to be considered.
  245. # region_var : Character. (Optional) Column name in 'result' that contains region identifiers. Default is "region".
  246. # cores : Numeric. Number of CPU cores to use for parallel processing. Default is 30.
  247. #
  248. # Returns:
  249. # A data frame with additional columns summarizing UCSC CpG annotation information.
  251. annotate_region <- function(result, array, genome, dir.data.aux, cpg, region_var = "region", cores = 30) {
  253. anno.gr <- get_anno_gr(array = array, genome = genome, dir.data.aux = dir.data.aux)
  254. anno_df <- as.data.frame(anno.gr)
  256. if(genome == "hg38" & array %in% c("EPIC","EPICv1")) island_var <- "Relation_to_UCSC_CpG_Island"
  257. else island_var <- "Relation_to_Island"
  259. # Create a GRanges object from the result data frame
  260. result.gr <- result %>%
  261. makeGRangesFromDataFrame(start.field = "start",
  262. end.field = "end",
  263. seqnames.field = "seqnames")
  265. doParallel::registerDoParallel(cores)
  267. anno <- plyr::ldply(
  268. 1:length(result.gr),
  269. .fun = function (i) {
  270. rg <- result.gr[i]
  271. overlap <- findOverlaps(rg, anno.gr)
  272. hit <- subjectHits(overlap)
  273. anno_df_sub <- anno_df[hit,]
  274. cpgs <- intersect(anno_df_sub$Name, cpg)
  275. anno_df_sub <- anno_df_sub[match(cpgs, anno_df_sub$Name),]
  277. # cpg in region
  278. cpgs_in_region <- paste(cpgs, collapse = ",")
  280. # UCSC_RefGene_Name
  281. UCSC_RefGene_Name <- paste0(unique(anno_df_sub[,"UCSC_RefGene_Name"]), collapse = ";")
  283. # UCSC_RefGene_Accession
  284. UCSC_RefGene_Accession <- paste0(unique(anno_df_sub[,"UCSC_RefGene_Accession"]), collapse = ";")
  286. # UCSC_RefGene_Group
  287. UCSC_RefGene_Group <- paste0(unique(anno_df_sub[,"UCSC_RefGene_Group"]), collapse = ";")
  289. # Relation_to_Island
  290. Relation_to_Island <- paste0(unique(anno_df_sub[,island_var]), collapse = ";")
  292. df <- data.frame(
  293. num_probes = length(cpgs),
  294. UCSC_RefGene_Name = UCSC_RefGene_Name,
  295. UCSC_RefGene_Accession = UCSC_RefGene_Accession,
  296. UCSC_RefGene_Group = UCSC_RefGene_Group,
  297. Relation_to_Island = Relation_to_Island,
  298. cpgs_in_region = cpgs_in_region
  299. )
  301. df[df == "NA"] <- NA
  303. df
  304. }, .parallel = T
  305. )
  307. cbind(result, anno)
  308. }
  309. # ------------------------------------------------------------------------------
  310. # Function: annotate_enhancer
  311. # ------------------------------------------------------------------------------
  312. # Description:
  313. # Annotates regions with enhancer overlap information using an external enhancer dataset.
  314. #
  315. # Parameters:
  316. # result : Data frame. Regions to be annotated.
  317. # nasser.enhancer.gr: GRanges object. Enhancer regions with associated cell type information.
  318. # cpg : Vector. List of CpG identifiers used in the analysis.
  319. # genome : Character. The genome version ("hg19" or "hg38").
  320. # array : Character. The methylation array type.
  321. #
  322. # Returns:
  323. # A data frame with two additional columns:
  324. # - nasser_is_enhancer: Logical indicating enhancer overlap.
  325. # - nasser_is_enhancer_cell_types: Character string of associated cell types.
  328. annotate_enhancer <- function(result, nasser.enhancer.gr, cpg, genome, array) {
  329. # Create a GRanges object from the result data frame
  330. result.gr <- result %>%
  331. makeGRangesFromDataFrame(start.field = "start",
  332. end.field = "end",
  333. seqnames.field = "seqnames")
  335. # Find overlaps between the result regions and the enhancer regions
  336. hits <- findOverlaps(result.gr, nasser.enhancer.gr) %>% as.data.frame()
  338. # Initialize the enhancer annotation columns
  339. result$nasser_is_enhancer <- FALSE
  340. result$nasser_is_enhancer[unique(hits$queryHits)] <- TRUE
  342. result$nasser_is_enhancer_cell_types <- NA
  343. result$nasser_is_enhancer_cell_types[unique(hits$queryHits)] <- sapply(unique(hits$queryHits), function(x) {
  344. paste(unique(nasser.enhancer.gr$CellType[hits$subjectHits[hits$queryHits %in% x]]), collapse = ",")
  345. })
  347. return(result)
  348. }
  350. annotate_chmm <- function(result, dir.data.aux = dir.data.aux) {
  352. load(file.path(dir.data.aux,"E073_15_coreMarks_segments.rda"))
  354. message("Annotating E073_15_coreMarks_segments")
  356. result$region <- paste0(result$seqnames,":",result$start,"-", result$end)
  357. result$start <- as.numeric(result$start)
  358. result$end <- as.numeric(result$end)
  359. result.gr <- result %>% makeGRangesFromDataFrame(
  360. start.field = "start",
  361. end.field = "end",
  362. seqnames.field = "seqnames"
  363. )
  364. hits <- findOverlaps(result.gr, ChmmModels.gr) %>% as.data.frame()
  365. hits$state <- ChmmModels.gr$state[hits$subjectHits]
  366. hits$region <- result$region[hits$queryHits]
  367. result$E073_15_coreMarks_segments_state <- hits$state[match(result$region,hits$region)]
  368. result$region <- NULL
  370. result
  371. }
  372. # ------------------------------------------------------------------------------
  373. # Function: add_dmr_annotation
  374. # ------------------------------------------------------------------------------
  375. # Description:
  376. # Integrates multiple annotation methods (region, enhancer, coreMarks segmentation, GREAT) for DMRs.
  377. #
  378. # Parameters:
  379. # result : Data frame. DMRs (differentially methylated regions) to be annotated.
  380. # cpg : Vector. List of CpG identifiers used for regional annotation.
  381. # dir.data.aux : Character. Directory path for auxiliary data files (e.g., segmentation, enhancer data).
  382. # region_var : (Optional) Character. Column name in 'result' containing region identifiers.
  383. # If NULL, a region identifier is created using chromosome, start, and end.
  384. # array : Character. The methylation array type (default "EPIC").
  385. # genome : Character. Genome version ("hg19" or "hg38").
  386. # annotate : Vector of characters. Specifies which annotation types to apply. Options include:
  387. # "region", "enhancer", "E073_15_coreMarks_segments", "GREAT".
  388. #
  389. # Returns:
  390. # A data frame with multiple annotation columns added.
  392. add_dmr_annotation <- function(result,
  393. cpg,
  394. dir.data.aux,
  395. region_var = NULL,
  396. array = "EPIC",
  397. genome = "hg19",
  398. annotate = c("region", "enhancer", "E073_15_coreMarks_segments", "GREAT")) {
  400. # Prepare the result data frame for genomic range creation
  402. if(is.null(region_var)) {
  403. result$seqnames <- paste0("chr", result$`#chrom`)
  404. result$region <- paste0(result$seqnames, ":", result$start, "-", result$end)
  405. region_var <- "region"
  406. } else {
  407. region <- str_split(result[[region_var]], ":|-", simplify = T)
  408. result$seqnames <- region[,1]
  409. result$start <- region[,2]
  410. result$end <- region[,3]
  411. }
  412. result$start <- as.numeric(result$start)
  413. result$end <- as.numeric(result$end)
  416. # Annotate with E073_15_coreMarks_segments if selected.
  417. if ("E073_15_coreMarks_segments" %in% annotate) {
  418. load(file.path(dir.data.aux, "E073_15_coreMarks_segments.rda"))
  419. message("Annotating E073_15_coreMarks_segments")
  420. result <- annotate_coreMarks_segments(result, ChmmModels.gr, region_var)
  421. }
  423. # Annotate with GREAT if selected.
  424. if ("GREAT" %in% annotate) {
  425. message("Annotating GREAT")
  426. result <- annotate_great(result, genome)
  427. }
  429. # Annotate with enhancer if selected.
  430. if ("enhancer" %in% annotate) {
  431. # Load enhancer data and filter
  432. data <- readr::read_tsv(file.path(dir.data.aux, "AllPredictions.AvgHiC.ABC0.015.minus150.ForABCPaperV3.txt.gz"),
  433. show_col_types = F)
  434. CellType.selected <- readxl::read_xlsx(file.path(dir.data.aux, "Nassser study selected biosamples.xlsx"),
  435. col_names = FALSE) %>% dplyr::pull(1)
  436. data.filtered <- data %>%
  437. dplyr::filter(CellType %in% CellType.selected) %>%
  438. dplyr::filter(!isSelfPromoter) %>%
  439. dplyr::filter(class != "promoter")
  440. nasser.enhancer.gr <- data.filtered %>%
  441. makeGRangesFromDataFrame(start.field = "start",
  442. end.field = "end",
  443. seqnames.field = "chr",
  444. keep.extra.columns = TRUE)
  446. message("Annotating enhancer")
  447. result <- annotate_enhancer(result, nasser.enhancer.gr, cpg, genome, array)
  448. }
  450. # Annotate with UCSC (island annotation) if selected.
  451. if ("region" %in% annotate) {
  452. message("Annotating region")
  453. result <- annotate_region(result, array, genome, dir.data.aux, cpg, region_var = region_var)
  454. }
  456. return(result)
  457. }
  458. # ===================================================================================================
  459. # Bacon correction
  460. # ===================================================================================================
  461. # Function: bacon_adj
  462. # ------------------------------------------------------------------------------
  463. # Description:
  464. # Adjusts for bias and genomic inflation in EWAS data using the bacon method.
  465. # Supports correction using either z-scores or effect sizes with standard errors.
  466. #
  467. # Parameters:
  468. # data : Data frame. The EWAS results containing statistics to be adjusted.
  469. # est_var : Character. Column name for effect size estimates.
  470. # z_var : Character. Column name for z-scores.
  471. # std_var : Character. Column name for standard errors.
  472. # use_z : Logical. If TRUE, bacon correction is performed using z-scores only; otherwise, effect sizes and SEs are used.
  473. # save : Logical. Whether to save the bacon-corrected data and inflation statistics to files.
  474. # dir.save: Character. Directory path where the output files will be saved.
  475. # prefix : Character. A prefix used for naming the output files.
  476. #
  477. # Returns:
  478. # A list containing:
  479. # - data.with.inflation: The corrected EWAS data frame.
  480. # - bacon.obj : The bacon object from the initial analysis.
  481. # - inflation.stat : A data frame summarizing inflation and bias statistics.
  483. bacon_adj <- function(data, est_var, z_var, std_var,
  484. use_z = F,
  485. save = F,
  486. dir.save = NULL,
  487. prefix = "Framingham"){
  489. ### 1. Compute genomic inflation factor before bacon adjustment
  490. data <- data %>% mutate(
  491. chisq = get(z_var)^2
  492. )
  494. # inflation factor - last term is median from chisq distrn with 1 df
  495. inflationFactor <- median(data$chisq,na.rm = TRUE) / qchisq(0.5, 1)
  496. print("lambda")
  497. print(inflationFactor)
  499. ### 2. bacon analysis
  500. if(use_z){
  501. ### 2. bacon analysis
  502. z_scores <- data[[z_var]]
  503. bc <- bacon(
  504. teststatistics = z_scores,
  505. na.exclude = TRUE,
  506. verbose = F
  507. )
  508. # inflation factor
  509. print("lambda.bacon")
  510. print(inflation(bc))
  511. # bias
  512. print("estimate bias")
  513. print(bias(bc))
  514. print("estimates")
  515. print(bacon::estimates(bc))
  517. ### 3. Create final dataset
  518. data.with.inflation <- data %>% mutate(
  519. zScore.bacon = tstat(bc)[,1],
  520. pValue.bacon.z = pval(bc)[,1],
  521. fdr.bacon.z = p.adjust(pval(bc), method = "fdr"),
  522. ) %>% mutate(z.value = z_scores)
  524. print("o After bacon correction")
  525. print("Conventional lambda")
  526. lambda.con <- median((data.with.inflation$zScore.bacon) ^ 2,na.rm = TRUE)/qchisq(0.5, 1)
  527. print(lambda.con)
  529. # percent_null <- trunc ( bacon::estimates(bc)[1]*100, digits = 0)
  530. # percent_1 <- trunc ( bacon::estimates(bc)[2]*100, digits = 0 )
  531. # percent_2 <- 100 - percent_null - percent_1
  532. bc2 <- bacon(
  533. teststatistics = data.with.inflation$zScore.bacon,
  534. na.exclude = TRUE,
  535. priors = list(
  536. sigma = list(alpha = 1.28, beta = 0.36),
  537. mu = list(lambda = c(0, 3, -3), tau = c(1000, 100, 100)),
  538. epsilon = list(gamma = c(90, 5, 5)))
  539. )
  540. } else {
  541. est <- data[[est_var]]
  542. se <- data[[std_var]]
  545. bc <- bacon(
  546. teststatistics = NULL,
  547. effectsizes = est,
  548. standarderrors = se,
  549. na.exclude = TRUE,
  550. verbose = F
  551. )
  552. # posteriors(bc)
  553. # inflation factor
  554. print("lambda.bacon")
  555. print(inflation(bc))
  556. # bias
  557. print("estimate bias")
  558. print(bias(bc))
  559. print("estimates")
  560. print(bacon::estimates(bc))
  562. ### 3. Create final dataset
  563. data.with.inflation <- data.frame(
  564. data,
  565. Estimate.bacon = bacon::es(bc),
  566. StdErr.bacon = bacon::se(bc),
  567. pValue.bacon = pval(bc),
  568. fdr.bacon = p.adjust(pval(bc), method = "fdr"),
  569. stringsAsFactors = FALSE
  570. )
  572. print("o After bacon correction")
  573. print("Conventional lambda")
  574. lambda.con <- median((data.with.inflation$Estimate.bacon/data.with.inflation$StdErr.bacon) ^ 2,na.rm = TRUE)/qchisq(0.5, 1)
  575. print(lambda.con)
  577. # percent_null <- trunc ( estimates(bc)[1]*100, digits = 0)
  578. # percent_1 <- trunc ( estimates(bc)[2]*100, digits = 0 )
  579. # percent_2 <- 100 - percent_null - percent_1
  580. bc2 <- bacon(
  581. teststatistics = NULL,
  582. effectsizes = data.with.inflation$Estimate.bacon,
  583. standarderrors = data.with.inflation$StdErr.bacon,
  584. na.exclude = TRUE,
  585. priors = list(
  586. sigma = list(alpha = 1.28, beta = 0.36),
  587. mu = list(lambda = c(0, 3, -3), tau = c(1000, 100, 100)),
  588. epsilon = list(gamma = c(99, .5, .5)))
  589. )
  590. }
  592. print("inflation")
  593. print(inflation(bc2))
  594. print("estimates")
  595. print(bacon::estimates(bc2))
  596. data.with.inflation$chisq <- NULL
  598. inflation.stat <- data.frame(
  599. "Inflation.org" = inflationFactor,
  600. "Inflation.bacon" = inflation(bc),
  601. "Bias.bacon" = bias(bc),
  602. "Inflation.after.correction" = lambda.con,
  603. "Inflation.bacon.after.correction" = inflation(bc2),
  604. "Bias.bacon.after.correction" = bias(bc2)
  605. )
  608. if(save){
  609. readr::write_csv(
  610. data.with.inflation,
  611. file.path(dir.save, paste0(prefix, "_bacon_correction.csv"))
  612. )
  614. writexl::write_xlsx(
  615. inflation.stat,
  616. file.path(dir.save, paste0(prefix, "_inflation_stats.xlsx"))
  617. )
  618. }
  619. return(
  620. list(
  621. "data.with.inflation" = data.with.inflation,
  622. "bacon.obj" = bc,
  623. "inflation.stat" = inflation.stat
  624. )
  625. )
  626. }

annotation_and_bacon.R at commit d2cc5c5, under MIT · at the source

Overview

Authors: Wei Zhang1, David Lukacsovich1, Juan I Young2,3, Lissette Gomez3, Michael A Schmidt2,3, Brian W Kunkle2,3, Xi Chen1,4, Eden R Martin2,3, Lily Wang1,2,3,4,5
ORCID iDs: Lily Wang
  1. Division of Biostatistics, Department of Public Health Sciences, University of Miami, Miller School of Medicine, Miami, FL 33136 USA
  2. Dr. John T Macdonald Foundation Department of Human Genetics, University of Miami, Miller School of Medicine, Miami, FL 33136 USA
  3. John P. Hussman Institute for Human Genomics, University of Miami Miller School of Medicine, Miami, FL 33136 USA
  4. Sylvester Comprehensive Cancer Center, University of Miami, Miller School of Medicine, Miami, FL 33136 USA
  5. Soffer Clinical Research Ctr, University of Miami Miller School of Medicine, 1120 NW 14Th St, Miami, FL 33136 USA
Institutions: University of Miami (United States); Sylvester Comprehensive Cancer Center (United States)
Journal: GeroScience, volume 48, issue 3, pages 3185-3203
Dates: received 8 September 2025; accepted 7 March 2026; published online 14 April 2026; in print June 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1007/s11357-026-02195-x · PMID 41975027 · PMCID PMC13356189 · OpenAlex W7154153233
Open access: hybrid, a free copy (OpenAlex)
Status: code verified
Categories: genetics / omics (modality), human (organism), Alzheimer's / dementia (population)
Methods: Statistics, Smoothing, state filtering, decompositions
Keywords: Aging, Alzheimer’s disease, Epigenetics, DNA methylation
MeSH: Aging*, Alzheimer Disease*, DNA Methylation*, Epigenome*, Aged, Aged, 80 and over, Epigenesis, Genetic, Female, Genome-Wide Association Study, Humans, Male (* major topic)
Topic: Epigenetics and DNA Methylation (Molecular Biology, Biochemistry, Genetics and Molecular Biology), according to OpenAlex
Funding: NIA NIH HHS (R01AG062634, R01 AG062634, U01 AG024904); NINDS NIH HHS (R61NS135587, R01 NS128145, R61 NS135587, RF1 NS128145, R01NS128145)
Citations: not cited yet (Europe PMC); 77 references in the paper

Abstract

Aging is the strongest risk factor for Alzheimer’s disease (AD), yet the role of age-associated DNA methylation (DNAm) changes in blood and their relevance to AD remains poorly understood. We performed a meta-analysis of blood DNAm samples from 475 dementia-free subjects aged over 65 years across two independent cohorts, the Framingham Heart Study (FHS) at Exam 9 and the Alzheimer’s Disease Neuroimaging Initiative (ADNI). We adjusted for sex and immune cell-type proportions and corrected batch effects and genomic inflation. Integrative analyses included pathway enrichment, mQTL analysis, colocalization with Alzheimer’s disease and related dementia (ADRD) GWAS summary statistics, brain-blood DNAm correlations, and comparison to independent AD methylation studies. We identified 3758 CpGs and 556 differentially methylated regions (DMRs) consistently associated with chronological age in both cohorts at a 5% false discovery rate. Our pathway enrichment analyses highlighted metabolic regulation and synaptic signaling, processes previously implicated in Alzheimer’s disease. Colocalization with ADRD GWAS summary statistics identified 32 genomic regions consistent with shared genetic signals for DNAm and ADRD risk. Roughly one-third of aging-associated CpGs overlapped CpGs associated with AD or AD neuropathology in external studies. Finally, we prioritized nine promoter CpGs (including those located in PDE1B, ELOVL2, and PODXL2) showing strong positive blood-to-brain methylation concordance and external AD associations, nominating them as candidate blood-based biomarkers. Our study demonstrated that late-life aging signatures in blood DNAm converge on processes implicated in AD and intersect with dementia genetics. A small set of CpGs with blood-brain concordance and external AD support offers promising candidate blood-based biomarkers for future validation.

Supplementary Information: The online version contains supplementary material available at 10.1007/s11357-026-02195-x.

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 6 matches between paragraphs and lines of code.

TransBioInfoLab/AD-Aging-blood-sample-analysis

License: MIT
State: the link answers, verified on 29 September 2026
Evidence: files inventoried
Commit: d2cc5c5d37b2512bdcdc11a402a38a9bd6c7fcdc, 5 June 2025
Languages: R (14)
Size: 18 files, 14 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README, license file, 5 notebooks
Not found: CITATION.cff, environment file, tests, continuous integration, documentation
Tools: tidyverse (11 files), survival (4 files), data.table (2 files), cowplot (1 file), ggplot2 (1 file), ggpubr (1 file), lme4 (1 file), lmerTest (1 file), metafor (1 file)
Availability: 1 check, the latest on 29 September 2026: the link answers
  • 29 September 2026: the link answers
16 files
At the source: github.com/TransBioInfoLab/AD-Aging-blood-sample-analysis

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;
  • 6 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

No dataset and no data link were found in the paper.

Data availability

The ADNI and Framingham Heart Study datasets can be accessed from http://adni.loni.usc.edu and the dbGap database (accession: phs000974.v5.p4). The scripts for the analyses performed in this study are at https://github.com/TransBioInfoLab/AD-Aging-blood-sample-analysis.

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, 29 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 9 authors, 4 keywords, 11 MeSH terms, 2 funders, 75 references.

Cite

This paper

Zhang, W., Lukacsovich, D., Young, J. I., Gomez, L., Schmidt, M. A., Kunkle, B. W., Chen, X., Martin, E. R., & Wang, L. (2026). The aging epigenome: integrative analyses reveal intersection with Alzheimer's disease. GeroScience, 48(3), 3185-3203. https://doi.org/10.1007/s11357-026-02195-x

BibTeX

@article{zhang2026aging,
author = {Zhang, Wei and Lukacsovich, David and Young, Juan I and Gomez, Lissette and Schmidt, Michael A and Kunkle, Brian W and Chen, Xi and Martin, Eden R and Wang, Lily},
title = {{The aging epigenome: integrative analyses reveal intersection with Alzheimer's disease}},
journal = {GeroScience},
year = {2026},
month = apr,
volume = {48},
number = {3},
pages = {3185--3203},
publisher = {Springer},
issn = {2509-2715},
doi = {10.1007/s11357-026-02195-x},
url = {https://doi.org/10.1007/s11357-026-02195-x},
pmid = {41975027},
pmcid = {PMC13356189}
}

RIS

TY - JOUR
AU - Zhang, Wei
AU - Lukacsovich, David
AU - Young, Juan I
AU - Gomez, Lissette
AU - Schmidt, Michael A
AU - Kunkle, Brian W
AU - Chen, Xi
AU - Martin, Eden R
AU - Wang, Lily
TI - The aging epigenome: integrative analyses reveal intersection with Alzheimer's disease
T2 - GeroScience
J2 - Geroscience
PY - 2026
DA - 2026/04/14
VL - 48
IS - 3
SP - 3185
EP - 3203
SN - 2509-2715
PB - Springer
DO - 10.1007/s11357-026-02195-x
UR - https://doi.org/10.1007/s11357-026-02195-x
LA - en
ER -

CSL-JSON

{
"id": "10.1007/s11357-026-02195-x",
"type": "article-journal",
"title": "The aging epigenome: integrative analyses reveal intersection with Alzheimer's disease",
"container-title": "GeroScience",
"author": [
{
"family": "Zhang",
"given": "Wei"
},
{
"family": "Lukacsovich",
"given": "David"
},
{
"family": "Young",
"given": "Juan I"
},
{
"family": "Gomez",
"given": "Lissette"
},
{
"family": "Schmidt",
"given": "Michael A"
},
{
"family": "Kunkle",
"given": "Brian W"
},
{
"family": "Chen",
"given": "Xi"
},
{
"family": "Martin",
"given": "Eden R"
},
{
"family": "Wang",
"given": "Lily"
}
],
"container-title-short": "Geroscience",
"volume": "48",
"issue": "3",
"page": "3185-3203",
"DOI": "10.1007/s11357-026-02195-x",
"PMID": "41975027",
"PMCID": "PMC13356189",
"ISSN": "2509-2715",
"publisher": "Springer",
"URL": "https://doi.org/10.1007/s11357-026-02195-x",
"language": "en",
"issued": {
"date-parts": [
[
2026,
4,
14
]
]
}
}

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.1186/s13073-026-01698-8 [code]
From aging to Alzheimer's disease: concordant brain DNA methylation changes in late life.
Journal: Genome medicine
In common: survival, ggpubr, ggplot2, 1 other tool, Alzheimer's / dementia, genetics / omics, 25 references, author Lily Wang
[2] doi:10.1002/trc2.70257 [code]
Blood DNA methylation signature of cognitive reserve moderates the association between CSF tau pathology and memory in prodromal Alzheimer's disease.
Journal: Alzheimer's & dementia (New York, N. Y.)
In common: lmerTest, lme4, ggpubr, 2 other tools, Alzheimer's / dementia, genetics / omics, 2 references, author Lily Wang
[3] doi:10.1038/s41586-026-10877-x [code]
Human brain organoids record the passage of time over multiple years.
Journal: Nature
In common: lmerTest, lme4, cowplot, 4 other tools, genetics / omics, 2 references
[4] doi:10.1038/s41467-026-73865-9 [code]
Histamine shapes the neurocomputational dynamics of human learning.
Journal: Nature communications
In common: metafor, lmerTest, lme4, 5 other tools
[5] doi:10.1038/s44400-026-00100-z [code]
Peripheral blood microarray-based transcriptomic and epigenetic analyses identify immune, inflammation, and metabolic dysregulation in Alzheimer's disease.
Journal: NPJ dementia
In common: survival, lme4, ggpubr, 2 other tools, Alzheimer's / dementia, genetics / omics, 2 references
[6] doi:10.1186/s13195-026-02036-1 [code]
Genetic drivers of progression in Alzheimer's disease are distinct from disease risk.
Journal: Alzheimer's research & therapy
In common: metafor, lmerTest, lme4, 4 other tools, Alzheimer's / dementia, genetics / omics
[7] doi:10.1002/ece3.73881 [code]
Epigenetic Aging in Brain Tissue of the Self-Fertilizing Vertebrate, &lt;i&gt;Kryptolebias marmoratus&lt;/i&gt;.
Journal: Ecology and evolution
In common: ggpubr, ggplot2, tidyverse, genetics / omics, 5 references
[8] doi:10.1038/s41588-026-02722-8 [code]
A multiancestry polygenic risk score for Alzheimer's disease is associated with cognitive decline and neuropathological hallmarks in diverse populations.
Journal: Nature genetics
In common: survival, lmerTest, lme4, 3 other tools, Alzheimer's / dementia, genetics / omics, 1 reference
[9] doi:10.1038/s41467-026-77170-3 [code]
DNA methylation profiling identifies long-range epigenetic silencing of clustered protocadherins as a key determinant of meningioma progression.
Journal: Nature communications
In common: survival, lmerTest, lme4, 4 other tools, genetics / omics
[10] doi:10.1093/bioinformatics/btag592 [code]
Network-based stratification of allele-specific expression reveals patient subgroups in Huntington's disease.
Journal: Bioinformatics (Oxford, England)
In common: survival, lme4, cowplot, 4 other tools, genetics / omics

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.