Construction of an immune-related prognostic model and potential drugs screening for esophageal cancer based on bioinformatics analyses and network pharmacology

Abstract

Background: Esophageal cancer (ESCA) is a highly invasive malignant tumor with poor prognosis. This study aimed to discover a generalized and high-sensitivity immune prognostic signature that could stratify ESCA patients and predict their overall survival, and to discover potential therapeutic drugs by the connectivity map.

Methods: The key gene modules significantly related to clinical traits (survival time and state) of ESCA patients were selected by weighted gene coexpression network analysis (WCGNA), then the univariate and least absolute shrinkage and selection operator (LASSO) Cox regression analyses were used to construct a 15-immune-related gene prognostic signature.

Results: The immune-related risk model was related to clinical and pathologic factors and remained an effective independent prognostic factor. Enrichment analyses revealed that the differentially expressed genes (DEGs) of the high- and low-risk groups were associated with tumor cell proliferation and immune mechanisms. Based on the gathered data, a small molecule drug named perphenazine (PPZ) was elected. The pharmacological analysis indicates that PPZ could help in adjuvant therapy of ESCA through regulation of metabolic process and cellular proliferation, enhancement of immunologic functions, and inhibition of inflammatory reactions. Furthermore, molecular docking was performed to explore and verify the PPZ-core target interactions.

Conclusion: We succeed in structuring the immune-related prognostic model, which could be used to distinguish and predict patients' survival outcome, and screening a small molecule drug named PPZ. Prospective studies also are needed to further validate its analytical accuracy for estimating prognoses and confirm the potential use of PPZ for treating ESCA.

Keywords: esophageal cancer; immune‐related gene; long noncoding RNA; network pharmacology; perphenazine; prognostic model.

Figure 1

Flow diagram of this study. CMAP, connectivity map; DEGs, differentially expressed genes; ESCA, esophageal cancer; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; KM, Kaplan–Meier; LASSO, least absolute shrinkage and selection operator; lncRNAs, long noncoding RNAs; PPI, protein–protein interaction; ROC, receiver operating characteristic; TCGA, The Cancer Genome Atlas; WGCNA, weighted gene coexpression network analysis.

Figure 2

Construction of the immune‐related risk model. (A) Univariate Cox analysis identified survival‐related 34 mRNAs and 24 lncRNAs. p < .05 indicates a significant correlation between genes and prognosis, hazard ratio (HR) value > 1 means that the gene is a risk gene, and HR < 1 means a protective gene. (B) The best criteria to build the model based on LASSO regression. (C) Sankey diagram showing the coexpression relationship between the mRNAs and lncRNAs. LASSO, least absolute shrinkage and selection operator; lncRNAs, long noncoding RNAs; mRNA, messenger RNA.

Figure 3

Verification between the immune‐related risk model. (A) and (B‐a to c). The distributions of risk scores (a), survival times and status (b), and the heatmaps of gene expression levels (c) in the training set and validation set. The black dotted lines represent the median risk score cut‐off dividing patients into high‐ and low‐risk groups. The red dots represent the ESCA patients in high‐risk group, and the green dots represent the low‐risk. (A) and (B‐d) Kaplan–Meier survival analyses of the ESCA patients. (A) and (B‐e) ROC curves of the training set and validation set. (C) Multi‐index ROC curves of age, gender, grade, stage, TNM, and risk score. Time‐dependent ROC curves for 1–3 years. AUC, area under the receiver operating characteristic curve; ESCA, esophageal cancer; ROC, receiver operating characteristic.

Figure 4

Quantitative polymerase chain reaction verification of CDK9, ELNF1, FABP9, and STC2.

Figure 5

The heatmap of the clinical phenotypes (A) and distribution characteristic of pathological stage‐N (B) in the high‐ and low‐risk groups. (C) Differences in risk scores among different clinical features. IRGPI, immune‐related gene prognostic index; IRS, immune‐related risk score; TCGA, The Cancer Genome Atlas.

Figure 6

GO analysis of Gene Set Enrichment Analysis (GSEA) in the high‐ (A) and low‐risk (B) groups. GO, gene ontology.

Figure 7

The enrichment analysis of DEGs. (A) GO functional enrichment analysis. (B) KEGG pathway enrichment analysis. Circle color indicates P. Circle size represents the gene number enriched in the pathway. DEGs, differentially expressed genes; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 8

Structures and genetic networks of PPZ. (A) Two‐ and three‐dimensional chemical structure diagrams. (B) Protein–protein interaction (PPI) network. (C) Core targets of PPI network. PPZ, perphenazine.

Figure 9

The enrichment analysis of drug targets. (A) GO functional enrichment analysis. (B) KEGG pathway enrichment analysis. Circle color indicates P. Size represents the gene number enriched in the pathway. GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; TGF, transforming growth factor.

Figure 10

Docked complexes of PPZ with VAV3 (A, PDB ID: 2D86), BTK (B, PDB ID: 3GEN), and AR (C, PDB ID: 1E3G). The left views depict three‐dimensional docked pose of ligands with the target proteins and the right views are specific pics, showing the interaction for hydrogen bonds formed between PPZ and amino acids. AR, androgen receptor; Btk, Bruton tyrosine kinase; PPZ, perphenazine; Vav3, Vav guanine nucleotide exchange factor 3.