# Identification of Signature Genes and Construction of an Artificial Neural Network Model of Prostate Cancer

Authors: Hongye Dong; Xu Wang
Journal: Journal of Healthcare Engineering
 (2022)

---

## Abstract

This study aimed to establish an artificial neural network (ANN) model based on prostate cancer signature genes (PCaSGs) to predict the patients with prostate cancer (PCa). In the present study, 270 differentially expressed genes (DEGs) were identified between PCa and normal prostate (NP) groups by differential gene expression analysis. Next, we performed Metascape gene annotation, pathway and process enrichment analysis, and PPI enrichment analysis on all 270 DEGs. Then, we identified and screened out 30 PCaSGs based on the random forest analysis and constructed an ANN model based on the gene score matrix consisting of 30 PCaSGs. Lastly, analysis of microarray dataset GSE46602 showed that the accuracy of this model for predicating PCa and NP samples was 88.9 and 78.6%, respectively. Our results suggested that the ANN model based on PCaSGs can be used for effectively predicting the patients with PCa and will be helpful for early PCa diagnosis and treatment.


---

## Body

## 1. Introduction

Prostate cancer (PCa) is a tumor caused by malignant hyperplasia of prostate epithelial cells. It has a very high incidence in elderly men, with 80% of cases occurring in men over 65 years old [1, 2]. In the early stage of PCa, most patients have no obvious symptoms due to the insidious onset and slow growth of the tumor [3]. Once PCa is advanced, it can cause symptoms such as abnormal urination, pelvic discomfort, erectile dysfunction, and even bone pain and spinal cord compression, which can greatly affect the quality of life of patients [4, 5]. Accordingly, there is an urgent need to develop effective biological approaches to improve diagnosis and prognosis of PCa.Over the past few decades, various computer-aided diagnostic models have been used to predict the risk of various cancers, such as logistic regression, Cox proportional risk models, and decision trees [6–8]. Artificial neural network (ANN) is a mathematical or computational model that uses structures similar to synaptic connections in the brain to process information [9]. ANN models have been applied to risk assessment of many diseases, including colon cancer, lung cancer, hepatocellular carcinoma, meningioma, and so on and have shown reliable and accurate performance in disease prediction and evaluation [10–13]. However, no studies have been reported on predicting prostate cancer risk based on ANN models.In this study, we downloaded RNA-Seq data from PCa and normal prostate (NP) samples from Gene Expression Omnibus (GEO) database, identified differentially expressed genes (DEGs), followed by Metascape gene list analysis and random forest analysis. An ANN model was established according to gene score calculation for PCa signature genes (PCaSGs) in samples. In addition, the reliability of ANN model prediction was validated by drawing a ROC curve and an independent microarray dataset of PCa, GSE46602. Microarray dataset GSE46602 has been utilized to calculate gene scores for further testing the accuracy of the ANN model. Our study results could provide new insights for identifying those patients with PCa.


## 2. Materials and Methods

### 2.1. Data Downloaded and Collated from the GEO Database

Firstly, we selected and downloaded three independent datasets (GSE60329, GSE71016, and GSE46602) and corresponding clinical information from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) [14]. Only datasets with effective sample size greater than 50, PCa group and control group, complete clinical follow-up information, and complete transcriptome expression matrix were accepted to ensure reliability of the findings. Then, each GEO dataset was annotated according to the platform annotation file, and the probe IDs were converted into gene symbols to obtain the whole gene expression matrix. Finally, two microarray gene expression datasets were merged as the training set (including 102 PCa and 61 normal prostate (NP) samples) and the other one as the testing set (including 36 PCa and 14 NP samples).


### 2.2. Differential Expression Analysis of the Training Set

To compare genes differentially expressed in PCa and NP groups, we performed differential expression analysis between groups by statistical tests using theR packages limma, pheatmap, and ggplot [15, 16] based on the gene expression matrix of the training set. Results were presented in a heatmap and a volcano plot. The cut-off criteria of DEGs were adjusted to P value < 0.05 and | logFC | > 1.


### 2.3. Metascape Gene List Analysis on DEGs

Metascape is a powerful gene function annotation analysis tool that enables researchers to apply currently popular bioinformatics analysis methods to batch gene and protein analysis to achieve knowledge of gene or protein function [17, 18]. We choose the Metascape because the database is updated monthly to ensure the reliability of data. Gene annotation, pathway and process enrichment analysis, and protein-protein interaction (PPI) enrichment analysis were performed on DEGs using Metascape. The list of DEGs was entered, and ‘Homo sapiens’ was selected as the organism.


### 2.4. Random Forest Analysis on DEGs

Random forest analysis is an analysis method that uses decision tree algorithm to evaluate the importance of variables [19]. With this algorithm, we could filter the DEGs to find the disease signature genes. We first constructed a random-forest model using 500 trees on the training set using the R package randomForest [20]. Then, we calculated the point with the minimum cross-validation error to find the optimal number of trees for random forest [21]. Next, we then ranked the importance of DEGs and selected the Top30 DEGs with the highest importance score and named them as PCaSGs. Finally, the expression of important PCaSGs was output and visualized with a heat map using R packages limma and pheatmap [22].


