LUAD HTSeq transcriptome profiles were acquired from TCGA data portal (http://cancergenome.nih.gov/, Data Release 9.0) (Tomczak et al., 2015). The GENCODE V22 (GRCh38) (Frankish et al., 2019) human genome annotation file was used to annotate mRNA and lncRNA transcripts. The gene types of “processed_transcript,” “lincRNA,” “3prime_overlapping_ncrna,” “antisense,” “non_coding,” “sense_intronic,” and “sense_overlapping” with length greater than 200 nt were defined as lncRNAs. The transcripts per kilobase million (TPM) value was converted from fragments per kilobase million (FPKM) according to the formula: T P M = F P K M ∑ a l l   g e n e s F P K M ∗ 10 6 (Li and Dewey, 2011). TPM transcriptome profiles were further normalized by log2(x+1) transformed, and the genes covering >30% missing values were filtered. We adopted clinic information from the work of Liu et al., which provided the usage recommendations for corresponding TCGA LUAD cohort in the meanwhile (Liu et al., 2018). Gene expression profiles and clinical information from two GEO datasets were collected as the external validation datasets (Supplementary Table S1).

Experimentally validated miRNA–mRNA interaction information was collected from miRTarBase 7.0 (Chou et al., 2018), miRecords 2013 (Xiao et al., 2009), and TarBase 6.0 (Vergoulis et al., 2012). A total of 388,895 unique interactions that contained 2,846 miRNAs and 18,936 mRNAs were obtained. We integrated starBase 2.0 (Li J.-H. et al., 2014), DIANA-LncBase V2 (Paraskevopoulou et al., 2016), and lncRNASNP2 (Miao et al., 2018) databases and obtained 10,318 experimentally validated nonredundant miRNA–lncRNA interactions involving 290 miRNAs and 1,162 lncRNAs.

Chromosomal instability (CIN) risk scores of samples were calculated by summing of the CIN70 genes’ TPM value, whose expression level was consistently correlated with total functional aneuploidy in several cancer types (Carter et al., 2006). A high CIN score indicates a high extent of chromosomal instability.

1) CeRNA Identification in LUAD CIN-Low and CIN-High Patients

Differentially expressed genes between tumor and normal adjacent samples were reserved for identifying ceRNA pairs in LUAD CIN-Low and CIN-High patients (Zhou et al., 2014; Xu et al., 2015). First, we screened mRNA–lncRNA pairs which shared significantly common upstream miRNA regulation according to the hypergeometric test FDR-adjusted p value <0.05. The p value was calculated as 1 − ∑ i = 0 r − 1 ( K i ) ( N − K M − i ) ( N M ) , where N was the number of background human miRNAs, M and K represented the number of miRNAs regulating the current mRNA and lncRNA, respectively, and r was the common miRNA number between the two RNAs. Also, Pearson correlation coefficients (PCC) > 0.1 as well as correlation FDR-adjusted p value <0.05 were considered significantly co-expressed pairs for having a similar expression trend. For CIN-Low and CIN-High LUAD patients’ groups separately, mRNA–lncRNA pairs that met the aforementioned conditions of co-regulated and co-expressed RNAs were recognized as ceRNAs which served as nodes in the dysregulated network.

2) CIN-Related Dysregulated CeRNA Network

In the dysregulated ceRNA network, alteration of ceRNA PCC between the CIN-High group and CIN-Low group in the LUAD cohort was assigned to measure the weight of edges and the node weight quantifying the prognostic role of the gene for LUAD DSS (Li et al., 2012; Wang et al., 2015). Given a ceRNA pair, r C I N _ H i g h and r C I N _ L o w represented the PCC in two groups, and n C I N _ H i g h and n C I N _ L o w represented the sample size. If the ceRNA identification conditions were not satisfied in one of the CIN-Low and CIN-High groups, the PCC was set to 0. The edge weight (E) was calculated as φ − 1 [ 1 − 2 ∗ ( 1 − φ ( | F ( r C I N _ H i g h ) − F ( r C I N _ L o w ) 1 n C I N _ H i g h + 3 + 1 n C I N _ L o w + 3 | ) ) ] , where F ( x ) = 1 2 l n 1 + x 1 − x , and φ denoted the standard normal distribution function. Meanwhile, the node weight (N) of the network was defined as φ − 1 ( 1 − p ) , where p was the univariate Cox proportional hazards analysis significant p value, adopting the samples without postoperative treatment for purpose of eliminating its effect on clinical outcome.

