Identification of neuronal synapse-related signatures and potential therapeutic drugs in colorectal cancer based on machine learning algorithms and molecular docking

Introduction

Colorectal cancer (CRC) is characterized by significant heterogeneity and aggressiveness (1). It ranks third in incidence but second in mortality globally (2). Currently, decisions regarding adjuvant chemotherapy (ACT) are predominantly based on clinical pathological staging, with little consideration of molecular biological characteristics (3). Traditional clinical staging is based on the tumor-node-metastasis (TNM) staging system, which focuses on tumor invasion (T), lymph node involvement (N), and the presence of metastasis (M) (4). This incomplete approach fails to account for individual variability in tumor biology and immune characteristics, which may lead to potential overtreatment or undertreatment, ultimately hindering the ability to provide patients with optimal clinical care. Moreover, despite significant advances in tumor immunotherapy, particularly in targeting immune checkpoints, challenges such as treatment resistance, tumor spatiotemporal heterogeneity, and immune-related adverse events persist (5-7). Therefore, in the era of personalized medicine, there is an urgent need for new and reliable biomarkers to predict the prognosis of CRC patients, as well as for developing novel therapeutic strategies and drugs.

Growing evidence highlights the crucial biological and pathological roles of interactions between the nervous system and cancer (8,9). Both the central nervous system (CNS) and the peripheral nervous system (PNS) are implicated in various cancer types throughout the body (10). The field of cancer neuroscience is rapidly developing, focusing on deciphering the key signaling pathways in cancer-nervous system communication and exploring these pathways as potential targets for enhancing cancer therapies (11). It has been reported that nerve presence in breast cancer is associated with a poorer prognosis. Additionally, sensory nerves have been identified as more abundant in triple-negative breast cancer (TNBC) tumors (8). Breast cancer cells form ‘pseudo-tripartite synapses’ with neurons that secrete glutamate as a neurotransmitter, allowing the released glutamate to activate N-methyl-d-aspartate receptors (NMDARs), thereby promoting brain metastasis in breast cancer (12). Furthermore, direct, functional bona fide chemical synapses have been identified between presynaptic neurons and postsynaptic glioma cells (13). In fact, the nervous system can directly influence immune cells within the tumor microenvironment (TME), thereby significantly impacting cancer progression (14). It is well known that the gastrointestinal tract has its own intrinsic nervous system—the enteric nervous system (15). A previous study found that enteric serotonergic neurons release 5-hydroxytryptamine, which binds to its receptors and activates the Wnt/β-catenin signaling pathway, thereby driving the self-renewal of CRC stem cells (16). The above suggests that developing biomarkers and therapies targeting the nervous system in CRC tumors is a promising new direction for further exploration.

In this study, we conducted a comprehensive analysis of transcriptomic data from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO). By applying a range of machine learning techniques, we identified seven prognostic neuronal synapse genes (NSGs) and developed a highly effective neuronal synapse-related signature (NSRS) to predict the prognosis of CRC patients. Additionally, molecular docking was utilized to screen for candidate drugs targeting the identified NSGs. This work contributes to the advancement of precision medicine and aims to further improve the clinical outcomes for CRC patients. We present this article in accordance with the TRIPOD reporting checklist (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1988/rc).

Methods

Data acquisition and preprocessing