### 2.5. Gene Score Calculation for PCaSGs in Samples

Batch effects, in simple terms, are incidental deviations in data that have nothing to do with the results of an experiment [23]. Therefore, in order to remove the batch effects of samples from different sources, we calculated gene scores for each PCaSGs in each sample [24]. Firstly, the expression matrix of PCaSGs and corresponding lgFC values were input into R software. Then, the relative expression quantities of PCaSGs were compared with the median expression value, if the quantity of up-regulated gene was higher than the median value, the gene score was marked as 1, otherwise marked as 0; if the quantity of downregulated gene was lower than the median value, the gene score was marked as 1, otherwise marked as 0. Finally, the results of all gene scores were output.


### 2.6. Construction of an ANN Model

The ANN model is a simplified model that mimics the way the human brain processes information. The model works by simulating a large number of abstract interconnection processing units similar to neurons [25, 26]. To test the reliability and accuracy of gene scoring results, we constructed a neural network model based on 30 PCaSGs using R packages neuralnet and NeuralNetTools [27, 28]. We imported the gene score data of 30 PCaSGs as the input layer, and set 5 nodes as the middle hidden layer. These units received training feedback through variable connection strength (or weight), next output results from the output layer [29]. The gene score of each sample was compared between the PCa group and the NP group to predict which group the sample belonged to. Finally, we draw ROC curve to verify the reliability of ANN model prediction using R package pROC [30].


### 2.7. Gene Score Calculation in the Testing Set

Firstly, the transcriptome expression matrix of the testing set and corresponding lgFC values of DEGs were input intoR software. Then, the relative expression quantities of DEGs were compared with the median expression value, if the quantity of upregulated gene was higher than the median value, the gene score was marked as 1, otherwise marked as 0; if the quantity of downregulated gene was lower than the median value, the gene score was marked as 1, otherwise marked as 0. Finally, the result of all gene scores were output.


### 2.8. Prediction Performance of the ANN Model in the Testing Set

In order to further test the accuracy of the ANN model constructed based on gene scores, we used the ANN model based on 30 PCaSGs to calculate the scores of all samples in the testing set and predicted which group the samples belonged to by comparing the scores of the PCa group and the NP group. Then, we combined the prediction results of the ANN model with the real grouping information to calculate the accuracy of the model prediction. Finally, we draw the ROC curve to verify the reliability of ANN model prediction usingR package pROC [30].


## 2.1. Data Downloaded and Collated from the GEO Database

Firstly, we selected and downloaded three independent datasets (GSE60329, GSE71016, and GSE46602) and corresponding clinical information from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/) [14]. Only datasets with effective sample size greater than 50, PCa group and control group, complete clinical follow-up information, and complete transcriptome expression matrix were accepted to ensure reliability of the findings. Then, each GEO dataset was annotated according to the platform annotation file, and the probe IDs were converted into gene symbols to obtain the whole gene expression matrix. Finally, two microarray gene expression datasets were merged as the training set (including 102 PCa and 61 normal prostate (NP) samples) and the other one as the testing set (including 36 PCa and 14 NP samples).


## 2.2. Differential Expression Analysis of the Training Set

To compare genes differentially expressed in PCa and NP groups, we performed differential expression analysis between groups by statistical tests using theR packages limma, pheatmap, and ggplot [15, 16] based on the gene expression matrix of the training set. Results were presented in a heatmap and a volcano plot. The cut-off criteria of DEGs were adjusted to P value < 0.05 and | logFC | > 1.


## 2.3. Metascape Gene List Analysis on DEGs

Metascape is a powerful gene function annotation analysis tool that enables researchers to apply currently popular bioinformatics analysis methods to batch gene and protein analysis to achieve knowledge of gene or protein function [17, 18]. We choose the Metascape because the database is updated monthly to ensure the reliability of data. Gene annotation, pathway and process enrichment analysis, and protein-protein interaction (PPI) enrichment analysis were performed on DEGs using Metascape. The list of DEGs was entered, and ‘Homo sapiens’ was selected as the organism.


## 2.4. Random Forest Analysis on DEGs

Random forest analysis is an analysis method that uses decision tree algorithm to evaluate the importance of variables [19]. With this algorithm, we could filter the DEGs to find the disease signature genes. We first constructed a random-forest model using 500 trees on the training set using the R package randomForest [20]. Then, we calculated the point with the minimum cross-validation error to find the optimal number of trees for random forest [21]. Next, we then ranked the importance of DEGs and selected the Top30 DEGs with the highest importance score and named them as PCaSGs. Finally, the expression of important PCaSGs was output and visualized with a heat map using R packages limma and pheatmap [22].


## 2.5. Gene Score Calculation for PCaSGs in Samples

