A novel risk model of three gefitinib-related genes FBP1, SBK1 and AURKA is related to the immune microenvironment and is predicting prognosis of lung adenocarcinoma patients

Gefitinib could inhibit cell growth and migration in LUAD

After adding an appropriate concentration of Gefitinib, we found that with increased Gefitinib concentration, the proliferative ability of cells decreased during the CCK-8 assay (Figure 1A–1C), and the apoptosis rate significantly raised during flow cytometry (Figure 1D). In addition, Gefitinib significantly inhibited cell migration and invasion during wound healing and Transwell assays (Figure 1E, 1F).

Figure 1. Gefitinib could inhibit cancer cell growth and migration. (A–C) Cell viability; (D) Cell apoptosis; (E) Cell Migration; (F) Cell Invasion.

The expression levels of Gefitinib target genes in LUAD

Gefitinib targets approximately 100 genes, primarily kinases (Table 1). Compared with normal lung tissues, 84 DEGs were identified in LUAD tissues (Table 2). The top 10 DEGs were visualized in a heatmap and histogram based on the fold change (Figure 2).

Figure 2. The expression levels of Gefitinib target genes in LUAD tissues. (A) High expression; (B) Low expression. Abbreviations: LUAD: Lung adenocarcinoma; T: LUAD tissues; N: normal tissues.

Biological functions, signaling mechanisms and protein-protein interaction (PPI) networks of the differentially expressed Gefitinib target genes (DEGTGs)

Gene ontology (GO) enrichment analysis revealed that the DEGTGs were significantly enriched in the intracellular signaling transduction, positive regulation of phosphatidylinositol 3-kinase activity, positive regulation of phosphatidylinositol 3-kinase signaling, receptor signaling protein tyrosine kinase activity, kinase activity, signal transduction, vascular endothelial growth factor-activated receptor activity, positive regulation of cell migration, vascular endothelial growth factor receptor signaling pathway, receptor complex, positive regulation of mitogen-activated protein kinase (MAPK) cascade, cell migration, growth factor binding, regulation of cell proliferation, negative regulation of apoptotic process, cell cycle, JNK cascade and others (Figure 3A–3C and Supplementary Table 1). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis revealed that the DEGTGs were enriched in the Rap1, Ras, MAPK, ERBB, PI3K-AKT, HIF-1, pathways in cancer, calcium, insulin, GnRH, and other signaling pathways (Figure 3D and Table 3). A PPI network constructed based on the DEGTGs is shown in Figure 4A, and the key genes in the PPI network are revealed in Figure 4B.

Figure 3. Functions and mechanisms of the DEGTGs. (A) BP; (B) MF; (C) CC; (D) KEGG. Abbreviations: DEGTGs: differentially expressed Gefitinib target genes; GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; BP: biological process; MF: molecular function; CC: cell component.

Figure 4. PPI network of the DEGTGs. (A) PPI network; (B) Hub DEGTGs in PPI network. Abbreviations: DEGTGs: differentially expressed Gefitinib target genes; PPI: protein-protein interaction.

Construction of a nomogram based on the DEGTGs

Survival analysis revealed that AURKA, AURKB, FBP1, FGR, HCK, LRRK2, and SBK1 levels were associated with poor prognosis in LUAD patients (Table 4). Receiver operating characteristic (ROC) analysis revealed that AURKA, AURKB, FBP1, FGR, HCK, LRRK2, and SBK1 had diagnostic values for LUAD, with the area under the curve (AUCs) of 0.961, 0.973, 0.910, 0.967, 0.895, 0.922, and 0.873, respectively. Based on the results of survival and ROC analysis, we constructed a nomogram based on the AURKA, AURKB, FBP1, FGR, HCK, LRRK2, and SBK1 (Figure 5).

Figure 5. Nomogram of the DEGTGs in LUAD. Abbreviations: DEGTGs: differentially expressed Gefitinib target genes; LUAD: lung adenocarcinoma.

Functions and clinical values involved in subgroups of the DEGTGs

