Integrated Analysis of miRNA and mRNA Expression Profiles Reveals the Molecular Mechanism of Posttraumatic Stress Disorder and Therapeutic Drugs

Introduction

Posttraumatic stress disorder (PTSD) often occurs after exposure to trauma. The PTSD patients often re-experience traumatic events in the form of nightmares, flashbacks, or intrusive memories. Patients also develop persistent avoidance of talking or thinking about the traumatic event, show negative cognitive and mood changes, as well as emotional and physical reactions.¹ Around 20% of military members suffer from PTSD, which is also common among the general population, affecting approximately 7% of adult Americans.² In China, the prevalence rate of PTSD was reported to be as high as 62.8% (257 of 430) one month after the 2008 Wenchuan earthquake, and the symptoms of 3 years of PTSD were also reported, ranging from 8.8% to 22.7%.³ In addition to the exacerbating psychological symptoms and negative effects on the quality of life, multiple studies have shown that patients with PTSD have an increased risk of hypertension and cardiovascular disease.^4–6 Genetic and molecular studies on PTSD have discovered some potential genes, which have also been confirmed by other related studies. Studies using candidate-based approaches have focused more on biologically indicative genes. However, high-throughput technology has greatly facilitated the screening of potential genome-wide genes, which have already been widely used in PTSD research.⁷ The molecular mechanism for the development of PTSD is still unknown, but studies from multi-omics data have revealed that changes are mainly involved in the disruption of glucocorticoid signaling and inflammatory systems, and these alterations occur at the level of gene expression.⁸

Genome-wide transcriptomic studies provide enormous evidence of peripheral blood transcriptional markers for PTSD; however, many limitations still exist. Previous studies based on transcriptome data have shown the important role of genes such as ATP6APIL and DSCAM,⁹ FKBP5 and STAT5B^10–12 and the biological pathways involved, as well as the subsequent cellular processes and immune-related pathway¹³ were significantly enriched in differentially expressed genes. Investigations at the system-level regarding the behaviors and properties of these genes revealed important roles of not only the modules of innate immune and interferon signaling transduction but also the modules of hemostasis and wound responsiveness.¹⁴

As the amount of publicly available high-throughput data in global databases continues to grow, an open question arises: How do we appropriately use this large-scale data to achieve a comprehensive understanding of disease at the molecular level? The Random Forest algorithm is a common machine learning algorithm often used as a routine bioinformatics protocol for disease and biology research and plays a key role in various areas of human life, especially it provides methods for the diagnosis and prognosis of human diseases.^15–17 Several studies have used a variety of biomarkers to establish diagnostic or prognostic models for clinical patients.^18,19 However, there are few diagnosis and prognostic analyses of PTSD.

In the current study, we integrated transcriptomic data of PTSD from public databases to identify differentially expressed genes at the transcriptome level. A miRNA-mRNA regulatory network was established based on the miRNA-mRNA interaction database. Drugs interacting with PTSD DEGs were further identified by searching the drug-gene interaction database, which may provide helpful information for the treatment of PTSD. Based on PTSD DEGs, we constructed a diagnostic model that showed good performance for the classification of PTSD and control samples. The design of this study is illustrated in Figure 1.

Figure 1 Workflow of the study. Abbreviations: PTSD, posttraumatic stress disorder group; DEmRNA, differentially expressed mRNAs; DEmiRNA, differentially expressed miRNAs; GO, Gene Ontology analysis; KEGG pathway analysis, Kyoto Encyclopedia of Genes and Genomes pathway analysis.

Materials and Methods

Data Preparation

The post-traumatic brain injury gene expression profiles and the corresponding clinical datasets were obtained from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/). Dataset with next-generation high-throughput sequencing platform were only included to avoid the effect of different platforms on data analysis, resulting in a total of 182 PTSD patients and 151 control samples from five datasets. Raw read counts of the datasets GSE64813,¹⁴ GSE89866, GSE87768²⁰ and GSE109409²¹ were downloaded and preprocessed before differential expression analysis. For GSE63878,¹⁴ the raw CEL data were downloaded and normalized by the robust multichip average (RMA) method. The criteria for selecting these datasets are as follows: 1) all blood samples; 2) all case-control data sets; 3) for differential expression analysis, the three datasets GSE64813, GSE89866, and GSE87768 are second-generation sequencing data and are the same platform.