Batch effects, in simple terms, are incidental deviations in data that have nothing to do with the results of an experiment [23]. Therefore, in order to remove the batch effects of samples from different sources, we calculated gene scores for each PCaSGs in each sample [24]. Firstly, the expression matrix of PCaSGs and corresponding lgFC values were input into R software. Then, the relative expression quantities of PCaSGs were compared with the median expression value, if the quantity of up-regulated gene was higher than the median value, the gene score was marked as 1, otherwise marked as 0; if the quantity of downregulated gene was lower than the median value, the gene score was marked as 1, otherwise marked as 0. Finally, the results of all gene scores were output.


## 2.6. Construction of an ANN Model

The ANN model is a simplified model that mimics the way the human brain processes information. The model works by simulating a large number of abstract interconnection processing units similar to neurons [25, 26]. To test the reliability and accuracy of gene scoring results, we constructed a neural network model based on 30 PCaSGs using R packages neuralnet and NeuralNetTools [27, 28]. We imported the gene score data of 30 PCaSGs as the input layer, and set 5 nodes as the middle hidden layer. These units received training feedback through variable connection strength (or weight), next output results from the output layer [29]. The gene score of each sample was compared between the PCa group and the NP group to predict which group the sample belonged to. Finally, we draw ROC curve to verify the reliability of ANN model prediction using R package pROC [30].


## 2.7. Gene Score Calculation in the Testing Set

Firstly, the transcriptome expression matrix of the testing set and corresponding lgFC values of DEGs were input intoR software. Then, the relative expression quantities of DEGs were compared with the median expression value, if the quantity of upregulated gene was higher than the median value, the gene score was marked as 1, otherwise marked as 0; if the quantity of downregulated gene was lower than the median value, the gene score was marked as 1, otherwise marked as 0. Finally, the result of all gene scores were output.


## 2.8. Prediction Performance of the ANN Model in the Testing Set

In order to further test the accuracy of the ANN model constructed based on gene scores, we used the ANN model based on 30 PCaSGs to calculate the scores of all samples in the testing set and predicted which group the samples belonged to by comparing the scores of the PCa group and the NP group. Then, we combined the prediction results of the ANN model with the real grouping information to calculate the accuracy of the model prediction. Finally, we draw the ROC curve to verify the reliability of ANN model prediction usingR package pROC [30].


## 3. Results

### 3.1. Identification of DEGs between PCa and NP Groups

Firstly, we obtained a gene expression matrix containing 22,014 genes by merging and cleaning of the datasets GSE60329 and GSE71016. Then, 270 DEGs were identified between PCa and NP groups by differential gene expression analysis, with 155 downregulated and 115 upregulated. The details of the expression matrix of DEGs were given in Supplementary File S1 (diff.xls and diffGeneExp.xls), and the result was represented by the heatmap (Figure1(a)) and the volcano plot (Figure 1(b)).Figure 1
The heatmap (a) and volcano plot (b) of the DEGs. Red, upregulated DEGs; blue or green, downregulated DEGs. Cut-off criteria: adj.P < 0.05 and | logFC | > 1.
(a)(b)


#### 3.1.1. Gene Annotation and Enrichment Analysis on DEGs

We performed a series of Metascape gene annotation, pathway and process enrichment analysis, and PPI enrichment analysis on all 270 DEGs. The 270 DEGs annotation and enrichment information were detailed in the Supplementary File S2 (metascape_result.xls). Figure2(a) summarizes the enrichment of DEGs functions or pathways. Terms with a P value < 0.01, a minimum count of 3, and an enrichment factor >1.5 (the enrichment factor is the ratio between the observed counts and the counts expected by chance) were collected and grouped into clusters based on their membership similarities (see Table 1). To further capture the relationships between the terms, a subset of enriched terms had been selected and rendered as a network plot, where terms with a similarity >0.3 were connected by edges. The networks were visualized using Cytoscape, where each node represents an enriched term and was colored first by its cluster ID (Figure 2(b)) and then by its P value (Figure 2(c)). For the gene list composed of 270 DEGs, PPI enrichment analysis was carried out with STRING and BioGrid databases. The PPI network and MCODE components identified for the gene list was gathered and shown in Figure 2(d), 2(e). And pathway and process enrichment analysis had been applied to each MCODE component independently, and the three best-scoring terms by P value were retained as the functional description of the corresponding components (see Figure 2(b)).Figure 2
(a) Bar graph of enriched terms across input gene lists, colored byp values. Network of enriched terms: (b) colored by cluster ID, where nodes that share the same cluster ID are typically close to each other; (c) colored by p value, where terms containing more genes tend to have a more significant p value. Protein-protein interaction network (d) and MCODE components (e) identified in the DEGs list.
(a)(b)(c)(d)(e)


#### 3.1.2. Construction of the Random Forest Tree and Identification of PCaSGs

