Blood DNA markers and exposure risk scores accurately predict PTSD in military and civilian groups

In order to maximize the available data from which to develop risk scores using machine learning approaches, we created a discovery cohort comprised of 1226 individuals drawn from five cohorts (Table 1). Two of these cohorts are civilian— Detroit Neighborhood Health Study (DNHS) and Grady Trauma Project (GTP), and three cohorts are military— Army Study to Assess Risk and Resilience in Servicemembers (Army STARRS), Marine Resilience Study (MRS I&II), and Prospective Research in Stress-related Military Operations (PRISMO). Details about each cohort are given in the supplementary file. The overall workflow of the pre-processing and methods combining data from the five cohorts is shown in supplementary file (Figure S1).

DNAm from whole blood was measured using the Illumina MethylationEPIC BeadChip following the manufacturer's recommended protocol. Raw DNAm β values were obtained, and a sex check was conducted using the minfi R package [20] to eliminate any sex-discordant samples. Quality control (QC) was performed on each cohort separately, using a standardized pipeline as previously described [21]. A total of 818,691 probes passed QC. Normalization was carried out using the single-sample Noob (ssNoob) method in the minfi R package [20]. Furthermore, ComBat adjustment was performed, using an empirical Bayesian framework implemented in the SVA R package [22, 23] to reduce the likelihood of bias due to known batch effects (chip and position), while preserving the variation for age, sex (if applicable), and PTSD. The resulting QC'd data was used in subsequent analyses.

All the discovery cohorts had a small percentage of missing values (ranging from 0.002 to 0.03%). As the machine learning models require complete data, we used the mean method—a common and simple imputation technique to impute the missing data while maintaining the distribution of the data [44, 45]. We then adjusted the DNAm data for potential confounding factors, including cell composition, ancestry, smoking score, sex (if applicable), and age, for models 1 and 2 (described below). The adjustment was made for each CpG by regressing out all the covariates using linear regression and then replacing the values of CpG with the corresponding residuals [46]. For model 3 (described below), we also regressed out the two exposure variables of interest, cumulative trauma and childhood trauma, in addition to the covariates used in models 1 and 2. This was done separately for each cohort to account for any differences related to exposure variables in individual cohorts.

We used SelectKBest in Scikit-learn [51], a univariate feature selection approach. This method computes ANOVA F-values based on univariate statistical tests to identify the best features in relation to a particular phenotype. We identified the most important features from DNAm and exposure variables (in cases of Models 1 & 2) based on the rank order of the features’ association with PTSD. For Model 3, we selected features solely from DNAm data. The feature selection process was repeated 500 times, ranging from 10 to 5000 features with a 10-feature increment each time to determine the optimal feature set for the Elastic Net model best accuracy. As different studies/cohorts used different instruments to measure cumulative trauma and childhood trauma, we normalized the data using a min–max scale that ranged from [1].

In order to identify the best model to classify PTSD and determine risk scores, we trained three popular machine learning models —Random Forest, Lasso, and Elastic Net on 75% of the data, and then tested them on the remaining 25% using the Scikit-learn [51] framework. We also conducted a tenfold cross-validation on training and testing data to evaluate the effectiveness of the models (Figure S1.1). After selecting the most accurate machine learning model, which was evaluated based on the methylation and exposure variable dataset, we used important features identified during the feature selection process to classify PTSD. Following covariate adjustment of the two additional exposure variables, we re-ran the feature selection process to identify important features for Model 3 (described below). Performance of the models was assessed using accuracy, precision, recall, f1-score and area under the curve (AUC) metrics.

Risk scores are the weighted sum of the important features. Using feature weights (i.e. effect sizes) from training data (75%), we created risk scores using discovery cohort test data (25%), in order to test for an association between risk scores and PTSD. Model 1 contributed to the development of exposure and methylation risk scores (eMRS), whereas Model 2 provided methylation-only risk scores (MoRS). Finally, Model 3 led to the creation of methylation-only risk scores with adjusted exposure variables (MoRSAE).

A logistic model was employed to test for an association between risk scores (eMRS, MoRS and MoRSAE) and PTSD, and the Nagelkerke approach was used to assess the models' resulting R-Squared (R²) values. For all analyses, a Wilcoxon rank-sum test was used to assess differences in risk scores between cases and controls. To assess the direction of effect and strength of association among study variables in the both discovery and independent cohorts, Pearson's and point-biserial correlation was used, as appropriate.

To validate the risk scores, we tested their ability to distinguish those with vs. without PTSD in four independent, external cohorts using the same pre-processing and covariate adjustment pipeline as in the discovery cohort. Brief descriptions of the external cohorts (NCPTSD-TRACTS, BEAR, DCHS and PROGrESS) are provided in Supplementary File 1. We utilized weights from significant features identified in models 1, 2, and 3 of the discovery cohort to generate risk scores (i.e., eMRS, MoRS, and MoRSAE) in the external cohorts. Similar to the discovery cohorts, we conducted Pearson and Point-Biserial correlation tests, association tests using logistic regression model, and Wilcoxon rank-sum tests on external cohorts.