When LUAD patients were classified based on the AURKA, AURKB, FBP1, FGR, HCK, LRRK2, and SBK1 expression levels and K = 2, PCA analysis revealed significant differences between cluster1 and cluster2 (Figure 6). According to our screening criteria analysis, 93 DEGs were found in LUAD tissues of cluster 2 compared with cluster 1 (Table 5). 93 DEGs were significantly enriched in neutrophil-mediated cytotoxicity, regulation of cytokine biosynthetic process, killing of cells of other organisms, disruption of cells of other organisms, regulation of neurotransmitter levels, cytokine biosynthetic process, cytokine metabolic process, positive regulation of response to drug, drug transport, cytolysis, hormone activity and others (Figure 7 and Supplementary Table 2). Cox regression analysis revealed that the expression levels of OLFM4, CDH7, CRABP1, RAB3B, SLC13A5, IGF2BP1, INHA, ABCC2, DKK1, KRT6A, MAGEC1, FGA, SPX, HOXA13, KRT6C, UPK1B, CPS1, and EVX1 were significantly associated with the OS in LUAD (Figure 8A). The levels of OLFM4, CDH7, CRABP1, RAB3B, SLC13A5, IGF2BP1, INHA, ABCC2, DKK1, KRT6A, MAGEC1, FGA, SPX, HOXA13, KRT6C, UPK1B, CPS1, and EVX1 in normal tissues and LUAD tissues are shown in a heatmap and histogram (Figure 8B, 8C).

Figure 6. The clinical roles involved in subgroups of the DEGTGs in LUAD. Abbreviations: DEGTGs: differentially expressed Gefitinib target genes; LUAD: lung adenocarcinoma.

Figure 7. Functions involved in subgroups of the DEGTGs in LUAD. (A) BP; (B) CC; (C) MF. Abbreviations: DEGTGs: differentially expressed Gefitinib target genes; GO: gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; BP: biological process; MF: molecular function; CC: cell component; LUAD: lung adenocarcinoma.

Figure 8. The prognosis values and levels of the DEGTGs in LUAD. (A) Prognosis related genes; (B, C) The levels of DEGTGs in LUAD. Abbreviations: DEGTGs: differentially expressed Gefitinib target genes; LUAD: lung adenocarcinoma.

Risk model of the DEGTGs

Univariate Cox regression analysis revealed that the expression levels of FBP1, FGR, AURKB, SBK1, LRRK2 and AURKA were significantly correlated with the OS of LUAD (Figure 9A). Multivariate COX regression analysis and the AIC method revealed that FBP1, SBK1, and AURKA were independent predictors of OS and were the influencing factors of the risk model (Figure 9B). Moreover, we found that Gefitinib promoted FBP1 mRNA and protein expression and inhibited the SBK1 and AURKA mRNA and protein expression in PC9 cells (Supplementary Figure 1). In addition, high-risk patients had a poor prognosis with the risk model based on FBP1, SBK1, and AURKA (Figure 9C–9E). Therefore, we constructed a nomogram for Gefitinib target genes FBP1, SBK1, and AURKA (Supplementary Figure 2). Univariate and multivariate Cox regression analysis revealed that the risk score and clinical stage correlated with poor prognosis in LUAD patients (Figure 10A, 10B). Therefore, we constructed a risk model-related nomogram (Figure 10C).

Figure 9. Risk model of the DEGTGs in LUAD. (A) Risk genes that affect patient prognosis in LUAD; (B–E) Patients with higher risk scores have the poor prognosis. Abbreviations: DEGTGs: differentially expressed Gefitinib target genes; LUAD: lung adenocarcinoma.

Figure 10. Prognosis-related risk scores and nomogram in LUAD. (A, B) Cox regression analysis that risk score is an independent risk factor for poor prognosis in LUAD patients; (C) Nomogram. Abbreviations: LUAD: lung adenocarcinoma.

Signaling mechanisms involved in risk model

