We downloaded the gene expression, clinical data, Pan-Cancer Atlas Hub from GEO database¹ (Chen C. et al., 2021), UCSC Xena² (Li F. et al., 2021), GDSC³.

Four hundred and thirty three LUAD tissues and 19 normal tissues were selected from the GSE68465 (GPL96). By GEO2R (Liu et al., 2019), P < 0.01, Log FC > 2 or Log FC < -2 were defined as screening criteria. Next, Metascape⁴ (Liu et al., 2020), KEGG⁵ (Zheng et al., 2021), and DAVID⁶ (Li J. et al., 2021) databases were used to analyze biological process (BP), cellular component (CC), molecular function (MF) and related pathways of differentially expressed mRNAs (DEmRNAs).

UALCAN⁷ (Liu H. et al., 2021) is a comprehensive, user-friendly, and interactive web resource for analyzing cancer OMICS data. Combing with clinical data, we explored the correlation between CDK2 expression and clinical indicators, such as age, grade, sex, smoking habits, stage, TP-53 mutation, weight. Heat maps were showed the positive and negative related genes with CDK2. We analyzed the overall survival (OS) and Relapse-Free Survival (RFS) about CDK2 for LUAD patients by PrognoScan⁸ (Liu X. et al., 2021).

Compared with normal lung tissues, CDK2 expression was highly expressed (P < 0.05) in LUAD by Oncomine⁹ (Dai et al., 2021) and TIMER¹⁰ (Wu R. et al., 2021) databases. As for 33-cancers, we studied the relationships between the expression of CDK2 and clinical stage, survival conditions through different “R” packages.

To research the major molecule CDK2, we showed 50 positive and negative genes in LUAD by heat maps (Person test method). The 1, 3, 5, and 8 years’ ROC curves were constructed by ‘‘ROC’’ R package. The abscissa of the prediction model was False Positive (FP) and the ordinate was True Positive (TP). Next, we screened out the top 300 genes with the most significant positive correlation with CDK2 for enrichment analysis. Bar and bubble charts showed the classical functions of 300 genes by ‘‘cluster-profiler’’ R package. We downloaded the GSEA4.1.0¹¹ (Cao et al., 2021a) to investigate Kyoto Encyclopedia of Genes and Genomes (KEGG) of CDK2-related 300 genes. The top of 50 terms were sorted by P-value in circle graph.

CIBERSORT¹² (Tan et al., 2021) is an analysis tool to provide an estimation of the abundances of member cell types in a mixed cell population using gene expression data. We selected immune cells associated with CDK2 (Pearson correlation coefficient r > 0.15). Four CDK2-related immune cells were shown in circle graph.

R package was used to analyze the relation of CDK2 and immune cell infiltration score (P ≤ 0.05). A total of 33-cancers, according to the median of CDK2 gene expression, the samples were divided into high expression groups and low expression groups. Then, we discussed the different expression of immune cells (P < 0.05). The above results were shown by scatter plots and box plots.

TISIDB¹³ (Huang X. Y. et al., 2021) is a web portal for tumor and immune system interaction, which integrates multiple heterogeneous data types. First, we found CDK2-related Immunomodulators with P-value < 0.05. Next, cBioPortal¹⁴ (Lu et al., 2021) was used to explore the genes related to Immunomodulators. DAVID and STRING¹⁵ (Zhuang et al., 2021) were used to analyze the Gene Ontology (GO), protein interaction of 49 genes. GSEA 4.1.0 further analyzed the NOM p-val and FDR q-val of these immune related genes.

Through LUAD clinical data, “survival,” “survminer” R packages were used to construct the COX regression model of age, gender, and stage. Forest plot (Xiang et al., 2021) was shown related Hazard ratio (HR), Log-Rank, Concordance index of LUAD. “rms” R package was used to structure the Nomogram model (Liu et al., 1976) of age, gender, stage, T, N, M.

Using “Survival” R package, we explored the univariate Cox analysis of 49 genes (P < 0.05). Different P-value and HR of statistically significant genes were shown in forest map. According to the results of univariate regression analysis, the significant genes were further analyzed by multivariate analysis. According to the risk score of genes selected by multiple factors, the LUAD samples were divided into lower-risk and higher-risk groups. The risk score was as follows: Risk score = Expgene1 × Coefgene1+Expgene2 × Coefgene2+ Expgene3 × Coefgene3+ Expgene4 × Coefgene4. “survival,” “survminer,” “survival ROC” R packages were performed survival probability curve and ROC curve. Combining the patient’s high and low risk situation and the state of life, death, we used “pheatmap” R package to show the heat map.

We use Human Protein Atlas (HPA) database to verify CDK2 expression in LUAD tissues and normal lung tissues.