1) Searching for Dense Modules in Dysregulated CeRNA Network

The objective of building a dysregulated network was to identify dense modules with relatively high weight of node and edge. We utilized the algorithm from the dmGWAS_3.0 R package (Wang et al., 2015) to define the module score and conduct module mining. The module score was λ ∑ E n E + ( 1 − λ ) ∑ N n N , where n E and n N were the number of edges and nodes in the network, and λ was between 0 and 1. The method of module search are summarized into the following three parts: 1) For each gene in the network, it was assigned as a seed module at the beginning; 2) add a neighbor node whose shortest path to any node in the module was shorter or equal to d into the module once at a time only when the increment of the current module score was greater than r multiply by the score of the previous module; 3) repeat these steps till no more neighbors could be added. The parameter d was decided as 1, r = 0.1, and λ was set to 0.44.

2) Screening and Integrating Dense Modules

In order to seek valuable messages from a large number of dense modules, we prioritized network modules by following three qualification rules. First, the top 10% modules with the highest scores were selected for further analysis. Second, a random perturbation strategy was used to screen out the non-randomized prognostic modules by multivariate Cox regression analysis for DSS. For each given module, random multivariate Cox regression p values were calculated after random acquisition of the same number of genes in the background 10,000 times. Then, sorting the 10,000 p values in ascending order and only the modules ranked in top 1% (random p value <0.01) were left. Furthermore, we merged the modules via genes in common and acquired the subnetwork named CIN-related dysregulated ceRNAs (CRDC).

The CRDC risk score of each LUAD patient was constructed by a linear combination of gene expression values and log-transformed hazard ratios (HR) from the univariate Cox regression model, in which the univariate Cox regression analysis was conducted using the DSS without undergoing postoperative treatment. The CRDC risk score was calculated as ∑ i = 1 N E x p r e s s i o n _ i *log ( H R _ i ) , where N was the number of the CRDC genes.