We ultimately identified and screened out 30 PCaSGs based on the random forest analysis. As shown in Figure3(a), the cross-validation error was minimized when the tree number reached 333. Figure 3(b) shows the top 30 DEGs on the importance scale, called PCaSGs. The expression of 30 PCaSGs in each PCa or NP sample was visualized (Figure 3(c)), and the further details were provided in the supplementary file S3 (rfGeneExp.xls). It could be seen from Figure 3(c) that the PCaSGs in the two groups have relatively obvious hierarchical clustering, indicating that the PCaSGs expression levels obtained through random forest tree analysis can distinguish whether a sample is in the PCa group or not.Figure 3
(a) The random forest tree. The abscissa represents the number of trees, and the ordinate represents the error of the cross-validation. The red, green, and black curves represent the error of PCa, NP, and all sample groups, respectively. (b) The bubble plot of PCaSGs. The abscissa represents the important score of genes, and the ordinate represents PCaSGs. (c) The heatmap of the Top30 PCaSGs. Red, upregulated genes; blue, downregulated genes. Cut-off criteria: adj.P < 0.05 and | logFC | > 1.
(a)(b)(c)


### 3.2. Gene Score for 30 PCaSGs in 163 Samples

After batch effect correction of 163 samples from different sources in multiple datasets, we obtained a gene score matrix consisting of 30 PCaSGs. The further details of the matrix were presented in the supplementary file S4 (geneScore.xls).


#### 3.2.1. Construction of the ANN Model Based 30 PCaSGs

We constructed the ANN model based the gene score matrix consisting of 30 PCaSGs in 163 samples usingR packages (see Figure 4(a)). Figure 3(a) shows the weights of the input layer (a gene score matrix consisting of 30 PCaSGs) to the hidden layer (consisting of 5 nodes), and Figure 3(b) shows the weights of the hidden layer to the output layer (representing the grouping of samples). As could be seen from Figure 3(c) and Figure 4(b), the accuracy of prediction of NP group by neural network model was 98.4% and that of PCa group was 97.1%, and the area under the ROC curve (AUC) of the training set was 0.998. The above results indicated that the ANN model constructed had high accuracy and reliability.Figure 4
(a) The ANN model based on the gene score of 30 PCaSGs. The circles in the left column are the input layer units (gene scores of 30 PCaSGs), the circles in the middle column are the hidden layer units (consisting of 5 nodes), and the circles in the right column are the output layer units (two groups). (b) The ROC curve of the ANN model. AUC: area under curve. 95% CI : the 95% confidence interval. Abscissa: 1-specificity (false positive rate); ordinate: sensitivity (true positive rate).
(a)(b)


### 3.3. Gene Score for 241 DEGs in 50 Samples

After batch effect correction of 50 samples of the datasets GSE46602, we obtained a gene score matrix consisting of 241 DEGs. The further details of the matrix were presented in the supplementary file S5 (testGeneScore.xls).


#### 3.3.1. Verification of the ANN Model by Testing Set

We used the ANN model based on 30 PCaSGs to calculate the scores of all 50 samples in the testing set. If the score of a sample in the PCa group was higher than that in the NP group, the sample was predicted to belong to the PCa group, otherwise it belonged to the NP group. The scoring matrix for each sample was detailed in the supplementary file S6 (test.neuralPredict.xls). As shown in Figures4 and 5, the accuracy of prediction of NP group by the ANN model in testing set was 78.6% and that of PCa group was 88.9%, and the AUC of the testing set was 0.869. The above results indicated that the prediction model constructed were credible after verification of the testing set.Figure 5
The ROC curve of the ANN model in testing set. AUC: area under curve. 95% CI : the 95% confidence interval. Abscissa: 1-specificity (false positive rate); ordinate: sensitivity (true positive rate).Table 1
Top 20 clusters with their representative enriched terms (one per cluster).
GOCategoryDescriptionCount%Log 10 (P)Log 10 (q)GO: 0006954GO Biological Processesinflammatory response218.02−8.4−4.06R-HSA-9031628Reactome Gene SetsNGF-stimulated transcription72.67−7.38−3.34R-HSA-8953897Reactome Gene SetsCellular responses to stimuli238.78−6.49−2.62GO: 0009636GO Biological Processesresponse to toxic substance124.58−6.06−2.56GO: 0045638GO Biological Processesnegative regulation of myeloid cell differentiation83.05−5.88−2.49GO: 0008015GO Biological Processesblood circulation155.73−5.59−2.27M5885Canonical PathwaysNABA MATRISOME ASSOCIATED218.02−5.53−2.27R-HSA-195258Reactome Gene SetsRHO GTPase Effectors134.96−5.19−2.12GO: 0050900GO Biological Processesleukocyte migration114.2−5.19−2.12GO: 0008202GO Biological Processessteroid metabolic process114.2−4.84−1.85GO: 1903034GO Biological Processesregulation of response to wounding93.44−4.81−1.84WP383WikiPathwaysStriated muscle contraction pathway51.91−4.73−1.8GO: 0045229GO Biological Processesexternal encapsulating structure organization114.2−4.59−1.68GO: 0009725GO Biological Processesresponse to hormone197.25−4.53−1.65WP465WikiPathwaysTryptophan metabolism51.91−4.51−1.65GO: 0015840GO Biological Processesurea transport31.15−4.46−1.65GO: 0002526GO Biological Processesacute inflammatory response62.29−4.31−1.58GO: 0032200GO Biological Processestelomere organization72.67−4.3−1.57GO: 2000147GO Biological Processespositive regulation of cell motility166.11−4.24−1.54GO: 0030856GO Biological Processesregulation of epithelial cell differentiation83.05−4.14−1.49“Count” is the number of genes in the user-provided lists with membership in the given ontology term. “%” is the percentage of all of the user-provided genes that are found in the given ontology term (only input genes with at least one ontology term annotation are included in the calculation). “Log 10 (P)” is the p-value in log base 10. “Log 10 (q)” is the multi-test adjusted p-value in log base 10.