Information on the platforms and number of samples of each dataset is provided in Table 1.

Table 1 Dataset Used in This Study

Differential Analysis of Gene Expression

Datasets GSE64813 and GSE89866 performed on the same platform of Illumina high-throughput were used to screen differentially expressed mRNAs (DEmRNAs). Dataset GSE87768 was also performed on the same platform, including six case – control samples matched in age and gender, and were used to identify differentially expressed miRNAs (DEmiRNAs). DESeq2 package²² (v1.18.1) was applied to screen the DEGs (DEmRNAs and DEmiRNAs) using the raw read count data, according to the R code shown in the manual. miRNAs and mRNAs with a mean of reads < 1 were excluded. DEGs referred to miRNAs and mRNAs with an adjusted P<0.05, using Benjamini-Hochberg false discovery rate.

miRNA-mRNA Interaction Analysis

The miRNA targets and interactions between selected DEmiRNAs and DEmRNAs were obtained using the R package multMiR²³ (v1.2.0) and the reference manual. We considered the genes detected as miRNA targets only in the validated database available in the package. STRING database²⁴ was used to predict interactive proteins using DEmRNAs as a query. The PPI network with a combined score > 0.4 was considered to be a significant result. The biological network of the protein interactions and miRNA-mRNA regulatory network was visualized using the open source Cytoscape²⁵ (v3.7.1).

Gene Ontology and Pathway Enrichment Analysis

Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) pathway enrichment analyses were carried out using GeneCodis,²⁶ a web-based tool that allows for the simultaneous evaluation of annotations from various sources to depict the potential biological functions of the DEmRNAs. For this purpose, prior to computing p-values, the a-priori algorithm was used to extract annotations frequently co-occurring in the analyzed genes.

Construction of PTSD Diagnostic Risk Model

The GSE64813 and GSE89866 datasets were merged because of their small number of samples. Batch effect was removed by using Bioconductor “combat”²⁷ and no difference is seen in both PTSD and control samples between these 2 datasets after batch correction. The least absolute shrinkage and selection operator (LASSO) regression model, and random forest algorithm were used to further identify the important DEmRNAs with the ability to differentiate control samples and PTSD.²⁸ For the random forest algorithm, we first set mtry = 1–50 and ntree = 500, to search for the best mtry parameter with the lowest error rate. The best ntree parameter is then selected by fixing the mtry parameter. Finally, all DEmRNs were ordered according to their importance, and only DEmRNs with cumulative importance > 95% were kept. For the LASSO regression model, 10-fold cross-validation was performed to determine the optimal lambda parameter, which corresponds to the lowest error rate. The overlapping DEmRNAs screened by the above two methods were used to establish a logistic diagnostic risk model for PTSD. The merged dataset was randomly grouped into a validation and training set with equal proportions. Additionally, in order to avoid the bias of a single random grouping, we repeated this process 100 times and calculated the mean AUC, specificity, and sensitivity. This research was conducted strictly following the statement of transparent reporting of a multivariable prediction model for individual prognosis or diagnosis (TRIPOD).²⁹

Drug Target Analysis

Based on searches of the Drug-Gene Interaction database (DGIdb), potential therapeutic targets of DEmRNAs were identified. The DGIdb³⁰ database contains the interaction information of more than 10,000 drugs and 40,000 specific genes. Currently, the latest version of the DGIdb consists of more than 56,309 interaction claims collected from 30 sources. Such storage of information greatly expands on existing anti-neoplastic drug–gene interactions and the catalogue of druggable genes in the DGIdb.

Statistical Analysis

Statistical analyses in this study were conducted using the R software³¹ (version 3.6.1). Group comparisons were performed using the Mann-Whitney U-test for variables showing an abnormal distribution. For continuous variables, an independent t-test was conducted on normally distributed variables. The model specificity and sensitivity were analyzed with a receiver operating characteristic (ROC) curve and quantified by the area under the ROC curve³² (AUC) 27. All statistical analyses were two-sided. Statistical significance was set at p <0.05. We used “*” to stand for p < 1e-5, “**” for p < 0.01, and “***” for p < 0.05 in this study.