The LUAD patients’ immune score, stromal score, and estimate score were calculated by the ESTIMATE R package (version 1.0.13); meanwhile, tumor purity was obtained by cos(0.6049872018 + 0.0001467884 × ESTIMATE score) (Yoshihara et al., 2013). Tumor mutation load (TMB) was performed as non-synonymous mutation counts in coding sequences (CDSs) divided by total number of CDSs then multiplied by 10⁶. The geometric mean of GZMA and PRF1 TPM expression values estimated were defined as cytolytic activity (CYT), which reflected the local immune cytolytic T-cell activity (Rooney et al., 2015). TIMER2.0 (Li T. et al., 2020) was applied to estimate immune cell infiltration levels in the tumor immune microenvironment for LUAD patients. The tumor immune dysfunction and exclusion (TIDE) score, a signature consistent with tumor immune evasion, which could predict immune checkpoint blockade (ICB) response, was performed by TIDE command-line interface with tumor TPM expression normalized by adjacent normal samples (Jiang et al., 2018; Fu et al., 2020). The immunophenoscore (IPS), a superior predictor of response to anti-CTLA4 and anti-PD-1 antibodies, was used to evaluate tumor immunogenicity (Charoentong et al., 2017). Also, we predicted the possible response to immunotherapy by the SubMap method from GenePattern (https://cloud.genepattern.org/gp/pages/index.jsf).

Differentially expressed genes (DEGs) were identified by the DESeq2 R package (Love et al., 2014) with the threshold of |log2FC| > 1 and FDR <0.01. CIN-High and CIN-Low as well as CRDC-High and CRDC-Low patients were divided by the StepMiner method, which could find an appropriate threshold for the dichotomizing numeric vector (Sahoo et al., 2007; Liu et al., 2013). For example, the CIN scores were sorted in ascending order and were used for finding an optimal location t which maximized the signal-to-noise ratio (SNR). S N R = ∑ i = 1 n ( μ 1 I ( i ≤ t ) + μ 2 I ( i > t ) − μ ) 2 ∑ i = 1 n ( μ 1 I ( i ≤ t ) + μ 2 I ( i > t ) − C I N i ) 2 , where I is the indicator function and is set to 0 when the condition in brackets is not satisfied, μ is the mean CIN score for all patients, and μ 1 and μ 2 represent the mean CIN score of two groups separated by the location t. Then, the CIN score of position t was regarded as the classification threshold. The conventional ROC curve and time-dependent ROC curve estimation was performed using the pROC (version 1.18.0) (Robin et al., 2011) and survivalROC (version 1.0.3) (Heagerty and Zheng, 2005) R package separately. The Kaplan–Meier survival curves between different groups divided by the CIN score or CRDC score was tested by the log-rank test using the R package survminer (version 0.4.9). The multivariable Cox proportional hazard model, stepwise regression method, and nomogram were constructed by the rms R package (version 6.2–0). The decision curve was built to evaluate benefit of CRDC at different prognostic time points with the DCA R package (version 2.0). Biological progress (BP) enrichment analysis and visualization were implemented by Metascape (Zhou Y. et al., 2019), clusterProfiler R package (version 4.0.5) (Wu et al., 2021), and fgsea R package (version 1.20.0). The hallmark, KEGG pathway, and immunologic gene sets enrichment analysis were performed by the clusterProfiler R package. BP, hallmark, and immunologic gene sets were obtained from MSigDB v7.4 (Subramanian et al., 2005; Liberzon et al., 2015). The R version adopted in this article was 4.0.3. Representation of significance of p value in figures were described as follows: *: p < 0.05, **: p < 0.01, ***: p < 0.001, and ns for nonsignificant.

For quantifying chromosomal instability of samples, we calculated CIN scores of all samples by summing up genes’ expression in CIN70 signature. With a paired-sample comparison analysis in 57 LUAD patients, CIN scores of lung adenocarcinoma samples (350.44 ± 60.00) were significantly higher than those of paired adjacent normal samples (214.1201 ± 22.98) (Figure 2A, Wilcoxon signed-rank test, p = 5.3e-11) as well as all samples in TCGA LUAD data cohort (Supplementary Figures S1A,B). In 513 LUAD and 59 adjacent normal samples, receiver-operating characteristic (ROC) analysis showed that the CIN score had a high performance in discriminating LUAD patients from adjacent normal samples with the area under the curve (AUC) reached up to 0.989 (Figure 2B). Meanwhile, we found that high CIN scores were related with the worse clinical stage and TNM stage generally (Figure 2C). Also, patients with high CIN scores had worse overall survival (OS), DSS, disease-free interval (DFI), and progression-free interval (PFI) (Supplementary Figures S1C–J). This indicated that LUAD samples possessed higher CIN than adjacent normal samples, which is consistent with the general characteristics of genomic instability in carcinoma (Grady and Carethers, 2008; Martin et al., 2010).

We identified significantly DEGs in LUAD samples by the DESeq2 R package. Within 3,400 upregulated and 1,975 downregulated mRNAs, 52 genes of the CIN70 signature were upregulated, and none of them were downregulated. Similarly, 42 of CIN70 genes showed upregulation in 289 CIN-High patients compared to the CIN-Low group that contained 224 patients. It was observed that these highly altered expressed CIN70 genes in two comparison groups possessed greater extent differential expression than other genes (Figure 2D). Functional enrichment analysis indicated that the DEGs of these two comparable groups shared common biological processes and pathways. For instance, both downregulated gene sets participated in microtubule-based movement and microtubule bundle formation biological progress as well as neuroactive ligand–receptor interaction and cAMP signaling pathways. On the contrary, the upregulated genes enriched nuclear division, organelle fission, chromosomal segregation, mitotic nuclear division, and meiotic cell cycle process, in addition to cell cycle and nicotine addiction pathways (Figures 2E,F).

Due to the complex indirect regulation of lncRNAs for target mRNAs, we constructed a CIN-related dysregulated ceRNA network. The mRNAs and lncRNAs markedly changed on the expression level between LUAD versus adjacent normal samples were considered LUAD-associated genes and used to identify ceRNA pairs in CIN-Low and CIN-High samples, respectively. The ceRNA network contained 2,195 edges and 908 nodes composed of 787 mRNAs and 121 lncRNAs. We further constructed the double-weighted CIN-related dysregulated ceRNA network, in which the edge weight and node weight indicated the extent of dysregulation of ceRNA in CIN-High samples contrast with CIN-Low and the quantized value for RNA’s prognostic role in DSS, respectively (Figure 3A, Supplementary Table S2, details were prescribed in Materials and Methods). The network had a typical biological network property, a scale-free characteristic, and followed power-law distribution with R² being 0.80 (Figure 3B). Network topological analysis suggested that lncRNAs may play an important role in regulating mRNA expression indirectly by the ceRNA crosstalk mechanism to exert biological functions in spite of low coding potential (Figure 3C and Supplementary Figures S2A–C).

In this work, the dysregulation type of ceRNA was defined as gain (PCC in the CIN-High group were higher than that in CIN-Low) or loss (PCC in CIN-High patients was smaller than that in CIN-Low). Among 2,195 ceRNAs, 828 of them were type of gain; however, the majority were loss (accounted for 62.3%), which suggested that most of ceRNAs no longer maintained previous correlation in the CIN-High group. PCC of gain pairs raised from 0.071 ± 0.115 to 0.224 ± 0.078 and in loss pairs reduced from 0.262 ± 0.097 to 0.090 ± 0.124. The extent of alteration in loss pairs had a greater alteration extent than that in the gain pairs (Supplementary Figures S2D,E and Figure 3D). BP functions were enriched by Metascape, in which gain and loss ceRNAs showed relevance to mitotic cell cycle, tube morphogenesis, regulation of cell cycle process, negative regulation of cell differentiation, epithelial cell differentiation, regulation of kinase activity, and cell adhesion, especially immune system development (Supplementary Figure S2F). It is implied that CIN-related dysregulated ceRNAs were involved with functions related to tumorigenesis and immunological functions.

We wish to seek a collection of ceRNAs possessing a high prognostic value and a great level of variation in samples with high and low CIN; therefore, we adopted the method of dmGWAS_3.0. The method balanced the weight of nodes and edges for scoring the modules and implemented a greedy algorithm to search for dense modules. The parameters were set as follows: λ for balancing weight of nodes and edges was set to “default” calculated as 0.44; the new nodes search range shortest path d = 1; r was decided as 0.1 which decides the magnitude of increment. A total of 651 modules were obtained by this method. These modules covered 776 nodes in the network, which demonstrated that these modules had the characteristic of high repeat coverage. Therefore, we screened and integrated these modules. First, 65 modules contained 83 non-repetitive RNAs with the top 10% highest module scores were selected, which displayed significant overlap with CIN70 genes on account of having four common mRNAs (HDGF, OIP5, CDK1 and CEP55, hypergeometric enrichment test p = 1.86e-05). Second, a random perturbation strategy was used to screen out the non-randomized prognostic modules, and only four modules with random p < 0.01 were left. It has been observed that they shared common RNAs, indicating that these modules can be linked to be a subnet in the dysregulated network. Therewith, the subnet consisting of these modules was extracted containing eight mRNAs (AMOTL1, EFNB2, FGF2, FURIN, CCND2, IFNG, ITGB4, and RHOV) and 4 lncRNAs (LINC00473, LINC00707, MIR497HG, and RP11-16E12.2) (Figure 3E). We further named the subnet as CRDC (CIN-related dysregulated ceRNAs).

The CRDC risk score for each sample was calculated as cumulative sum of expression values of 12 RNAs multiplied by log-transformed hazard ratios from univariate Cox regression analysis in non-postoperative treatment patients (the univariate Cox regression results of RNAs are shown in Figure 4A and Supplementary Table S3). Most of the CRDC genes tended to be risk factors of LUAD prognosis, except for MIR497HG, CCND2, RP11-16E12.2, and IFNG. Studies have shown that CRDC genes were certified related with tumor suppression or progression which is consistent with our research (Supplementary Table S4). In addition, the protein expression levels of partial CRDC genes were further explored by representative immunohistochemistry (IHC) images from the Human Protein Atlas (HPA) database (Sjostedt et al., 2020). In total, five CRDC proteins, FURIN, ITGB4, RHOV, EFNB2, and AMOTL1, were not detected in normal lung tissues, however, expressed in lung cancer samples. Also, IFNG had a reverse situation (Supplementary Figure S3).

The CRDC risk scores in LUAD were significantly higher than adjacent normal in all (Wilcoxon signed-rank test, p = 1.1e-15, Supplementary Figure S4A) and paired samples (Wilcoxon signed-rank test, p = 8e-9, Figure 4B), and it can distinguish the two types of samples clearly (AUC = 0.835, Figure 4C). The higher CRDC score related to the worse clinical stage (Kruskal–Wallis test, p = 0.0012, Supplementary Figure S4B). Spearman rank correlation analysis was performed on the CRDC risk score, CIN score, clinical stage, and smoking exposure in pack-years. There was no significant correlation between smoking and CRDC risk score or clinical stage (p > 0.05), but a margin significant positive correlation with CIN risk scores, which was consistent with the findings that smoking was probably associated with chromosomal instability (Saletta et al., 2007; Li E. et al., 2014). What was noteworthy is that the CIN score, clinical stage, and CRDC risk score were significantly positively correlated, which indicated that the higher the CRDC risk score, the higher the chromosomal instability and the worse clinical stage (Figure 4D). Also, the high CRDC score was associated with poor OS, DSS, and PFI in all LUAD patients (Figures 4E–H). Because postoperative treatment may have some effect on clinical outcomes which was taken into account in the assessment of the impact of CRDC on prognosis. Among samples without postoperative treatment, the OS, DSS, DFI, and PFI of the CRDC high risk group were significantly worse than those of the low risk group (Figures 4I–L). These results have shown that LUAD maintained higher CRDC scores than adjacent normal tissues, such as CIN. Also, the higher CRDC score was concordance with higher CIN and worse clinical stage. The impact of CRDC on prognosis was rarely affected by the postoperative treatment in the TCGA LUAD cohort.

In this part, prognostic independence and robustness of the CRDC score will be tested. In the TCGA LUAD cohort, multivariate Cox proportional hazard analysis indicated that CRDC was an independent risk factor of DSS considered with factors of age, gender, smoking, clinical stage, and CIN score (p < 0.001, HR = 1.27, 95%CI: 1.12–1.4, Figure 5A). On account of the multicollinearity between variables (Figure 4D), we utilized stepwise multivariate Cox hazard regression to overcome this issue. The result showed that age, gender, smoking, and CIN had no significant effect on DSS. Only the clinical stage (p = 0.017, HR = 1.3, 95%CI: 1.1–1.7) and CRDC risk score (p < 0.001, HR = 1.3, 95%CI: 1.1–1.5) tended to be the independent risk factors for DSS, with the CRDC being the most significant factor (Figure 5B). A 1, 3, and 5 years DSS nomogram was constructed by the clinical stage and CRDC risk score showed that CRDC remained certain reference and prediction value for clinicians (corrected C-index = 0.661, Figure 5C). The CRDC score had a great benefit for predicting DSS outcome, which was observed from the decision curve analysis (DCA) of the nomogram (Figure 5D). Time-dependent AUCs of the CRDC score in predicting 1, 3, and 5 years DSS were, respectively, 0.585, 0.654, and 0.678 (Figure 5E). Moreover, in two external GEO LUAD datasets (GSE31210 and GSE72094), the CRDC score still performed a stable and robust prognostic value in predicting patients’ overall survival and acted as an independent risk factor for prognosis (Figures 5F–I). In short, we demonstrated the robust performance of CRDC for LUAD outcome prediction in TCGA and two GEO cohorts.

For investigating CRDC-related functions, differential gene analysis between CRDC-High and CRDC-Low patients was performed by DESeq2 at first. We identified that the CRDC had positive connections with multiple hallmarks such as G2M checkpoint, E2F targets, and MYC targets. Targeting WEE1, which was crucial in the G2M cell-cycle checkpoint arrest for DNA repair before mitotic entry, for inhibition and compromising the G2M checkpoint presents an opportunity to potentiate therapy (Matheson et al., 2016). Alterations in one or more key components of the core transcriptional machinery formed by the cyclin-dependent kinase (CDK)-RB-E2F axis result in heightened oncogenic E2F activity, leading to uncontrolled proliferation in cancer (Kent and Leone, 2019). Genetic deregulation of MYC expression and loss of checkpoint components, such as TP53, permit MYC to drive malignant transformation in cancer (Stine et al., 2015). On the contrary, Interferon gamma is a member of the type II interferon class, secreted by cells of both the innate and adaptive immune systems, and is crucial in antitumor response (Shankaran et al., 2001). The downregulation of interferon gamma response activity indicated that there was some connection between CRDC and immunity (Figure 6A). Also, we found high activity in the cell cycle and chemical carcinogenesis pathways, conversely, low activity in the cAMP signaling pathway, cell adhesion molecules, intestinal immune network for IgA production, B cell receptor signaling pathway, and T cell receptor signaling KEGG items (Figures 6B–D). CAMP signaling increases histone deacetylase 8 (HDAC8s) expression through the Epac–Rap1–Akt pathway leading to augmenting cisplatin-induced apoptosis (Park and Juhnn, 2017) and inhibits radiation-induced ATM phosphorylation promoting apoptosis in lung cancer (Cho et al., 2014). Some cell adhesion molecules are now considered clinical biomarkers in multiple tumor types, contributing to carcinoma progression and metastasis (Beauchemin and Arabzadeh, 2013; Smart et al., 2021). Researchers have found that the intestinal immune network for IgA production was significantly dysregulated in lung metastases from colorectal cancer (Shen et al., 2021). Decreasing activity of B cell and T cell receptor signaling implied that high CRDC risk may be related to reduced humoral immunity and cell-mediated immunity. We observed that CRDC displayed relevance to plenty of immunologic signatures (Supplementary Figures S4D,E); likewise, the suppression of immunity function appeared in the BP GSEA results (Figure 6E), which provided further evidence that CRDC was related to immunologic function. These function enrichment analyses showed that genes divided by the CRDC risk score enriched a number of known carcinoma-associated biological functions and reminded us of the immunological correlation with CRDC at the same time.

The functional analysis mentioned previously indicated that CRDC may correlate with immunological function changes; therefore, we implemented further exploration of immune characteristics in LUAD samples. CRDC correlated with tumor purity positively and had negative correlation with CYT, ESTIMATEScore, ImmuneScore, and StromalScore (Spearman correlation test, p < 0.05, Figure 7A). LUAD patients with lower CRDC scores exhibited higher PD-L1, PD-1, and CTLA4 expression levels than CRDC-High patients (Wilcoxon signed-rank test, p < 0.05, Figures 7B–D). The correlation heatmap showed that multiple CRDC genes have a significant correlation with immune cell proportion in the microenvironment (Figure 7E). Also, patients with lower CRDC scores had more B cell, T cell CD4⁺, T cell CD8⁺, neutrophil, macrophage, and myeloid dendritic cell infiltration than patients with higher CRDC scores (Figures 7F,G). The TIDE and IPS scores exhibited differences in patients distinguished by the CRDC score, which displayed that LUAD patients with decreased CRDC risk scores had lower immune dysfunction and exclusion score as well as more possibility to benefit from immune checkpoint therapy (Figures 7H,I). In addition, based on the common subtype comparison analysis between two CRDC groups and melanoma samples with the information of immunotherapy reaction (Roh et al., 2017; Chen et al., 2021), we could infer that patients with low CRDC risk were more likely to respond to anti-PD1 therapy (Bonferroni-corrected p = 0.044), yet high CRDC group may be resistant to this form of treatment (Nominal p = 0.042, Figure 7J). In brief, we revealed that CRDC had relevance to immune characteristics, and patients with low CRDC risk exhibited greater immunocompetence than those with high risk.