Network Pharmacology Integrated Molecular Docking to Explore the Mechanism of Blister Beetle Therapy for Lung Adenocarcinoma

Abstract

Lung adenocarcinoma (LUAD) is one of the major causes of cancer death in the world. Studies show that the effective anticancer component in blister beetles is cantharidin, which can improve chemotherapy efficacy, median survival, and prognosis of LUAD. However, the antitumor mechanism of blister beetles has not been fully clarified. This study aimed to identify the key targets of the treatment of LUAD by blister beetles based on the principle of network pharmacology. An integrated approach including network pharmacology and a molecular docking technique was conducted, which mainly comprises target prediction, weighted gene correlation network analysis (WGCNA) analysis, network construction, gene ontology, and pathway enrichment analysis. 35 key targets were obtained and significantly associated with response to external stimuli, collagen binding, cyclin binding, organic acid binding, pyruvate metabolism, glycolysis, and amino acid biosynthesis pathways. Both LASSO regression and the RF model had a high predictive ability, and 9 candidate genes were screened, among which BIRC5 and PLK1 were the key targets for the treatment of LUAD by using blister beetles and showed significant survival significance. Cantharidin exerts its antitumor effects through 8 targets in 32 pathways, while BIRC5 and PLK1 have obvious survival significance.

Figure 1

WGCNA analysis. (a) Sample clustering and phenotypic heat map. (b) Filtering of soft thresholds. (c) The merge module. (d) Heat map of the correlation between modules and LUAD. (e) Correlation between modular genes and LUAD. (f) A volcano map of differentially expressed genes between LUAD/normal samples. (g) Screening of key genes in modules. Blue represents the differentially expressed genes between LUAD/control and pink represents the key module genes associated with lung adenocarcinoma acquired by WGCNA.

Figure 2

Enrichment analysis results of GO and KEGG of target genes. (a) Venn diagram of the intersection of genes, drug targets, and disease targets. (b) GO enrichment bars of key target genes. (c) KEGG pathway enrichment bubble map of key target genes.

Figure 3

Analysis results of the key target-function/-pathway regulatory network diagram and PPI regulatory network diagram. (a) Regulatory network diagram of the key targets and GO function. The yellow hexagon represents GO_MF terms, the purple hexagon represents GO_BP terms, the green hexagon represents GO_CC Terms, and the blue dot represents the key target. (b) Regulatory network diagram of key targets and pathways. The pink hexagons represent KEGG pathways and the blue dots represent key targets. (c) Protein interaction networks of key targets. The lines represent the interaction between them, the color indicates their degree value; the darker the color, the higher the degree value and the higher the core position.

Figure 4

Construction of the LUAD diagnostic model and screening of biomarkers. (a) The characteristic genes were screened by LASSO regression analysis. Deviance on the horizontal axis represents the proportion of residual explained by the model, showing the relationship between the number of characteristic genes and the proportion of residual explained by the model (Dev), and the coefficient of genes on the vertical axis (left). The x-coordinate is log (Lambda), and the y-coordinate represents the error of cross-validation (right). (b) Evaluation and validation of the diagnostic model by the ROC curve. (c, d) The cumulative residual distribution diagram of the sample (left) and the boxplot of the sample residual (right). The curve area indicates the cumulative residual value of the whole sample. The smaller the curve area is, the smaller the cumulative residual value of the sample is. The red dot represents the root mean square of the residual. (e) Importance of genetic variables in the RF models. (f) Evaluation of the RF model by the ROC curve.

Figure 5

Screening and expression analysis of characteristic genes. (a) Intersection of the LASSO diagnostic model gene and RF model gene. Pink represents the LASSO diagnostic model genes and blue represents the RF model genes. (b) ROC curves for candidate biomarkers. (c) Expression of biomarkers. (d) Correlation of biomarkers. Blue is negative and red is positive.

Figure 6

Construction of the pharmacological regulatory network map of Chinese medicine with biomarkers and results of molecular docking. (a) Network diagram of pharmacological regulation of biomarkers in traditional Chinese medicine. (b) 3D conformer structure schematic of cantharidin. (c) Docking result diagram of 3UEH and cantharidin. The green double-ring stick model is the active molecule cantharidin. (d) 3D conformer structure schematic of oleic acid. (e) Docking result diagram of 5TRM and oleic acid. The green double-ring stick model is the active molecule oleic acid. (f) 3D conformer structure schematic of 3-phenyl-4-azafluorene. (g) Docking result diagram of 1Q40 and 3-phenyl-4-azafluorene. The green double-ring stick model is the active molecule 3-phenyl-4-azafluorene. (h) 3D conformer structure schematic of 3-phenyl-4-azafluorene. (i) Docking result diagram of 3UEH and 3-phenyl-4-azafluorene. The green double-ring stick model is the active molecule 3-phenyl-4-azafluorene. Dotted gray/yellow lines represent hydrophobic bonds formed between active ingredients and amino acid residues. Each dotted gray/yellow line represents a hydrophobic bond.

Figure 7

Potential function, prognostic significance, and clinical relevance of biomarkers. (a) Bar chart of GO enrichment results for biomarkers. (b) Bar chart of KEGG enrichment results for biomarkers. (c) Kaplan–Meier survival profile of biomarkers. (d) Heat maps of biomarkers associated with clinical factors.