We explored CDK2 related miRNAs from Targetscan¹⁶ (Barnett-Itzhaki et al., 2021), miRWalk¹⁷ (Sheng L. P. et al., 2021), mirDB¹⁸ (Ren et al., 2021), Starbase¹⁹ (Li et al., 2014). All miRNAs from four databases were analyzed by Venn map²⁰. Kaplan–Meier Plotter software was used to explore the prognostic value of miRNAs with LUAD. The lncRNAs that regulate common miRNAs were screened by Starbase database (Li et al., 2014). The Cytoscape (Han et al., 2021) was used to build ceRNA network and Cytohubba selected the key node genes according to the degree in the network.

We did the basic experimental verification of PCR about CDK2 in LUAD cell lines (A549, H1299, H1975) and normal bronchial epithelial cell line (BEAS-2B). GraphPad Prism7²¹ software was used to count the differences between cancer cell lines and normal cell line.

We downloaded the response data of 192 anti-tumor drugs in more than 1000 cancer cell lines. “Ggplot2” R package was used to explore the correlation between the CDK2 expression and IC50 of 192 anti-tumor drugs by box diagrams and point diagrams.

A total of 250 genes were chosen by GEO2R according to the P-value and LogFC. 161 mRNAs were highly expressed and 78 mRNAs were expressed at lower levels in 19 adjacent non-LUAD tissues and 433 LUAD tissues. The study design of this project was shown in Figure 1. The volcano map, scatter plots, peak plot of GSE68465 samples were displayed in Figure 2. A total of 15 highly expressed genes and 15 low expressing genes were shown in Table 1.

Through Metascape software, we identified a number of pathways significantly enriched, including activation of immune response, regulation of cell adhesion, T cell activation involved in immune response, regulation of cell adhesion, PID-HNF3B PATHWAY (Figure 3A). The interactions between different pathways were shown in Figure 3B. Figures 3C,D drew the P-value of different pathways. According to the relationship of DE mRNAs, several protein analyses were performed in the Figure 3E. CDK2, MYB, GATA3 related to each other in some tumor pathway and the results were listed in circle graph (Figure 3F). DAVID was used to analyze the KEGG results of genes and we found 9 mRNAs participated in the classic P13-AKT signaling pathway, such as MYB, COL3A1, COL4A1, CSF1R, CDK2, ITGB7, LAMA2, TLR2, and VWF. Combining with the Overall Survival (OS) of LUAD patients, we chose CDK2 as targeted molecule (P < 0.05). KEGG software further was used to identify the upstream molecules of CDK2 in P13-AKT signaling pathway. Finally, the CDKN1A might regulate downstream molecule CDK2 to influence cell cycle progression in LUAD. Through the scatter plot, there was a positive correlation between CDK2 and CDKN1A using GEPIA²² (P = 0, R = 0.59) (Figure 3G).

Through the Ualcan database, we studied the expression of CDK2 in LUAD tissues and normal tissues. CDK2 expression was high in 515 LUAD tissues in comparison with 59 normal tissues (P = 1.624E-12) (Figure 4A). Different ages had differential expression levels of CDK2 (Normal-vs.-Age 61–80, P = 1.044E-2) (Figure 4B). The expression of CDK2 was associated with various clinical features, such as Grade (Normal-vs.-Grade2, P = 5.088808E-03) (Figure 4C), the types of adenocarcinoma (Normal-vs.- Adenocarcinoma, P = 2.195E-02) (Figure 4D), gender (Male-vs.-Female, P = 1.338E-03) (Figure 4E), smoking habits (Normal-vs.-Smoker, P = 1.319E-11) (Figure 4F), stage (Normal-vs.-stage1, P = 3.458E-02) (Figure 4G), TP-53 mutation state (Normal-vs.-TP53 mutation, P = 1.624E-12) (Figure 4H), Weight (Normal-vs.-weight, P = 2.978E-2) (Figure 4I). The expression of CDK2 was related to the survival and prognosis of patients with LUAD (HR = 1.66, P = 5.8e-15) (Figure 4J). The higher the expression of CDK2, the shorter the survival time. Through this database, we discovered the CDK2-related genes and the heat maps were shown in Figures 4K,L. Many molecules involved in tumor classical signaling pathways were related to CDK2 expression. PrognoScan database: A new database for meta-analysis of the prognostic value of genes. We further found the expression of CDK2 influenced the OS and RFS of LUAD patients in different GSE datasets (GSE13213: COX P-value = 0.027245; GSE31210: COX P-value = 0.017966; GSE31210: COX P-value = 0.028818; GSE31210: COX P-value = 0.004197) (Figures 4M–P).

