Association of epigenetic age acceleration with MRI biomarkers of aging and Alzheimer's disease neurodegeneration.

1. #' Generate the Dunedin Methylation Pace of Aging Scores!
2. #'
3. #' \code{PACEProjector} returns the Dunedin Pace of Aging Methylation Scores
4. #'
5. #' @param betas A numeric matrix containing the percent-methylation for each probe. Missing data should be 'NA's. The rows should be probes, with the probe ID as the row name, and the columns should be samples, with sample names as the column name.
6. #' @param proportionOfProbesRequired (default: 0.8). This value specificies the threshold for missing data (see description for more details on how missing data is handled). Note: that if the function detects EPICv2 data based on rownames that contain the replicate suffix, the proportion of probes required will automatically be lowered to 0.7.
7. #' @return A list of mPACE values. There will be one element in the list for each mPACE model. Each element will consist of a numeric vector with mPACE values. The names of the values in the vector will be the sample names from the 'betas' matrix. The output is equivalent to the pace of aging (in years) expected over a 1-year period.
8. #' @details This function returns the Dunedin Methylation Pace of Aging scores for methylation data generated from either the Illumina 450K array or the Illumina EPIC array. The Age45 score is one that has been trained on data based on 3 waves of collection (26, 38, and 45). The manuscript is currently in preparation, but has been shown to be more accurate than the Age38 score.
9. #' Missing data handled in two different ways (and the threshold for both is set by the 'proportionOfProbesRequired' parameter). First, if a sample is missing data for more probes than the threshold, the sample will get an NA back for a score. If a particular probe is missing fewer samples than the threshold, then missing data is set to the mean in the provided 'betas' matrix. If a probe is missing more samples than the threshold, then all samples in the 'betas' matrix have their value replaced with the mean of the training data for that particular model.
10. #' Because of how we handle missing data, it is recommended that entire cohorts be run at once as a large 'betas' matrix.
11. #' @examples
12. #' PACEProjector(betas)
15. PACEProjector = function( betas, proportionOfProbesRequired=0.8 ) {
16. requireNamespace("preprocessCore")
18. if( any(grepl("TC|BC", rownames(betas))) )
19. {
20. # Print message to user
21. print("This looks like EPICv2 array data. If EPICv2, DunedinPACE will lower the proportionOfProbesRequired to 0.7 and proceed with missing probes present on 450k and EPICv1. Averaging the opposite strand replicates may take a bit more time.")
23. # Proceed with averaging: Code taken from ENmix package.
24. cgid = sapply(strsplit(rownames(betas), split = "_"), unlist)[1, ]
25. dupcg = unique(cgid[duplicated(cgid)])
26. betas2 = betas[cgid %in% dupcg, ]
27. cid = sapply(strsplit(rownames(betas2), split = "_"), unlist)[1, ]
28. betas2 = aggregate(betas2, by = list(cid), FUN = function(x) mean(x, na.rm = TRUE))
29. rownames(betas2) = betas2[, 1]
30. betas2 = as.matrix(betas2[, -1])
31. betas = betas[!(cgid %in% dupcg), ]
32. rownames(betas) = sapply(strsplit(rownames(betas), split = "_"), unlist)[1, ]
33. betas <-rbind(betas, betas2)
35. # Set proportionOfProbesRequired to 0.7 if condition is TRUE
36. proportionOfProbesRequired <- 0.7
37. }
39. # loop through models
40. model_results <- lapply(mPACE_Models$model_names, function(model_name) {
41. # make sure it has been converted to a matrix
42. if( !is.numeric(as.matrix(betas)) ) { stop("betas matrix/data.frame is not numeric!") }
43. probeOverlap <- length(which(rownames(betas) %in% mPACE_Models$model_probes[[model_name]])) / length(mPACE_Models$model_probes[[model_name]])
44. probeOverlap_background <- length(which(rownames(betas) %in% mPACE_Models$gold_standard_probes[[model_name]])) / length(mPACE_Models$gold_standard_probes[[model_name]])
45. # make sure enough of the probes are present in the data file
46. if( probeOverlap < proportionOfProbesRequired | probeOverlap_background < proportionOfProbesRequired ) {
47. result <- rep(NA, ncol(betas))
48. names(result) <- colnames(betas)
49. result
50. } else {
51. # Work with a numeric matrix of betas
52. betas.mat <- as.matrix(betas[which(rownames(betas) %in% mPACE_Models$gold_standard_probes[[model_name]]),])
53. # If probes don't exist, we'll add them as rows of values based on their mean in the gold standard dataset
54. probesNotInMatrix <- mPACE_Models$gold_standard_probes[[model_name]][which(mPACE_Models$gold_standard_probes[[model_name]] %in% rownames(betas.mat) == F)]
55. if( length(probesNotInMatrix) > 0 ) {
56. for( probe in probesNotInMatrix ) {
57. tmp.mat <- matrix(0, nrow=1, ncol=ncol(betas.mat))
58. rownames(tmp.mat) <- probe
59. colnames(tmp.mat) <- colnames(betas.mat)
60. tmp.mat[probe,] <- rep(mPACE_Models$gold_standard_means[[model_name]][probe], ncol(tmp.mat))
61. betas.mat <- rbind(betas.mat, tmp.mat)
62. }
63. }
65. # Identify samples with too many missing probes and remove them from the matrix
66. samplesToRemove <- colnames(betas.mat)[which(apply(betas.mat, 2, function(x) { 1 - ( length(which(is.na(x))) / length(x) ) < proportionOfProbesRequired}))]
67. if( length(samplesToRemove) > 0 ) {
68. betas.mat <- betas.mat[,-which(colnames(betas.mat) %in% samplesToRemove)]
69. }
70. if(ncol(betas.mat) > 0) {
71. # Identify missingness on a probe level
72. pctValuesPresent <- apply( betas.mat, 1, function(x) { 1 - (length(which(is.na(x))) / length(x)) } )
73. # If they're missing values, but less than the proportion required, we impute to the cohort mean
74. probesToAdjust <- which(pctValuesPresent < 1 & pctValuesPresent >= proportionOfProbesRequired)
75. if( length(probesToAdjust) > 0 ) {
76. if( length(probesToAdjust) > 1 ) {
77. betas.mat[probesToAdjust,] <- t(apply( betas.mat[probesToAdjust,], 1 , function(x) {
78. x[is.na(x)] = mean( x, na.rm = TRUE )
79. x
80. }))
81. } else {
82. betas.mat[probesToAdjust,which(is.na(betas.mat[probesToAdjust,]))] <- mean(betas.mat[probesToAdjust,], na.rm=T)
83. }
84. }
85. # If they're missing too many values, everyones value gets replaced with the mean from the Dunedin cohort
86. if( length(which(pctValuesPresent < proportionOfProbesRequired)) > 0 ) {
87. probesToReplaceWithMean <- rownames(betas.mat)[which(pctValuesPresent < proportionOfProbesRequired)]
88. for( probe in probesToReplaceWithMean ) {
89. betas.mat[probe,] <- rep(mPACE_Models$model_means[[model_name]][probe], ncol(betas.mat))
90. }
91. }
93. # Normalize the matrix to the gold standard dataset
94. betas.norm <- preprocessCore::normalize.quantiles.use.target(betas.mat, target=mPACE_Models$gold_standard_means[[model_name]])
95. rownames(betas.norm) <- rownames(betas.mat)
96. colnames(betas.norm) <- colnames(betas.mat)
97. # Calculate score:
98. score = mPACE_Models$model_intercept[[model_name]] + rowSums(t(betas.norm[mPACE_Models$model_probes[[model_name]],]) %*% diag(mPACE_Models$model_weights[[model_name]]))
99. names(score) <- colnames(betas.norm)
100. if( length(samplesToRemove) > 0 ) {
101. score.tmp <- rep(NA, length(samplesToRemove))
102. names(score.tmp) <- samplesToRemove
103. score <- c(score, score.tmp)
104. }
105. score <- score[colnames(betas)]
106. score
107. } else {
108. result <- rep(NA, ncol(betas.mat))
109. names(result) <- colnames(betas.mat)
110. result
111. }
112. }
113. })
114. names(model_results) <- mPACE_Models$model_names
115. model_results
116. }

PACEProjector.R at commit 4b56998, under GPL-3.0 · at the source