We used ${\mathrm{X}}_{\mathrm{i}}\in {\mathrm{R}}^{\mathrm{n}\times {\mathrm{p}}_{\mathrm{i}}}(i=1,2,3)$ to represent the region of interest (ROI), SNP, and gene expression data, respectively. $n$ is the total number of samples and ${\mathrm{p}}_{\mathrm{i}}(i=1,2,3)$ represents the number of features of the three data, respectively. ${\mathrm{X}}_{\mathrm{i}}$ can be decomposed into a basis matrix $\mathrm{W}\in {\mathrm{R}}^{\mathrm{n}\times \mathrm{k}}$ and three feature matrices ${\mathrm{H}}_{\mathrm{i}}\in {\mathrm{R}}^{\mathrm{k}\times {\mathrm{p}}_{\mathrm{i}}}(i=1,2,3)$. $\mathrm{k}$ is the number of dimensionality reduction. JDSNMF uses ${\mathrm{H}}_{\mathrm{i}}$ multilayer neural networks to capture nonlinear features based on the JDSNMF algorithm. The JDSNMF algorithm also controls the growth of the decomposed matrices by the F-parameter with the objective function shown below.

(1)$$\begin{array}{c}\hfill \min \sum _{i=1}^I{\parallel X_i-W{H}_{i_0}\parallel }_F^2+\lambda {\parallel M\parallel }_F^{\prime }\hfill \\ \hfill \text{s.t.}{H}_{i_0}=s\left(Z_{i_1}s\left(Z_{i_2}\mathrm{\dots }s\left(Z_{i_N}{H}_{i_N}\right)\right)\right)\hfill \\ \hfill H_{i_{n-1}}=s\left(Z_{i_n}{H}_{i_n}\right),H_{i_0}\mathrm{\dots }{H}_{i_{N-1}}\ge 0,\hfill \\ \hfill M\in \{U,Z_{i_1},\mathrm{\dots },Z_{i_N},H_{i_N}\},n=1,\mathrm{\dots },N\hfill \end{array}$$

$\mathrm{W}\in {ℝ}^{\mathrm{n}\times {\mathrm{k}}_0}$ is known as the sample potential matrix. ${\mathrm{H}}_{{\mathrm{i}}_0}\in {ℝ}^{{\mathrm{k}}_0\times {\mathrm{p}}_{\mathrm{i}}}$ is referred to as the feature potential matrix, which incorporates the first layer of the network. ${\mathrm{H}}_{{\mathrm{i}}_{\mathrm{n}}}\in {ℝ}^{{\mathrm{k}}_{\mathrm{n}}\times {\mathrm{p}}_{\mathrm{i}}}$ is called the characteristic potential feature of the $\mathrm{n}+1$ layer. $s(x)$ represents the activation function. The expression is as follows if it is a sigmoid activation function.

(2)$$s(x)=\frac{1}{1+e^{-x}}$$

If $s(x)$ is a tanh activation function, the expression is as follows.

(3)$$\tanh (x)=\frac{e^z-e^{-z}}{e^z+e^{-z}}$$

If $s(x)$ is the Rectified Linear Unit (ReLU) activation function, the expression can be defined as Eqn. 4.

(4)

According to the value of $k$, the algorithm will generate $k$ co-expression modules by using $H_{i_0}(i=1,2,3)$. We used z score to standardize each line of $H_{i_0}$, and then compared it with the artificially set threshold value T. If the normalized value was greater than the T value, it was considered that the index feature corresponding to this value was qualified to enter the co-expression module. In this study, we set the value of T to 1.

In the feature selection, we used the scikit-learn package for Python (v3.7, Python Software Foundation, Portland, OR, USA) to achieve the weight assignment to ROI, SNP, and genes in the co-expression module using the random forest (RF) algorithm. The parameters were explicitly set: ‘n_estimators’ was selected between 100 and 600, and ‘criterion’ was identified between ‘gini’ and ‘entropy’.