## 3.1. Identification of DEGs between PCa and NP Groups

Firstly, we obtained a gene expression matrix containing 22,014 genes by merging and cleaning of the datasets GSE60329 and GSE71016. Then, 270 DEGs were identified between PCa and NP groups by differential gene expression analysis, with 155 downregulated and 115 upregulated. The details of the expression matrix of DEGs were given in Supplementary File S1 (diff.xls and diffGeneExp.xls), and the result was represented by the heatmap (Figure1(a)) and the volcano plot (Figure 1(b)).Figure 1
The heatmap (a) and volcano plot (b) of the DEGs. Red, upregulated DEGs; blue or green, downregulated DEGs. Cut-off criteria: adj.P < 0.05 and | logFC | > 1.
(a)(b)


### 3.1.1. Gene Annotation and Enrichment Analysis on DEGs

We performed a series of Metascape gene annotation, pathway and process enrichment analysis, and PPI enrichment analysis on all 270 DEGs. The 270 DEGs annotation and enrichment information were detailed in the Supplementary File S2 (metascape_result.xls). Figure2(a) summarizes the enrichment of DEGs functions or pathways. Terms with a P value < 0.01, a minimum count of 3, and an enrichment factor >1.5 (the enrichment factor is the ratio between the observed counts and the counts expected by chance) were collected and grouped into clusters based on their membership similarities (see Table 1). To further capture the relationships between the terms, a subset of enriched terms had been selected and rendered as a network plot, where terms with a similarity >0.3 were connected by edges. The networks were visualized using Cytoscape, where each node represents an enriched term and was colored first by its cluster ID (Figure 2(b)) and then by its P value (Figure 2(c)). For the gene list composed of 270 DEGs, PPI enrichment analysis was carried out with STRING and BioGrid databases. The PPI network and MCODE components identified for the gene list was gathered and shown in Figure 2(d), 2(e). And pathway and process enrichment analysis had been applied to each MCODE component independently, and the three best-scoring terms by P value were retained as the functional description of the corresponding components (see Figure 2(b)).Figure 2
(a) Bar graph of enriched terms across input gene lists, colored byp values. Network of enriched terms: (b) colored by cluster ID, where nodes that share the same cluster ID are typically close to each other; (c) colored by p value, where terms containing more genes tend to have a more significant p value. Protein-protein interaction network (d) and MCODE components (e) identified in the DEGs list.
(a)(b)(c)(d)(e)


### 3.1.2. Construction of the Random Forest Tree and Identification of PCaSGs

We ultimately identified and screened out 30 PCaSGs based on the random forest analysis. As shown in Figure3(a), the cross-validation error was minimized when the tree number reached 333. Figure 3(b) shows the top 30 DEGs on the importance scale, called PCaSGs. The expression of 30 PCaSGs in each PCa or NP sample was visualized (Figure 3(c)), and the further details were provided in the supplementary file S3 (rfGeneExp.xls). It could be seen from Figure 3(c) that the PCaSGs in the two groups have relatively obvious hierarchical clustering, indicating that the PCaSGs expression levels obtained through random forest tree analysis can distinguish whether a sample is in the PCa group or not.Figure 3
(a) The random forest tree. The abscissa represents the number of trees, and the ordinate represents the error of the cross-validation. The red, green, and black curves represent the error of PCa, NP, and all sample groups, respectively. (b) The bubble plot of PCaSGs. The abscissa represents the important score of genes, and the ordinate represents PCaSGs. (c) The heatmap of the Top30 PCaSGs. Red, upregulated genes; blue, downregulated genes. Cut-off criteria: adj.P < 0.05 and | logFC | > 1.
(a)(b)(c)


## 3.1.1. Gene Annotation and Enrichment Analysis on DEGs

