Machine Learning Analysis of Blood microRNA Data in Major Depression: A Case-Control Study for Biomarker Discovery.

BACKGROUND: There is a lack of reliable biomarkers for major depressive disorder (MDD) in clinical practice. However, several studies have shown an association between alterations in microRNA levels and MDD, albeit none of them has taken advantage of machine learning (ML).
METHOD: Supervised and unsupervised ML were applied to blood microRNA expression profiles from a MDD case-control dataset (n = 168) to distinguish between (1) case vs control status, (2) MDD severity levels defined based on the Montgomery-Asberg Depression Rating Scale, and (3) antidepressant responders vs nonresponders.
RESULTS: MDD cases were distinguishable from healthy controls with an area-under-the receiver-operating characteristic curve (AUC) of 0.97 on testing data. High- vs low-severity cases were distinguishable with an AUC of 0.63. Unsupervised clustering of patients, before supervised ML analysis of each cluster for MDD severity, improved the performance of the classifiers (AUC of 0.70 for cluster 1 and 0.76 for cluster 2). Antidepressant responders could not be successfully separated from nonresponders, even after patient stratification by unsupervised clustering. However, permutation testing of the top microRNA, identified by the ML model trained to distinguish responders vs nonresponders in each of the 2 clusters, showed an association with antidepressant response. Each of these microRNA markers was only significant when comparing responders vs nonresponders of the corresponding cluster, but not using the heterogeneous unclustered patient set.
CONCLUSIONS: Supervised and unsupervised ML analysis of microRNA may lead to robust biomarkers for monitoring clinical evolution and for more timely assessment of treatment in MDD patients.

The identification of biomarkers for complex disorders such as major depression is an important challenge to overcome to enable precision medicine. In this study, we demonstrate how machine learning could be an effective approach to address this challenge, in terms of diagnosis, disease severity or response to medication.

Introduction