Results

Identification of DEmRNAs and miRNAs

The GSE64813 and GSE89866 datasets identified 669 and 4856 differentially expressed mRNAs in the PTSD case and control samples, respectively (Figure 2A and B). A total of 65 DEmRNAs were obtained in both datasets, which accounted for 9.72% and 1.17% of the total DEmRNAs, respectively (Figure 2C). Further analysis of the common DEmRNAs revealed that 19 genes showed a significant upregulation, and 46 genes were sharply downregulated in PTSD samples (Figure 2D and E), indicating that the activation and inhibition of some biological processes regulated by these DEmRNAs may be important molecular mechanisms for the formation of PTSD. The Genotype-Tissue Expression (GTEx) project allows users to access the transcriptomic data of healthy tissues collected from autopsies.³³ The utility of GTEx databases revealed a different expression pattern of these overlapped DEmRNAs between healthy brain tissues and other healthy tissues (Figure 2F), which implied a high degree of similarity between gene expression in peripheral blood and that in brain tissue. The expression of five miRNAs expressions were significantly upregulated in the PTSD group compared to the control group (Figure 3A). By searching the multiMiR database, a total of 8348 miRNA-targets connections were obtained, including four DEmiRNAs and 10 DEmRNAs screened (Figure 3B).

Figure 2 Expression profiles of differentially expressed genes (DEmRNAs). (A and B) Volcano map of DEmRNAs expression levels between control samples and PTSD in GSE64813 and GSE89866. The red points represent up-regulated DEmRNAs with p-value < 0.05; the green points are down-regulated DEmRNAs with p < 0.05. (C) Overlapped DEmRNAs in GSE64813 and GSE89866. (D and E) Hierarchical clustering of the overlapped DEmRNAs in GSE64813 and GSE89866. Each row in the heatmap represents a gene, and each column represents a sample. The color scale at the right of the heatmap represents the raw Z-score ranging from blue (low expression) to red (high expression). (F) Expression heatmap of overlapped DEmRNAs in GTEx normal sample tissues.

Figure 3 Paired differential expression analyses between PTSD and control. (A) Volcano map of DEmiRNAs expression levels between PTSD and control samples in GSE87768. The red points represent up-regulated DEmRNAs with p < 0.05. (B) Barplot of DEmiRNAs and targets interaction frequency annotated by multiMR tool.

GO Term and KEGG Pathway Enrichment Analyses on DEmRNAs

KEGG pathway and gene ontology analysis results indicated that DEmRNAs from GSE64813 were significantly enriched in the immune system process, innate immune response, cell adhesion, and defense response to virus (Figure 4A). It was found from the results of functional enrichment analysis that DEmRNAs from GSE89866 were particularly enriched in herpes simplex virus 1 infection, shigellosis, ribosome, endocytosis, and Salmonella infection (Figure 4B). Although the existence of difference between the function enrichments of two datasets DEmRNAs, interestingly, we found the GO term and KEGG pathway of the overlapping DEmRNAs were significantly enriched in Platelet activation, cell adhesion, and ECM-receptor interaction (Figure 4C). Also, we added functional enrichment analysis of 19 upregulated genes and 46 downregulated genes respectively. The analysis results (Tables S1–S4) showed that KEGG of downregulated genes also enriched in Platelet activation, cell adhesion, and ECM-receptor interaction pathways.

Figure 4 GO and KEGG enrichment analysis for DEmRNAs. (A) GO and KEGG enrichment analysis for DEmRNAs in GSE64813 dataset. (B) GO and KEGG enrichment analysis for DEmRNAs in GSE89866 dataset. (C) GO and KEGG enrichment analysis for DEmRNAs in GSE64813 and GSE89866 datasets. X-axis represents the ratio of mapped genes between DEmRNAs and each GO or KEGG term.

miRNA Interactome of PTSD

