The GEOquery package was utilized to download the dataset from the GSE76826 GEO database (https://www.ncbi.nlm.nih.gov/) (Miyata et al., 2016). This dataset examined biomarkers in older outpatients and inpatients with Major Depressive Disorder (MDD), specifically those aged 50 years and above. The depressive status of participants was assessed using the Hamilton Depression Structured Interview Guide (SIGH-D) rating scale. Subsequently, the normalizeBetweenArrays function from the limma package was employed for further normalization of the data. Following this, differential expression analysis of genes between 12 normal samples and 10 patients diagnosed with depression was conducted using the limma package. Differentially expressed genes were identified based on criteria of |log2FC| > 1.2 and P < 0.05. The results of this differential analysis were visualized through a volcano plot. During annotation processing, probes corresponding to multiple molecules were removed, ensuring that only unique probes remained for each molecule encountered; among these, only the probe exhibiting the highest signal value was retained. Finally, KEGG enrichment analysis was performed on the selected differentially expressed genes. Additionally, GWAS data on depression (dataset finn-b-F5_DEPRESSIO) comprising23,424 European individuals and covering 8,281,749 Sips was obtainedfrom the IEU Open GWAS database (https://gwas.mrcieu.ac.uk/datasets/finn-b-F5_DEPRESSIO/).

The genes associated with OS were screened in the GSEA database. Subsequently, an intersection analysis was conducted between depression-related target genes and oxidative stress-related target genes. To optimize the classification of these intersecting genes, two machine learning algorithms, namely, support vector machine (SVM) and random forest, were employed for comprehensive analysis and processing. The R software [4.2.1] was utilized for data analysis, employing the random Forest package (version 4.7.1.1) to construct a random forest model aimed at evaluating feature importance. Initially, a random seed (seed = 2024) was established to ensure the reproducibility of results. Following the framework of the random forest algorithm proposed by Breiman (Breiman, 2001), an ensemble model comprising 100 decision trees was developed using default parameter settings for data analysis. The Mean Decrease Accuracy of each feature was computed to quantify its importance score during model training. Simultaneously, SVM (Meyer et al., 2015) analysis was conducted utilizing the e1071 package [1.7.13] for classification modeling purposes. A random seed (seed = 2024) was preset prior to data processing, ensuring repeatability in the calculation process through the set.seed() function. When constructing the classification model with the svm() function, a 5-fold cross-validation strategy was employed to assess model stability effectively. The Gaussian radial basis kernel was selected as the algorithm’s kernel, facilitating nonlinear mapping to handle complex data structures efficiently. Results were visualized using the ggplot2 package [3.4.4].

To assess the diagnostic value of a gene, we employ receiver operating characteristic (ROC) analysis to construct a diagnostic model and subsequently validate its performance using a Diagnostic Nomogram. The gene-based nomogram was constructed using the “rms' package. Subsequently, a binary logistic model was built using the glm function, and model calibration analysis and visualization were performed using the rms [6.4.0] package. The rmda package was utilized to calculate the corresponding net return rate and facilitate visual display, with results visualized through ggplot2 [3.3.6].

In the two-way MR analysis, a screening threshold of P < 1 × 10⁻⁵ was applied to ensure an adequate number of SNPs. Subsequently, instrumental variables (IVs) exhibiting strong linkage disequilibrium (LD) were excluded based on r2 < 0.001 and kb = 10,000 criteria. The F statistic was then calculated as F = [R2/(1-R2)] × (N-K-1)/K, where R2 represents the proportion of exposure factor change explained by IVs, N denotes the sample size associated with the exposure factor, and K represents the number of instrumental variables used. If F < 10 is considered indicative of weak instrumentality, it should be eliminated.

In order to investigate the causal relationship between depression and key genes, this study employed a two-sample MR analysis. The instrumental variable single nucleotide polymorphism (SNP) used in this analysis should satisfy the assumptions of “correlation with exposure factors', which include being correlated with exposure factors, unaffected by confounding factors, and only affecting outcome factors through exposure factors. Key genes related to the phenotype were screened from the National Center for Biotechnology’s database (https://www.ncbi.nlm.nih.gov/). A GWAS aggregate statistical dataset available on the MR Platform was utilized, consisting of 23,424 individuals of European descent as phenotypes for depression. Various methods including MR Egger, weighted median, inverse variance weighting (IVW), simple model, and weighted model were employed to conduct two-sample MR analyses. The IVW results served as the primary evaluation indicators while other method results estimating causal effects were considered sensitivity analyses. Statistical significance was determined if P < 0.05. To facilitate result interpretation, Beta values obtained in this study were converted into odds ratios (OR), accompanied by calculation of their corresponding 95% confidence intervals (CI).

If the instrumental variable affects the outcome through factors other than the exposure factor, it indicates the presence of horizontal pleiotropy, which undermines the validity of the assumptions of independence and exclusivity. To assess horizontal pleiotropy, we employ MR-Egger intercept analysis. A P-value greater than 0.05 suggests no statistically significant deviation from zero intercept, indicating that the results are not influenced by horizontal pleiotropy and confirming that the instrumental variable solely impacts the outcome through exposure.

Numerous studies have demonstrated that depression is a psychiatric disorder characterized by aberrant functionality in specific cerebral regions. In this investigation, the therapeutic potential of medicinal gene intervention for depression was explored using the Human eFP (“electronic fluorescent hieroglyphs') browser (http://bar.utoronto.ca/efp_human/), which facilitated rapid examination of gene expression profiles and intuitively revealed the diverse brain regions expressing IL1R1.

The GSEA database was utilized to screen for genes associated with lipid metabolism, enabling the analysis of differential gene expression related to lipid metabolism in depression. Additionally, the STRING data facilitated network interaction analysis of these differentially expressed genes. GSE76826 data set was used for Spearman correlation analysis of differential genes.

We employed DrugnomeAI (https://astrazeneca-cgr-publications.github.io/DrugnomeAI/) to assess the medicinal potential of IL1R1, considering a target with a percentile score >60 in DrugnomeAI gene rankings as having promising medicinal properties. Subsequently, we utilized the dgidb follow-up database (https://www.dgidb.org/), ChEMBL database (https://www.ebi.ac.uk/chembl/), and relevant literature searches on gene-drug interactions to evaluate the effects of potential drug candidates on IL1R1. Furthermore, stringent criteria for safety, efficacy, and ethical approval were applied to select suitable drug candidates.

We utilized AutoDockTools-1.5.7 software to conduct molecular docking of the medicinal gene with the potential drug candidate and determine the binding energy, which served as an indicator for evaluating the binding activity between pharmaceutically available genes and drug candidates. A higher absolute value of binding energy indicates better binding activity and stability. Specifically, a binding energy ≤ −5.0 kcal/mol suggests a more favorable interaction between pharmaceutically available genes and drug candidates. The results demonstrated that the combination of the medicinal gene and candidate drug exhibited significant efficacy.

Molecular dynamics simulations were performed using Gromacs 2023Gromacs 2023.2 with the receptor topology using the AMBER99SB-ILDN force field and the ligand topology constructed from GAFF2 force field parameters. Pre-processing such as energy minimization, isothermal isovolumic (NVT) and isothermal isobaric (NPT) equilibrium were performed with a run time of 100 ns and a time step of 2 fs.

The analysis identified a total of 1,472 differentially expressed genes, with 713 genes upregulated and 759 genes downregulated (Figure 1A). These differential gene expression patterns were visualized using volcano plot analysis. To gain insights into the potential biological functions associated with these differentially expressed genes, KEGG enrichment analysis was conducted. Notably, the enriched pathways included glutathione metabolism, arachidonic acid metabolism, MAPK signaling pathway, and NF-kappa B signaling pathway. Of particular interest is glutathione (GSH), which plays a crucial role as an antioxidant in cellular defense against reactive oxygen and nitrogen species while maintaining redox homeostasis. Collectively, these findings suggest that the observed differential gene expression may be implicated in OS and inflammation (Figure 1B).

A total of 655 genes related to OS were screened (Supplementary Table S1), out of which 21 genes were found to be associated with depression (Figure 1C; Supplementary TABLE S2). Correlation analysis was performed on these 21 genes, revealing a significant correlation among them as depicted in Figure 1D. To achieve more accurate classification between the disease group and normal group, two analytical methods, namely, support vector machine (SVM) and random forest were employed for processing. The impact of each gene on the heterogeneity of observations at each node in the classification tree is assessed through random forest analysis. The findings reveal that 13 genes significantly contribute to the classification of depression (Figure 1E; Supplementary Table S3). In support vector machine analysis, a set of 20 feature variables achieved the lowest error rate of 0.08, indicating their importance. Notably, IL1R1 exhibited the lowest score in cross-validation, suggesting its highest significance in classification (Figure 1F; Supplementary Table S4).

The two machine learning techniques identified a set of twelve genes (IL1R1, EPAS1, HSPB1, ETS1, IL6, SRXN1, PAX2, ANKRD2, GPX7, CD36, PRKD1 and CCR7) that exhibited significant potential for classification (Figure 2A). To further evaluate the diagnostic value of these genes we employed a ROC diagnostic model. Our findings revealed that both ILIRI (AUC = 0.950) and EPASl (AUC = 0.917) demonstrated robust diagnostic power for depression (Figures 2B,C; Supplementary Table S5). Furthermore, we validated this using a Diagnostic Nomogram, where we assessed the probability of a patient being diagnosed with depression by integrating the scores of the IL1R1 and EPAS1 genes (Figure 2D). The calibration curve also confirmed the excellent predictive ability of our ROC model (Figure 2E). Additionally, The DCA visually demonstrates that these tools and models can improve the efficiency of the diagnosis and treatment of depression (Figure 2F). Collectively, the results suggest that both EPASl and ILIRI hold promise as biomarkers for depression with ILIRI exhibiting superior diagnostic capability.

Prospective MR results demonstrated that the key genes IL1R1 (ENSG00000115594, OR = 0.928; 95% CI, 0.865 to 0.995; P = 0.035) and LPO (ENSG00000167419, OR = 0.937; 95% CI, 0.893 to 0.983; P = 0.008) were significantly associated with a decreased risk of depression, suggesting their potential as protective factors against depression development (Figure 3). However, reverse MR analysis did not reveal any evidence of a causal relationship between depression and these key genes. To assess the stability of our findings, we conducted an evaluation for horizontal pleiotropy using the MR Egger intercept test which yielded a non-significant result (P > 0.05), indicating no presence of horizontal pleiotropy bias in our study. Furthermore, through leave-one-out analysis, we found that none of the individual SNPs included in this study had a substantial impact on the robustness of our results, thus confirming their stability (Supplementary Figures S1, S2).

The analysis of the GSE76826 dataset revealed a significant upregulation of IL1R1 expression. (Figure 4A). To further validate this finding, we examined the expression pattern of IL1R1 in the prefrontal cortex using the GSE12654 dataset (Figure 4B). Consistently, our results demonstrated a pronounced elevation of IL1R1 expression in the prefrontal cortex. Similarly, the GSE98793(64 normal individuals and 128 individuals with depression) dataset found a high expression pattern of IL1R1 in peripheral blood (Figure 4C). Moreover, utilizing the Human eFP Browser tool, we investigated the expression profile of IL1R1 in various brain tissues and observed its prominent enrichment in the olfactory bulb, corpus callosum, and hippocampus (Figure 4D). Intriguingly, cardiac muscle cells and trigeminal ganglia exhibited particularly elevated levels of IL1R1 (Figure 4E).

Dyslipidemia is closely associated with depression, a common condition in patients with dyslipidemia that affects brain lipid metabolism and potentially impacts mood and mental health. In patients with depression, the downregulation of CCR7, GNB2L1, PLA2G6, PIK3IP1, PDGFB, CD81, HTR2C, ORMDL1, and PIK3R1 genes was observed. Conversely, upregulation of NOD2, SAMD, ORMDL2, FGR, KIT, XBP, NSMAF, TNFRSF1, PRKC, LYN, CDC42, VAV3, RAC, SLC27A1, IDH PRKAA genes related to lipid metabolism were detected (Figure 5A). IL-1 serves as a crucial pro-inflammatory cytokine that activates intracellular signaling pathways by binding to IL-IRI receptor and triggers inflammation. Prolonged or excessive inflammatory responses may negatively impact lipid metabolism. Therefore, we identified associations between IL-IRI and TNFRSFIA, NODZ, PDGFB, LVN, and PIK3Rl through the STRING database (Figure 5B). Negative correlations were found between PDGFB and PIK3Rl while positive correlations were observed for TNFRSFIA, NODZ, and LYN (Figure 5C).

The DrugnomeAI analysis revealed a percentile score of 95.87 for IL1R1, indicating its potential medicinal properties. To identify potential drug candidates targeting IL1R1, we conducted a comprehensive gene-drug search. Table 1 presents the small molecule drugs that were selected after rigorous screening and evaluation. These include Epoetin alfa, Anakinra, Tretinoin, Vitamin D, Genistein, BAY-1797, BAY-474, BI 665915, BAY-7598 and BAY-386. In addition to this, Venlafaxine and Paroxetine, which are marketed for the treatment of depression, were selected as positive drugs for affinity comparison. Molecular docking studies were performed using IL1R1 (PDB ID:4gaf) and the aforementioned small molecules to assess their binding energies with IL1R1. The results demonstrated strong interactions between these drug candidates and IL1R1 as evidenced by binding energies below −7; thus suggesting their potential efficacy in treating depression. Whereas, the docking scores with both Venlafaxine and Paroxetine were greater than −7.0, indicating that these drug candidates showed better affinity. Among them, BI 665915, BAY-7598, BAY-386, and BI 639667 exhibited the lowest binding energies and were further visualized through docking simulations with ILIRI (Figures 5D–G). To verify the binding stability, we performed 100 ns molecular dynamics simulations of BI 639667 with IL1R1. The root mean square deviation (RMSD) results showed that the fluctuations of the system were in the range of 0.1 nm, which indicated that the system had strong stability (Figure 5H). The root mean square fluctuation (RMSF) could reflect the overall amino acid residue fluctuation of the complex was small (Figure 5I). The radius of gyration (Rg) suggested that the complexes remained relatively stable during the simulation (Figure 5J). The solvent accessible surface area (SASA) indicated that the binding of BI 639667 had less effect on the structure of IL1R1 (Figure 5K). Further calculation of the hydrogen bonding contacts between BI 639667 and IL1R1 was consistent with the docking results, and the complexes were able to form stable hydrogen bonds (Figure 5L).