We performed a five-fold cross-validation on the training set to construct a diagnostic model using the GridSearchCV function. Finally, the optimal parameters were ‘entropy’ for ‘criterion’ and 500 for ‘n_estimators’. In addition, we constructed diagnostic models for ROI, SNP, and genes based on the logistic regression (LR) algorithm utilizing IBM SPSS Statistics 26 software (IBM SPSS Statistics, Chicago, IL, USA).

The BrainNet Viewer package of Matlab 2018a software (https://www.nitrc.org/projects/bnv/) was used to visualize the important brain regions selected by the JDSNMF algorithm.

The R package “clusterProfiler” (https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html) was applied to perform the Kyoto Encyclopedia of Genomes (KEGG) and Gene Ontology (GO) enrichment analysis of the important genes in the module. Bubble plots were then visualized using the R package “ggplot2” (v3.4.4, https://cran.r-project.org/web/packages/ggplot2/index.html).

AD, mild cognitive impairment (MCI), and healthy control (HC) samples were downloaded from the The Alzheimer’s Disease Neuroimaging Initiative (ADNI) database. The information on their sample set is shown in Table 1. For sMRI, we aligned the data to the Montreal Neurological Institute (MNI) standard space using the statistical parametric mapping (SPM) toolkit of Matlab software. Correction, segmentation, and alignment were then performed and, finally, the gray matter density of 90 brain regions with the cerebellar regions removed were extracted as the sMRI ROIs. Regarding the SNP data, we used the PLINK tool (https://www.cog-genomics.org/plink2/) to remove SNPs that did not meet the criteria for sex detection, Hardy-Weinberg equilibrium, and minor allele frequencies less than 0.05. The SNPs were then genetically annotated using ANNOVAR (http://annovar.openbioinformatics.org/), and 2378 SNPs within $\pm$5000 base pairs of multiple AD risk genes were extracted as input genetic data for the model. For the gene expression data, we obtained 414 genes differentially expressed in the HC and diseased groups using the limma algorithm. In addition, the differences in age and sex among the samples were explored based on the Wilcoxon test and Chi-squared test, respectively. The statistical results are presented in Table 2.

Since the JDSNMF algorithm was unsupervised, only the AD and MCI samples were used as input to the algorithm. The hyperparameters to be selected for the JDSNMF algorithm included the activation function, learning rate, and the number of dimensionality reduction. We set the number of iterations to 10,000. After fixing the other parameters, we first selected the activation function. The losses using the tanh, sigmoid, and rectified linear unit (ReLU) functions were 21,679.373, 4781.157, and 5563.1367, respectively. Therefore, we chose the sigmoid function for the subsequent analysis. Next, we selected the learning rate from the range of [0.1, 0.01, 0.001], with losses of 30,049, 4835, and 4781, respectively. Consequently, we set the learning rate to 0.001. Finally, we selected $k$. We selected it to be from 40 to 60. The changing loss with value trend is shown in Fig. 1A. This shows that loss takes the minimum value when it equals 60. Fig. 1B shows the line plot of the loss variation for the best combination of parameters. It can be seen that the loss tends to stabilize after iterating to 2000 for loss variation.

Fig. 1.

The process of parameter selection. (A) Histogram showing the loss at values between 40 and 60. (B) Line graph showing the variation of the loss for the best combination of parameters.

Since we set $k$ to 60, we obtained 60 co-expression modules. Module 42 did not contain any ROIs and was excluded. We calculated the reconstruction error for the other modules (Fig. 2). Among them, module 35 had the smallest reconstruction error. This module had 15 ROIs, 62 genes, and 490 SNPs.

Fig. 2.

Histogram showing the reconstruction error of co-expression modules.

In this section, we ranked the feature importance of ROIs, genes, and SNPs in module 35 based on the RF algorithm. The histograms of the top 10 feature weights are given in Fig. 3A–C, respectively. The diagnostic performance of these top markers was then explored, and diagnostic models were constructed using these ROIs, genes, and SNPs, respectively.

Fig. 3.

Heat maps showing weights of top features. (A–C) Weighted histograms of the top 10 regions of interest, top 10 single nucleotide polymorphisms, and top 10 genes, respectively.

We analyzed the correlation of the top 10 ROIs as shown in Fig. 4. Fig. 4A is a heat map drawn from the Pearson correlation coefficients among the 10 ROIs. Fig. 4B shows the connectivity of these ROIs in the brain template. There was a maximum positive correlation between rInfPar and rAng (corr = 0.9533), except for the correlation between each brain region and itself, which is 1. The maximum negative correlation can be found between rPal and rPal (corr = –0.9679). For the top 10 genes, Fig. 5 displays the results of their GO and KEGG enrichment analysis. We aimed to explore the biological significance of the top 10 ROIs, the top 10 genes, and the biological pathways in which they are involved. In addition, to confirm the algorithm’s performance for association analysis, we plotted the correlation heat map of the top ROIs and top SNPs (Fig. 6A) and the correlation heat map of the top ROIs and top genes (Fig. 6B). For Fig. 6A, rs4844384 and lPal had the highest positive correlation (corr = 0.5754). rs957191 and lPut had the highest negative correlation (corr = –0.4794). As shown in Fig. 6B, sulfiredoxin 1 (SRXN1) and lTha had the highest positive correlation (corr = 0.4962). collagen type III alpha 1 chain (COL3A1) and Iput had the highest negative correlation.

Fig. 4.

Correlation analysis and visualization of the top 10 ROIs. (A) Heat map showing the correlation between the top 10 ROIs. (B) Visualization of the top 10 ROIs and the relationship between them on the brain template. ROIs, region of interests.

Fig. 5.

Results of Gene Ontology and Kyoto Encyclopedia of Genomes enrichment analysis of the top 10 genes.

Fig. 6.

Heat maps showing the correlation between the top 10 ROIs and the top 10 SNPs and top 10 genes. (A) Heat map showing the correlation between the top 10 ROIs and the top 10 SNPs. (B) Heat map showing the correlation between the top 10 ROIs and the top 10 genes. SNPs, single nucleotide polymorphism.

In order to construct a diagnostic model on top markers, we constructed a diagnostic model for AD using the top 10 ROIs, top 10 SNPs, and top 10 genes based on the LR algorithm, respectively. Fig. 7A–C shows the receiver operating characteristic (ROC) curves of the three diagnostic models. The top 10 genes in the test set area under the curve (AUC) can reach the maximum (AUC = 0.947). The respective ROC curves of the top 10 ROIs, SNPs, and genes are presented in Figs. 8,9,10. As can be seen from the figures, the majority of top markers have AUCs greater than 0.5.

Fig. 7.

Diagnostic model construction based on top markers. (A–C) ROC curves of the top 10 ROIs, top 10 SNPs, and top 10 genes, respectively. ROC, receiver operating characteristic; TPR, true positive rate; FPR, false positive rate.

Fig. 8.

Diagnostic performance validation of the top 10 ROIs. (A–J) ROC curves of the top 10 ROIs, respectively.

Fig. 9.

Diagnostic performance validation of the top 10 SNPs. (A–J) ROC curves of the top 10 SNPs, respectively.

Fig. 10.

Diagnostic performance validation of the top 10 genes. (A–J) ROC curves of the top 10 genes, respectively.

We compared the performance of the JDSNMF algorithm with the JNMF algorithm and the JCB-SNMF algorithm for correlation analysis (Table 3). Specifically, we introduced the Pearson correlation coefficients of the original and reconstructed matrices for comparison to achieve similarity between the original matrix and the reconstructed matrix after decomposition.

$Corr(X_i,W{H}_i)(i=1,2,3)$ in the table represents the Pearson correlation coefficients between the original expression matrix and the reconstructed matrix for ROIs, genes, and SNPs, respectively. We found that the JDSNMF algorithm had the best performance.