The integrated approach outlined in Figure 1 was applied using DEmRNAs and DEmiRNAs to generate a PTSD-specific miRNA interactome, with potential miRNA–mRNA target pairs. Eighteen prioritized miRNA–mRNA pairs, including four DEmiRNAs and ten DEmRNAs validated from the multiMiR database, were determined. The finding that 10 DE target genes interacted with four DEmiRNAs indicated that in a complex structure, miRNAs bound to target many genes. The network was composed of a large cluster of inter-connected miRNA–mRNA pairs in which the targets were downregulated, and miRNAs were upregulated, which contains 81 nodes and 175 edges (Figure 5A). Within the upregulated miRNA cluster, previously unknown PTSD-related CDKN1A showed two connected miRNAs. A total of 286 potential drugs interacting with 20 DEmRNAs were identified when drug-gene interactions were explored using DGIdb (Table S5). Most potential drugs might interact with prostaglandin-endoperoxide synthase 1 (PTGS1) (129/286), beta-tubulin gene (TUBB1) (50/286), and cyclin-dependent kinase inhibitor 1A (CDKN1A) (37/286) either in an unknown manner or as an inhibitor (Figure 5B, 310/337). Additionally, the drugs that interact with the CLIP2 gene were all antagonists and three up-regulated genes shared the same drugs that were all inhibitors (Figure 5B), indicating that these genes may be important potential therapeutic targets for PTSD.

Figure 5 PTSD miRNA–mRNA interactome and interacting drugs. (A) Network of differentially expressed miRNAs targeting differentially expressed genes between PTSD and control samples. Triangles are DEmiRNAs and circles are DEmRNAs; edges denote that a miRNA targets the connected gene. The size of the node labels is proportional to fold change. Node color characterizes the direction of the expression change between controls and PTSDs. (B) Drugs interacted with overlapped DEmRNAs. Gene and potential interacted drugs are connected with lines. The interaction types are represented by different arrows.

DEmRNAs for Diagnostic Prediction of PTSD

After merging the two datasets (GSE64813 and GSE89866 datasets), the control and PTSD samples are randomly (seed=12345) divided into the same proportion of training set (control/PTSD: 38/38) and validation set (control/PTSD: 38/38). Batch correction analysis showed that batch effects were eliminated after batch correction when merging the GSE64813 and GSE89866 datasets (Figure S1A and B). The random forest (Figure 6A and B) and LASSO (Figure 6C) methods yielded 48 and 17 DEmRNAs, of which all the DEmRNAs were common to both methods. Using the logistic regression method, we further used these 17 DEmRNAs to construct a risk score model for PTSD diagnosis. The Violin diagram (Figure 6D) shows that the risk scores of the PTSD samples of the two datasets were significantly higher than those of the control samples. In addition, we also evaluated the ability of the diagnostic model’s risk score to classify PTSD and control samples, and found that on the training and validation sets, the diagnostic model’s risk score classified the two samples well, and the AUC values reached 0.867 and 0.845 (Figure 6E). Two independent datasets, GSE109409 and GSE63878, were then used to assess the model performance and showed an AUC of 0.767 and 0.795, respectively, indicating a good discrimination between PTSD and controls (Figure 6F).

Figure 6 PTSD risk score model construction and validation. (A) Error rate generated from random forest analysis in the training cohort. (B) The importance of DEmRNAs evaluated by random forest algorithm. Top 10 important gene symbols were presented. (C) Misclassification error for different numbers of variables revealed by the LASSO regression model. The red dots represent the value of the misclassification error, the grey lines represent the standard error (SE), the two vertical dotted lines on the left and right, respectively, represent optimal values by the minimum criteria and 1-SE criteria. “Lambda” is the tuning parameter. (D) Distribution of risk score in different datasets. The box plots inside the violin indicate the median value and interquartile range of PTSD risk score. (E) ROC of the PTSD risk score model in the training and validation cohorts. (F) ROC analysis for the 17 DEmRNA based diagnostic prediction model in two independent datasets GSE109409 and GSE63878. Abbreviations: AUC, area under curve; CI, confidence interval; Speci, specificity; Sensi, sensitivity.

Discussion