In the United States, the lifetime prevalence for major depressive disorder (MDD) is 20.6% among individuals aged 18 years or older. Almost one-half (49%) of the cases have severe and 39.7% moderate depression (Hasin et al., 2018). Without early treatment, there can be permanent consequences on the patient’s brain function that increase their risk of experiencing additional depressive episodes (Moylan et al., 2013). Overall, the economic burden of MDD is more than US$ 170 billion per year and appears to be increasing over time (Greenberg et al. (2015). However, there is still a lack of reliable biomarkers that can guide patient monitoring and timely assessment of treatment efficacy.

Increasing evidence suggests that molecular signaling for depression is linked with microRNA expression and that the dysregulation of microRNA signaling can initiate or exacerbate depressive pathophysiology (Hansen and Obrietan, 2013). MicroRNAs are small noncoding RNA molecules that play a role in the regulation of gene expression and neuronal physiology. Smalheiser et al. (2012) found that the expression of several microRNAs was significantly downregulated in the prefrontal cortex of depressed suicide individuals compared with matched psychiatric control participants. Bocchio-Chiavetto et al. (2013) measured the expression of microRNA in 10 depressed individuals before and after treatment with antidepressants. After the treatment with antidepressants, 2 microRNAs were significantly downregulated and 28 were upregulated. In a recent randomized placebo-controlled trial, we identified several microRNA markers of duloxetine treatment response that were replicated in 2 independent clinical trials, an animal model, and postmortem brain samples (Lopez et al., 2017). The findings suggest that there is a strong possibility that microRNAs are involved in the pathophysiology of depression and affect the mechanism of action of antidepressants.

Machine learning (ML) algorithms have been created for analyzing complex multivariate data with a focus on empirical predictive power and generalizability. ML has demonstrated success in clinical psychiatry in terms of diagnosis, prognosis, treatment decisions, and biomarker detection (Dwyer et al., 2018). A review of the literature on ML and MDD shows a shortage of studies that apply ML methods to analyze microRNA data (Gao et al., 2018). However, a recent paper has demonstrated the effective use of ML in identifying a serum microRNA signature for Alzheimer’s disease that could predict disease status with 85.7% accuracy (Zhao et al., 2019).

Given the important role of microRNAs in MDD and the effectiveness of ML in taking advantage of complex data, we aimed to explore whether ML analysis of blood microRNA profiles can serve as a new approach for biomarker discovery in MDD.

METHODS

Participant Recruitment

The study protocol was approved by the Research Ethics Board of the Douglas Mental Health University Institute. Informed written consent was obtained from all participants. All participants were recruited from an outpatient clinic at the Douglas Mental Health University Institute in Montréal, Canada, and assessed by an experienced psychiatrist using the SCID-I (First et al., 2012) following DSM-IV criteria. Patients were all suffering from a current major depressive episode (MDE) as part of an MDD. Exclusion criteria included comorbidity with other major psychiatric disorders, bipolar disorder, alcohol or substance abuse over the last 6 months, or a severe medical condition. None of the participants were medicated at baseline. None had received fluoxetine or lithium over the last month or any psychotropic medication over the last week. Depression severity was determined using the Montgomery-Asberg Depression Rating Scale (MADRS). MADRS measures are based on 10 different symptoms, including (1) apparent sadness, (2) reported sadness, (3) inner tension, (4) reduced sleep, (5) reduced appetite, (6) concentration difficulties, (7) lassitude, (8) inability to feel, (9) pessimistic thoughts, and (10) suicidal thoughts. The MADRS scores were collected at baseline and again after 8 weeks of antidepressant treatment.

Sample Processing

Peripheral blood samples were collected at baseline and after 8 weeks, and tubes were frozen using a sequential freezing process. Whole blood for RNA was collected in ethylenediaminetetraacetic acid tubes and filtered using LeukoLOCK filters (Life Technologies). Total RNA was extracted using a modified version of the LeuEkoLOCK Total RNA Isolation System protocol and included DNase treatment to remove genomic DNA. The RNA quality was assessed using the Agilent 2200 Tapestation, and only samples with an RNA integrity number ≥6.0 were used.

Small RNA-seq

All libraries were prepared using the Illumina TruSeq small RNA Library preparation protocol following the manufacturer’s instructions. Samples were sequenced at the McGill University and Genome Quebec Innovation Centre (Montreal, Canada) using the Illumina HiSeq2000 with 50-nucleotide (nt) single-end reads. All sequencing data were processed using CASAVA 1.8+ (Illumina) and extracted from FASTQ files. The Fastx_toolkit was used to trim the Illumina adapter sequences. Additional filtering based on defined cutoffs was applied, including (1) Phred quality (Q) mean scores higher than 30, (2) reads between 15 and 40 nt in length, (3) adapter detection based on perfect-10 nt match, and (4) removal of reads without a detected adapter. Additionally, we used Bowtie (Song et al., 2014) to align reads to the human genome (GRCh37) and ncPRO-seq (Chen et al., 2012) in combination with miRBase (V20) (Kozomara and Griffiths-Jones, 2013) to match them to known microRNA sequences. Furthermore, all sequencing data were normalized with the Bioconductor-DESeq2 package (Love et al., 2014), using a detection threshold of 10 counts per miRNA.

The number of microRNA samples included for subsequent analyses includes the baseline (T0) and week 8 (T8) of 140 MDD cases and 28 healthy controls. The total number of microRNA features is 285.

ML Analysis

Many powerful ML algorithms render themselves uninterpretable, making it difficult to understand their decision-making process. For our ML analysis of the data, we decided to use a state-of-the-art yet interpretable regularized gradient boosted machines approach (XGBoost implementation, Chen and Guestrin, 2016), which we also demonstrated as an effective algorithm in our previous study of schizophrenia (Trakadis et al., 2019).

Datasets are split into 70% and 30% for training and testing. A model selection procedure based on 5-fold cross-validation with 2500 iterations of parameter search was used to obtain the best training parameters using only the training dataset (n = 122). After obtaining the best training parameters, we retrained the model without 5-fold cross-validation (i.e., using at once the entire training set) and evaluated the model on the testing set. The model performance metric we used is the area under the receiver-operating characteristics curve (AUC).

Classification Analyses

With regards to discriminating cases from healthy controls, we trained the ML model using only the T0 microRNA data to ensure that the medication effect would not act as a confounder in the analysis.

For the severity class classification, we used the MADRS cutoff scores suggested by Snaith et al. (1986). Individuals’ MADRS scores were classified as “normal-mild” (MADRS scores from 0 to 19) or “moderate-severe” (MADRS scores 20 and above). Using the class-labeled dataset, we identified the best classification model for classifying samples into these 2 MDD grades. However, for this analysis, we used the T8 microRNA data and MADRS scores, because at T0 all but 2 cases had MADRS scores ranging from 0 to 19 (thus, using T0, almost all samples would be labeled as “moderate-severe”).

We then repeated severity classification after unsupervised clustering of the T8 microRNA data, which was done to factor in the heterogeneity of MDD. More specifically, 500 iterations of a consensus k-means clustering method (Monti et al., 2003) were applied to the entire case-control dataset (n = 174). The model selection and evaluation procedure were then performed separately for each cluster under the assumption that the patients in each cluster are less heterogeneous at the pathophysiology and microRNA levels. If our assumption is correct, training the ML algorithm to identify signatures specific for the “normal-mild” vs “moderate-severe” class would be more efficient after unsupervised clustering, and thus the classification based on supervised ML analysis of the microRNA data would improve with this approach.

Lastly, we explored the relationship between microRNA and antidepressant response among MDD patients using the difference between the T8 and T0 microRNA values (n = 138, because 2 MDD cases were missing MADRS scores at T0). With regards to severity levels, the scores used were based on previous definitions: “normal” (0–6), “mild” (7–19), “moderate” (20–34), and “severe” (>34) (Snaith et al., 1986). Antidepressant response in our study was defined as a decrease of 2 severity levels when comparing the patient’s T8 and T0 MADRS scores. For example, a patient with a change from “severe” to “mild” or a change from “moderate” to “normal” would be labeled as a responder. We obtained a split of 46 RES and 92 nonresponders. We then repeated the same ML procedure described above to obtain a classifier for responders vs nonresponders.

Since each patient was taking a mixture of multiple antidepressants, to address this heterogeneity and improve the performance of the classifier we performed unsupervised clustering on the “T8-T0” dataset using the consensus k-means algorithm described above in the MDD severity classification section. Samples were split into 2 clusters. We then performed ML classification analysis for antidepressant response separately in each individual cluster. To explore if the top microRNA for each cluster (i.e., the microRNA with the maximum importance in the ML classification model) was associated with antidepressant response, a permutation test was performed. Specifically, 500 000 iterations were performed to derive the empirical P value of a difference in mean between responders and nonresponders. The significance threshold was set at .05. In the case of multiple top microRNAs (multiple top microRNAs having equal maximum importance), the P values were adjusted using the Bonferroni correction method. To explore if the top microRNA(s) identified in each cluster were specific to that cluster, we performed permutation testing for microRNAs extracted from the first cluster using samples from the second cluster, and vice versa. We also performed permutation tests for the top microRNAs identified from the clusters using all MDD samples (unclustered). Finally, we extracted the top microRNA from the antidepressant classification model trained on all samples and performed permutation testing to explore if ML analysis of the data before stratification was helpful in the identification of a marker for treatment response.

Clinical History Analysis

To explore how clinical history factors into antidepressant response, we examined whether patients with a prior history of treatment with antidepressants responded differently compared with antidepressant-naïve patients. We also examined whether patients who present with their first MDEs (i.e., no prior episodes besides the current one) responded differently compared with patients with recurring MDEs (collected from SCID-IA [DSM-IV], question A29). Antidepressant response here is defined as a ratio of the T8 to T0 MADRS score (T8/T0) in order to capture more precise differences in antidepressant response between groups using the permutation tests. Permutation tests were performed for 500 000 iterations to derive the P value for a significant difference in antidepressant response between the groups compared. Multiple testing was adjusted using the Bonferroni correction method.

Bioinformatic Analysis

We extracted the microRNA features used by the best case-control classification model and performed pathway analysis using the DIANA-miRPath v3.0 pathway analysis webserver (Vlachos et al., 2015) to obtain KEGG pathway terms significantly related to the set of microRNA features. Pathways with a false-discovery rate (FDR) < 0.05 were selected.

Software

The ML model was implemented using the Python (v.3.7.1) programming language (https://www.python.org/) with the “xgboost” (v.0.81) library (https://xgboost.readthedocs.io/). The consensus clustering procedure was implemented using the “scikit-learn” (v.0.21.2) (https://scikit-learn.org/) and “scipy” (v.1.3.0) (https://www.scipy.org/) Python libraries.

RESULTS

The demographics of patients and controls are summarized in supplemental Table 1. A total 65% of MDD cases and 46% of controls were female. Moreover, 80% of MDD cases and 82% of controls were Caucasians. The mean MADRS score at T0 was 33 (SD = 6.2) for cases and 0.6 (SD = 1.1) for controls. The mean MADRS score at T8 was 17.4 (SD = 10.9) for cases and 1.1 (SD = 1.7) for controls. A total 56% (n = 79) of patients reported presenting with their first MDE, and 14% (n = 19) of MDD patients were antidepressant naïve prior to current treatment.

As summarized in Table 1, for classification of cases and controls, the best trained model achieved an average cross-validation AUC of 0.93 (SD 0.06) and testing set AUC of 0.97. The best trained model trained to distinguish cases from controls utilized 33 out of 285 total microRNAs measured (Table 2). Pathway analysis for the 33 microRNAs found the following significantly enriched pathways with FDR < 0.05: (1) prion diseases, (2) transforming growth factor beta (TGF-beta) signaling pathway, (3) morphine addiction, (4) signaling pathways regulating pluripotency of stem cells, (5) mucin type O-glycan biosynthesis, and (6) proteoglycans in cancer.

Table 1.

Model Cross-Validation and Testing Set AUC Scores

Analysis	Mean AUC (SD) of trained model from cross-validation	Testing set AUC for final retrained model
Classification of cases and controls	0.93 (0.06)	0.97
Classification of MDD severity grades	0.76 (0.11)	0.63
Classification of MDD severity grades: cluster 1	0.75 (0.18)	0.76
Classification of MDD severity grades: cluster 2	0.72 (0.15)	0.70
Classification of antidepressant response	0.622 (0.13)	0.57
Classification of antidepressant response: cluster 1	0.652 (0.085)	0.54
Classification of antidepressant response: cluster 2	0.670 (0.16)	0.49

Abbreviations: AUC, area under the receiver-operating characteristics curve; MDD, major depressive disorder.

Model selection and evaluation were performed for each of the analyses listed in the table. The mean AUC across 5-folds of cross-validation during model training for the best model is presented, as well as the AUC from the evaluation on the testing set for the final retrained model.

Table 2.

Most Important microRNA Features Used by the Case-Control Classification Model

MicroRNA features ordered by decreasing importance (n = 33)

hsa-miR-27a-3p, hsa-miR-197-3p, hsa-miR-22-5p, hsa-miR-221-3p, hsa-miR-126-3p, hsa-miR-128-1-5p, hsa-miR-30b-5p, hsa-miR-339-3p, hsa-miR-301a-3p, hsa-miR-345-5p, hsa-miR-505-3p, hsa-miR-1249, hsa-miR-132-3p, hsa-miR-550a-5p, hsa-miR-589-5p, hsa-miR-769-5p, hsa-miR-10b-5p, hsa-miR-210-3p, hsa-miR-628-3p, hsa-let-7d-3p, hsa-miR-148a-5p, hsa-miR-155-5p, hsa-miR-140-3p, hsa-miR-150-3p, hsa-miR-181a-5p, hsa-miR-24-3p, hsa-miR-629-5p, hsa-let-7a-3p, hsa-miR-194-5p, hsa-miR-28-3p, hsa-miR-378a-3p, hsa-miR-6852-5p, hsa-miR-7706

The most important features used by the best performing machine learning model from the classification of cases and controls analysis, which could distinguish between cases and controls with an AUC of 0.97, were extracted and listed in order of decreasing importance.

Classification of individuals as normal-mild vs moderate-severe MADRS grades using their microRNA data based on best trained model showed an average cross-validation AUC of 0.76 (SD = 0.11). After retraining the best model on the full dataset and evaluating on the testing set, we obtained an AUC of 0.63.

For the clustering approach, we obtained 2 clusters (cluster 1: 89 participants; cluster 2: 79 participants) of similar sample size, which did not show differences in terms of MDD severity. The best model for cluster 1 samples achieved an average cross-validation AUC of 0.75 (SD = 0.18), while the best model for cluster 2 samples achieved an average cross-validation AUC of 0.72 (SD = 0.15). When evaluated on the testing sets, the cluster 1 model achieved an AUC of 0.76, while the cluster 2 model achieved an AUC of 0.70. Table 1 summarizes the results for each of the analyses.

For antidepressant response classification, we obtained an average cross-validation AUC of 0.62 (SD = 0.13) and an AUC of 0.57 on the testing set. After clustering, we again obtained 2 balanced clusters (cluster 1: 69 participants; cluster 2: 69 participants). We did not notice a separation of responders from nonresponders based on clustering. The best model for cluster 1 samples achieved an average cross-validation AUC of 0.65 (SD = 0.085), while for cluster 2 the average cross-validation AUC was 0.67 (SD = 0.16). On testing set evaluation, the cluster 1 model achieved an AUC of 0.54, while the cluster 2 model achieved an AUC of 0.49. For cluster 1, after supervised ML for classification of treatment response, the top and only microRNA utilized by the ML model was hsa-miR-5701. Following permutation testing, this microRNA was found to be significantly different between responders and nonresponders in cluster 1 (P = .021), but not in cluster 2 or the original (unclustered) dataset. For cluster 2, there were 4 microRNAs, all with equal importance, including hsa-let-7b-3p, hsa-let-7g-5p, hsa-miR-130b-3p, and hsa-miR-30d-3p. Following permutation testing, the only nominally significant microRNAs were hsa-let-7b-3p (P = .021) and hsa-miR-130b-3p (P = .045), albeit neither was significant after Bonferroni correction (P = .082 and P = .18, respectively). None of these 4 markers significantly differed between responders and nonresponders in cluster 1 or the original (unclustered) dataset. Finally, when extracting the top microRNA from the antidepressant classification model trained on all samples, the top microRNA was not found to be associated with treatment response following permutation testing (P = .12).

Of note, we observed that antidepressant-naïve patients responded significantly better than those who had taken antidepressants in the past (P = .00058, with Bonferroni correction) but did not observe a significant difference in response between patients who presented with their first vs recurring MDEs (P = .59, with Bonferroni correction).

Discussion

In this paper, we demonstrate how ML analysis of blood microRNA data could lead to biomarkers with potential clinical utility. Our assumption was that if this is true, ML analysis of microRNA data should not only lead to the successful classification of cases from controls but also to the efficient separation of individuals with mild vs severe depression.

First, we showed that microRNA data could be used to discriminate baseline medication-free MDD cases from controls (AUC = 0.97 using the test dataset). Of note, this result is not expected to be confounded by medication effects since we used only the T0 pretreatment trial microRNA data. To show that the microRNA signals are relevant to MDD, we conducted a pathway analysis using the microRNAs identified by the ML model (FDR < 0.05). We identified 6 pathways and highlighted the evidence in the literature for a link with MDD.

For example, there is evidence that endogenous prion protein (PrP(C)) is associated with MDD. PrP(C) was reduced in the white matter (Weis et al., 2008) and Brodmann’s areas 6 and 10 (Dean et al., 2019) in patients with MDD. PrP(C) has also been shown to modulate depressive-like behavior in mice (Gadotti et al., 2012).

The TGF-beta family of cytokines may also play a role in MDD. TGF-beta has been observed to be significantly elevated in the peripheral blood of MDD patients (Davami et al., 2016). Furthermore, a study found a significant decrease in TGF-beta1 in MDD patients after 6 weeks of treatment with an antidepressant (Kim et al., 2007).

The dopamine and reward systems are major parts of the morphine addiction pathway (Kim et al., 2016), and there is a link between dopamine neurons and depression (Knowland and Lim, 2018). There is also evidence linking stem cell and cell renewal capacity to MDD. In mice with interferon-α–induced depression, neural stem cell proliferation was found to be suppressed (Zheng et al., 2014). Furthermore, shorter telomere length (Verhoeven et al., 2014) is also associated with a higher severity of depression.

Although no direct link exists between mucin type O-glycans and MDD, a study showed that the p75 neurotrophin receptor, a heavily glycosylated protein, had a polymorphism, Ser205Leu, for a predicted O-glycosylation site, which had a protective effect for MDD (Fujii et al., 2011).

Next, we showed that microRNAs could be leveraged to distinguish participants with normal-mild from moderate-severe MDD (AUC = 0.63). We also demonstrated that the use of unsupervised clustering, aimed at reducing MDD heterogeneity, can improve model performance in our MDD grade classification task (AUC of 0.76 for cluster 2 and 0.70 for cluster 1). This supported our assumption that the individuals in each cluster were less heterogeneous after unsupervised clustering. This lead to a more efficient training of the ML algorithm to identify signatures specific for the “normal-mild” vs “moderate-severe” class. The sample size of our dataset was relatively small. However, given our results, we expect that performance estimates would improve and become more precise with ML models trained on larger samples.

We found that the differences between the T8 and T0 microRNAs were not strongly predictive of response status (AUC = 0.57 on the testing set). This came as no surprise given that the patients were undergoing treatment with different antidepressants, thus leading to heterogeneity negatively impacting the performance of the ML model. Patient stratification partially addressed this, as we saw a slight boost to the 5-fold cross-validation performance for each cluster compared with the unclustered analysis. However, we did not see any improvement in classifying response status on the testing set. We believe that the poor performance of the ML models on the testing set is likely due to the small sample size of each cluster but that there may still be intelligence derived from the approach. This is supported by the significant association of the top microRNA within each cluster with antidepressant response, which was specific (i.e., was not observed when analyzing the data of the other cluster or the unclustered data). Furthermore, no marker was found for antidepressant response when extracting the top microRNA from the antidepressant classification model trained on all samples (i.e., unclustered dataset). Putting everything together, we take this as evidence that a clustering approach, combined with supervised ML, could be useful to identify biomarkers in subgroups of patients that would otherwise be missed when analyzing heterogeneous populations.

Of note, our approach, using regularized ML with empirical cross-validation and testing as a method to prioritize features rather than multiple univariate testing, facilitates finding relevant biomarkers with minor effects or complex interactions that would otherwise be filtered out by multiple testing correction. This is very important given the complex relationships of different factors contributing to MDD such as duration of the depressive episode, duration of illness, and recurrence. For example, we found that patients who have no history of taking antidepressants responded significantly more to treatment compared with those with past history. At first sight, this should come as no surprise, since the usage of more antidepressant in the past indicates that the patient did not respond to the previous antidepressants and thus that they are harder to treat. However, we did not observe a difference between patients experiencing their first MDE vs patients with recurring MDEs, which contradicts this line of thinking and underlines the complexity of the different factors and their interaction in MDD.

Changes at the microRNA level are downstream to the different contributing clinical factors, thus explaining why we were able to distinguish cases from controls (AUC = 0.97) successfully. However, to better understand the contribution of each factor in MDD, further studies with larger sample size and more optimal patient stratification, with the inclusion of genetic, functional genetic, and detailed clinical data, would be recommended. Moreover, future studies should not be focused on examining changes between binary time points for antidepressant treatment but rather serial (i.e., at multiple time points) MADRS evaluation and collection of microRNA data. With this design, we could explore if early changes at the microRNA level after treatment initiation could predict treatment response at a later point, which would have major clinical implications in treatment optimization.

CONCLUSION

Our manuscript provides preliminary evidence that ML analysis of blood microRNA profiles may constitute a reliable approach for biomarker discovery for MDD (affected vs unaffected) clinical status, but also for clinical evolution (severity and treatment response), thus facilitating a more personalized approach in treating patients with MDD.

Click here for additional data file.