RNA sequencing (RNA-seq) data and corresponding clinical data for colon adenocarcinoma (COAD) from the TCGA-COAD (n=437) were downloaded from the University of California Santa Cruz (UCSC) Xena database (http://xena.ucsc.edu/), including 471 tumor samples and 41 adjacent normal samples. The corresponding somatic mutation data for TCGA-COAD were downloaded from the TCGA database via the Genomic Data Commons Data Portal (https://portal.gdc.cancer.gov/). In the TCGA-COAD dataset, we excluded samples that lacked clinical follow-up information. A total of 5,423 differentially expressed genes were identified between tumor and normal samples using the DESeq2 R package (17), including 2,863 upregulated and 2560 downregulated genes (with thresholds of P<0.05 and |log2 fold change| >1). The transcriptomic expression profiles and corresponding clinical data of the CRC cohorts GSE39582 (n=562) and GSE161158 (n=250) were downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/) for validation. In the GEO datasets, we excluded samples that lacked clinical follow-up information. The clinical information for all datasets can be found in Table S1. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Identification of prognostic NSGs

To identify neuron synapse-related prognostic genes, we first performed weighted gene co-expression network analysis (WGCNA) (18) on the TCGA-COAD dataset to search for co-expression modules associated with tumor events. We constructed a similarity matrix based on expression data and calculated the Pearson correlation coefficients between gene pairs. We further converted the similarity matrix into an adjacency matrix using a signed network, with the soft threshold set to b =12. Next, the topological overlap matrix (TOM) was used to convert it into a topological matrix to describe the strength of relationships between genes. A 1−TOM distance was applied to cluster the genes, and dynamic tree cutting was used to identify the modules. Each network module contained at least 50 genes, in accordance with the mixed dynamic tree-cutting criteria. In total, 13 co-expression modules were identified. The greenyellow module, which contains 550 tumor-related genes, showed the strongest positive correlation with tumor events. Subsequently, we collected 3027 neuron synapse-encoding genes from the Molecular Signatures Database (MSigDB) (https://www.gsea-msigdb.org/gsea/msigdb) and a previous study (19), which can be found in table available at https://cdn.amegroups.cn/static/public/tcr-24-1988-1.pdf. We intersected the neuron synapse genes, the greenyellow module genes, and differentially expressed genes (DEGs), identifying a total of 31 prognostic NSGs. The chromosomal locations of the 31 NSGs were visualized using the RCircos (v1.2.2) R package (20). Additionally, the protein interaction network of these NSGs was analyzed using the STRING database (https://string-db.org/). To elucidate their unique molecular characteristics and biological functions, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed using the clusterProfiler (v4.8.3) R package (21).

Construction and validation of the NSRS

To construct the NSRS prognostic model, we collected three bulk RNA-seq datasets: TCGA-COAD (training cohort), GSE39582 (testing cohort), and GSE161158 (testing cohort). First, we intersected the neuron synapse genes, the greenyellow module genes, and DEGs, identifying 31 overlapping NSGs. Then, we further utilized least absolute shrinkage and selection operator Cox regression (LASSO-Cox) algorithm to screen seven prognostic NSGs from the 31 candidate NSGs. The NSRS was constructed based on normalized gene expression levels weighted by corresponding coefficients from the multivariate stepwise Cox regression. The NSRS was calculated by the following formula:

riskScore=∑i=1nki∗Ei[1]

In this formula, E_i represents the normalized gene expression level of the NSGs, and K_i represents the corresponding coefficient. Patients in the training and testing cohort were divided into high- and low-risk groups according to the median cut-off of the risk score. We further assessed the predictive effectiveness of NSRS for clinical outcomes using Kaplan-Meier survival curves and log-rank tests on the training dataset TCGA-COAD, as well as the testing datasets GSE39582 and GSE161158. Additionally, we performed univariate and multivariate Cox regression analyses to determine whether NSRS could serve as an independent prognostic indicator for CRC, accounting for variables such as age, gender, grade, and TNM clinical stage. The “timeROC” R package (v0.4) was utilized to generate receiver operating characteristic (ROC) curves and the predictive ability of the NSRS prognostic model was evaluated by the area under the curve (AUC) value of the ROC curve.

Verification of hub prognostic NSGs

We employed the Kaplan-Meier Plotter database (https://kmplot.com/analysis/) to investigate the relationship between the expression of prognostic NSGs and the prognosis of CRC patients. Furthermore, we analyzed the association between the expression of NSGs and CRC prognosis in the TCGA-COAD dataset by categorizing patients into high and low expression groups based on the median expression levels of each gene. The protein levels of prognostic NSGs were analyzed using the cProSite database (https://cprosite.ccr.cancer.gov/). To validate the expression of the seven hub prognostic NSGs in CRC and normal tissues, immunohistochemical staining of CRC and normal colorectal tissues were obtained from The Human Protein Atlas database (https://www.proteinatlas.org/).

Somatic mutation analysis and gene set enrichment analysis (GSEA)

To explore the relevant genetic variations, we utilized the Maftools (v2.16.0) R package (22) to analyze gene mutations in the two NSRS subgroups. We performed GSEA functional enrichment analysis using GSEA software (v4.3.3) to identify key biological pathways and potential molecular mechanisms between the high and low NSRS subgroups (23). The reference molecular signature database used for this analysis was ‘h.all.v2024.1.Hs.symbols.gmt’.

The landscape of the TME and the potential significance of NSRS in patient immunotherapy

To identify the TME characteristics of CRC samples, the xCell (v1.1.0) R package (24) was used to evaluate the abundance of 54 cell types, including lymphocytes and non-lymphocytes, in the GSE39582 dataset. Table S2 provides a list of the cell types included in the xCell analysis. The relative proportions of these cell types were then compared between the two NSRS subgroups. The relationship between the NSRS score and tumor immune microenvironment (TIME) characteristics (stromal score, immune score, estimate score, and tumor purity) was analyzed using the ESTIMATE algorithm (25) and the Wilcoxon rank-sum test.

Target gene drug prediction and molecular docking

In this study, the Drug Signatures Database (DSigDB) was used to screen for candidate drugs targeting the identified NSGs. The DSigDB database can be accessed through the Enrichr platform (https://maayanlab.cloud/Enrichr/). To gain deeper insights into the impact of candidate drugs on their target genes and assess the druggability of these targets, molecular docking was performed at the atomic level in this study. We evaluated the binding affinities and interaction patterns between the candidate drugs and their targets. In this study, the protein-ligand docking software AutoDockTools (v1.5.7) (https://autodock.scripps.edu/) and AutoDock Vina (v4.2.6) (https://vina.scripps.edu/) were used to perform molecular docking between the identified significant and feasible drugs and the proteins encoded by their corresponding target genes (26). The structural data of the small molecule drugs were sourced from the PubChem Database (https://pubchem.ncbi.nlm.nih.gov/). The 3D structural data of the target proteins were obtained from the Protein Data Bank (PDB, https://www.rcsb.org/). Initially, all water molecules were removed from the protein and ligand files, followed by the addition of polar hydrogen atoms. The grid boxes were centered around the structural domain of each protein. Finally, the entire molecular docking process was visualized using AutoDockTools (v1.5.7) for analysis.

Statistical analysis

The correlation between two continuous variables was evaluated using the Pearson correlation coefficient. Categorical variables were compared using the Chi-squared test, while continuous variables were compared using the Wilcoxon rank-sum test or the t-test. The optimal cut-off value was identified using the survminer package (v0.4.9). Cox regression and Kaplan-Meier analyses were conducted using the survival package (v3.5.5). In this study, all statistical analyses were conducted using R (v4.3.1) and GraphPad Prism (v9.5.0). The statistical significance was set at P<0.05.

Results

Identification of prognostic NSGs in CRC patients

Initially, transcriptomic and clinical data of CRC patients were retrieved from the UCSC Xena database. Subsequently, differential expression analysis was performed on the TCGA-COAD dataset, which included 471 tumor samples and 41 adjacent normal samples, identifying 2863 upregulated genes and 2560 downregulated genes with P<0.05 and |log₂ fold change| >1 (Figure 1A).

Figure 1 Identification of prognostic NSGs in CRC patients. (A) Volcano plot exhibiting 5,423 differentially expressed genes between tumor and normal samples. (B) The selection of soft-thresholding powers. (C) Gene dendrogram for identifying gene co-expression modules. (D) Correlations of 13 modules with clinical phenotypes. The color scale represents the range of correlation coefficients. Red indicates positive correlation, while blue indicates negative correlation, with deeper colors signifying stronger correlations. Each cell contains the correlation coefficient, with the P-value in parentheses, which assesses the statistical significance of the correlation. (E) Hierarchical clustering and correlation heatmap of modules with tumor trait. The scale bar on the right indicates correlation coefficient values. (F) Correlation scatter plot between the green and yellow module and tumor phenotype. (G) Venn diagram to identify 31 overlapping prognostic NSGs. (H) Chromosomal localization of the 31 prognostic NSGs. (I) Protein-protein interactions analysis of the 31 prognostic NSGs. (J,K) The bubble plots showing the results of GO enrichment analysis and KEGG enrichment analysis for the 31 prognostic NSGs. BP, biological process; CC, cellular component; CRC, colorectal cancer; DEGs, differentially expressed genes; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function; N, node; NSGs, neuronal synapse genes; OS, overall survival; T, tumor; WGCNA, weighted gene co-expression network analysis.

We applied WGCNA to construct a gene co-expression network using the TCGA-COAD dataset to identify gene modules associated with clinical traits. The relationships between different soft thresholds and both scale independence and mean connectivity are shown in Figure 1B. The soft thresholding power was set to 12 to ensure the network was scale-free. According to the mixed dynamic shear tree criteria, each network module contained a minimum of 50 genes. Gene co-expression modules were identified using the dynamic tree cut package, and similar modules were merged, resulting in a final total of 13 co-expression modules (Figure 1C). Pearson’s test was used to analyze the correlation between each gene module and clinical traits. The greenyellow module, which contains 550 tumor-related genes, showed the strongest positive correlation with tumor events (Figure 1D). The heatmap of topological overlap further highlighted the strongest correlation between the greenyellow module and tumor events (Figure 1E). The correlation between the greenyellow module genes and clinical traits was analyzed and visualized in Figure 1F (cor =0.69, P=5.8e−79).

To identify prognostic NSGs in CRC patients, we collected 3,027 neuron synapse-encoding genes from MSigDB database and a previous study (table available at https://cdn.amegroups.cn/static/public/tcr-24-1988-1.pdf) (19). The Venn diagram illustrated the overlap between neuron synapse genes, the greenyellow module, and DEGs (Figure 1G). This analysis identified 31 overlapping prognostic NSGs. We further visualized the chromosomal locations and protein-protein interactions of the 31 prognostic NSGs (Figure 1H,1I). Finally, to investigate the potential molecular functions of NSGs, we conducted GO and KEGG pathway enrichment analyses. GO analysis revealed that these NSGs were significantly enriched in the ribose phosphate biosynthetic process, postsynaptic density, and transferase activity, transferring one-carbon groups (Figure 1J). KEGG analysis revealed that these NSGs were significantly enriched in the purine metabolism and biosynthesis of cofactors pathways (Figure 1K), suggesting that these prognostic NSGs may influence CRC tumor progression primarily by affecting enteric neuronal secretion and metabolic activities.

Construction of the NSRS prognostic risk model

We applied LASSO-Cox regression with 10-fold cross-validation to further identify 9 key prognostic NSGs from the 31 candidate NSGs (Figure 2A). Finally, seven core prognostic NSGs (ATIC, TRAP1, AGAP3, NME1, RUVBL1, WDR77 and LRFN4) were selected and used to construct the NSRS prognostic model through multivariate stepwise Cox regression analysis, using the lowest Akaike information criterion (AIC) for model selection. The hazard ratios (HRs) of these seven prognostic NSGs are shown in Figure 2B, indicating that AGAP3, RUVBL1, and LRFN4 may be significant risk factors for CRC, while TRAP1 and NME1 may be significant protective factors for CRC. The detailed steps for constructing the NSRS prognostic model could be found in the Methods section. Notably, Kaplan-Meier survival analysis revealed in both the training cohort (TCGA-COAD) and the validation cohorts (GSE39582 and GSE161158), the high-risk group consistently exhibited a poorer prognosis compared to the low-risk group (Figure 2C-2E). Moreover, ROC curve analysis revealed that NSRS had a great predictive accuracy, with AUC values of 0.716, 0.714, and 0.691 at 1, 3, and 5 years (Figure 2C). Additionally, we analyzed and visualized the relationship between RiskScore, patient follow-up time, events, and the expression changes of the seven NSGs. We observed that as the RiskScore increased, the survival rate of CRC patients significantly decreased (Figure 2C-2E).

Figure 2 Construction and validation of neuronal synapse-related signature in CRC patients. (A) LASSO Cox regression analysis for 31 prognostic NSGs. The tuning parameter (λ) was calculated based on the partial likelihood deviance with 10-fold cross-validation. (B) Forest plot displaying 7 prognostic NSGs identified through multivariate stepwise Cox regression analysis. (C) Kaplan-Meier curve of NSRS subgroups in the TCGA-COAD training cohort. Log-rank P values are indicated. The shaded areas around each curve indicate the 95% confidence intervals for the survival estimates. ROC curves analysis of NSRS at 1, 3, and 5 years. Distribution of risk score and survival status of NSRS in TCGA-COAD training cohort. (D,E) Kaplan-Meier curve of NSRS subgroups in the GSE161158 and GSE39582 testing cohort. Log-rank P values are indicated. The shaded areas around each curve indicate the 95% confidence intervals for the survival estimates. ROC curves analysis of NSRS at 1, 3, and 5 years. Distribution of risk score and survival status of NSRS in GSE161158 and GSE39582 testing cohort. AUC, area under the curve; CI, confidence interval; COAD, colon adenocarcinoma; CRC, colorectal cancer; LASSO, least absolute shrinkage and selection operator; NSGs, neuronal synapse genes; NSRS, neuronal synapse-related signature; ROC, receiver operating characteristic; TCGA, the Cancer Genome Atlas.

Verification of hub prognostic NSGs in CRC

ROC curve analysis revealed that the AUC for all seven prognostic NSGs was greater than 0.68, indicating their diagnostic efficiency (Figure 3A). We also investigated the relationship between the expression of prognostic NSGs and the prognosis of CRC patients. In the Kaplan-Meier Plotter database, we observed that CRC patients with high expression levels of AGAP3 and LRFN4 exhibited poorer survival outcomes, whereas those with high expression levels of ATIC, TRAP1, NME1, and RUVBL1 had better survival outcomes. In contrast, the WDR77 expression was not significantly associated with the prognosis of CRC patients (Figure 3B-3H). In the TCGA-COAD dataset, we demonstrated that patients with high expression levels of AGAP3, LRFN4, and RUVBL1 exhibited poor prognosis, whereas those with high expression levels of TRAP1 and NME1 had favorable prognosis. In contrast, expression levels of ATIC and WDR77 were not significantly associated with CRC prognosis (Figure S1A-S1G). The prognostic significance of ATIC and RUVBL1 expression levels showed significant differences between the two analysis methods, potentially due to substantial variations in patient cohort sizes. Next, we analyzed the mRNA expression levels of these seven prognostic NSGs and found that their expression was consistently higher in tumor samples compared to normal samples (Figure 3I), suggesting that these genes might play a potential role in the development of CRC. Additionally, we analyzed the protein levels of the seven prognostic NSGs in CRC patients using the cProSite database. We found that ATIC, TRAP1, NME1, and RUVBL1 were significantly upregulated in tumor samples from CRC patients, while AGAP3 and WDR77 showed no significant difference (Figure 3J). However, no corresponding protein data for LRFN4 was available in the cProSite database. This indicated that the protein and mRNA levels of these NSGs were not entirely consistent, potentially due to mRNA alternative splicing and post-translational modifications of the proteins. Immunohistochemistry staining further revealed the protein levels of AGAP3, ATIC, TRAP1, NME1, LRFN4, RUVBL1 and WDR77 in both CRC and normal tissues (Figure 3K).

Figure 3 Validation of the mRNA and protein levels of the 7 prognostic NSGs. (A) ROC analysis of 7 prognostic NSGs. (B-H) Prognosis analysis between high and low expression levels of NSGs in colorectal cancer patients. (I) mRNA expression levels of prognostic NSGs in CRC cancer. (J) Protein levels of prognostic NSGs in CRC cancer. (K) Representative immunohistochemistry staining results of AGAP3 (HPA012883, https://www.proteinatlas.org/ENSG00000133612-AGAP3/cancer/colorectal+cancer#img and https://www.proteinatlas.org/ENSG00000133612-AGAP3/tissue/colon#img), ATIC (CAB013462, https://www.proteinatlas.org/ENSG00000138363-ATIC/cancer/colorectal+cancer#img and https://www.proteinatlas.org/ENSG00000138363-ATIC/tissue/colon#img), TRAP1 (HPA041082, https://www.proteinatlas.org/ENSG00000126602-TRAP1/cancer/colorectal+cancer#img and https://www.proteinatlas.org/ENSG00000126602-TRAP1/tissue/colon#img), NME1 (HPA008467, https://www.proteinatlas.org/ENSG00000239672-NME1/cancer/colorectal+cancer#img and https://www.proteinatlas.org/ENSG00000239672-NME1/tissue/colon#img), LRFN4 (HPA077570, https://www.proteinatlas.org/ENSG00000173621-LRFN4/cancer/colorectal+cancer#img and https://www.proteinatlas.org/ENSG00000173621-LRFN4/tissue/colon#img), RUVBL1 (HPA019947, https://www.proteinatlas.org/ENSG00000175792-RUVBL1/cancer/colorectal+cancer#img and https://www.proteinatlas.org/ENSG00000175792-RUVBL1/tissue/colon#img) and WDR77 (HPA027271, https://www.proteinatlas.org/ENSG00000116455-WDR77/cancer/colorectal+cancer#img and https://www.proteinatlas.org/ENSG00000116455-WDR77/tissue/colon#img) in normal colon tissue and colorectal cancer tissue, obtained from the Human Protein Atlas database, with permission from its website. Scale bar: 100 μm. **, P<0.01; ***, P<0.001; ns, not significant. AUC, area under the curve; CI, confidence interval; CRC, colorectal cancer; HR, hazard ratio; NSGs, neuronal synapse genes; ROC, receiver operating characteristic; TMT, tandem mass tag; TPM, transcripts per million.

Clinical significance of the NSRS model

To evaluate the clinical significance and predictive accuracy of the NSRS model, we analyzed risk scores in combination with other clinical variables and constructed a new nomogram (Figure 4A). The C-index of the nomogram model was 0.804. We observed significant differences in riskScore among CRC patients at different tumor stages (Figure 4B), suggesting that the NSRS model offers robust diagnostic value in assessing tumor progression. The calibration curves demonstrated excellent concordance between the predicted and observed survival probabilities at the 1-, 3-, and 5-year intervals (Figure 4C-4E). Univariate and multivariate Cox regression analyses demonstrated that the riskScore was an independent and effective prognostic factor (univariate Cox model: HR =1.452, P<0.001; multivariate Cox model: HR =1.375, P<0.001) (Figure 4F,4G). The above results further confirm the reliability and clinical applicability of the NSRS predictive model that we constructed.

Figure 4 Clinical value of NSRS model by independent prognostic analysis. (A) Nomogram was established in TCGA-COAD cohort with risk score and clinical factors. (B) The box plot showing the risk score levels across tumor stages. (C-E) The calibration plots of the nomogram illustrate the concordance between predicted and observed survival rates at 1, 3, and 5 years. The error bars represent the 95% confidence intervals. (F) Univariate and (G) multivariate cox regression analysis of clinicopathologic factors and risk score. CI, confidence interval; COAD, colon adenocarcinoma; N, node; NSRS, neuronal synapse-related signature; OS, overall survival; T, tumor; TCGA, The Cancer Genome Atlas.

Exploration of the molecular functions of NSRS subgroups

To investigate the potential biological functions of different NSRS subgroups, we employed the “DESeq2” R package to perform differential expression analysis on the TCGA-COAD dataset, identifying transcriptomic variations between NSRS subgroups. Next, we conducted GO enrichment analysis and GSEA to identify gene sets enriched in distinct NSRS subgroups. GO analysis revealed that upregulated genes in the high NSRS group were primarily enriched in biological processes associated with peripheral sympathetic nerve activity, including “temperature homeostasis”, “negative regulation of secretion”, “cellular response to salt”, “adaptive thermogenesis”, “regulation of amine transport”, and “regulation of glutamate secretion” (Figure 5A). Conversely, downregulated genes in the high NSRS group were predominantly enriched in immune-related functions, such as “human immune response”, “antimicrobial humoral response”, and “cell killing” (Figure 5B). This suggests that peripheral neuronal signaling in the TME may influence immune cell function, thereby promoting immune evasion and cancer progression (27).

Figure 5 Function enrichment analysis of different NSRS subgroups. (A,B) Bubble plot illustrating Gene Ontology enrichment analysis of the up-regulated and down-regulated genes. (C-H) GSEA analysis of Hallmark Gene Sets results showing the enriched pathways in the high-NSRS group. (I-N) GSEA analysis of Hallmark Gene Sets results showing the enriched pathways in the low-NSRS group. GSEA, gene set enrichment analysis; NSRS, neuronal synapse-related signature.

GSEA enrichment analysis revealed that the “epithelial mesenchymal transition”, “inflammatory response”, “tnfa signaling via nfkb”, “hypoxia” and “kras signaling” pathways were significantly enriched in the high NSRS group (Figure 5C-5H), while the “E2F targets”, “MYC targets v2”, “oxidative phosphorylation” and “fatty acid metabolism” pathways were primarily enriched in the low NSRS group (Figure 5I-5N). This suggests that the high NSRS group exhibits enhanced epithelial-mesenchymal transition (EMT), hypoxia, and inflammatory features, all of which are known to contribute to tumor evasion and progression (28-30). In contrast, the low NSRS group exhibits higher cell cycle and metabolic activity.

Somatic mutation characteristics of NSRS subgroups

We further investigated the relationship between somatic mutations, NSRS, and survival outcomes in CRC patients, as well as the differences in somatic mutation profiles across the different NSRS subgroups. First, through Kaplan-Meier survival curve analysis, we first observed that CRC patients with both high tumor mutation burden (TMB) and high NSRS had significantly poorer survival outcomes compared to those with low TMB and low NSRS (Figure 6A). Additionally, we analyzed the distribution of somatic mutations between the low and high NSRS subgroups in CRC to uncover their genetic signatures. In both the high and low NSRS groups, APC, TP53, KRAS, and TTN genes exhibited the highest mutation rates (Figure 6B,6C), with mutations in these genes closely linked to tumorigenesis. Furthermore, we observed differences in somatic mutation frequencies between the high and low NSRS groups (Figure 6B,6C).

Figure 6 Tumor mutation characteristics in different NSRS subgroups. (A) The prognosis of the combination of NSRS and TMB. (B,C) The difference of most frequent mutation genes in low NSRS and high NSRS subgroups. NSRS, neuronal synapse-related score; TMB, tumor mutation burden.

The landscape of the TIME and the immune characteristics of NSRS subgroups

In order to comprehensively investigate the TIME characteristics in high and low- NSRS samples, the xCell algorithm (24) was utilized to evaluate 54 cell types, including 21 types of lymphocytes and 33 other cell types, in GSE39582 dataset (table available at https://cdn.amegroups.cn/static/public/tcr-24-1988-2.xlsx). We observed that most immune cells and stromal cells were significantly higher in the high NSRS group, particularly dendritic cells, macrophages, pericytes, and fibroblasts, indicating a more complex TIME. Conversely, cell types such as Th1 cells, Th2 cells, and megakaryocyte-erythroid progenitors (MEP) were significantly more abundant in the low NSRS group (Figure 7A). Notably, the high NSRS group exhibited a higher abundance of neuronal cells (Figure 7B), consistent with our classification results, suggesting that increased neural system enrichment may be associated with greater malignancy in CRC tumors (31). We analyzed the relationship between the seven prognostic NSGs and 54 cell types, revealing that all seven NSGs were significantly correlated with various immune and stromal cells (Figure 7C). Additionally, we used the ESTIMATE algorithm to investigate the relationship between NSRS and tumor purity in CRC patients. Our analysis revealed that the high NSRS group had higher ESTIMATE score, immune score and stromal score compared to the low NSRS group, while tumor purity was lower, suggesting a more complex tumor microenvironment in the high NSRS group (Figure 7D-7G). It has been reported that tumors with lower purity have more somatic mutations and are associated with poorer prognosis (32). We further analyzed and visualized the correlation between riskScore and immune score (P<0.0001) (Figure 7H). Finally, we investigated the differences in immune checkpoints, major histocompatibility complex (MHC), and cytokines between the different NSRS subgroups. We found that the high NSRS group exhibited higher levels of immune checkpoints, MHC, and cytokines, indicating a more pronounced immune phenotype (Figure 7I).

Figure 7 The landscape of the TME in CRC and the immune characteristics of different NSRS subgroups. (A) Boxplots illustrating the enrichment levels of the 29 immune signatures in the high-score and low-score groups. (B) The boxplot illustrating the enrichment abundance of neurons in two groups. (C) Visualization of immune signatures and prognostic genes correlation results. (D-G) The estimate score, immune score, stromal score and tumor purity in different NSRS groups. (H) Scatter plot showing the correlation between the NSRS risk score and the immune score. (I) Specific immune checkpoints, MHC and chemokines of the NSRS subgroups. *, P<0.05; **, P<0.01; ***, P<0.001; ns, not significant. CRC, colorectal cancer; MHC, major histocompatibility complex; NSRS, neuronal synapse-related signature; TME, tumor microenvironment.

Candidate drug screening and molecular docking

We used the DSigDB database to identify potentially effective intervention drugs targeting the prognostic NSGs. We listed the potential small-molecule drugs with adjusted P values less than 0.05 (Table 1). The results identified alprostadil as the most significant drug, associated with ATIC, RUVBL1 and WDR77. Nocodazole, a well-known anti-tumor agent, was found to significantly target RUVBL1 and NME1 (Table 1). Additionally, RUVBL1 was identified as a target for all candidate drugs, highlighting its potential as a critical therapeutic target for CRC.

To further evaluate the binding affinity of candidate drugs to their targets, we utilized AutoDockTools (v1.5.7) to conduct molecular docking between selected feasible candidate drugs and their corresponding gene-encoded proteins, obtaining the respective binding energies of the interactions. NME1 and dihydroergocristine showed the lowest binding energy (−8.9 kcal/mol), indicating a highly stable interaction (Table 2). Figure 8 visualized the hydrogen bonds and key amino acid residues involved in the binding between the candidate drugs and their target proteins.

Figure 8 The docking conformation analysis of prognostic NSGs with their docked compounds. (A) ATIC docking alprostadil. (B) RUVBL1 docking alprostadil. (C) TRAP1 docking clindamycin. (D) NME1 docking clindamycin. (E) NME1 docking dihydroergocristine. (F) NME1 docking nocodazole. NSGs, neuronal synapse genes.

Discussion

The American Joint Committee on Cancer (AJCC) staging system is a traditional tool employed in clinical management to guide treatment decisions and develop surveillance strategies. Nevertheless, its utility is restricted by the variability in clinical outcomes observed among patients at the same stage (33-35). Hence, there is an urgent need for more advanced personalized assessment methods to enhance clinical decision-making and customize treatment strategies for patients. With advancements in artificial intelligence and bioinformatics, several prognostic biomarkers, including long noncoding RNAs (lncRNAs), microsatellite instability (MSI), TMB, and immune activation, have been developed for CRC patients (36-39). However, these biomarkers are still constrained by spatiotemporal heterogeneity, moderate accuracy, and limited personalization. To our knowledge, this is the first study to employ machine learning for developing neuron synapse-related prognostic markers from the perspective of neuron-tumor interactions. Additionally, we predicted candidate drugs targeting prognostic NSGs through molecular docking analysis.

Based on the seven prognostic NSGs identified using WGCNA and LASSO-Cox regression algorithms, we constructed an NSRS risk model for CRC patients, providing new insights into precise patient stratification and personalized treatment. The two NSRS risk levels demonstrated significantly different prognostic outcomes. Moreover, ROC and C-index analyses revealed that NSRS consistently exhibited high accuracy and stable performance across a training cohort and two testing cohorts, highlighting its significant potential for clinical applications. Univariate and multivariate Cox regression analyses demonstrated that NSRS is an independent prognostic factor, with higher NSRS scores observed in patients with more advanced tumor stages. GO analysis demonstrated that the high NSRS group was significantly enriched in functions related to neural activity, including neurotransmitter secretion, while the low NSRS group was enriched in functions associated with immune response. This finding suggests that the abundance of the nervous system in the TME may be linked to immunosuppression. Prior evidence indicates that the CNS plays a role in regulating cancer progression and modulating immune system activity (40-43). Xu et al. demonstrated that the hypothalamic-pituitary unit promotes tumor immunosuppression via α-melanocyte-stimulating hormone (α-MSH) and its receptor, MC5R (44). Overall, targeting the tumor nervous system may represent a promising therapeutic strategy for the treatment of CRC. Notably, GSEA analysis revealed significant enrichment of pathways related to EMT, inflammation, and hypoxia in the high NSRS group, while pathways associated with oxidative phosphorylation, cell cycle, and differentiation were enriched in the low NSRS group. These findings suggest that patients in the high NSRS group exhibit more aggressive tumor characteristics. EMT is involved in various aspects of malignancy, including metastasis formation, tumor invasion, and treatment resistance (45-48). Hypoxia or inadequate tissue oxygenation is a key factor in the difficulty of curing most cancers and also contributes to increased resistance of cancer cells to treatment (49,50). Chronic inflammation plays a critical role in tumor progression and therapy resistance (51). Based on these findings, targeting the EMT, inflammation, and hypoxia pathways may offer a more effective treatment strategy for high NSRS CRC patients, while targeting the cell cycle and differentiation pathways may be more beneficial for low NSRS CRC patients.

Additionally, to further investigate the relationship between NSRS and tumor immune characteristics in CRC, the xCell algorithm was used to evaluate 54 cell types in both high and low NSRS groups. We observed that most immune cells and stromal cells, including dendritic cells, macrophages, mast cells, and fibroblasts, were significantly higher in the high NSRS group. Specifically, macrophages and mast cells produce factors that sustain chronic inflammation and promote tumor growth (52-54). The above indicates more intricate TIME in the high NSRS group. It is possible that more malignant tumors, due to their rapid proliferation, the formation of an inflammatory microenvironment, and the activation of immune evasion mechanisms, may exhibit increased immune infiltration and responses. Immune checkpoints are a class of immunosuppressive molecules that play a key role in preventing autoimmune reactions (55-57). We found that the high NSRS group expressed higher levels of immune checkpoints, indicating increased sensitivity to immune checkpoint inhibitors (ICIs). Previous studies have defined infiltrated-inflamed TIMEs as ‘hot’ tumors, characterized by high levels of immune cell infiltration and a stronger response to ICIs (58,59), though this is not always associated with a favorable prognosis (60). Additionally, we found that the high NSRS group had higher immune and stromal scores, while tumor purity was lower, which is consistent with our xCell analysis results. This indicates that low-purity CRC tumors may be indicative of a poor prognosis (61). Previous study has shown that low-purity gliomas exhibit significantly heightened local immune responses (32), which are unfavorable for glioma prognosis (62).

Finally, using the DSigDB database and molecular docking, we identified several candidate drugs targeting prognostic NSGs, such as alprostadil, dihydroergocristine, and nocodazole. Among these, alprostadil, which targets ATIC and RUVBL1, is a synthetic analogue of prostaglandin E1 (PGE1) and is widely used for the treatment of male erectile dysfunction (63). Li et al. reported that PGE1 and its analog misoprostol inhibit the self-renewal of chronic myelogenous leukemia stem cells by activating the EP4 receptor (64). Additionally, PGE1 effectively inhibits the growth of drug-refractory human medulloblastoma (65). Dihydroergocristine, a drug approved by the US Food and Drug Administration (FDA), has been recently reported to inhibit the proliferation of CRC tumor cells (66). In addition, it has been proven to be a non-toxic and relatively safe drug (67).

The NSRS model can be easily reproduced using a polymerase chain reaction (PCR)-based assay, making it highly appealing for clinical translation and application. These results offer valuable insights into targeting the tumor nervous system for diagnosis and treatment, and they lay a solid foundation for future experimental research. However, the model also has certain limitations. Firstly, further clinical studies and practice are needed to validate our findings, given that the current conclusions are derived from a relatively limited patient cohort. The reliability of the NSRS model should be further verified in large-sample retrospective and prospective studies with long-term follow-up data to ensure its robustness and generalizability across diverse patient populations. Secondly, we need to further elucidate the specific mechanisms underlying the prognostic significance of NSGs in CRC and validate these mechanisms through in vitro and in vivo experiments. Additionally, further experimental research and clinical validation are essential for candidate drugs targeting prognostic NSGs to facilitate their translation into clinical treatment pathways for CRC patients.

Conclusions

In conclusion, this is the first study to develop neuron synapse-related prognostic biomarkers from the perspective of neuron-tumor interactions using machine learning. These findings have the potential to enable more targeted treatment for CRC patients and improve their prognosis.

Acknowledgments

We would like to thank the researchers for providing open access to the raw bulk RNA-seq data, as well as all other individuals who contributed to this study.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1988/rc

Peer Review File: Available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1988/prf

Funding: This research was supported by grants from the Young Scientists Fund of National Natural Science Foundation of China (No. 81702428).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tcr.amegroups.com/article/view/10.21037/tcr-24-1988/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Punt CJ, Koopman M, Vermeulen L. From tumour heterogeneity to advances in precision treatment of colorectal cancer. Nat Rev Clin Oncol 2017;14:235-46. [Crossref] [PubMed]
2. Bray F, Laversanne M, Sung H, et al. Global cancer statistics 2022: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2024;74:229-63. [Crossref] [PubMed]
3. Weiser MR. AJCC 8th Edition: Colorectal Cancer. Ann Surg Oncol 2018;25:1454-5.
4. Puppa G, Sonzogni A, Colombari R, et al. TNM staging system of colorectal carcinoma: a critical appraisal of challenging issues. Arch Pathol Lab Med 2010;134:837-52. [Crossref] [PubMed]
5. Vitale I, Shema E, Loi S, et al. Intratumoral heterogeneity in cancer progression and response to immunotherapy. Nat Med 2021;27:212-24. [Crossref] [PubMed]
6. Jin Z, Sinicrope FA. Mismatch Repair-Deficient Colorectal Cancer: Building on Checkpoint Blockade. J Clin Oncol 2022;40:2735-50. [Crossref] [PubMed]
7. Weng J, Li S, Zhu Z, et al. Exploring immunotherapy in colorectal cancer. J Hematol Oncol 2022;15:95. [Crossref] [PubMed]
8. Le TT, Payne SL, Buckwald MN, et al. Sensory nerves enhance triple-negative breast cancer invasion and metastasis via the axon guidance molecule PlexinB3. NPJ Breast Cancer 2022;8:116. [Crossref] [PubMed]
9. Krishna S, Choudhury A, Keough MB, et al. Glioblastoma remodelling of human neural circuits decreases survival. Nature 2023;617:599-607. [Crossref] [PubMed]
10. Hanahan D, Monje M. Cancer hallmarks intersect with neuroscience in the tumor microenvironment. Cancer Cell 2023;41:573-80. [Crossref] [PubMed]
11. Pan C, Winkler F. Insights and opportunities at the crossroads of cancer and neuroscience. Nat Cell Biol 2022;24:1454-60. [Crossref] [PubMed]
12. Zeng Q, Michael IP, Zhang P, et al. Synaptic proximity enables NMDAR signalling to promote brain metastasis. Nature 2019;573:526-31. [Crossref] [PubMed]
13. Venkataramani V, Tanev DI, Strahle C, et al. Glutamatergic synaptic input to glioma cells drives brain tumour progression. Nature 2019;573:532-8. [Crossref] [PubMed]
14. Mancusi R, Monje M. The neuroscience of cancer. Nature 2023;618:467-79. [Crossref] [PubMed]
15. Marklund U. Diversity, development and immunoregulation of enteric neurons. Nat Rev Gastroenterol Hepatol 2022;19:85-6. [Crossref] [PubMed]
16. Zhu P, Lu T, Chen Z, et al. 5-hydroxytryptamine produced by enteric serotonergic neurons initiates colorectal cancer stem cell self-renewal and tumorigenesis. Neuron 2022;110:2268-2282.e4. [Crossref] [PubMed]
17. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 2014;15:550. [Crossref] [PubMed]
18. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 2008;9:559. [Crossref] [PubMed]
19. Zhang W, Fu Y, Peng L, et al. Immunoproximity biotinylation reveals the axon initial segment proteome. Nat Commun 2023;14:8201. [Crossref] [PubMed]
20. Zhang H, Meltzer P, Davis S. RCircos: an R package for Circos 2D track plots. BMC Bioinformatics 2013;14:244. [Crossref] [PubMed]
21. Yu G, Wang LG, Han Y, et al. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 2012;16:284-7. [Crossref] [PubMed]
22. Mayakonda A, Lin DC, Assenov Y, et al. Maftools: efficient and comprehensive analysis of somatic variants in cancer. Genome Res 2018;28:1747-56. [Crossref] [PubMed]
23. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
24. Aran D, Hu Z, Butte AJ. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol 2017;18:220. [Crossref] [PubMed]
25. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
26. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function, efficient optimization, and multithreading. J Comput Chem 2010;31:455-61. [Crossref] [PubMed]
27. Khanmammadova N, Islam S, Sharma P, et al. Neuro-immune interactions and immuno-oncology. Trends Cancer 2023;9:636-49. [Crossref] [PubMed]
28. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell 2011;144:646-74. [Crossref] [PubMed]
29. Chen Z, Han F, Du Y, et al. Hypoxic microenvironment in cancer: molecular mechanisms and therapeutic interventions. Signal Transduct Target Ther 2023;8:70. [Crossref] [PubMed]
30. de Visser KE, Joyce JA. The evolving tumor microenvironment: From cancer initiation to metastatic outgrowth. Cancer Cell 2023;41:374-403. [Crossref] [PubMed]
31. Li X, Ye C, Wang M, et al. Crosstalk Between the Nervous System and Colorectal Cancer. Neurosci Bull 2025;41:93-106. [Crossref] [PubMed]
32. Zhang C, Cheng W, Ren X, et al. Tumor Purity as an Underlying Key Factor in Glioma. Clin Cancer Res 2017;23:6279-91. [Crossref] [PubMed]
33. Poston GJ, Figueras J, Giuliante F, et al. Urgent need for a new staging system in advanced colorectal cancer. J Clin Oncol 2008;26:4828-33. [Crossref] [PubMed]
34. Amin MB, Greene FL, Edge SB, et al. The Eighth Edition AJCC Cancer Staging Manual: Continuing to build a bridge from a population-based to a more "personalized" approach to cancer staging. CA Cancer J Clin 2017;67:93-9.
35. Fan Z, Edelmann D, Yuan T, et al. Developing survival prediction models in colorectal cancer using epigenome-wide DNA methylation data from whole blood. NPJ Precis Oncol 2024;8:191. [Crossref] [PubMed]
36. Diao Z, Han Y, Chen Y, et al. The clinical utility of microsatellite instability in colorectal cancer. Crit Rev Oncol Hematol 2021;157:103171. [Crossref] [PubMed]
37. Liu Z, Liu L, Weng S, et al. Machine learning-based integration develops an immune-derived lncRNA signature for improving outcomes in colorectal cancer. Nat Commun 2022;13:816. [Crossref] [PubMed]
38. Mezheyeuski A, Backman M, Mattsson J, et al. An immune score reflecting pro- and anti-tumoural balance of tumour microenvironment has major prognostic impact and predicts immunotherapy response in solid cancers. EBioMedicine 2023;88:104452. [Crossref] [PubMed]
39. Budczies J, Kazdal D, Menzel M, et al. Tumour mutational burden: clinical utility, challenges and emerging improvements. Nat Rev Clin Oncol 2024;21:725-42. [Crossref] [PubMed]
40. Gillespie S, Monje M. The Neural Regulation of Cancer. Annual Review of Cancer Biology 2020;4:371-90.
41. Schiller M, Ben-Shaanan TL, Rolls A. Neuronal regulation of immunity: why, how and where? Nat Rev Immunol 2021;21:20-36. [Crossref] [PubMed]
42. Liu SQ, Li B, Li JJ, et al. Neuroendocrine regulations in tissue-specific immunity: From mechanism to applications in tumor. Front Cell Dev Biol 2022;10:896147. [Crossref] [PubMed]
43. Winkler F, Venkatesh HS, Amit M, et al. Cancer neuroscience: State of the field, emerging directions. Cell 2023;186:1689-707. [Crossref] [PubMed]
44. Xu Y, Yan J, Tao Y, et al. Pituitary hormone α-MSH promotes tumor-induced myelopoiesis and immunosuppression. Science 2022;377:1085-91. [Crossref] [PubMed]
45. Yang J, Antin P, Berx G, et al. Guidelines and definitions for research on epithelial-mesenchymal transition. Nat Rev Mol Cell Biol 2020;21:341-52. [Crossref] [PubMed]
46. Lüönd F, Sugiyama N, Bill R, et al. Distinct contributions of partial and full EMT to breast cancer malignancy. Dev Cell 2021;56:3203-3221.e11. [Crossref] [PubMed]
47. Huang Y, Hong W, Wei X. The molecular mechanisms and therapeutic strategies of EMT in tumor progression and metastasis. J Hematol Oncol 2022;15:129. [Crossref] [PubMed]
48. Fontana R, Mestre-Farrera A, Yang J. Update on Epithelial-Mesenchymal Plasticity in Cancer Progression. Annu Rev Pathol 2024;19:133-56. [Crossref] [PubMed]
49. Semenza GL. Targeting HIF-1 for cancer therapy. Nat Rev Cancer 2003;3:721-32. [Crossref] [PubMed]
50. Lequeux A, Noman MZ, Xiao M, et al. Targeting HIF-1 alpha transcriptional activity drives cytotoxic immune effector cells into melanoma and improves combination immunotherapy. Oncogene 2021;40:4725-35. [Crossref] [PubMed]
51. Zhao H, Wu L, Yan G, et al. Inflammation and tumor progression: signaling pathways and targeted intervention. Signal Transduct Target Ther 2021;6:263. [Crossref] [PubMed]
52. Lichterman JN, Reddy SM. Mast Cells: A New Frontier for Cancer Immunotherapy. Cells 2021;10:1270. [Crossref] [PubMed]
53. Christofides A, Strauss L, Yeo A, et al. The complex role of tumor-infiltrating macrophages. Nat Immunol 2022;23:1148-56. [Crossref] [PubMed]
54. Toledo B, Zhu Chen L, Paniagua-Sancho M, et al. Deciphering the performance of macrophages in tumour microenvironment: a call for precision immunotherapy. J Hematol Oncol 2024;17:44. [Crossref] [PubMed]
55. de Miguel M, Calvo E. Clinical Challenges of Immune Checkpoint Inhibitors. Cancer Cell 2020;38:326-33. [Crossref] [PubMed]
56. Morad G, Helmink BA, Sharma P, et al. Hallmarks of response, resistance, and toxicity to immune checkpoint blockade. Cell 2021;184:5309-37. [Crossref] [PubMed]
57. Korman AJ, Garrett-Thomson SC, Lonberg N. The foundations of immune checkpoint blockade and the ipilimumab approval decennial. Nat Rev Drug Discov 2022;21:509-28. [Crossref] [PubMed]
58. Zhou Z, Wang J, Wang J, et al. Deciphering the tumor immune microenvironment from a multidimensional omics perspective: insight into next-generation CAR-T cell immunotherapy and beyond. Mol Cancer 2024;23:131. [Crossref] [PubMed]
59. Wu B, Zhang B, Li B, et al. Cold and hot tumors: from molecular mechanisms to targeted therapy. Signal Transduct Target Ther 2024;9:274. [Crossref] [PubMed]
60. Finkin S, Yuan D, Stein I, et al. Ectopic lymphoid structures function as microniches for tumor progenitor cells in hepatocellular carcinoma. Nat Immunol 2015;16:1235-44. [Crossref] [PubMed]
61. Mao Y, Feng Q, Zheng P, et al. Low tumor purity is associated with poor prognosis, heavy mutation burden, and intense immune phenotype in colon cancer. Cancer Manag Res 2018;10:3569-77. [Crossref] [PubMed]
62. Cheng W, Ren X, Zhang C, et al. Bioinformatic profiling identifies an immune-related risk signature for glioblastoma. Neurology 2016;86:2226-34. [Crossref] [PubMed]
63. Hamzehnejadi M, Tavakoli MR, Homayouni F, et al. Prostaglandins as a Topical Therapy for Erectile Dysfunction: A Comprehensive Review. Sex Med Rev 2022;10:764-81. [Crossref] [PubMed]
64. Li F, He B, Ma X, et al. Prostaglandin E1 and Its Analog Misoprostol Inhibit Human CML Stem Cell Self-Renewal via EP4 Receptor Activation and Repression of AP-1. Cell Stem Cell 2017;21:359-373.e5. [Crossref] [PubMed]
65. Wu F, Zhang C, Zhao C, et al. Prostaglandin E1 Inhibits GLI2 Amplification-Associated Activation of the Hedgehog Pathway and Drug Refractory Tumor Growth. Cancer Res 2020;80:2818-32. [Crossref] [PubMed]
66. Patrício RPS, Videira PA, Pereira F. A computer-aided drug design approach to discover tumour suppressor p53 protein activators for colorectal cancer therapy. Bioorg Med Chem 2022;53:116530. [Crossref] [PubMed]
67. Malík M, Tlustoš P. Nootropics as Cognitive Enhancers: Types, Dosage and Side Effects of Smart Drugs. Nutrients 2022;14:3367. [Crossref] [PubMed]