We performed a series of Metascape gene annotation, pathway and process enrichment analysis, and PPI enrichment analysis on all 270 DEGs. The 270 DEGs annotation and enrichment information were detailed in the Supplementary File S2 (metascape_result.xls). Figure2(a) summarizes the enrichment of DEGs functions or pathways. Terms with a P value < 0.01, a minimum count of 3, and an enrichment factor >1.5 (the enrichment factor is the ratio between the observed counts and the counts expected by chance) were collected and grouped into clusters based on their membership similarities (see Table 1). To further capture the relationships between the terms, a subset of enriched terms had been selected and rendered as a network plot, where terms with a similarity >0.3 were connected by edges. The networks were visualized using Cytoscape, where each node represents an enriched term and was colored first by its cluster ID (Figure 2(b)) and then by its P value (Figure 2(c)). For the gene list composed of 270 DEGs, PPI enrichment analysis was carried out with STRING and BioGrid databases. The PPI network and MCODE components identified for the gene list was gathered and shown in Figure 2(d), 2(e). And pathway and process enrichment analysis had been applied to each MCODE component independently, and the three best-scoring terms by P value were retained as the functional description of the corresponding components (see Figure 2(b)).Figure 2
(a) Bar graph of enriched terms across input gene lists, colored byp values. Network of enriched terms: (b) colored by cluster ID, where nodes that share the same cluster ID are typically close to each other; (c) colored by p value, where terms containing more genes tend to have a more significant p value. Protein-protein interaction network (d) and MCODE components (e) identified in the DEGs list.
(a)(b)(c)(d)(e)


## 3.1.2. Construction of the Random Forest Tree and Identification of PCaSGs

We ultimately identified and screened out 30 PCaSGs based on the random forest analysis. As shown in Figure3(a), the cross-validation error was minimized when the tree number reached 333. Figure 3(b) shows the top 30 DEGs on the importance scale, called PCaSGs. The expression of 30 PCaSGs in each PCa or NP sample was visualized (Figure 3(c)), and the further details were provided in the supplementary file S3 (rfGeneExp.xls). It could be seen from Figure 3(c) that the PCaSGs in the two groups have relatively obvious hierarchical clustering, indicating that the PCaSGs expression levels obtained through random forest tree analysis can distinguish whether a sample is in the PCa group or not.Figure 3
(a) The random forest tree. The abscissa represents the number of trees, and the ordinate represents the error of the cross-validation. The red, green, and black curves represent the error of PCa, NP, and all sample groups, respectively. (b) The bubble plot of PCaSGs. The abscissa represents the important score of genes, and the ordinate represents PCaSGs. (c) The heatmap of the Top30 PCaSGs. Red, upregulated genes; blue, downregulated genes. Cut-off criteria: adj.P < 0.05 and | logFC | > 1.
(a)(b)(c)


## 3.2. Gene Score for 30 PCaSGs in 163 Samples

After batch effect correction of 163 samples from different sources in multiple datasets, we obtained a gene score matrix consisting of 30 PCaSGs. The further details of the matrix were presented in the supplementary file S4 (geneScore.xls).


### 3.2.1. Construction of the ANN Model Based 30 PCaSGs

We constructed the ANN model based the gene score matrix consisting of 30 PCaSGs in 163 samples usingR packages (see Figure 4(a)). Figure 3(a) shows the weights of the input layer (a gene score matrix consisting of 30 PCaSGs) to the hidden layer (consisting of 5 nodes), and Figure 3(b) shows the weights of the hidden layer to the output layer (representing the grouping of samples). As could be seen from Figure 3(c) and Figure 4(b), the accuracy of prediction of NP group by neural network model was 98.4% and that of PCa group was 97.1%, and the area under the ROC curve (AUC) of the training set was 0.998. The above results indicated that the ANN model constructed had high accuracy and reliability.Figure 4
(a) The ANN model based on the gene score of 30 PCaSGs. The circles in the left column are the input layer units (gene scores of 30 PCaSGs), the circles in the middle column are the hidden layer units (consisting of 5 nodes), and the circles in the right column are the output layer units (two groups). (b) The ROC curve of the ANN model. AUC: area under curve. 95% CI : the 95% confidence interval. Abscissa: 1-specificity (false positive rate); ordinate: sensitivity (true positive rate).
(a)(b)


## 3.2.1. Construction of the ANN Model Based 30 PCaSGs

We constructed the ANN model based the gene score matrix consisting of 30 PCaSGs in 163 samples usingR packages (see Figure 4(a)). Figure 3(a) shows the weights of the input layer (a gene score matrix consisting of 30 PCaSGs) to the hidden layer (consisting of 5 nodes), and Figure 3(b) shows the weights of the hidden layer to the output layer (representing the grouping of samples). As could be seen from Figure 3(c) and Figure 4(b), the accuracy of prediction of NP group by neural network model was 98.4% and that of PCa group was 97.1%, and the area under the ROC curve (AUC) of the training set was 0.998. The above results indicated that the ANN model constructed had high accuracy and reliability.Figure 4
(a) The ANN model based on the gene score of 30 PCaSGs. The circles in the left column are the input layer units (gene scores of 30 PCaSGs), the circles in the middle column are the hidden layer units (consisting of 5 nodes), and the circles in the right column are the output layer units (two groups). (b) The ROC curve of the ANN model. AUC: area under curve. 95% CI : the 95% confidence interval. Abscissa: 1-specificity (false positive rate); ordinate: sensitivity (true positive rate).
(a)(b)


## 3.3. Gene Score for 241 DEGs in 50 Samples