In the risk model based on FBP1, SBK1 and AURKA, the high-risk score group was enriched in the cell cycle, oocyte meiosis, RNA degradation, basal transcription factors, spliceosome, ubiquitin-mediated proteolysis, homologous recombination, mismatch repair, p53 signaling pathway, regulation of autophagy, DNA replication, base excision repair and others (Table 6).

The risk score correlated with the immune microenvironment in LUAD

As shown in Figure 11, the risk score was significantly correlated with the immune score, estimate score, stromal score, and immune cells, including memory B cells, plasma cells, CD8 T cells, resting CD4 memory T cells, activated CD4 memory T cells, follicular helper T cells, regulatory T cells (Tregs), gamma delta T cells, resting NK cells, activated NK cells, monocytes, M1 macrophages, M2 macrophages, resting dendritic cells, resting mast cells, and activated mast cells (Figure 12). The expression levels of naïve B cells, memory B cells, resting dendritic cells, monocytes, M1 macrophages, plasma cells, resting mast cells, gamma delta T cells, activated CD4 memory T cells, resting NK cells, resting low-risk groups were significantly different (Supplementary Figure 3).

Figure 11. Risk score is related to the immune, estimate, and stromal scores in LUAD. (A) Immune score; (B) Estimate score; (C) Stromal score. Abbreviation: LUAD: lung adenocarcinoma.

Figure 12. Risk model is related to the immune cells in LUAD. (A) B cells memory; (B) Plasma cells; (C) T cells CD8; (D) T cells CD4 memory resting; (E) T cells CD4 memory activated; (F) T cells follicular helper; (G) T cells regulatory (Tregs); (H) T cells gamma delta; (I) NK cells resting; (J) NK cells activated; (K) Monocytes; (L) Macrophages M1; (M) Macrophages M2; (N) Dendritic cells resting; (O) Mast cells resting; (P) Mast cells activated. Abbreviation: LUAD: lung adenocarcinoma.

The risk score was significantly correlated with immune cell markers such as CSF1R, STAT5A, STAT5B, ITGAM, BCL6, CD274 (PD-L1), ITGAX, CD8B, and others (Figure 13 and Supplementary Figure 4). In addition, immune cell markers such as CSF1R, STAT5A, STAT5B, ITGAM, BCL6, CD274, ITGAX, CD8B, CD1C, HLA-DPB1, CD8A, NRP1, GZMB, CEACAM8, HLA-DRA, KIR3DL3, IRF5, TGFB1, HLA-DPA1, IL21, HLA-DQB1, KIR2DL3, and TNF were significantly different between the high- and low-risk groups (Supplementary Figure 5).

Figure 13. Risk model related to the immune cell markers in LUAD. (A) CSF1R; (B) STAT5A; (C) STAT5B; (D) ITGAM; (E) BCL6; (F) CD274; (G) ITGAX; (H) CD8B; (I) HLA-DPB1; (J) CD8A; (K) TGFB1; (L) HLA-DPA1; (M) IL21; (N) HLA-DQB1; (O) S4A4A; (P) CCR7; (Q) CD1C; (R) STAT3; (S) NRP1; (T)GZMB. Abbreviation: LUAD: lung adenocarcinoma.

PC9 cell culture and construction of cell models

The PC9 cell line was cultured in the RPMI-1640 with 10% fetal bovine serum at 37°C in an incubator with 5% CO₂. The PC9 cells were cultured in 6-well plates for 24 h in an incubator, then 0, 0.002, 0.02, 0.2, and 2 μM of Gefitinib were added.

Cell growth

The effects of Gefitinib on the proliferation and apoptosis of PC9 cells were detected by CCK-8 and flow cytometry. The experimental steps of CCK-8 were as follows: cell counting and 96-well plate plating, and 0, 0.002, 0.02, 0.2 and 2 μM Gefitinib were added to incubate cells. The CCK-8 solution at 24, 48, and 72 h incubation, and reading after 1 h incubation to appraisal the absorbance values of each group. The flow cytometry experiment steps were as follows: cell counting, and 6-well plate plating. 0, 0.002, 0.02, 0.2, and 2 μM Gefitinib were added to incubate cells, and stained according to the manufacturer’s instructions. The apoptotic PC9 cells were immediately detected using the FACS Caliber II Sorter and the Cell Quest FACS system (BD Biosciences, USA). The data were analyzed using FlowJo software (Version: 7.6.5).

