The data of this study, including simple nucleotide variation, copy number variation (CNV), RNA-seq gene expression, miRNA expression, survival, and clinical information, were retrieved from The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/). 368 patients with intact information of simple nucleotide variation, CNV, RNA-seq gene expression, miRNA expression, survival, and clinical features were included in this study for further analysis in R software (version 4.0.3).

We also collected 34 LUAD samples from patients who underwent lobectomy and systematic lymph node resection at the Department of Thoracic Surgery, Zhongshan Hospital, Fudan University between 2016 and 2017. All pulmonary resections were performed by professional thoracic surgeons in our institution, and the diagnoses of LUAD were confirmed by at least two qualified pathologists. All patients have signed informed consent to conduct genomic studies consistent with the ethical principles of the Declaration of Helsinki. Patients with a history of chemotherapy, radiotherapy, and immunotherapy or who had evidence of metastasis were excluded. RNA sequencing for all tumor samples was performed using Illumina Hiseq 2500 and BGI-500RNAseq platforms ( Supplementary Material 1 ). Besides, CDKN2A, CCND1, CDK4, CCNE1, and RB1 mutations in our samples were identified by whole exon sequencing. GATK4 software was used to detect single nucleotide polymorphisms (SNPs), insertions, and deletions (Indels). Only SNPs and Indels with quality/depth ratio >=2.0 were considered for subsequent analysis and were annotated using Annovar software. The study was approved by the ethical committees of Zhongshan Hospital (No. 201986122).

Patients who harbored any of the designated cell cycle progression pathway-related gene mutation (CDKN2A, CCND1, CDK4, CCNE1, and RB1) were defined as the mutation group, and the remaining were classified as the wild group. Clinical characteristics including age, gender, race, anatomic location of the lesion, smoking history, and TNM stage were compared between the two groups. Next, Kaplan-Meier survival analysis was performed, and log-rank test was used to show the difference between the two groups.

The simple nucleotide variation data were stratified into two groups, and their mutational patterns were investigated separately. DMGs were detected using Fisher exact test and visualized by maftools package. Besides, mutational load and CNV were calculated for every patient, and their differences between the two groups were explored.

After removing the genes whose average expression levels in all patients were less than 10, we used edgeR package to explore the protein-coding DEGs based on the RNA-seq gene expression data (count format). The absolute value of Log₂FoldChange (|Log FC |) >1 and p-value <0.05 were considered as significantly different.

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed based on the DEGs using org.Hs.eg.db and Clusterprofiler package. Subsequently, we performed gene set enrichment analysis (GSEA) based on the 16684 protein-coding genes. Adjusted p-value less than 0.05 was considered as significantly different. Finally, we used gene set variation analysis (GSVA) to investigate the difference of cell-cycle-related events ( Supplementary material 2 ) between the two groups. Rather than ssGSEA, we harnessed the GSVA method in this study because GSVA includes the normalization of gene expression to reduce the noise of the data and has been shown to outperform ssGSEA when measuring the signal-to-noise ratio in differential gene expression and differential pathway activity identification analyses (13).