Through the Oncomine database, CDK2 expression was higher in 15 cancer tissues in comparison with normal tissues (Figure 5A). There was a great difference between cancer tissue and normal tissue. The statistical significance between normal and tumor tissues was further found in the TIMER database (the more “^∗” symbol, the greater the difference) (Figure 5B). Mata analysis of 15 published studies on LUAD showed that expression of CDK2 was higher in LUAD (Median Rank = 5384.0, P = 4.16E-6) (Figure 5C). By t-test, box plot and peak plot of CDK2 in LUAD were shown in Figures 5D,E (t-Test = 3.249, Fold Change = 1.332, P = 0.003). Using the R package, we ranked CDK2 by its expression in 33 cancers (Figure 5F). Through different R packages, the expression of CDK2 was relatively higher in 17 cancers (P < 0.05) (Figures 6A–Q). The different expression levels of CDK2 had statistical significance on the stage of patients (P < 0.05) (Figures 6R–Y). There was significant difference in the expression of CDK2 between stage I and stage III, I and IV, II, and IV LUAD (^∗/^∗∗). For different types of cancer, we conducted paired differential expression of CDK2. Simultaneously, we found that there was statistical significance in 14 cancer types (P < 0.05) (Figures 7A–N). In 10 cancer types, high CDK2 expression was associated with poor prognosis of the patients (P < 0.05) (Figures 7O–X).

We compared CDK2 high expression groups with low expression groups in LUAD, and Figures 8A,B showed the top differentially expressed genes in two groups of cancer specimens. We took FP as the abscissa and TP as the ordinate and showed the AUC curve of 1, 3, 5, and 8 years to forecast the survival of patients. The areas under the curve were 0.465, 0.524, 0.612, and 0.575 (Figure 8C). “Clusterprofiler” R package was used to show the function analysis. And CDK2 was related to the DNA replication, regulation of cell cycle, cell cycle checkpoint, P53 signaling pathway. The bubble and bar charts were shown in Figures 8D–G. Wave charts about GO, KEGG, Reactome showed that CDK2 involved in classic tumor signaling pathways, such as regulation of TP53 activity, PTEN regulation, apoptosis, P13K-AKT signaling pathway (Figures 8H–J). The circle graph was shown in Figure 8K. For different immune cells, CDK2 was positively correlated with CD4 T cells (Cor = 0.1679, P = 0.0003), macrophage M1 (Cor = 0.2437, P = 1.40E-07) and negatively correlated with Mast cells (Cor = –0.1545, P = 0.0009) by CIBERSOPT algorithm (Figure 8L).

In LUAD, we further investigated different expression of CDK2 in different immune cell types. We found that 14 immune cells were closely related to CDK2 expression, such as T cells Follicular (P = 0.05, r = 0.09), T cells Regulatory Tregs (P = 0.01, r = –0.13), Macrophages.M0 (P = 0, r = 0.14), Macrophages (P = 0.02, r = 0.11), Eosinophils (P = 0.01, r = 0.13), Mast cells activated (P = 0, r = 0.13), Macrophages M1 (P = 0, r = 0.24), Mast cells (P = 0, r = –0.15), Monocytes (P = 0.02, r = –0.11), Mast cells Resting (P = 0, r = –0.2), Plasma cells (P = 0, r = –0.15), Tcells CD8 (P = 0, r = 0.14), Neutrophils (P = 0, r = 0.15), T cells CD4 (P = 0, r = 0.17) (Figures 9A–N). According to the median value of CDK2 expression, patients with LUAD were divided into the high expression groups and low expression groups. In different groups, the expression of 10 immune cells had systematic differences (P < 0.05) (Figures 9O–X).