In cohorts with available pre-deployment data, a logistic model was used to predict post-deployment PTSD using risk scores calculated from pre-deployment DNAm data and exposure data (Army STARRS, MRS Iⅈ n = 276). Note that these participants had their post-deployment DNAm data included in the discovery cohort analyses described above.

To investigate the biological significance of the important CpGs identified in the feature selection step, we performed Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis using missMethyl [52]. Gene ontologies and KEGG pathways that reached a nominal significance level of p < 0.05 were considered important.

Table 1 provides a summary of the demographic characteristics and clinical information of all participants (n = 1226) in the discovery cohort with current PTSD. More information about cumulative and childhood trauma is provided in Table S1. A slight majority of participants were male (n = 629). Two cohorts, DNHS and GTP, were comprised mostly of African Americans, while the remaining three cohorts were predominantly of European ancestry. In all cohorts, a significant difference in PTSD symptom severity was observed between cases and controls (p < 0.05). With the exception of Army STARRS, childhood trauma also demonstrated a significant difference between PTSD cases and controls (p < 0.05) in all cohorts. Finally, a significant difference was observed in cumulative trauma between cases and controls in DNHS and GTP (p < 0.001).

We developed three different risk scores with the goal of distinguishing those with vs. without PTSD using machine learning approaches. Our first model, eMRS, included both exposure and DNAm variables and identified 3730 features (3728 CpGs, cumulative trauma, and childhood trauma) as important in the discovery cohort. Using these 3730 features, Elastic Net approaches were employed to achieve the best accuracy (92%; Fig. 1), precision (91%), recall (87%), and f1-score (89%); Table 2 (See Fig. S2 for AUCs with Lasso and Random Forest approaches). The eMRS significantly predicted PTSD (beta = 2.64, p < 0.001), R² = 0 0.70), with higher eMRS values in PTSD cases than controls (p < 0.001; Fig. 2A, left plot). Our second MoRS model, based solely on the 3728 methylation features in model 1, accurately classified PTSD with 89% accuracy and had an AUC of 95% (Fig. 3; Table 2). Additionally, the precision, recall, and f1-score were at 86%, 83%, and 84%, respectively, as shown in Table 2. As with eMRS, the MoRS significantly predicted PTSD (beta = 2, p < 0.001, R² = 0.54) and had higher MoRS values in cases vs controls (p < 0.001) (Fig. 2A, middle plot). Our third and final model (i.e., MoRSAE), which used DNAm data adjusted for the two exposure variables as well as the other covariates in models 1 and 2, identified 4150 significant features that classified PTSD with 84% accuracy and an AUC of 89% (Fig. 4, with precision, recall, and f1-score at 80%, 77%, and 78%, respectively (Table 2). As with models eMRS and MoRS, MoRSAE significantly predicted PTSD (beta = 1.20, p < 0.001, R² = 0.36) and had significantly (p < 0.001) different, and higher, MoRSAE in PTSD cases vs. controls (Fig. 2A, right plot). In summary, while all three models produced risk scores that significantly predicted PTSD in the test dataset, and showed higher scores in aggregate between cases and controls, there was a decline in effect size (b) and explanatory power (R²) such that eMRS > MoRS > MoRSAE.

A significant positive point-biserial correlation between eMRS and current PTSD was observed (ρ = 0.72, p < 0.001; Figure S3). Cumulative trauma (ρ = 0.40, p < 0.001) and childhood trauma (ρ = 0.57, p < 0.001) also showed a positive and significant correlation with eMRS. Notably, there was also a significant and positive point-biserial correlation (ρ = 0.62, p < 0.001) between MoRS and PTSD, significant and positive correlation between cumulative trauma and MoRS (ρ = 0.16, p < 0.01) and childhood trauma and MoRS (ρ = 0.169, p < 0.01) (Figure S3). In contrast, while we observed a significant (p < 0.001) and positive point-biserial correlation (ρ = 0.49) between MoRSAE and PTSD (Figure S3), we observed a negative correlation between MoRSAE and cumulative trauma (ρ = -0.13, p = 0.02) and childhood trauma (ρ = -0.12, p = 0.03), respectively.

We conducted external validation on risk scores from the three different models across four external cohorts— NCPTSD-TRACTS, BEAR, DCHS and PROGrESS. The NCPTSD-TRACTS cohort demonstrated a noticeable distinction (p < 0.05) in childhood trauma, but not in cumulative trauma (Table S1) between cases and controls. Similar to the discovery cohorts, the BEAR cohort exhibited a significant difference in both cumulative trauma and childhood trauma when comparing cases and controls. The DCHS cohort, on the other hand, only showed a significant difference in cumulative trauma, while the PROGrESS cohort did not display any significant difference in trauma variables between cases and controls.

The eMRS significantly predicted PTSD in one external cohort, BEAR (beta = 0.6839, p = 0.006) (Table S2); in this cohort, there was also a significant correlation (ρ = 0.24, p = 0.003) between eMRS and PTSD (Figure S4) and a significant difference in eMRS between PTSD cases and controls (p = 0.02, Figure S5). The eMRS did not significantly predict PTSD in any of the other three independent cohorts; however, the correlation between eMRS and PTSD showed the same (i.e., positive) direction of effect in the NCPTSD-TRACTS (beta = 0.0598, p = 0.35), PROGrESS (beta = 0.1141, p = 0.53) and DCHS (beta = 0.0631, p = 0.81) cohorts (Figures S6-S11). For model 2, the MoRS did not significantly predict PTSD in any external cohort (NCPTSD-TRACTS: beta = -0.0977, p = 0.28; BEAR: beta = 0.0239, p = 0.93; PROGrESS: beta = 0.2156, p = 0.52; DCHS: beta = 0.3739, p = 0.37). On the other hand, for model 3, the MoRSAE approached significance in association with PTSD in the NCPTSD-TRACTS cohort (beta = -0.1707, p = 0.05) and had significant difference in risk scores between cases and controls (p = 0.018) (Figure S7); however, the direction of effect was opposite to that observed in the discovery cohort.

A compelling feature of risk scores is their ability to predict future disease risk. In our data, we were able to test the predictive ability of the MRS derived from our diagnostic/classification models on prospective risk of PTSD in two of our pre-deployment military cohorts, MRS I&II and Army STARRS (with data from the two cohorts analyzed together). MRS were calculated using “unseen” DNAm data from a pre-deployment timepoint, i.e. using DNAm data not included in the discovery cohort. All three models significantly predicted future PTSD based on risk scores calculated with pre-deployment data (eMRS beta = 1.92, p < 0.001, R² = 0.53; MoRS beta = 1.99, p < 0.001, R² = 0.46; and MoRSAE beta = 1.77, p < 0.001, R² = 0.47) and had significant difference in risk scores between individuals who developed PTSD and those who did not (Figs. 5, 6, 7).

Gene ontology (GO) analysis on the set of 3728 CpGs from models 1 and 2 revealed 403 nominally significant GO terms; among the 4150 important CpGs from Model 3, 382 nominally significant GO terms were identified. There were 115 GO terms common between models, including regulation of muscle adaptation, positive regulation of autophagy of mitochondrion, and sucrose metabolic process. Additionally, the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis identified 47 pathways for models 1 and 2 and 25 pathways for model 3 at p < 0.05 (list of GO and KEGG terms are provided in Supplementary File 2). Further, 14 pathways were common in models 1 and 2, and model 3, including, HIF-1 signaling pathway, mTOR signaling pathway, Insulin signaling pathway and Galactose metabolism. None of the GO terms or KEGG pathways passed the multiple hypothesis correction test.

All participants included in analyses of the discovery and external validation cohorts provided informed consent for their participation. Details on the specific IRB approval associated with each study are provided in the Supplementary Material.

Not applicable.

Murray B. Stein has in the past 3 years received consulting income from Acadia Pharmaceuticals, Aptinyx, atai Life Sciences, BigHealth, Biogen, Bionomics, BioXcel Therapeutics, Boehringer Ingelheim, Clexio, Delix Therapeutics, Eisai, EmpowerPharm, Engrail Therapeutics, Janssen, Jazz Pharmaceuticals, NeuroTrauma Sciences, PureTech Health, Sage Therapeutics, Sumitomo Pharma, and Roche/Genentech. Dr. Stein has stock options in Oxeia Biopharmaceuticals and EpiVario. He has been paid for his editorial work on Depression and Anxiety (Editor-in-Chief), Biological Psychiatry (Deputy Editor), and UpToDate (Co-Editor-in-Chief for Psychiatry). He has also received research support from NIH, Department of Veterans Affairs, and the Department of Defense. He is on the scientific advisory board for the Brain and Behavior Research Foundation and the Anxiety and Depression Association of America. Dr. Chia-Yen Chen is an employee of Biogen. Dr. Nikolaos P. Daskalakis has served on scientific advisory boards for BioVie Pharma, Circular Genomics and Sentio Solutions for unrelated work. Dr. Nicole R. Nugent is a member of the scientific advisory board for Ilumivu. Dr. Sheila Rauch support from Wounded Warrior Project (WWP), Department of Veterans Affairs (VA), National Institute of Health (NIH), McCormick Foundation, Tonix Pharmaceuticals, Woodruff Foundation, and Department of Defense (DOD). Dr. Rauch also receives royalties from Oxford University Press and American Psychological Association Press. Dr Ressler reported receiving personal consulting fees from Sage Therapeutics, Senseye, Boerhinger Ingelheim, Jazz Pharmaceuticals, and Acer, Inc. and a sponsored research grant from Alto Neuroscience outside the submitted work.

Owing to military cohort data sharing restrictions, data from MRS I& II, Army STARRS, PRISMO, and NCPTSD-TRACTS cannot be publicly posted. For other cohorts, individual-level data from the cohorts or cohort-level summary statistics may be made available to researchers following an approved analysis proposal through the PGC Post-traumatic Stress Disorder EWAS group with agreement of the cohort PIs. For additional information on access to these data, including PI contact information for the contributing cohorts, please contact the corresponding author.