To further investigate the inner correlation of the DEGs, we mapped the DEGs in STRING (14) (https://cn.string-db.org, version 11.0b) to build a PPI network. Interaction with a confidence score greater than 0.4 was included in the network. Subsequently, the Molecular Complex Detection (MCODE) plugin in Cytoscape (https://cytoscape.org, version 3.8.2) was harnessed to identify modules of the PPI network (degree cut-off=2, node score cut-off=0.2, k-score=2, and max depth=100).

Differentially expressed miRNAs were detected by limma package using the miRNA expression data (TPM normalized). The criteria for defining differentially expressed miRNAs were set as follows: |Log FC| >0.5 and p-value <0.05. Next, we investigated differentially expressed long non-coding RNAs (lncRNAs) between the two groups, and the thresholds for |Log FC| and p-value were the same as the protein-coding DEGs.

After matching the differentially expressed miRNAs with DEGs and differentially expressed lncRNAs using miRWalk (http://mirwalk.umm.uni-heidelberg.de) and miRcode (http://www.mircode.org) webtools, a ceRNA network of DEGs - differentially expressed miRNAs -differentially expressed lncRNAs was constructed in Cytoscope.

The profiling of 22 immune cell types was analyzed using RNA-seq gene expression data (FPKM format) and the leucocyte signature matrix (LM22) in CIBERSORT (https://cibersort.stanford.edu). LM22, composed of 547 genes, could distinguish 22 human hematopoietic cells, including various subtypes of T cells, B cells, plasma cells, NK cells, and myeloid cells (15). Subsequently, we compared the immune and stromal scores between the two groups by estimate package. At last, we detected the expression levels of 15 immune checkpoints and 20 co-stimulators between the two groups. The relationships between these immune molecules were detected by Spearman correlation analysis and visualized by ggcorrplot package.

After taking the intersection of the DEGs and 1872 cell-cycle-related genes ( Supplementary material 3 ), we obtained 34 cell-cycle-related DEGs (ccDEGs). Based on the ccDEGs, we performed least absolute shrinkage and selection operator (LASSO) and subsequent binary logistic regression analysis to build a model which could predict if a LUAD patient harbored specific cell cycle progression pathway-related gene (CDKN2A, CCND1, CDK4, CCNE1, and RB1) mutation. The LASSO regression analysis is a penalized method to select data with high dimensions and reduce the impact of overfitting (16, 17). Ten-fold cross-validation was adopted using the glmnet package to determine the optimal parameter λ and corresponding genes. The cut-off value of this model was determined by the receiver operating characteristic curve (ROC). Next, we explored the effect of the candidate genes on proliferation. The CERES dependency scores of these candidate genes in 51 LUAD cell lines were retrieved from the Cancer Dependency Map (DEPMAP, https://depmap.org/portal/), which represented the effect on cell viability after knocking out corresponding genes by CRISPR-Cas9 genetic perturbation reagents (18). A score of zero indicated that a gene was not essential. Scores less than zero indicated knockout of the corresponding genes could inhibit cell proliferation; the smaller the score, the more pronounced the effect. Scores greater than zero showed the opposite effect.

All statistical analyses were conducted in R software (version 4.0.3). The comparison of the clinical features of the LUAD patients between the two groups was carried out, of which categorical variables were compared by Chi-square test or Fisher exact test when appropriate and continuous variables were compared by Student’s t-test. The Student’s t-test was also used to compare continuous variables such as the mutational load, immune and stromal scores. Log-rank test was used to compare overall survival between the two groups in Kaplan-Meier survival analysis. All the tests used in our study were two-sided, and the significance threshold of p-value was set as 0.05.

The mutational rates of RB1, CDKN2A, CCNE1, CDK4, and CCND1 in all LUAD patients included in our study were 8.2%, 3.5%, 1.4%, 0.5%, and 0.3%, respectively. After dividing the patients into mutation (n=46) and wild (n=322) groups, we found that more than half of the patients in the mutation group harbored RB1 mutation (65.2%, the common variants in lung cancer were C706F, G748K, and R661W), followed by CDKN2A (28.3%, the common variants in LUAD were R80Q and H83Y), CCNE1 (10.9%), CDK4 (4.3%, the common variants in LUAD were R24L and R24C) and CCND1 (2.2%). Besides, 5 of 46 patients in the mutation group harbored more than one mutant gene mentioned above. Comparison of the baseline information showed no significant differences between the two groups in age, gender, anatomical location of the lesion, smoking history, and TNM stage ( Table 1 ). Next, we performed the Kaplan-Meier survival analysis between the two groups, but log-rank test showed no significance ( Figure 1 ). As the number of patients in the mutation group was much smaller than the wild group, this result needs to be further verified in larger cohorts.

Different mutational patterns were detected between the mutation and wild groups ( Figure S1 ). The top 20 mutant genes in the mutation and wild groups were shown separately ( Figures 2A, B ), among which fourteen ones were shared by the two groups, and all of them were observed with higher mutational rates in the mutation group. Next, we used Fisher exact test to detect the DMGs between the two groups. 286 genes were defined as DMGs ( Figure 2C , Supplementary material 4 ), 30 of which such as PDE3A (22% mutation group vs. 3% wild group, P<0.0001) and TP53 (76% mutation group vs. 46% wild group, P<0.0001) were related to cell cycle events.

At last, we compared the mutational load and CNV between the two groups. The mutation group had a higher level of mutational load ( Figure 2D ) and CNV ( Figure 2E ), which indicated that the mutation group harbored more mutational events and thus is more likely to benefit from immunotherapies.

Among the 16684 protein-coding genes analyzed, 302 were upregulated and 354 were downregulated in the mutation group ( Figure 3A , Supplementary Material 5) . Among the DEGs, VGLL2 (logFC=6.41, P<0.0001), GCG (logFC=5.64, P<0.0001), and SOST (logFC=4.04, P<0.0001) were the most significantly upregulated genes; DEFA5 (logFC= -8.11, P<0.0001), PRB4 (logFC= -6.87, P<0.0001), and SPAG11B (logFC= -6.66, P<0.0001) were the most significantly downregulated genes.

To explore the biological functions affected by DEGs, we performed the GO and KEGG analysis. In GO analysis ( Figures 3B–D ), the top 3 of the 630 significantly different biological processes were cornification, antimicrobial humoral response, and signal release. KEGG analysis showed the DEGs were related to fat digestion and absorption, drug metabolism − cytochrome P450, and retinol metabolism ( Figure 3E ).

Apart from GO and KEGG analysis, we further investigated the influence of cell cycle progression pathway status using GSEA and GSVA analysis. Processes like cell cycle, DNA replication, and ribosome were significantly different between the two groups by GSEA analysis ( Figures 4A–C ). GSVA analysis showed cell-cycle-related events such as DNA replication, E2F targets, mitotic spindle, and G2M checkpoints were significantly different between the two groups ( Figures 4D–I ).

PPI clusters were investigated by MCODE (19) in Cytoscape. The top four modules (clustering scores were 9.391, 7.000, 6.526, and 4.414, respectively) are shown in Figures 5A–D . Next, we carried out the GO enrichment analysis of the four clusters separately. The result demonstrated that these clusters were related to biological processes such as cell-cell signaling, cornification, arachidonic acid secretion, and nephrogenic mesenchyme development ( Figure 5E ).

Both miRNAs and lncRNAs could act as important modulators which may influence the expression and function of mRNAs. Here, we analyzed the differentially expressed miRNAs and lncRNAs between the two groups ( Figure S2 , Supplementary material 5) . The results showed that 31 miRNAs were upregulated in the mutation group, the most significantly different ones of which were miR−3689e (logFC=2.30, P<0.01), miR−3689f (logFC=1.95, P<0.01), and miR−1197 (logFC=1.70, P<0.001). 5 miRNAs were downregulated in the mutation group, and the most significantly different one was miR−203b−3p (logFC=-0.86, P<0.001). For lncRNAs, LINC01305 (logFC=3.21, P<0.0001) was the most significantly different one among 31 upregulated lncRNAs in the mutation group; AF003626.1 (logFC=-4.63, P<0.0001) was the most significantly different one among 64 downregulated lncRNAs.

After matching DEGs and differentially expressed lncRNAs with differentially expressed miRNAs, we built a ceRNA network ( Figure 6 , S3 , Supplementary material 6 ) to elucidate their relationships (Here, only genes negatively regulated by miRNAs were shown in this network). We found that both miR-3131 and miR-185-3p could regulate more than 40 genes, suggesting the two miRNAs and the genes they regulate might play a significant role in the differences between the mutation and wild group.

Tumor immune cell infiltration refers to the migration of immune cells from the peripheral blood to the tumor tissue, where they exert their function (20). Different immune cell infiltration patterns were detected between the two groups ( Figure 7A , S4 )using CIBERSORT (15) (permutations set =1000). Among all the significantly different cell types, we found that the proportion of naïve B cells (P<0.05) and M1 macrophages (P<0.01) were higher in the mutation group, while CD4⁺ memory resting T cells (P<0.01), monocytes (P<0.001), activated mast cells (P<0.05), resting (P<0.05) and activated (P<0.01) dendritic cells (DCs) were higher in the wild group. Next, we further explored the tumor purity between the mutation and wild groups, while we failed to find any significant difference between the two groups ( Figures 7B, C ).

Then we analyzed the expression levels of 15 immune checkpoints and 20 co-stimulators between the two groups ( Figures 7D, E ). Among co-stimulators, TNFSF13 showed a significant difference (P<0.05) and was higher in the wild groups. As for immune checkpoints, LAG3 was significantly higher (P<0.05) in the mutation group. The Spearman correlation test demonstrated that most of these immune molecules were positively correlated ( Figure 7F ).

The optimal parameter λ of LASSO analysis was set as the smallest partial likelihood deviance ( Figure S5A, B ). A panel of fifteen ccDEGs including SOX2, IGF2, CDKN2A, EREG, BMP7, IL1A, HEPACAM2, MAPK4, SOX11, FBXO43, TDRD12, INSC, MOV10L1, PIWIL2, and NLRP5 were employed ( Figure S6 ). Subsequently, we conducted the binary Logistic regression analysis to get the coefficients of each item, and the predictive model could be demonstrated as the following formula: score = 0.042*expression level of SOX2+ (-0.2086)*expression level of IGF2+ 0.6522*expression level of CDKN2A+ (-0.1066)*expression level of EREG+ 0.1315*expression level of BMP7+ (-0.3672)*expression level of IL1A+ 0.0585*expression level of HEPACAM2+ 0.1757*expression level of MAPK4+ 0.3392*expression level of SOX11+ 0.7282*expression level of FBXO43+ 0.2677*expression level of TDRD12+ (-1.9734)*expression level of INSC+ 0.1685*expression level of MOV10L1+ 0.7578*expression level of PIWIL2+ (-1.2477)*expression level of NLRP5 ( Figure S5C ). ROC demonstrated the cut-off value of this model was -2.466, and the area under the curve (AUC) is 0.845 ( Figure S5 D). When the score is less than -2.466, the LUAD patient is more likely to harbor CDKN2A/CCND1/CDK4/CCNE1/RB1 mutations. In the validation cohort from our institution, 28 were found to harbor CDKN2A/CCND1/CDK4/CCNE1/RB1 mutations by whole exon sequencing. After testing the model in our samples, we found the sensitivity and specificity of this model were 0.667 and 0.929, separately.

After obtaining the CERES scores of the fifteen ccDEGs ( Supplementary material 7 ), we found the CERES scores of five ccDEGs (BMP7, IGF2, MAPK4, TDRD12, and NLRP5) were greater than zero in more than 3/4 LUAD cell lines, suggesting they’re more likely to exert an anti-proliferation effect in LUAD patients. On the contrary, four ccDEGs (FBXO43, IL1A, INSC, and SOX2) were thought to promote proliferation, given that their CERES scores were less than zero in more than 3/4 LUAD cell lines.