Through the TISIDB database, we found CDK2-related immunomodulators and the heat maps of immunopotentiators and immunosuppressants were shown in Figures 10A,B. By sorting the P-value (P < 0.05), we identified 13 immunosuppressants (ADORA2A, BTLA, CD274, CSF1R, IL10, KDR, LAG3, LGALS9, PDCD1, PDCD1LG2, TGFB1, TIGIT, VTCN1) and 21 immunopotentiators (CD27, CD276, CD28, CD40LG, CD48, CD70, CD80, CXCL12, CXCR4, ENTPD1, HHLA2, ICOSLG, IL2RA, IL6, IL6R, KLRC1, MICB, PVR, TMEM173, TNFRSF13B, ULBP1) with a high correlation of CDK2. In the cBioProtal, we explored forty-nine genes associated with 34 immunomodulators. Clinical data and gene expression data were downloaded from TCGA database. Forty-nine genes were mixed from the TCGA using “perl” and “R” package (Figure 10C). Through the Metascape database, the most abundant pathway of 49 genes was the immune system process, other function terms were biological adhesion, biological regulation, cell proliferation (Figure 10D). The BP, CC, MF was showed specifically in Table 2. The protein network interaction (PPI) of these 49 immune related genes were shown in the Figure 10E (P < 1.0e-16). GSEA analyzed the function, ES, NES, NOM p-val, FDR q-val of 49 genes. The statistically significant items were PUJANA_ATM_PCC_NETWORK and INTRACELLULAR_SIGNAL_TRANSDUCTION (P < 0.01) (Figure 10F and Table 3). Combing with clinical data and expression matrix, univariate and multivariate regression analyses were performed on age, gender and tumor stage (Total P = 1.399e-08) (Table 4). The forest map showed HR and concordance index (0.7) of LUAD. As predicted, the tumor stage was an independent risk factor for the prognosis of patients with LUAD (Figure 10G). Nomogram model combined several clinical factors, such as age, gender, stage, T, N, M to intuitively analyze the prognosis of LUAD (Figure 10H). Each patient can be assessed by nomogram model based on baseline clinical data. Univariate regression analysis showed that 36 genes were associated with the prognosis of LUAD (P < 0.05). We showed the HR, HR95L, HR95H, P-value of these 36 significant genes. CDK2 was an independent prognostic gene in LUAD (HR > 1) (Table 5 and Figure 10I). Multivariate regression analysis was used to analyze the 36 genes and only four genes (SIT1, SNAI3, ASB2, and CDK2) were included in the prediction model (Table 6). Through the GEPIA database, SIT1 (P = 0.04), SNAI3 (P = 0.00035), ASB2 (P = 0.00027) were lower expressed in LUAD (P < 0.05) (Figures 10J–L). To verify the prognostic model, the risk curve showed that the high-risk group was more lethal to lung cancer patients (P = 6.223E-04) (Figure 10M). Based on different risk scores, patients were divided into high and low risk groups using “R” package (Figures 10N,O). The ROC curve showed the different AUC of 4 interesting genes and CDK2 had more predictive value in the prognostic model (AUC = 0.579) (Figure 10P).

Immunohistochemical result showed that the CDK2 expression in LUAD tissues was significantly higher than that in normal lung tissues (Figures 11A,B). We selected 455 miRNAs related to CDK2 from TargetScan database. 10,018 miRNAs were found in miRWalk database. And we explored 76 miRNAs from mirDB and 28 miRNAs from Starbase. Combining with differentially expressed miRNAs in LUAD, 7 miRNAs were joined in the network using Venn map, such as hsa-miR-302b-3p, hsa-miR-372-3p, hsa-miR-302a-3p, hsa-miR-373-3p, hsa-miR-520a-3p, hsa-miR-520d-3p, and hsa-miR-302d-3p (Figure 11C). These miRNAs had obvious influence on the prognosis of LUAD patients (P < 0.05) (Figures 11D–J). We used the Starbase database to find the potential lncRNAs that regulate seven miRNAs. These coding genes and non-coding genes were interacted with each other by using Cytoscape software (Figure 11L). According to the degree in the network, we screened out the top 15 genes using CytoHubba (six lncRNAs: XIST, SNHG16, RP11-145M9.4, MAP3K14, MIR4720, and RP11-379K17.11) (Figure 11K).

Compared with the expression of CDK2 in BEAS-2B cell line, the expression of CDK2 in A549 cell line was increased, but there was no statistical difference (P = 0.0808). There were significant differences between H1299 cell line (P = 0.0006) and H1975 cell line (P = 0.0030) (Figure 11M).

A total of 192 anti-tumor drugs were included in the study. The IC50 level of 89 anti-tumor drugs were related to the expression of CDK2. According to the size of P value (P < 0.05), we screened out the top 20 anti-tumor drugs with positive or negative correlation, such as Camptothecin (r = –0.074, P = 0.000018), Vinblastine (r = –0.085, P = 0.0000243), Cisplatin (r = –0.099, P = 0.0000843), Cytarabine (r = –0.0746, P = 0.0000975), Navitoclax (r = –0.106, P = 0.000158), Vorinostat (r = –0.113, P = 0.0002), Nilotinib (r = –0.127, P = 0.000258), Olaparib (r = –0.0966, P = 0.000302), Axitinib (r = 0.343, P = 0.000381), AZD7762 (r = –0.0902, P = 0.000382), SB216763 (r = 0.284, P = 0.000403), KU-55933 (r = 0.315, P = 0.000404), PLX-4720 (r = –0.0857, P = 0.000428), Wee1 Inhibitor (r = –0.141, P = 0.000780), PD173074 (r = –0.0963, P = 0.000794), Obatoclax Mesylate (r = –0.08009, P = 0.000872), Sorafenib (r = –0.076, P = 0.0009), Irinotecan (r = –0.0801, P = 0.00120), BMS-536924 (r = 0.0765, P = 0.0012), and GSK1904529A (r = –0.113, P = 0.0012) (Table 7 and Figures 12A–T).