After batch effect correction of 50 samples of the datasets GSE46602, we obtained a gene score matrix consisting of 241 DEGs. The further details of the matrix were presented in the supplementary file S5 (testGeneScore.xls).


### 3.3.1. Verification of the ANN Model by Testing Set

We used the ANN model based on 30 PCaSGs to calculate the scores of all 50 samples in the testing set. If the score of a sample in the PCa group was higher than that in the NP group, the sample was predicted to belong to the PCa group, otherwise it belonged to the NP group. The scoring matrix for each sample was detailed in the supplementary file S6 (test.neuralPredict.xls). As shown in Figures4 and 5, the accuracy of prediction of NP group by the ANN model in testing set was 78.6% and that of PCa group was 88.9%, and the AUC of the testing set was 0.869. The above results indicated that the prediction model constructed were credible after verification of the testing set.Figure 5
The ROC curve of the ANN model in testing set. AUC: area under curve. 95% CI : the 95% confidence interval. Abscissa: 1-specificity (false positive rate); ordinate: sensitivity (true positive rate).Table 1
Top 20 clusters with their representative enriched terms (one per cluster).
GOCategoryDescriptionCount%Log 10 (P)Log 10 (q)GO: 0006954GO Biological Processesinflammatory response218.02−8.4−4.06R-HSA-9031628Reactome Gene SetsNGF-stimulated transcription72.67−7.38−3.34R-HSA-8953897Reactome Gene SetsCellular responses to stimuli238.78−6.49−2.62GO: 0009636GO Biological Processesresponse to toxic substance124.58−6.06−2.56GO: 0045638GO Biological Processesnegative regulation of myeloid cell differentiation83.05−5.88−2.49GO: 0008015GO Biological Processesblood circulation155.73−5.59−2.27M5885Canonical PathwaysNABA MATRISOME ASSOCIATED218.02−5.53−2.27R-HSA-195258Reactome Gene SetsRHO GTPase Effectors134.96−5.19−2.12GO: 0050900GO Biological Processesleukocyte migration114.2−5.19−2.12GO: 0008202GO Biological Processessteroid metabolic process114.2−4.84−1.85GO: 1903034GO Biological Processesregulation of response to wounding93.44−4.81−1.84WP383WikiPathwaysStriated muscle contraction pathway51.91−4.73−1.8GO: 0045229GO Biological Processesexternal encapsulating structure organization114.2−4.59−1.68GO: 0009725GO Biological Processesresponse to hormone197.25−4.53−1.65WP465WikiPathwaysTryptophan metabolism51.91−4.51−1.65GO: 0015840GO Biological Processesurea transport31.15−4.46−1.65GO: 0002526GO Biological Processesacute inflammatory response62.29−4.31−1.58GO: 0032200GO Biological Processestelomere organization72.67−4.3−1.57GO: 2000147GO Biological Processespositive regulation of cell motility166.11−4.24−1.54GO: 0030856GO Biological Processesregulation of epithelial cell differentiation83.05−4.14−1.49“Count” is the number of genes in the user-provided lists with membership in the given ontology term. “%” is the percentage of all of the user-provided genes that are found in the given ontology term (only input genes with at least one ontology term annotation are included in the calculation). “Log 10 (P)” is the p-value in log base 10. “Log 10 (q)” is the multi-test adjusted p-value in log base 10.


## 3.3.1. Verification of the ANN Model by Testing Set

We used the ANN model based on 30 PCaSGs to calculate the scores of all 50 samples in the testing set. If the score of a sample in the PCa group was higher than that in the NP group, the sample was predicted to belong to the PCa group, otherwise it belonged to the NP group. The scoring matrix for each sample was detailed in the supplementary file S6 (test.neuralPredict.xls). As shown in Figures4 and 5, the accuracy of prediction of NP group by the ANN model in testing set was 78.6% and that of PCa group was 88.9%, and the AUC of the testing set was 0.869. The above results indicated that the prediction model constructed were credible after verification of the testing set.Figure 5
The ROC curve of the ANN model in testing set. AUC: area under curve. 95% CI : the 95% confidence interval. Abscissa: 1-specificity (false positive rate); ordinate: sensitivity (true positive rate).Table 1
Top 20 clusters with their representative enriched terms (one per cluster).
GOCategoryDescriptionCount%Log 10 (P)Log 10 (q)GO: 0006954GO Biological Processesinflammatory response218.02−8.4−4.06R-HSA-9031628Reactome Gene SetsNGF-stimulated transcription72.67−7.38−3.34R-HSA-8953897Reactome Gene SetsCellular responses to stimuli238.78−6.49−2.62GO: 0009636GO Biological Processesresponse to toxic substance124.58−6.06−2.56GO: 0045638GO Biological Processesnegative regulation of myeloid cell differentiation83.05−5.88−2.49GO: 0008015GO Biological Processesblood circulation155.73−5.59−2.27M5885Canonical PathwaysNABA MATRISOME ASSOCIATED218.02−5.53−2.27R-HSA-195258Reactome Gene SetsRHO GTPase Effectors134.96−5.19−2.12GO: 0050900GO Biological Processesleukocyte migration114.2−5.19−2.12GO: 0008202GO Biological Processessteroid metabolic process114.2−4.84−1.85GO: 1903034GO Biological Processesregulation of response to wounding93.44−4.81−1.84WP383WikiPathwaysStriated muscle contraction pathway51.91−4.73−1.8GO: 0045229GO Biological Processesexternal encapsulating structure organization114.2−4.59−1.68GO: 0009725GO Biological Processesresponse to hormone197.25−4.53−1.65WP465WikiPathwaysTryptophan metabolism51.91−4.51−1.65GO: 0015840GO Biological Processesurea transport31.15−4.46−1.65GO: 0002526GO Biological Processesacute inflammatory response62.29−4.31−1.58GO: 0032200GO Biological Processestelomere organization72.67−4.3−1.57GO: 2000147GO Biological Processespositive regulation of cell motility166.11−4.24−1.54GO: 0030856GO Biological Processesregulation of epithelial cell differentiation83.05−4.14−1.49“Count” is the number of genes in the user-provided lists with membership in the given ontology term. “%” is the percentage of all of the user-provided genes that are found in the given ontology term (only input genes with at least one ontology term annotation are included in the calculation). “Log 10 (P)” is the p-value in log base 10. “Log 10 (q)” is the multi-test adjusted p-value in log base 10.


