Entrectinib in Lung Adenocarcinoma: Network Toxicology Identification of Prognostic Genes, Toxicity Targets and Combination Drugs ()
1. Introduction
ROS1 fusions are rare—just 1% to 2% of all lung adenocarcinoma (LUAD) cases. But what they lack in prevalence, they make up for in clinical importance. This tiny patient subgroup has rewritten how we treat lung cancer. These patients don’t fit the typical lung cancer profile: they’re young (median age 50), almost never smoked (80%), and get blood clots far more often than patients with EGFR or ALK mutations [1].
Entrectinib changed everything for them. It’s an oral TKI that actually gets into the brain—something many older drugs couldn’t do. It hits TRK, ROS1, and ALK hard, with IC50 values in the single-digit nanomolar range (0.1 to 1.7 nM) [2]. Three big trials—STARTRK-1, STARTRK-2, and ALKA-372-001—told the same story: 67% of patients responded, and responses lasted a median of 16.8 months. For patients whose cancer had spread to the brain, 49% saw their brain tumors shrink [3]. The BFAST study went even further: in patients found to have ROS1 fusions via liquid biopsy, the response rate jumped to 81.5% [4].
Let’s put this in perspective. LUAD makes up 40% of all lung cancers. GLOBOCAN 2020 numbers show lung cancer still kills more men worldwide than any other cancer. And in Asia, LUAD rates are going up faster than anywhere else [5]. For all the hype around immunotherapy and targeted drugs, the 5-year survival rate for LUAD is still only 20%. We need better options. Badly.
Entrectinib isn’t a cure, though. Two big problems stand in the way. First, resistance. The L2026M gatekeeper mutation breaks the drug’s binding site and turns on MEK/ERK signaling, causing tumors to grow back—sometimes in as little as 6 weeks. Other times, cancer cells find workarounds through AXL or PI3K/MAPK pathways [4]. Lab studies show combining Entrectinib with trametinib, a MEK inhibitor, works better than either drug alone (combination index <0.9) [6]. Second, toxicity. Most side effects are mild—taste changes, dizziness, tiredness—but some are serious: brain fog, balance problems, heart rhythm issues. And while rare, fatal heart inflammation has been reported [3] [7]. Just last year, researchers found Entrectinib damages nerve cells by turning down THBS1 and blocking PI3K-AKT/TGF-β signaling [8].
No one has used network toxicology to pull all this together—how Entrectinib works, who benefits most, why it causes side effects, and what drugs to pair it with. That’s exactly what we’re doing here. We will: 1) find the key genes linking Entrectinib to LUAD and build a prognostic model; 2) connect these genes to Entrectinib’s toxicity; 3) identify new drug combinations that could overcome resistance.
2. Materials and Methods
2.1. Data Sources and Database Applications
All public data used in this study were obtained from internationally authoritative open databases. The application scenarios, data acquisition content, and core analysis parameters of each database are uniformly described as follows: Transcriptome sequencing raw count data and matched clinical follow-up information of lung adenocarcinoma (LUAD) patients were downloaded from The Cancer Genome Atlas (TCGA) database [9]. After sample quality control and gene ID annotation, a total of 542 LUAD tumor tissue samples and 59 adjacent normal lung tissue samples were included for differential expression analysis. After excluding samples with missing survival time, survival status, and key clinicopathological features, 522 LUAD tumor samples with complete clinical follow-up information were finally obtained for prognosis-related analysis. Meanwhile, gene type classification was performed based on this transcriptome dataset, and the expression matrix of protein-coding mRNAs was extracted to distinguish mRNAs from long non-coding RNAs (lncRNAs). The standardized two-dimensional chemical structure and three-dimensional conformation file of Entrectinib (CAS No.: 1108743-60-7) were obtained from the PubChem database for subsequent molecular docking verification [10]. Experimentally verified Entrectinib-target interaction data were retrieved from the ChEMBL database to construct the potential drug target set [11]. Based on the SwissTargetPrediction platform [12], with the species limited to human proteins, potential targets of Entrectinib were predicted by the reverse pharmacophore matching algorithm. The Similarity Ensemble Approach (SEA) database was used to predict potential targets of Entrectinib based on the principle of ligand structure similarity [13]. Using “Lung Adenocarcinoma” as the search term, LUAD-related disease genes were collected from the GeneCards database [14], the Online Mendelian Inheritance in Man (OMIM) database [15], and the Comparative Toxicogenomics Database (CTD) database [16]. Genes with disease association scores higher than the median in the GeneCards database were included in subsequent analysis to jointly construct the LUAD disease-related gene set. Using the well-validated STRING database, we mapped all protein-protein interactions among our candidate genes, limiting our search strictly to human proteins and keeping only interactions with a confidence score of 0.7 or higher, as is standard practice in the field [17]. We pulled all available LUAD cell line IC50 data directly from the GDSC repository to run our drug sensitivity correlation analyses [18], and for the molecular docking work that followed, we pulled core gene protein sequences from UniProt and their corresponding 3D structures from RCSB PDB, which gave us everything we needed to prepare the proteins for docking simulations [19].
2.2. Differential Expression Analysis
Differential expression analysis of transcriptome data between LUAD tumor tissues and adjacent normal tissues was performed using the DESeq2 package (version 1.32.0) of R software, with the screening threshold set to |log2(fold change)| ≥ 1 and adjusted P-value < 0.05 [20]. The ggplot2 package (version 3.3.6) was used to draw the volcano plot of differentially expressed genes, and the pheatmap package (version 1.0.12) was used to draw the expression clustering heatmap of the top 50 most significantly upregulated and downregulated genes (100 genes in total).
2.3. Integration and Intersection Analysis of Drug Targets and Disease-Related Genes
The obtained Entrectinib-related target genes were subjected to unified gene ID conversion and deduplication, and the union was taken to construct the potential target set of Entrectinib. The obtained LUAD disease-related genes were deduplicated, and the union was taken to construct the LUAD disease-related gene set. The overlap of the two gene sets was visualized using a Venn diagram [21]. Furthermore, the intersection of the potential Entrectinib target set, LUAD disease-related gene set, and differentially expressed gene (DEG) set was obtained to identify the common genes, which were regarded as core candidate genes through which Entrectinib may affect the occurrence and development of LUAD by regulating their expression. A ternary Venn diagram was drawn using an online Venn analysis tool or the VennDiagram package (version 1.7.3) of R software for visualization.
2.4. Functional Enrichment Analysis
Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) signaling pathway enrichment analysis were performed on the above intersection genes using the clusterProfiler package (version 4.4.0) of R software [22]-[24], with adjusted P-value < 0.05 as the significance threshold for enrichment. GO analysis covered three dimensions: biological process (BP), cellular component (CC), and molecular function (MF). The ggplot2 package (version 3.3.6) was used to draw bar charts, bubble charts, and circular charts to visualize the enrichment results.
2.5. PPI Network Construction and Hub Gene Screening
We dumped our final intersection gene list straight into STRING to pull the raw PPI network. Isolated nodes got cut right away—they don’t just clutter the plot, they skew all topological parameter calculations. The cleaned network loaded directly into Cytoscape 3.10.0; we plotted interactions and crunched all topology metrics, scaling both node size and color depth to match each gene’s degree score [25]. No single hub algorithm is perfect, so we ran four separate ones through CytoHubba: Degree, MCC, MNC, EPC. We took the top 20 from each. Only genes that appeared in all four lists made our final hub set [26]. We mapped this core subnetwork and reran full GO/KEGG enrichment on these key genes using the exact same steps laid out in Section.
2.6. Construction and Interpretability Analysis of the Prognostic Model
The standardized expression data of 93 intersection genes were integrated with the survival time, survival status, and clinicopathological information of 522 LUAD patients to construct a prognostic analysis dataset. First, univariate Cox proportional hazards regression analysis was used to screen genes significantly associated with the overall survival (OS) of patients (P < 0.05), and a forest plot was drawn to distinguish high-risk genes from low-risk genes. The significantly associated genes screened by univariate Cox regression analysis were included in multivariate Cox regression analysis to calculate the regression coefficient (β) of each gene. The risk score of each patient was calculated based on the following formula: Risk score = Σ (gene expression level × β). Patients were divided into high-risk and low-risk groups using the median risk score as the cutoff value. The Shapley Additive exPlanations (SHAP) method was used to perform interpretability analysis on the multivariate Cox model, quantify the contribution of each gene to model prediction (mean |SHAP value|), and visualize the regulatory effects of genes at the population and single-sample levels through beeswarm plots, force plots, and waterfall plots.
2.7. Survival Analysis and Model Prediction Performance Evaluation
The Kaplan-Meier method was used to draw the OS survival curves of LUAD patients in the high- and low-risk groups, and the log-rank test was used to compare the survival differences between the two groups. A risk score distribution plot (with risk score on the ordinate and samples sorted in ascending order of risk score on the abscissa, divided into high/low-risk groups by the median risk score, with the low-risk group in blue-purple and the high-risk group in orange), a survival status scatter plot (with survival time in years on the ordinate, dead samples in orange and alive samples in blue-purple), and a hub gene expression heatmap were drawn to evaluate the rationality of risk grouping. Time-dependent receiver operating characteristic (ROC) curves were used to evaluate the prediction accuracy of the model for 1-, 3-, and 5-year OS of patients, and the area under the curve (AUC) was calculated. Correlation analysis was performed using the timeROC package (version 0.4) of R software [27].
2.8. Independent Prognostic Analysis and Nomogram Construction
Clinicopathological features of enrolled LUAD patients, including gender, age, tumor histological grade (Grade), and tumor clinical stage (Stage), were collected. First, the chi-square test or Fisher’s exact test was used to analyze the association between high/low-risk groups and each clinicopathological feature [28]. Subsequently, the risk score (high/low-risk group) and the above clinicopathological features were jointly included in the multivariate Cox proportional hazards regression model to analyze whether the risk score was an independent prognostic factor for LUAD patients. A forest plot was drawn, and the concordance index (C-index), Akaike Information Criterion (AIC), and global log-rank P-value of the model were calculated [29]. Based on the independent prognostic factors with statistical significance in multivariate Cox regression analysis, the rms package (version 6.3-0) of R software was used to construct a nomogram for predicting 1-, 3-, and 5-year OS of patients. A calibration curve was used to evaluate the prediction accuracy of the nomogram, with the predicted survival rate on the abscissa and the actual observed survival rate on the ordinate. The consistency between the predicted value and the actual value was reflected by the fitting degree of the curve to the 45˚ diagonal line.
2.9. Molecular Docking Verification
We ran all core gene protein 3D structures through Discovery Studio 2019 for standard cleanup and preparation: we stripped out crystallographic water molecules, added polar hydrogens, assigned formal charges, and patched any missing amino acid residues or side chains. We saved the fully processed proteins in PDB format, then used the same software to run energy minimization on the Entrectinib 3D structure and saved it separately. Next, we converted both protein and ligand PDB files to PDBQT format with AutoDock4.0, and ran all docking simulations using AutoDock Vina 1.2.6. We sized the Grid Box to completely cover each protein’s putative binding pocket and left all other parameters at software defaults [30]. For each run, we picked the lowest binding energy conformation as the optimal binding mode, then used PyMOL 3.1 and Discovery Studio 2019 to visualize and document key ligand-receptor interactions: hydrogen bonds, hydrophobic contacts, and π-π stacking.
2.10. Drug Sensitivity Analysis
The OncoPredict (formerly calcPhenotype) method was used to predict the drug response of LUAD cell lines based on the expression levels of core genes [31]. Spearman correlation analysis was used to calculate the correlation coefficient (r) and significance between the expression level of each core gene and drug IC50, with |r| > 0.3 and false discovery rate (FDR) < 0.05 as the significant association threshold. The pheatmap package was used to draw a correlation heatmap (drug-gene-heat) to show the positive/negative correlation between gene expression and drug sensitivity, as well as the clustering results of genes with similar expression patterns. Meanwhile, a drug-gene summary table (drug-summary) was generated to reflect the number and association strength of core genes covered by a single drug. The ggplot2 package was used to draw scatter plots of the top 20 significantly associated drug-gene pairs to reflect the effect of gene expression on drug sensitivity intuitively.
3. Results
3.1. Data Filtering and Differential Expression Analysis
A total of 542 lung adenocarcinoma (LUAD) tumor tissue samples and 59 adjacent normal lung tissue samples (601 samples in total) were included for differential expression analysis. Meanwhile, 522 LUAD tumor samples with complete follow-up information were obtained after clinical data cleaning for subsequent prognostic model construction and validation. A total of 3635 significantly differentially expressed genes (DEGs) were identified by differential expression analysis, including 1798 upregulated genes and 1837 downregulated genes.
The overall experimental results are shown in Figure 1. The volcano plot (Figure 1(A)) illustrated the magnitude of differential expression and statistical significance of genes at the whole-genome level. The top 50 most significantly upregulated and downregulated DEGs were selected to draw an expression clustering heatmap (Figure 1(B)). The results showed that these 100 genes could completely distinguish LUAD tumor tissues from normal lung tissues, and the expression patterns of samples within the same group were highly clustered.
Figure 1. Identification and visualization of differentially expressed genes (DEGs). (A) Volcano plot showing the distribution of DEGs. (B) Heatmap of DEGs based on expression profiles.
3.2. Entrectinib Target Prediction and Collection of LUAD Disease-Related Genes
A total of 475 potential targets of Entrectinib were obtained through joint prediction of multiple databases. Among them, 66 targets were from the ChEMBL database, 137 from the SwissTargetPrediction platform, and 371 from the SEA database. The common intersection target of the three databases was 5, and the SEA database had the largest number of unique targets (282, accounting for 59.4%) (Figure 2(A)).
Meanwhile, 19858, 204, and 490 LUAD-related genes were obtained from the GeneCards, OMIM, and CTD databases, respectively. After merging and deduplication, a total of 19912 LUAD disease-related genes were obtained. The common intersection gene of the three databases was 6, and the unique genes from the GeneCards database accounted for 96.5% (Figure 2(B)).
3.3. Identification of Compound-Disease-Differential Expression Intersection Genes
The intersection of Entrectinib target genes (475), LUAD disease-related genes (19912), and DEGs (3635) was taken, and a total of 93 common intersection genes were obtained (Figure 2(C)). These genes simultaneously possessed three major characteristics: being potential target binding genes of Entrectinib, participating in the occurrence and development of LUAD, and being significantly differentially expressed in tumor tissues. Therefore, they were identified as core candidate genes.
![]()
Figure 2. Venn diagram analysis for target screening and hub gene identification in Entrectinib against lung adenocarcinoma. (A) Intersection of potential targets of Entrectinib predicted by ChEMBL, SEA, and SwissTargetPrediction databases. (B) Intersection of LUAD-associated genes obtained from CTD, GeneCards, and OMIM databases. (C) Overlap of Entrectinib targets, lung adenocarcinoma (LUAD)-related genes, and genes from the TCGA-LUAD cohort. (D) Intersection of hub genes screened by four topological algorithms (Degree, EPC, MCC, MNC) in cytoHubba.
3.4. Functional Enrichment Analysis of Intersection Genes
3.4.1. Gene Ontology (GO) Functional Enrichment Analysis
The results of GO enrichment analysis showed that the 93 intersection genes were mainly enriched in biological process (BP) terms related to cell cycle and mitosis regulation, including protein autophosphorylation, mitotic cell cycle phase transition, positive regulation of cell cycle, nuclear division, and G1/S phase transition (adjusted P < 1 × 10−6). In terms of cellular component (CC), they were mainly enriched in structures such as chromosomal regions, kinase complexes, centromeres, and spindles. In terms of molecular function (MF), they were mainly enriched in protein serine/threonine kinase activity, protein tyrosine kinase activity, and transmembrane receptor protein kinase activity (adjusted P < 1 × 10−6). The bar chart, bubble chart, and circular chart are shown in (Figures 3(A)-(C)), respectively.
![]()
Figure 3. Functional enrichment and PPI network analysis of hub genes. (A) Bar plot of Gene Ontology (GO) biological process enrichment. (B) Gene Ontology (GO) cellular component enrichment results visualized by bubble plot. (C) Circos plot showing the annotation of hub genes in related pathways. (D) Bar plot of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment. (E) Bubble plot showing the top enriched KEGG pathways. (F) Heatmap of the correlation between hub genes and enriched pathways.
3.4.2. Kyoto Encyclopedia of Genes and Genomes (KEGG) Signaling Pathway Enrichment Analysis
KEGG enrichment analysis showed that the intersection genes were significantly enriched in pathways such as cell cycle, PI3K-Akt signaling pathway, MAPK signaling pathway, EGFR tyrosine kinase inhibitor resistance, and non-small cell lung cancer (adjusted P < 0.05). Among them, the cell cycle pathway was the core pathway with the largest number of enriched genes. The bubble chart showing the top 30 significantly enriched pathways is presented in (Figure 3(D), Figure 3(E)). KEGG enrichment analysis showed that they were significantly enriched in PI3K-Akt, FoxO, non-small cell lung cancer, and PD-L1/PD-1 pathways (Figure 3(F)).
3.5. Construction of Protein-Protein Interaction (PPI) Network and Screening of Core Genes
The constructed PPI network contained 93 nodes and 1872 protein interaction relationships, showing significant modular aggregation characteristics (Figure 4(A)). The redder the node color and the larger the node volume, the higher the degree value.
Four algorithms (Degree, MCC, MNC, and EPC) of the CytoHubba plugin were used to screen the top 20 genes, and the results were completely consistent. Finally, 20 core genes (Hub genes) were identified: CDK1, CCNA2, CCNB1, CHEK1, EZH2, AURKA, AURKB, PLK1, TOP2A, CCNB2, TTK, BUB1, KIF11, CDC7, CCNE2, PLK4, NEK2, MELK, CCNA1, and INCENP. Among them, CDK1 was the core node with the highest Degree value in the network (Figure 4(B)).
GO enrichment analysis of core genes further confirmed their involvement in mitosis regulation, PI3K-Akt/MAPK signal transduction, and tyrosine phosphorylation.
Figure 4. Protein-protein interaction (PPI) network construction and hub gene selection. (A) Global PPI network of candidate target genes. (B) The top 20 hub genes identified using the cytoHubba plugin in Cytoscape.
3.6. Construction of Prognostic Model Based on Core Genes
3.6.1. Screening of Prognosis-Related Genes
Univariate Cox regression analysis showed that the expression levels of 34 out of the 93 intersection genes were significantly associated with the overall survival (OS) of LUAD patients (P < 0.05), and all 20 PPI core genes were included. The forest plot (Figure 5) showed that there were 26 high-risk genes (HR > 1), whose high expression was associated with poor prognosis, and 8 low-risk genes (HR < 1), whose high expression was associated with favorable prognosis.
Figure 5. Forest plot of univariate Cox regression analysis for prognostic hub genes.
3.6.2. Construction of Prognostic Risk Score Model
Based on the multivariate Cox regression coefficient (β), the risk score formula was constructed: Risk score = Σ (gene expression level × β coefficient). The representative regression coefficients were as follows: MCM4 (−0.5113), TTK (−0.3258), and BUB1 (−0.4050) among low-risk genes; PRKDC (0.3641), CDK1 (0.3328), CCNA2 (0.3095), PLK1 (0.3104), and CHEK1 (0.2451) among high-risk genes. Taking the median risk score as the cutoff, 494 patients were divided into the high-risk group (232 cases) and the low-risk group (262 cases).
SHAP Interpretability Analysis of the Prognostic Model
SHAP (SHapley Additive exPlanations) analysis quantified the contribution of each gene to the model at both the population and single-sample levels. The results showed that MCM4 (mean|SHAP| = 0.396), BUB1 (0.359), and CDK1 (0.309) were the driver genes with the highest contribution. The SHAP direction was completely consistent with the Cox regression coefficient: the higher the expression of high-risk genes such as CDK1, and PLK1, the larger the SHAP value (positive contribution); the higher the expression of low-risk genes such as MCM4, BUB1, and TTK, the smaller the SHAP value (negative contribution). The swarm plot, force plot, and waterfall plot (Figures 6(A)-(D)) further verified the central role of core genes in prognostic regulation.
![]()
Figure 6. SHAP analysis results for the prognostic gene signature. (A) Mean SHAP values showing the relative importance of each gene. (B) Beeswarm plot illustrating the distribution of SHAP values for each gene. (C) Waterfall plot demonstrating the contribution of top genes to a single sample prediction. (D) Force plot visualizing the prediction decomposition for a representative sample.
3.7. Survival Analysis and Model Evaluation
3.7.1. Kaplan-Meier Overall Survival Analysis
The Kaplan-Meier survival curve showed that the overall survival rate of patients in the low-risk group was significantly higher than that in the high-risk group, and the difference between groups was extremely statistically significant (P < 0.001). With the extension of follow-up time, the decline rate of the number of samples in the high-risk group was significantly faster (Figure 7(A)). The risk score distribution plot showed that the scores of the high and low-risk groups were clearly demarcated without overlap (Figure 7(B)). The survival status plot showed that most death cases were concentrated in the high-risk score area, and most surviving cases were distributed in the low-risk area (Figure 7(C)). The core gene expression heatmap showed that genes such as CDK1, TOP2A, CCNA2, PLK1, CCNB1, and CHEK1 were highly expressed in the high-risk group and lowly expressed in the low-risk group; genes such as BTK, EPHB2, HTR1D, and showed the opposite pattern (Figure 7(E)).
![]()
Figure 7. Prognostic value validation of the constructed gene signature. (A) Kaplan-Meier overall survival curves for high- and low-risk groups. (B) Distribution of risk scores in all patients. (C) Scatter plot of patient survival status corresponding to risk scores. (D) Time-dependent ROC curves for 1-, 3- and 5-year overall survival prediction. (E) Heatmap of gene expression profiles between high- and low-risk groups. (F) Univariate and multivariate Cox regression analysis of independent prognostic factors.
3.7.2. ROC Curve Evaluation of Model Predictive Performance
Time-dependent receiver operating characteristic (ROC) curve analysis showed that the AUC values of the model for predicting 1-year, 3-year, and 5-year overall survival were 0.707, 0.708, and 0.678, respectively (Figure 7(D)), all significantly higher than 0.5, confirming that the model had stable predictive accuracy for short-term and medium-term prognosis.
3.8. Independent Prognostic Analysis
3.8.1. Correlation Analysis between Clinicopathological Features and Risk Score
Cross-analysis of high/low-risk groups with clinicopathological features such as gender, age, and tumor clinical stage showed that the risk score was significantly associated with tumor clinical stage (P < 0.05), but not with other features (Figure 7(F)).
3.8.2. Multivariate Cox Proportional Hazards Regression Analysis
To verify whether the prognostic risk score was an independent prognostic factor, multivariate Cox regression analysis was performed by incorporating risk score, tumor clinical stage, gender, and age (Figure 7(F)). A total of 476 LUAD patients with complete information were included in the analysis. The results showed that the hazard ratio (HR) of the risk score (high-risk group vs low-risk group) was 2.034 (95% CI: 1.479 - 2.80, P < 0.001). For tumor clinical stage, taking Stage I as the reference, the HR of Stage II was 2.148 (95% CI: 1.475 - 3.13, P < 0.001), Stage III was 2.912 (95% CI: 1.963 - 4.32, P < 0.001), and Stage IV was 3.772 (95% CI: 2.142 - 6.64, P < 0.001). The HR of gender (male vs female) was 0.973 (95% CI: 0.718 - 1.32, P = 0.8625), and the HR of age was 1.013 (95% CI: 0.997 - 1.03, P = 0.1129). The overall evaluation indicators of the model showed that the global Log-Rank P value was 8.64 × 10−13, the AIC was 1769.18, and the concordance index (C-index) was 0.71, indicating that the model fitted well.
The above results confirmed that the prognostic risk score was still an independent prognostic factor for the overall survival of LUAD patients after excluding the confounding effects of clinicopathological factors such as gender, age, and tumor clinical stage.
3.9. Molecular Docking Validation
To verify the direct binding potential of Entrectinib with the proteins encoded by core prognostic genes, four high-risk genes, CCNB1, PRKDC, F2, and CHEK1, were selected for molecular docking. The docking results of CCNB1 showed that Entrectinib stably embedded into the hydrophobic pocket of CCNB1, formed hydrogen bonds with ASP-86 (2.7 Å) and LYS-88 (2.3 Å), and also had hydrophobic and π-π stacking interactions (Figure 8(A)). The docking results of PRKDC showed that Entrectinib bound to the functional pocket, formed a hydrogen bond with GLY-491 (3.4 Å), and had various non-covalent interactions such as halogen bonds (Figure 8(B)). The docking results of F2 showed that Entrectinib embedded into the hydrophobic pocket, formed a hydrogen bond with LEU-40 (3.4 Å), and had π-π stacking and van der Waals forces (Figure 8(C)). The docking results of CHEK1 showed that Entrectinib precisely bound to the kinase domain, forming key hydrogen bonds (2.1 Å) and salt bridges (Figure 8(D)). All docking binding energies were less than −5.0 kcal/mol, indicating that Entrectinib formed stable molecular binding with the above proteins.
![]()
Figure 8. Molecular docking results of Entrectinib with four key hub genes. (A) Binding mode of Entrectinib with CCNB1. (B) Binding mode of Entrectinib with PRKDC. (C) Binding mode of Entrectinib with F2. (D) Binding mode of Entrectinib with CHEK1. Each panel shows the overall protein structure, detailed binding interactions, surface representation and 2D interaction diagram.
3.10. Correlation Analysis between Core Genes and Drug Sensitivity
3.10.1. Quantitative Characteristics of Significantly Associated Drugs and Genes
Spearman correlation analysis (|r| > 0.3, FDR < 0.05) showed that BMS-754807 was significantly associated with 8 core genes (the most), followed by SB216763 (7), Cisplatin (6), and Lapatinib (6). Doramapimod had the highest average correlation coefficient (mean|r|) (Figures 9 (A)-(C)).
Figure 9. Screening results of potential therapeutic drugs targeting hub genes. (A) Enrichment analysis of drug sensitivity with |NES|>1 and FDR<0.25. (B) Bubble plot showing the association between candidate drugs and hub genes. (C) Bar plot displaying the number of significantly associated hub genes for each drug.
3.10.2. Correlation between Core Genes and Chemotherapy and Targeted Drug Sensitivity
IRF1 was extensively positively correlated with various chemotherapeutic drugs such as Paclitaxel, Cisplatin, and Gemcitabine (r = 0.30~0.45), suggesting that it was a core drug resistance driver gene. IRF2 and IRF3 were significantly negatively correlated with multiple drugs (r = −0.30 - −0.41), indicating that they were core sensitivity driver genes. CCNB1 was negatively correlated with targeted drugs such as OSI-027 and Vorinostat (r = −0.31 - −0.36) (Figure 10).
Figure 10. Spearman correlation analysis of hub genes and drug sensitivity.
The correlation heatmap (Figure 10) showed that genes such as IRF1, IRF2, and IRF3 exhibited extensive drug sensitivity associations, while genes such as CCNB1, IKBKE, and IMAC1 showed specific associations (Figure 10).
Entrectinib sensitivity analysis showed that high expression of IRF1 (r = 0.37), LCK (r = 0.31), IKBKE (r = 0.33), and CCNB1 (r = 0.31) was significantly positively correlated with drug resistance; no significantly negatively correlated (sensitivity-related) core genes were observed (Figure 10).
3.10.3. Univariate Correlation Validation of Key Gene-Drug Pairs
The scatter plot validation results of key gene-drug pairs were highly consistent with the Spearman analysis (Figures 10 and Figures 11(A)-(T)). Specifically, ABCB1 was significantly negatively correlated with JQ1 (R = −0.44, P = 2.2e−16) (Figure 11(A)), suggesting that high expression of ABCB1 could enhance the sensitivity of LUAD cells to JQ1. CCND2 was significantly negatively correlated with BMS-754807, SB216763, and Doramapimod (R = −0.41 - −0.42), indicating that high expression of CCND2 was associated with increased sensitivity to the above drugs (Figures 11(B)-(D)). CDK2 was significantly positively correlated with Doramapimod (R = 0.65, P = 2.2e−16), suggesting that high expression of CDK2 mediated Doramapimod resistance. ITK was significantly negatively correlated
Figure 11. Spearman correlation validation scatter plots of crucial prognostic hub gene-drug pairs.
with JQ1 and SB216763 (R = −0.41, P = 2.2e−16) (Figure 11(M), Figure 11(N)), indicating that high expression of ITK enhanced drug sensitivity. NTRK3 was significantly negatively correlated with BMS-754807 and SB216763 (R = −0.45, P = 2.2e−16), suggesting that high expression of NTRK3 enhanced drug sensitivity (Figure 11(P), Figure 11(T)).
3.10.4. Summary of Correlation Analysis Results
In summary, the expression of core Hub genes was significantly and specifically associated with LUAD drug sensitivity. IRF1 was the core drug resistance driver gene, IRF2/3 were the core sensitivity driver genes, and CCNB1, IKBKE, IMAC1, etc., could be used as potential biomarkers for individualized medication.
4. Discussion
This study integrated network toxicology and bioinformatics to screen 20 core genes of Entrectinib in lung adenocarcinoma (LUAD), constructed a prognostic risk model, and explored its multidimensional application in ROS1-positive LUAD. Unlike previous studies, this work established a “core gene-prognosis-toxicity-combination therapy” framework.
Strikingly, every single one of the 20 core genes we pulled encodes a cell cycle regulator—some of these are sharply upregulated in LUAD tumors and all track with poorer patient survival. SHAP value analysis zeroed in on MCM4, BUB1, and CDK1 as the three strongest prognostic drivers. CDK1, the non-negotiable master regulator of G2/M transition, has long been tied to TKI treatment failure [32]; miR-6077-driven CDK1 overexpression bypasses G2/M arrest entirely and confers robust cisplatin-pemetrexed resistance. PLK1, too, drives chemoresistance: it is markedly upregulated in taxane-refractory LUAD cells, induces ABCB1 expression via c-Myc, and PLK1 inhibition completely restores drug sensitivity. Our docking data further shows CHEK1 forms a tight 2.1 Å hydrogen bond with Entrectinib, and high tumoral CHEK1 levels independently predict shorter overall survival in NSCLC patients [33] [34].
Regarding toxicity, Entrectinib binds PRKDC, F2, and CHEK1. PRKDC deficiency increases neurotoxicity in mice [35]; our docking (GLY-491, 3.4 Å) suggests Entrectinib may inhibit PRKDC, contributing to neurotoxicity (58% dizziness). F2 binding (LEU40, 3.4 Å) provides a molecular basis for thrombosis risk in ROS1positive patients [36]. CHEK1 binding may explain myelosuppression, as CHEK1 maintains hematopoietic stem cells [37].
Our drug sensitivity screens hit on a weird IRF family pattern: high IRF1 levels tracked perfectly with Entrectinib resistance (r = 0.37), while IRF2 and IRF3 consistently marked better treatment response. This dual nature of IRF1—it can either boost TKI sensitivity or drive immune escape depending on tumor context—goes a long way toward explaining its wildly conflicting prognostic reports, and no one has really connected this to Entrectinib before [38]. JQ1, a BET bromodomain inhibitor, had the cleanest negative correlation with ABCB1 and ITK (R = −0.44) and cranks down MYC signaling to reverse TKI resistance [39]. BMS-754807 (IGF-1R inhibitor) hits eight of our 20 core genes, which makes total sense for dual ROS1/IGF-1R pathway blockade [40]. SB216763 (GSK-3β inhibitor) targets seven core genes, and no one’s really explored this for ROS1-positive LUAD. Most tellingly, CHEK1 inhibitors already show robust synergy with TKIs in preclinical NSCLC models, which lines up perfectly with our core gene findings [41].
Our docking runs turned up something no one had reported before: Entrectinib binds four completely non-canonical targets—CCNB1, PRKDC, F2, and CHEK1—hinting at a far broader inhibition profile than the 8 out of 456 kinases measured in the original kinase panel [41]. Our work has clear limitations: all data comes exclusively from TCGA, drug sensitivity calls are cell line-only, we have no in vivo toxicity data, and most mechanistic links are correlative rather than causal. To solidify these findings, future work will need patient-derived organoids, prospective patient cohorts, CRISPR knockout screens, and direct SPR binding assays to confirm these novel interactions.
5. Conclusions
This study established the first integrated “core gene-prognosis-toxicity-combination therapy” framework for Entrectinib in ROS1-positive LUAD via network toxicology and bioinformatics.
We identified 20 hub genes (e.g., CDK1, PLK1, CHEK1) enriched in cell cycle pathways, and constructed a prognostic model with stable performance (1/3/5-year OS AUCs: 0.707/0.708/0.678) that was an independent prognostic factor (HR = 2.034, P < 0.001).
Molecular docking revealed Entrectinib stably binds to PRKDC, F2 and CHEK1, explaining its neurotoxicity, thrombotic risk and myelosuppression. IRF1 was identified as a core resistance biomarker, while IRF2/3 as sensitivity biomarkers. JQ1, BMS-754807 and SB216763 were screened as promising combination agents.
This study provides a theoretical basis for precise Entrectinib treatment, and the findings require further preclinical and clinical validation.
Authors’ Statement
Bin Ning: Writing—original draft, Methodology. Fangwei Xiong: Software, Data curation. Jiale Shao: Validation, Conceptualization. Liangcheng Huang: Validation, Data curation. Yanting Xiao: Writing—review & editing, Methodology, Supervision, Funding acquisition.
Funding
This research was supported by: Guangxi provincial-level project of university students’ innovation and entrepreneurship (Grant No. S202310599064).
Acknowledgements
I am sincerely thankful to my colleagues for their dedication and support throughout this research. I would also like to thank the WeChat Official Account “DoubleStrand Bioinformatics” for sharing practical bioinformatics tutorials, data analysis skills, and valuable research ideas, which helped me overcome many technical difficulties and complete this study smoothly.
NOTES
*These authors made the same contribution to this article.
#Corresponding author.