Cell migration and invasion

In 6 well plate, the clean tip of pipette is used to create “wound” when cells form confluent monolayer, and cells cultured in RPMI-1640 were supplemented with 0% fetal bovine serum. Transfected cells were used to capture the image of wound edge with optical microscope at 0 h. After 24 h of incubation, the cell images of the same area were captured for measurement, and the wound healing rate was calculated [19]. Transfected cells were investigated for transwell migration. The basal lateral chamber was balanced with 600 μL RPMI-1640 medium containing 10% fetal bovine serum, while the apical chamber was added with 200 μL PC9 cell suspension without FBS. After 24 h of incubation at 37°C and 5% CO₂. According to the manufacturer's agreement of transwell System, fix the transwell chamber with 5% paraformaldehyde for 15 min, and then dye it with 0.1% crystal violet for 10 min. Transmembrane cells were calculated under inverted fluorescence microscope [19].

Acquisition and identification of Gefitinib target genes

The secondary structure of Gefitinib was obtained from the official website of the PubChem database, and then the target genes of Gefitinib were obtained from the SwissTargetPrediction database using the structure of Gefitinib. The transcriptome data of 59 normal tissues and 535 LUAD tissues were downloaded from the official website of The Cancer Genome Atlas (TCGA) database, and the gene expression data of Gefitinib targets in normal tissues and cancer tissues were extracted and analyzed using the Limma package.

Biological functions, signaling mechanisms and PPI network of DEGTGs

A P-value < 0.05 was used to screen for significant DEGTGs in LUAD tissues. The DEGTGs underwent GO and KEGG analyses to investigate the biological functions and signaling mechanisms of Gefitinib target genes. A PPI network of the DEGTGs was constructed in the STRING database. The PPI network was visualized by Cytoscape software, and the key genes were selected by the degree method using the CytoHubba plugin [19, 20].

Construction of a nomogram of prognostic-related Gefitinib target genes

The data of DEGTGs were merged with the prognostic information of LUAD patients, and patients with incomplete prognostic information were excluded. The patients were stratified based on the median value of the DEGTGs, and the relationship between the expression level of DEGTGs and prognosis was identified using survival analysis. The diagnostic value of prognostic genes AURKA, AURKB, FBP1, FGR, HCK, LRRK2, and SBK1 for LUAD was assessed using the ROC analysis. A prognostic nomogram of the Gefitinib target genes was constructed based on survival and ROC analysis results.

The roles of Gefitinib target gene subgroups in LUAD using consensus clustering

The gene expression data of AURKA, AURKB, FBP1, FGR, HCK, LRRK2, and SBK1 in LUAD tissues from TCGA were extracted. The Consensus ClusterPlus package clustered LUAD patients and determined the optimal k value to group the LUAD patients. Principal component analysis (PCA) was used to explore the differences between subgroups constructed based on the AURKA, AURKB, FBP1, FGR, HCK, LRRK2, and SBK1. The differentially expressed genes (DEGs) in the tissues of the two groups of LUAD patients were identified using the screening criteria logFC (>2 or <−2) and FDR < 0.05. The biological functions and signaling mechanisms involved in DEGs in the subgroups of Gefitinib target genes were explored using GO and KEGG analyses. The relationship between DEGs in the subgroups of Gefitinib target genes and OS of LUAD patients was analyzed using COX regression analysis, and the expression of prognostic genes in normal lung tissues and LUAD tissues was visualized.

Construction of risk model based on the DEGTGs