## 4. Discussion

PCa is an epithelial malignant tumor occurring in the prostate and is the most common malignant tumor of male genitourinary system [31]. PCa is a very slow-progressing cancer. In the early stages of the disease, many patients do not know they have it. Once the cancer begins to grow rapidly or spread outside the prostate, it becomes more serious [32, 33]. PCa remains one of the major health challenges due to lacking reliable prognostic biomarkers and therapeutic targets [34]. In this paper, 270 DEGs were identified between PCa and NP groups by differential gene expression analysis. Next, we performed Metascape gene annotation, pathway and process enrichment analysis, and PPI enrichment analysis on all 270 DEGs. Then, we identified and screened out 30 PCaSGs based on the random forest analysis and constructed an ANN model based the gene score matrix consisting of 30 PCaSGs. Lastly, we successfully validated our ANN model by testing set.One of the important findings of this study was to identify the important functions, key pathways, and protein interactions of DEGs in PCa, among which inflammation response was more closely related to PCa. Many studies had shown that the occurrence and development of tumors were closely related to the microenvironment of chronic inflammation. It had been reported that “benign” prostatic hypertrophy was not benign, but might be a chronic inflammation of the prostate's lower reproductive tract, and this chronic inflammation could be a common precursor of PCa [35]. By comparing the seropositivity of PCa patients to trichomonas vaginalis with that of the normal control population, Kim J et al. [36] found that the seropositivity of the former (19.7%) was significantly higher than that of the latter (1.7%, P < 0.001). Kwon OJ et al. [37] constructed a mouse model of prostatitis and found that inflammation alters the tissue microenvironment of the normal prostatic epithelial differentiation process and, through this cellular process, accelerates the development of PCa originating from basal cells. There had been a lot of evidence that inflammatory response plays a key role in PCa development, so we speculated that the biological functions and pathways of these DEGs may be closely related to the risk of PCa.Moreover, the ANN model is a powerful tool for disease prediction, which has higher accuracy and reliability than logistic regression, Cox proportional risk models, and decision trees [38–40]. So far, there is no study report on predicting PCa risk based on neural network model. However, in other areas of tumour research, many studies have reported using ANN models to predict cancer risk. Cegla P et al. [41] used ANN model to evaluate the influence of semiquantitative PET derived parameters and hematological parameters on the overall survival of patients with head and neck squamous cell carcinoma (HNSCC), and the results showed that ANN can be used as a supplement to PET derived parameters, which was helpful to find the prognostic parameters of HNSCC overall survival. Guo W et al. [42] collected 80 patients with advanced lung cancer who needed palliative chemotherapy, established multiple prognostic prediction models by screening clinical variables, and verified the model by ROC curve. The results showed that ANN model had high accuracy in predicting pneumonia infection during chemotherapy in lung cancer patients. Similarly, potential CT-benefit ANN model constructed by Lu J et al. [43] could accurately predict the potential benefit and long-term prognosis of adjuvant chemotherapy in patients with advanced gastric cancer and showed good prognostic stratification ability.Consistent with this finding, through independent dataset and self-verification of samples, our ANN model constructed had strong prediction ability and identification accuracy of PCa (see Figure 4(b), 5 and Figures 3(c) and 4). However, the performance of the ANN model still needs to be verified by comparison with other reliable computer-based diagnostic models, and the application value of the ANN model should be comprehensively evaluated in combination with clinical imaging and pathological biopsy.


## 5. Conclusion

In summary, our results suggested that the ANN model based on PCaSGs can be used for effectively predicting the patients with PCa and will be helpful for clinicians in guiding early diagnosis and treatment of PCa patients.



---

*Source: 1562511-2022-04-07.xml*