Differential expression analysis of two independent datasets was performed to identify pathways and genes involved in PTSD etiology. This study found a small number of differentially expressed genes in peripheral blood collected from individuals with and without PTSD in the two datasets. To understand the role of DEmRNAs in PTSD, we detected an over-representation of genes associated with ECM-receptor interaction and platelet activation pathways. These results are consistent with other studies,^34–36 and show a robust relationship between the ECM-receptor interaction pathway and PTSD.

PTSD is a disorder that originates in the brain, despite the existence of a blood-brain barrier between the periphery and the brain, and may affect the gene expression patterns between them. However, many studies have confirmed that there is a good correspondence between gene expression levels in peripheral blood and brain tissue, as demonstrated by similar gene expression profiles observed in the postmortem brain.³⁷ The gene expression profiles in peripheral blood can be used as effective biomarkers for many psychiatric diseases.³⁸ Furthermore, given the persistence of psychosomatic symptoms of PTSD in many patients, gene expression in the peripheral blood is also associated with brain tissue expression.³⁹ Finally, we observed similar expression characteristics of these differentially expressed genes identified in peripheral blood samples from normal brain tissues (Figure 2F). These results suggest that peripheral blood gene expression profiles could also reflect changes in gene expression in the brains of patients with PTSD.

From the current miRNA-mRNA interaction network analyses, one of the novel findings was the identification of 4 DEmiRNAs associated with 12 DEmRNAs. Interestingly, gene CDKN1A network analysis also identified modules correlated with wound response and hemostasis with significant overlap.⁴⁰ We also identified several novel differentially expressed miRNAs in PTSD subjects and normal subjects, including hsa-miR-374b-5p, hsa-miR-629-5p, hsa-miR-101-3p, hsa-miR-20b-5p, and hsa-miR-485-3p. Three of the five DEmiRNAs, hsa-miR-20b-5p, hsa-miR-629-5p, and hsa-miR-101-3p, were also reported in previous studies.^20,41 The top hit was hsa-miR-20b-5p, which has 14 validated targets (Figure 3B). This indicates that hsa-miR-20b-5p is altered and is likely to be a candidate marker for PTSD, a stress-related disorder. It has been shown that hsa-miR-20b-5p is increased in the hippocampus of Appswe/PS^ΔE 9 mice, and inhibition of miR-20b-5p attenuates apoptosis induced by Aβ25-35 in Alzheimer’s disease by targeting RhoC.⁴² The expression of hsa-miR-485-3p was found to be upregulated in mice models that showed susceptibility to chronic unpredictable mild stress (CUMS) compared to CUMS- resilience.⁴³ In Cohen’s study, miR-485, in conjunction with synaptic vesicle protein SV2A, has been found to regulate dendritic spine number in an activity-dependent manner.⁴⁴ The expression of hsa-miR-374b-5p, which was newly identified in this study, was found to be a poor prognostic factor in human glioma.⁴⁵ These results suggested the reliability of our study.

In the current study, a diagnostic prediction model for PTSD was established using the screened DEmRNAs. The model risk score increased significantly from a control to PTSD, and a high AUC indicated good performance of the model in identifying patients with PTSD from individuals and healthy controls. This result may facilitate the development of a new diagnostic strategy from a molecular perspective. Nonetheless, this study also has the following limitations. First, more studies are warranted to confirm the robust performance of the model on a larger cohort, and the consistency between mRNAs and miRNAs in the circulation also needs to be established to determine whether the detection of DEmiRNAs from peripheral blood could serve as a novel tool for PTSD screening. Second, more efforts based on experimental data are required.

Conclusions

In summary, our data provide a new molecular aspect that ECM-receptor interaction and platelet activation biological process may be the potential molecular mechanism of PTSD, and the genes involved in this process may be promising therapeutic targets. The present study also provided a high-performance diagnostic model for PTSD. Specifically, DEmiRNAs and DEmRNAs associated with cell proliferation appear to be involved in PTSD. The data showed that modules related to wound healing and platelet activation may contribute to resilience against PTSD development. Future studies could confirm the role of platelets in the stress response. Due to limitations including small sample size and lack of epigenetic testing, in future studies, participants with and without PTSD will be studied more closely in relation to important potential confounding variables and multivariate data sets by performing a meta-analysis.