Univariate Cox regression analysis was used to reveal the relationship between the gene expression levels of DEGTGs AURKA, AURKB, FBP1, FGR, HCK, LRRK2, and SBK1 and OS in patients with LUAD, with P < 0.05 as the screening criteria. Then, the relationship between the expression levels of FBP1, FGR, AURKB, SBK1, LRRK2, and AURKA and the OS of LUAD patients was analyzed using multivariate COX regression analysis and the Akaike information criterion (AIC) method. The effects of Gefitinib on AURKA, FBP1, and SBK1 were assessed by qRT-PCR and western blotting, and a risk model was established based on these results.

qRT-PCR

Trizol (Thermo Fisher Scientific, China) was used to isolate total RNA from the PC9 cells after supplementing Gefitinib, and total RNAs were reverse-transcribed into cDNA using the Revertra Ace qPCR RT Kit (Toyobo life science, Japan). The UltraSYBR mixture (cwbio, China) was used to complete the PCR cycles using a StepOne Plus Real-Time PCR System (ABI, USA). The change in relative AURKA, FBP1 and SBK1 mRNA expression levels were calculated using the 2^ΔΔct method. The primer sequences of AURKA, FBP1 and SBK1, and β-actin used were as follows. AURKA forward: 5′-GTGGAGCATCAGCTCAGAAGA-3′ and reverse 5′-ACTGTTCCAAGTGGTGCATA-3′; FBP1 forward: 5′-ATCGATTGCCTTGTGTCCGT-3′ and reverse 5′-AAGCAGTTGACCCCACAGTC-3′; SBK1 forward: 5′-ACGTCACCAAGCA CTACGAA-3′ and reverse 5′-GTGATGCTCACCTCCCGTAG-3′; β-actin forward 5′-ACTCTTCCAGCCTTCCTTCC-3′ and reverse 5′-CGTCA TACTCCTGCTTGCTG-3′.

Western blotting

Total protein was isolated from PC9 cells by using cell lysate after Gefitinib supplementation. Protein gel electrophoresis was performed after quantitative analysis using the BCA method. After electroporation, the PVDF film was cut and incubated overnight at low temperature with the primary antibody. The second antibody was incubated, and the protein was exposed the following day. The primary antibodies used, along with their respective concentrations, were the 1:3000 FBP1 (Proteintech, China), 1:1000 SBK1 (HUABIO, China), 1:1000 AURKA (HUABIO, China), and 1:10000 GAPDH (Proteintech, China).

Construction of a nomogram of risk model

Univariate and multivariate COX regression analysis was used to reveal the relationship between the age, gender, clinical stage, T stage, N stage, M stage, risk score, and prognosis of cancer patients. Finally, a risk model related to the nomogram was constructed.

Signaling mechanisms of the risk model

Gene Set enrichment analysis (GSEA) is a commonly used method to reveal the functions and mechanisms involved in the risk models [21–23]. In our study, we divided the gene expression data of LUAD patients into high- and low-risk groups. The impact of high- and low-risk groups was analyzed on each gene set to reveal the signaling mechanism involved in the risk model in the GSEA software [21–23].

Analysis of the relationship between the risk model and the LUAD immune microenvironment

The gene expression data in 535 LUAD tissues were immunoscored by using the CIBERSORT and estimate methods. Correlation analysis was used to reveal the relationship between the risk score and the levels of immune cells, stromal score, immune score, estimate score and immune cell markers in LUAD. The levels of immune cells, stromal score, immune score, estimate score, and immune cell markers in the high- and low-risk groups were explored in LUAD.

Statistical analysis

The significance of the effect of Gefitinib on LUAD cell growth and migration and the difference in the expression of the immune cells, stromal score, immune score, estimate score and immune cell markers in high- and low-risk groups were identified using the Student’s t-test. Survival and Cox regression analyses were used to assess the prognostic value of Gefitinib target genes. The diagnostic values of Gefitinib target genes for LUAD were identified using the ROC analysis. The AUC close to 1 was associated with a good diagnostic value. A P-value < 0.05 was statistically significant.

Availability of data and materials

The data were acquired from publicly available databases, and additional data could be obtained from the corresponding author upon reasonable request.