Prognostic Implications and Immune Infiltration Characteristics of Chromosomal Instability-Related Dysregulated CeRNA in Lung Adenocarcinoma

Introduction

Lung cancer is one of the leading causes of cancer-related deaths in both men and women worldwide, with LUAD, a kind of non-small-cell lung cancer (NSCLC), accounting for a large proportion (Reck and Rabe, 2017). CIN, the major type of genomic instability including gain/loss of whole chromosomes or large segments (aneuploidy), structural rearrangements, and focal aberrations (Geigl et al., 2008), confers tumorigenesis, metastasis, drug resistance, and poor prognosis (Lee et al., 2011; Duijf and Benezra, 2013; Bakhoum and Cantley, 2018; Bakhoum et al., 2018). CIN can be effectively detected and tended to be a risk factor for poor prognosis in LUAD (Choi et al., 2009), which mediates intratumor heterogeneity, therefore increasing the risk of recurrence or death (Jamal-Hanjani et al., 2017). CIN may counteract the therapeutic effectiveness of oncogene withdrawal treatment and be responsible for tumor relapse in lung cancer (Sotillo et al., 2010).

LncRNAs indirectly regulate target mRNAs by sharing common microRNA response elements (MREs) formatting a posttranscriptional mechanism: competing endogenous RNA (ceRNA) (Salmena et al., 2011). CeRNAs play an indispensable role in the development of carcinogenesis and could be detected based on multiple computational methods (Le et al., 2017; Zhang et al., 2019) An oncogenic lncRNA HOTAIR and a protein-coding gene HER2 associated with gastric carcinogenesis regulated each other owing to inhibition from miR-331-3p, which provides a potential anticancer treatment scheme (Liu et al., 2014). LncRNA WDFY3-AS2 and mRNA RORA are both involved in suppression of ovarian cancer by the WDFY3-AS2/miR-18a/RORA axis, in which miR-18a is reported to be an oncogene (Li W. et al., 2020). Given the ceRNA crosstalk exhibits reciprocal and complexity features, researchers focused on the complicated ceRNA regulation network to exploit cancer-associated key molecules (Sumazin et al., 2011; Karreth and Pandolfi, 2013; Song et al., 2017). Perturbation of ceRNA interaction wildly exists in disease versus normal status. LncRNAs with ceRNA activity could be candidate epigenetic diagnostic biomarkers for early detection of osteoporosis by constructing an osteoporosis-related dysregulated ceRNA network (Zhang et al., 2020). Also, epigenetically related lncRNAs involved in dysregulated ceRNA–ceRNA networks offered novel potential molecular therapeutic targets across pan-cancer (Xu et al., 2021). However, the molecular function of dysregulated ceRNA in LUAD remains further elucidated.

Based on the molecular pattern, immunotherapy has become a new treatment option for NSCLC in recent years (Reck and Rabe, 2017). The blockade of immune checkpoints targeting cytotoxic T lymphocyte—associated antigen 4 (CTLA4) and the programmed cell death protein 1 pathway (PD-1/PD-L1) have demonstrated promise in stimulating antitumor immunity (Pardoll, 2012). Existing research studies have shown that clinical outcome and treatment response to immune checkpoint blockers are affected by the composition and proportion of various immune cells (Brahmer, 2013; Schalper et al., 2015; Kadara et al., 2017). A vital element in precision diagnosis and personalized treatment of LUAD is exploring novel molecular signatures, especially those associated with the tumor immune microenvironment. Moreover, the ceRNA regulation mechanism could affect immune cell infiltration in multiple cancers. For instance, LINC00301 could serve as a competing endogenous RNA (ceRNA) against miR-1276 to expedite the HIF1α pathway in the cytoplasm of NSCLC cells facilitating tumor progression and triggering an immune-suppressing microenvironment (Sun et al., 2020). SNHG16 serves as a ceRNA by sponging miR-16-5p, which led to the derepression of its target gene SMAD5 and resulted in potentiation of the TGF-β1/SMAD5 pathway to induce immunosuppressive CD73+γδ1 Treg cells (Ni et al., 2020). So far, no studies have systematically excavated the role of CIN-related ceRNAs in immune regulation for LUAD patients.

Here, we analyzed the CIN features in LUAD samples and systematically integrated CIN-related ceRNA characteristics for the first time by constructing a double-weighted CIN-related dysregulated ceRNA network, aiming to identify a potential biomarker for LUAD prognosis. The marker we found was assessed combined with multiple clinical factors. Furthermore, it was associated with immune infiltration which may help predict whether patients will benefit from immunotherapy (Figure 1). This research provided a perspective for identifying prognostic markers and exploring the immune microenvironment characteristics of LUAD patients.

FIGURE 1

FIGURE 1. Flowchart of this research. Exploration of CIN feature in expression level, clinical outcome relevance, and functional analysis; construction of a CIN-related double-weighted dysregulated ceRNA network; acquisition of a potential biomarker CRDC after dense module searching, filtering, and integration; analyzation and evaluation of CRDC in the prognostic value; inference for immune cell infiltration and benefit from immunotherapy.

Materials and Methods

Data Acquisition and Preprocessing

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: $TPM=\frac{FPKM}{{{\displaystyle \sum }}_{all genes}FPKM}\ast {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).

MiRNA–Target Regulation Relationships

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.

Estimation of the Chromosomal Instability Risk Score of Samples

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.

Construction of a Chromosomal Instability-Related Dysregulated CeRNA Network

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-{\displaystyle \sum }_{i=0}^{r-1}\frac{\left(\begin{array}{c}K\\ i\end{array}\right)\left(\begin{array}{c}N-K\\ M-i\end{array}\right)}{\left(\begin{array}{c}N\\ M\end{array}\right)}$, 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_{CIN\_High}$ and $r_{CIN\_Low}$ represented the PCC in two groups, and $n_{CIN\_High}$ and $n_{CIN\_Low}$ 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 ${\phi }^{-1}\left[1-2\ast \left(1-\phi \left(\left|\frac{F\left(r_{CIN\_High}\right)-F\left(r_{CIN\_Low}\right)}{\sqrt{\frac{1}{n_{CIN\_High}+3}+\frac{1}{n_{CIN\_Low}+3}}}\right|\right)\right)\right]$, where $F\left(x\right)=\frac{1}{2}ln\frac{1+x}{1-x}$, and $\phi$ denoted the standard normal distribution function. Meanwhile, the node weight (N) of the network was defined as ${\phi }^{-1}\left(1-p\right)$, 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.

Dense Modules Identification and Filtration

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 $\lambda \frac{{\displaystyle \sum E}}{\sqrt{n_E}}+\left(1-\lambda \right)\frac{{\displaystyle \sum N}}{\sqrt{n_N}}$, where $n_E$ and $n_N$ were the number of edges and nodes in the network, and $\lambda$ 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).

Calculation of the CIN-Related Dysregulated CeRNAs Risk Score

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 ${\displaystyle \sum }_{i=1}^{N}Expressio{n}_{\_i}\text{*log}\left(H{R}_{\_i}\right)$, where N was the number of the CRDC genes.

Immune Infiltration Analysis in LUAD Patients

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).

Statistical Analysis

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). $SNR=\frac{{{\displaystyle \sum }}_{i=1}^n{\left({\mu }_{1}I\left(i\le t\right)+{\mu }_{2}I\left(i>t\right)-\mu \right)}^2}{{{\displaystyle \sum }}_{i=1}^n{\left({\mu }_{1}I\left(i\le t\right)+{\mu }_{2}I\left(i>t\right)-CI{N}_i\right)}^2},$ where I is the indicator function and is set to 0 when the condition in brackets is not satisfied, $\mu$ is the mean CIN score for all patients, and ${\mu }_1$ and ${\mu }_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.

Results

Characteristics of Chromosomal Instability in Lung Adenocarcinoma

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).

FIGURE 2

FIGURE 2. Characteristics of chromosomal instability in lung adenocarcinoma. (A) Boxplot for CIN scores in LUAD and adjacent paired samples. The green column indicated adjacent normal samples, and the red column indicated LUAD tissues. (B) ROC curve of the CIN score in discriminating LUAD patients from adjacent normal samples. The green color indicated AUC. (C) Boxplots of patients’ CIN scores in different TNM stages and clinical stages. (D) Volcano plots for identifying significantly differentially expressed genes in (left) lung adenocarcinoma samples versus adjacent samples and (right) CIN-High compared to CIN-Low group. (E) Venn diagram and (F) BP function and KEGG pathway enrichment results for differentially expressed genes in comparative groups of LUAD vs. Adjacent and CIN-High vs. CIN-Low.

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).

Construction of the CIN-Related Dysregulated CeRNA Network

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).

FIGURE 3

FIGURE 3. Construction of a CIN-related dysregulated ceRNA network. (A) CIN-related dysregulated double-weighted ceRNA network, in which mRNAs and lncRNAs were colored green and red, respectively. Node degree was indicated by the node size, and node weight and edge weight were indicated by node transparency and edge width, respectively. (B) Degree distribution of the network. Boxplots of (C) mRNA and lncRNA degree and (D) alteration of ceRNA dysregulation extent in CIN-High compared with CIN-Low samples. (E) Subnetwork CRDC connected by dense modules after screening and integration.

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.

Filtration for Double-Weighted Dense Modules

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).

Formulation and Evaluation of the CRDC Risk Score and Prognostic Analysis

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).

FIGURE 4

FIGURE 4. Formulation of the CRDC risk score and prognostic analysis. (A) Barplot depicted the log-transformed hazard ratio of CRDC genes in the univariate Cox regression model. (B) Boxplot for CRDC of LUAD and adjacent paired samples. The green column indicated adjacent samples, and the red column indicated LUAD tissues. (C) ROC curve of CRDC in discriminating LUAD patients from adjacent control samples. The green color indicated AUC. (D) Spearman correlation analysis of smoking, clinical stage, CIN, and CRDC risk score. Kaplan–Meier survival curves of the CRDC risk score in (E–H) all LUAD patients (TCGA_ALL, for short) and (I–L) patients without postoperative treatment (TCGA_Notreat, for short).

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.

Prognostic Independence Evaluation and External Prognosis Validation of CRDC

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.

FIGURE 5

FIGURE 5. Prognostic independence evaluation and external validation. (A) Forest plot for the multiple Cox proportional hazards model in TCGA patients. (B) Forest plot for the multiple Cox proportional hazards model contained the clinical stage and CRDC risk score in TCGA patients (C) Nomogram constructed by the clinical stage and CRDC risk score for predicting 1, 3, and 5 years DSS. (D) Decision curve of the nomogram for 1, 3, and 5 years DSS. (E) Time-dependent ROC of the CRDC risk score for predicting 1, 3, and 5 years DSS outcome. Kaplan–Meier survival curves of CRDC-High and CRDC-Low patients and forest plots for the multiple Cox proportional hazards model in (F,G) GSE31210 and (H,I) GSE72094 datasets.

Exploration of CRDC-Related Biological Functions Indicating the Immunology Correlation

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.

FIGURE 6

FIGURE 6. CRDC relevant functional analysis. Ridgeline plots for log2FC (CRDC-High versus CRDC-Low) distribution of core-enriched genes in (A) hallmarks and (B) KEGG pathways based on the GSEA method. GSEA plots for (C) positive and (D) negative KEGG pathways associated with CRDC. (E) Table of the CRDC-related negative GO (BP) enrichment graph.

Immune Infiltration Microenvironment Characteristics Related With CRDC

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.

FIGURE 7

FIGURE 7. CRDC-related immune infiltration characteristics. (A) Spearman correlation analysis plot of the correlativity of the CRDC score, StromalScore, ImmuneScore, ESTIMATEScore, TumorPurity, TMB, and CYT. Boxplots described the expression difference of (B) PD-L1, (C) PD-1, and (D) CTLA4 between CRDC-High and CRDC-Low samples. (E) Spearman correlation heatmap of CRDC genes and immune microenvironment cell infiltration, StromalScore, ImmuneScore, ESTIMATEScore, TumorPurity, TMB, CIN, CYT, PD-L1, PD-1, and CTLA4. (F,G) Immune microenvironment cell infiltration level in CRDC-High and CRDC-Low samples. The violin plots displayed the TIDE (H) and IPS (I) scores between CRDC-High and -Low groups. (J) Heatmap of comparison and identification of common subtypes between groups divided by CRDC risk and response for ICB therapy. NR and R represented not respond and respond to immunotherapy, respectively.

Discussion

Previous studies have shown that CIN was closely related to tumorigenesis. Gain of chromosomes showed increased DNA damage and sensitivity to replication stress, thereby promoting genomic instability and possibly contributing to tumorigenesis (Passerini et al., 2016). Based on discussion of mathematical models of situations in which inactivation of one or two tumor suppressor genes is required for tumorigenesis; CIN will arise before inactivation of the first suppressor genes, therefore initiating the mutational sequence that leads to cancer (Michor et al., 2005). In restricted genetic and cellular contexts, the mitotic checkpoint protein insufficiency can cause whole chromosome instability and drive tumor formation through tumor suppressor gene loss of heterozygosity (Baker et al., 2009). In this work, we investigated the CIN characteristics in 513 TCGA LUAD patients and found that the extent of CIN in LUAD was significantly higher than that in adjacent normal samples, and it can be used to strictly separate the tumor from adjacent normal samples (Figures 2A,B). CIN 70 marker genes expressed a high level in LUAD (Figure 2D). Also, CIN was related to high clinical stage as well as poor survival in LUAD patients (Figure 2C and Supplementary Figure S1). Functional analysis results indicated that CIN may have an effect on tumor progression by promoting 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). These results were identical with the previous theory that tumors are generally characterized by genomic instability.

For the past few years, the molecular regulatory mechanism of ceRNA has been proved to play a crucial role in development of diseases, especially cancers and their dysregulation may conduce to tumor pathogenesis. In our previous work, we have constructed a clear cell kidney carcinoma dysregulated ceRNA–ceRNA network and identified two dysregulated patterns of ceRNAs interaction (gain and loss), which were demonstrated to be able to distinguish normal samples from cancer samples (Wang et al., 2018). The ceRNA dysregulation mechanism in LUAD needed expansion, and no research has explored characteristics of ceRNA alteration in the CIN level yet. In the present work, we explored CIN-related dysregulated ceRNAs in LUAD aiming for acquiring a prognosis indicator for the first time. Based on differential expression genes (mRNAs and lncRNAs) in LUAD compared to adjacent samples, we constructed a double-weighted CIN-related dysregulated ceRNA network. The edge weight in the network was the alteration of ceRNA PCC in CIN-High compared with CIN-Low patients representing the extent of ceRNA dysregulation. Also, the node weight represented the potency of the gene in predicting LUAD patient’s outcome. At this step, we adopted the samples without postoperative treatment eliminating the possible impact of treatment on prognosis. Also, we used DSS as a measure of prognosis in consideration of its definition, death from the diagnosed cancer type, has much greater relevance with cancer-associated clinical outcome than OS in which the endpoint may record non-cancer causing death (Liu et al., 2018). Based on the dysregulated network we built, 62.3% of the ceRNAs maintained a weaker co-expression relationship in the CIN-High than CIN-Low group (Figure 3D and Supplementary Figures S2D,E). Also, the greedy search algorithm was applied to detect modules with high double weight, whereafter we obtained a subnetwork CRDC formed by 12 genes after module screening and integrating (Figure 3E).

The CRDC risk score was higher in LUAD samples than that in adjacent normal samples and could separate the two types of samples, which was resembled with CIN (Figure 4B,C). Subsequently, we laid special stress on testing and evaluating the prognostic prediction efficiency of CRDC. The results implied that the high CRDC score can serve as an indicator for poor OSS, DSS, and PFI in patients with or without postoperative treatment and a key independent risk factor in TCGA and two GEO datasets (Figure 4E–L, Figure 5). In the multivariate Cox proportional hazard model, CIN did not show as an independent indicator for predicting DSS. It demonstrated that prognostic efficacy of CIN was not as stable as that of CRDC. We suspected that there was probably a connection with the dual natures and complicacy of CIN in cancer progressing. Extreme CIN was related to long-term survival in primary breast cancer (Roylance et al., 2011). In addition, the paradoxical and nonmonotonic relationship between CIN and prognosis was observed in ovarian, gastric, and non-small-cell lung tumors (Birkbak et al., 2011). CIN plays a multifaceted role in cancer, and its microenvironment, for instance, by introducing double-stranded DNA into the cytosol, CIN could engage the cGAS–STING antiviral pathway to facilitate inflammatory signaling (Bakhoum and Cantley, 2018).

The functional analysis results indicated that CRDC was probably relevant with biological function terms associated with tumor progression. For instance, high activity in G2M checkpoint, E2F targets, MYC targets, cell cycle, and chemical carcinogenesis, and low activity in the cAMP signaling pathway and cell adhesion molecules in opposite. At the same time, a number of immune-associated terms, such as interferon gamma response and B cell and T cell receptor signaling, exhibited a negative changing trend. These results suggested that the poor prognosis caused by CRDC may be achieved by promoting cell proliferation and migration as well as reducing antitumor immune response.

Based on the possibility of the connection between CRDC and immunity from the functional analysis, we consulted literatures and found that several CRDC genes have been proven to be associated with immunologic function. LINC00473 silencing enhanced miR-195-5p-targeted downregulation of PD-L1 in pancreatic cancer may block the cancer progression (Zhou W. Y. et al., 2019). The proprotein convertase furin was negatively correlated with immune cell infiltration in triple negative breast cancer fitting with our correlation analysis about the CRDC genes and immune features (Figure 7E). Also, furin deficiency in T cells decreased Tregs resulting in CD8⁺ T cell activation and IFN-γ upregulation (He et al., 2020). IFN-γ can augment immune function, however, induce expression of PD-L1 by which IFN-γ impairs antitumor immunity (Mandai et al., 2016). Also, tumor IFNG signaling blocking improves ICB response by CD8⁺ T cell and NK/ILC1-mediated killing (Benci et al., 2019). We found a higher expression level of IFNG and immune checkpoint genes in the CIN-Low group versus CIN-High group (Supplementary Figure S4C and Figures 7B–D), which suggested that CIN-Low patients had possibility to benefit from ICB. Also, Patients with low CRDC risk had lower immune dysfunction and exclusion score as well as higher CYT, immunophenoscore, and immune cell infiltration, which illustrated that these patients represented characteristics of “immune-hot” (Figures 7A,F–I). Immunotherapy response prediction offered evidence of the anti-PD1 therapy potential applicability for CIN-Low patients (Figure 7J).

There are two limitations in our work. We adopted the patients in TCGA cohort without postoperative treatment to evaluate prognosis for excluding the impact of treatment as much as possible. However, these untreated samples may be derived from information records missing, and we could not acquire the unambiguous treatment information from the external GEO verification datasets. In addition, the CRDC-related biological functions and immunological characteristics needed in-depth experimental verification.

Conclusion

In conclusion, we analyzed the CIN feature and constructed a double-weighted CIN-related dysregulated ceRNA network. Based on this network, a potential marker CRDC composed of 12 RNAs was acquired which could distinguish LUAD samples from adjacent normal tissues and correlated with multiple tumorigenic hallmarks and pathways. Also, CRDC probably serves as an indicator for evaluating prognosis and recommending immune checkpoint inhibitor therapy. This study will offer a novel perspective for understanding the molecular action in LUAD tumorigenesis and progression.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

Author Contributions

KL, YC, and QX conceived and designed the experiments. SG, TL, DX, JX, HW, JL, XB, MC, and ZX analyzed the data. KL and YC wrote the manuscript. All authors have read and approved the final manuscript.

Funding

This work was supported by the Key Research and Development Project of Hainan Province (ZDYF2021SHFZ097), National Natural Science Foundation of China (32160152), Major Science and Technology Program of Hainan Province (ZDKJ2021040, ZDKJ202003), Hainan Provincial Natural Science Foundation of China (820RC637), and Hainan Province Clinical Medical Center.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2022.843640/full#supplementary-material

References

Baker, D. J., Jin, F., Jeganathan, K. B., and van Deursen, J. M. (2009). Whole Chromosome Instability Caused by Bub1 Insufficiency Drives Tumorigenesis through Tumor Suppressor Gene Loss of Heterozygosity. Cancer Cell 16 (6), 475–486. doi:10.1016/j.ccr.2009.10.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Bakhoum, S. F., and Cantley, L. C. (2018). The Multifaceted Role of Chromosomal Instability in Cancer and its Microenvironment. Cell 174 (6), 1347–1360. doi:10.1016/j.cell.2018.08.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Bakhoum, S. F., Ngo, B., Laughney, A. M., Cavallo, J.-A., Murphy, C. J., Ly, P., et al. (2018). Chromosomal Instability Drives Metastasis through a Cytosolic DNA Response. Nature 553 (7689), 467–472. doi:10.1038/nature25432

PubMed Abstract | CrossRef Full Text | Google Scholar

Beauchemin, N., and Arabzadeh, A. (2013). Carcinoembryonic Antigen-Related Cell Adhesion Molecules (CEACAMs) in Cancer Progression and Metastasis. Cancer Metastasis Rev. 32 (3-4), 643–671. doi:10.1007/s10555-013-9444-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Benci, J. L., Johnson, L. R., Choa, R., Xu, Y., Qiu, J., Zhou, Z., et al. (2019). Opposing Functions of Interferon Coordinate Adaptive and Innate Immune Responses to Cancer Immune Checkpoint Blockade. Cell 178 (4), 933–948. doi:10.1016/j.cell.2019.07.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Birkbak, N. J., Eklund, A. C., Li, Q., McClelland, S. E., Endesfelder, D., Tan, P., et al. (2011). Paradoxical Relationship between Chromosomal Instability and Survival Outcome in Cancer. Cancer Res. 71 (10), 3447–3452. doi:10.1158/0008-5472.CAN-10-3667

PubMed Abstract | CrossRef Full Text | Google Scholar

Brahmer, J. R. (2013). Harnessing the Immune System for the Treatment of Non-small-cell Lung Cancer. J. Clin. Oncol. 31 (8), 1021–1028. doi:10.1200/JCO.2012.45.8703

PubMed Abstract | CrossRef Full Text | Google Scholar

Carter, S. L., Eklund, A. C., Kohane, I. S., Harris, L. N., and Szallasi, Z. (2006). A Signature of Chromosomal Instability Inferred from Gene Expression Profiles Predicts Clinical Outcome in Multiple Human Cancers. Nat. Genet. 38 (9), 1043–1048. doi:10.1038/ng1861

PubMed Abstract | CrossRef Full Text | Google Scholar

Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cel Rep. 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Zhu, Z., Li, X., Yao, X., and Luo, L. (2021). The Ferroptosis-Related Noncoding RNA Signature as a Novel Prognostic Biomarker in the Tumor Microenvironment, Immunotherapy, and Drug Screening of Gastric Adenocarcinoma. Front. Oncol. 11, 778557. doi:10.3389/fonc.2021.778557

PubMed Abstract | CrossRef Full Text | Google Scholar

Cho, E.-A., Kim, E.-J., Kwak, S.-J., and Juhnn, Y.-S. (2014). cAMP Signaling Inhibits Radiation-Induced ATM Phosphorylation Leading to the Augmentation of Apoptosis in Human Lung Cancer Cells. Mol. Cancer 13, 36. doi:10.1186/1476-4598-13-36

PubMed Abstract | CrossRef Full Text | Google Scholar

Choi, C.-M., Seo, K. W., Jang, S. J., Oh, Y.-M., Shim, T.-S., Kim, W. S., et al. (2009). Chromosomal Instability Is a Risk Factor for Poor Prognosis of Adenocarcinoma of the Lung: Fluorescence In Situ Hybridization Analysis of Paraffin-Embedded Tissue from Korean Patients. Lung Cancer 64 (1), 66–70. doi:10.1016/j.lungcan.2008.07.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Chou, C.-H., Shrestha, S., Yang, C.-D., Chang, N.-W., Lin, Y.-L., Liao, K.-W., et al. (2018). miRTarBase Update 2018: a Resource for Experimentally Validated microRNA-Target Interactions. Nucleic Acids Res. 46 (D1), D296–D302. doi:10.1093/nar/gkx1067

PubMed Abstract | CrossRef Full Text | Google Scholar

Duijf, P. H. G., and Benezra, R. (2013). The Cancer Biology of Whole-Chromosome Instability. Oncogene 32 (40), 4727–4736. doi:10.1038/onc.2012.616

PubMed Abstract | CrossRef Full Text | Google Scholar

Frankish, A., Diekhans, M., Ferreira, A.-M., Johnson, R., Jungreis, I., Loveland, J., et al. (2019). GENCODE Reference Annotation for the Human and Mouse Genomes. Nucleic Acids Res. 47 (D1), D766–D773. doi:10.1093/nar/gky955

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, J., Li, K., Zhang, W., Wan, C., Zhang, J., Jiang, P., et al. (2020). Large-scale Public Data Reuse to Model Immunotherapy Response and Resistance. Genome Med. 12 (1), 21. doi:10.1186/s13073-020-0721-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Geigl, J. B., Obenauf, A. C., Schwarzbraun, T., and Speicher, M. R. (2008). Defining 'chromosomal Instability'. Trends Genet. 24 (2), 64–69. doi:10.1016/j.tig.2007.11.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Grady, W. M., and Carethers, J. M. (2008). Genomic and Epigenetic Instability in Colorectal Cancer Pathogenesis. Gastroenterology 135 (4), 1079–1099. doi:10.1053/j.gastro.2008.07.076

PubMed Abstract | CrossRef Full Text | Google Scholar

He, Z., Khatib, A.-M., and Creemers, J. W. M. (2020). Loss of the Proprotein Convertase Furin in T Cells Represses Mammary Tumorigenesis in Oncogene-Driven Triple Negative Breast Cancer. Cancer Lett. 484, 40–49. doi:10.1016/j.canlet.2020.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Heagerty, P. J., and Zheng, Y. (2005). Survival Model Predictive Accuracy and ROC Curves. Biometrics 61 (1), 92–105. doi:10.1111/j.0006-341X.2005.030814.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Jamal-Hanjani, M., Wilson, G. A., McGranahan, N., Birkbak, N. J., Watkins, T. B. K., Veeriah, S., et al. (2017). Tracking the Evolution of Non-small-cell Lung Cancer. N. Engl. J. Med. 376 (22), 2109–2121. doi:10.1056/NEJMoa1616288

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T Cell Dysfunction and Exclusion Predict Cancer Immunotherapy Response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadara, H., Choi, M., Zhang, J., Parra, E. R., Rodriguez-Canales, J., Gaffney, S. G., et al. (2017). Whole-exome Sequencing and Immune Profiling of Early-Stage Lung Adenocarcinoma with Fully Annotated Clinical Follow-Up. Ann. Oncol. 28 (1), 75–82. doi:10.1093/annonc/mdw436

PubMed Abstract | CrossRef Full Text | Google Scholar

Karreth, F. A., and Pandolfi, P. P. (2013). ceRNA Cross-Talk in Cancer: when Ce-Bling Rivalries Go Awry. Cancer Discov. 3 (10), 1113–1121. doi:10.1158/2159-8290.CD-13-0202

PubMed Abstract | CrossRef Full Text | Google Scholar

Kent, L. N., and Leone, G. (2019). The Broken Cycle: E2F Dysfunction in Cancer. Nat. Rev. Cancer 19 (6), 326–338. doi:10.1038/s41568-019-0143-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Le, T. D., Zhang, J., Liu, L., and Li, J. (2017). Computational Methods for Identifying miRNA Sponge Interactions. Brief Bioinform 18 (4), 577–590. doi:10.1093/bib/bbw042

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, A. J. X., Endesfelder, D., Rowan, A. J., Walther, A., Birkbak, N. J., Futreal, P. A., et al. (2011). Chromosomal Instability Confers Intrinsic Multidrug Resistance. Cancer Res. 71 (5), 1858–1870. doi:10.1158/0008-5472.CAN-10-3604

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, B., and Dewey, C. N. (2011). RSEM: Accurate Transcript Quantification from RNA-Seq Data with or without a Reference Genome. BMC Bioinformatics 12, 323. doi:10.1186/1471-2105-12-323

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, E., Xu, Z., Liu, F., Wang, H., Wen, J., Shao, S., et al. (2014). Continual Exposure to Cigarette Smoke Extracts Induces Tumor-like Transformation of Human Nontumor Bronchial Epithelial Cells in a Microfluidic Chip. J. Thorac. Oncol. 9 (8), 1091–1100. doi:10.1097/JTO.0000000000000219

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, J., Roebuck, P., Grunewald, S., and Liang, H. (2012). SurvNet: a Web Server for Identifying Network-Based Biomarkers that Most Correlate with Patient Survival Data. Nucleic Acids Res. 40 (Web Server issue), W123–W126. doi:10.1093/nar/gks386

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, T., Fu, J., Zeng, Z., Cohen, D., Li, J., Chen, Q., et al. (2020). TIMER2.0 for Analysis of Tumor-Infiltrating Immune Cells. Nucleic Acids Res. 48 (W1), W509–W514. doi:10.1093/nar/gkaa407

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, W., Ma, S., Bai, X., Pan, W., Ai, L., and Tan, W. (2020). Long Noncoding RNA WDFY3‐AS2 Suppresses Tumor Progression by Acting as a Competing Endogenous RNA of microRNA‐18a in Ovarian Cancer. J. Cel Physiol 235 (2), 1141–1154. doi:10.1002/jcp.29028

CrossRef Full Text | Google Scholar

Li, J.-H., J.-H., Liu, S., Zhou, H., Qu, L.-H., and Yang, J.-H. (2014). starBase v2.0: Decoding miRNA-ceRNA, miRNA-ncRNA and Protein-RNA Interaction Networks from Large-Scale CLIP-Seq Data. Nucl. Acids Res. 42 (Database issue), D92–D97. doi:10.1093/nar/gkt1248

PubMed Abstract | CrossRef Full Text | Google Scholar

Liberzon, A., Birger, C., Thorvaldsdóttir, H., Ghandi, M., Mesirov, J. P., and Tamayo, P. (2015). The Molecular Signatures Database Hallmark Gene Set Collection. Cel Syst. 1 (6), 417–425. doi:10.1016/j.cels.2015.12.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Lichtenberg, T., Hoadley, K. A., Poisson, L. M., Lazar, A. J., Cherniack, A. D., et al. (2018). An Integrated TCGA Pan-Cancer Clinical Data Resource to Drive High-Quality Survival Outcome Analytics. Cell 173 (2), 400–e11. doi:10.1016/j.cell.2018.02.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, X.-h., Sun, M., Nie, F.-q., Ge, Y.-b., Zhang, E.-b., Yin, D.-d., et al. (2014). Lnc RNA HOTAIR Functions as a Competing Endogenous RNA to Regulate HER2 Expression by Sponging miR-331-3p in Gastric Cancer. Mol. Cancer 13, 92. doi:10.1186/1476-4598-13-92

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, Y., Ji, Y., and Qiu, P. (2013). Identification of Thresholds for Dichotomizing DNA Methylation Data. J. Bioinform Sys Biol. 2013 (1), 8. doi:10.1186/1687-4153-2013-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Love, M. I., Huber, W., and Anders, S. (2014). Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2. Genome Biol. 15 (12), 550. doi:10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Mandai, M., Hamanishi, J., Abiko, K., Matsumura, N., Baba, T., and Konishi, I. (2016). Dual Faces of IFNγ in Cancer Progression: A Role of PD-L1 Induction in the Determination of Pro- and Antitumor Immunity. Clin. Cancer Res. 22 (10), 2329–2334. doi:10.1158/1078-0432.CCR-16-0224

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, S. A., Hewish, M., Lord, C. J., and Ashworth, A. (2009). Genomic Instability and the Selection of Treatments for Cancer. J. Pathol. 220 (2), 281–289. doi:10.1002/path.2631

CrossRef Full Text | Google Scholar

Matheson, C. J., Backos, D. S., and Reigan, P. (2016). Targeting WEE1 Kinase in Cancer. Trends Pharmacol. Sci. 37 (10), 872–881. doi:10.1016/j.tips.2016.06.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, Y.-R., Liu, W., Zhang, Q., and Guo, A.-Y. (2018). lncRNASNP2: an Updated Database of Functional SNPs and Mutations in Human and Mouse lncRNAs. Nucleic Acids Res. 46 (D1), D276–D280. doi:10.1093/nar/gkx1004

PubMed Abstract | CrossRef Full Text | Google Scholar

Michor, F., Iwasa, Y., Vogelstein, B., Lengauer, C., and Nowak, M. A. (2005). Can Chromosomal Instability Initiate Tumorigenesis? Semin. Cancer Biol. 15 (1), 43–49. doi:10.1016/j.semcancer.2004.09.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Ni, C., Fang, Q.-Q., Chen, W.-Z., Jiang, J.-X., Jiang, Z., Ye, J., et al. (2020). Breast Cancer-Derived Exosomes Transmit lncRNA SNHG16 to Induce CD73+γδ1 Treg Cells. Sig Transduct Target. Ther. 5 (1), 41. doi:10.1038/s41392-020-0129-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Paraskevopoulou, M. D., Vlachos, I. S., Karagkouni, D., Georgakilas, G., Kanellos, I., Vergoulis, T., et al. (2016). DIANA-LncBase V2: Indexing microRNA Targets on Non-coding Transcripts. Nucleic Acids Res. 44 (D1), D231–D238. doi:10.1093/nar/gkv1270

PubMed Abstract | CrossRef Full Text | Google Scholar

Pardoll, D. M. (2012). The Blockade of Immune Checkpoints in Cancer Immunotherapy. Nat. Rev. Cancer 12 (4), 252–264. doi:10.1038/nrc3239

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, J.-Y., and Juhnn, Y.-S. (2017). cAMP Signaling Increases Histone Deacetylase 8 Expression via the Epac2-Rap1A-Akt Pathway in H1299 Lung Cancer Cells. Exp. Mol. Med. 49 (2), e297. doi:10.1038/emm.2016.152

PubMed Abstract | CrossRef Full Text | Google Scholar

Passerini, V., Ozeri-Galai, E., de Pagter, M. S., Donnelly, N., Schmalbrock, S., Kloosterman, W. P., et al. (2016). The Presence of Extra Chromosomes Leads to Genomic Instability. Nat. Commun. 7, 10754. doi:10.1038/ncomms10754

PubMed Abstract | CrossRef Full Text | Google Scholar

Reck, M., and Rabe, K. F. (2017). Precision Diagnosis and Treatment for Advanced Non-small-cell Lung Cancer. N. Engl. J. Med. 377 (9), 849–861. doi:10.1056/NEJMra1703413

PubMed Abstract | CrossRef Full Text | Google Scholar

Robin, X., Turck, N., Hainard, A., Tiberti, N., Lisacek, F., Sanchez, J.-C., et al. (2011). pROC: an Open-Source Package for R and S+ to Analyze and Compare ROC Curves. BMC Bioinformatics 12, 77. doi:10.1186/1471-2105-12-77

PubMed Abstract | CrossRef Full Text | Google Scholar

Roh, W., Chen, P.-L., Reuben, A., Spencer, C. N., Prieto, P. A., Miller, J. P., et al. (2017). Integrated Molecular Analysis of Tumor Biopsies on Sequential CTLA-4 and PD-1 Blockade Reveals Markers of Response and Resistance. Sci. Transl. Med. 9 (379), eaah3560. doi:10.1126/scitranslmed.aah3560

PubMed Abstract | CrossRef Full Text | Google Scholar

Rooney, M. S., Shukla, S. A., Wu, C. J., Getz, G., and Hacohen, N. (2015). Molecular and Genetic Properties of Tumors Associated with Local Immune Cytolytic Activity. Cell 160 (1-2), 48–61. doi:10.1016/j.cell.2014.12.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Roylance, R., Endesfelder, D., Gorman, P., Burrell, R. A., Sander, J., Tomlinson, I., et al. (2011). Relationship of Extreme Chromosomal Instability with Long-Term Survival in a Retrospective Analysis of Primary Breast Cancer. Cancer Epidemiol. Biomarkers Prev. 20 (10), 2183–2194. doi:10.1158/1055-9965.EPI-11-0343

PubMed Abstract | CrossRef Full Text | Google Scholar

Sahoo, D., Dill, D. L., Tibshirani, R., and Plevritis, S. K. (2007). Extracting Binary Signals from Microarray Time-Course Data. Nucleic Acids Res. 35 (11), 3705–3712. doi:10.1093/nar/gkm284

PubMed Abstract | CrossRef Full Text | Google Scholar

Saletta, F., Matullo, G., Manuguerra, M., Arena, S., Bardelli, A., and Vineis, P. (2007). Exposure to the Tobacco Smoke Constituent 4-aminobiphenyl Induces Chromosomal Instability in Human Cancer Cells. Cancer Res. 67 (15), 7088–7094. doi:10.1158/0008-5472.CAN-06-4420

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA Hypothesis: the Rosetta Stone of a Hidden RNA Language? Cell 146 (3), 353–358. doi:10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Schalper, K. A., Brown, J., Carvajal-Hausdorf, D., McLaughlin, J., Velcheti, V., Syrigos, K. N., et al. (2015). Objective Measurement and Clinical Significance of TILs in Non-small Cell Lung Cancer. J. Natl. Cancer Inst. 107 (3), dju435. doi:10.1093/jnci/dju435

PubMed Abstract | CrossRef Full Text | Google Scholar

Shankaran, V., Ikeda, H., Bruce, A. T., White, J. M., Swanson, P. E., Old, L. J., et al. (2001). IFNγ and Lymphocytes Prevent Primary Tumour Development and Shape Tumour Immunogenicity. Nature 410 (6832), 1107–1111. doi:10.1038/35074122

PubMed Abstract | CrossRef Full Text | Google Scholar

Shen, T., Liu, J.-L., Wang, C.-Y., Rixiati, Y., Li, S., Cai, L.-D., et al. (2021). Targeting Erbin in B Cells for Therapy of Lung Metastasis of Colorectal Cancer. Sig Transduct Target. Ther. 6 (1), 115. doi:10.1038/s41392-021-00501-x

CrossRef Full Text | Google Scholar

Sjöstedt, E., Zhong, W., Fagerberg, L., Karlsson, M., Mitsios, N., Adori, C., et al. (2020). An Atlas of the Protein-Coding Genes in the Human, Pig, and Mouse Brain. Science 367 (6482), eaay5947. doi:10.1126/science.aay5947

PubMed Abstract | CrossRef Full Text | Google Scholar

Smart, J. A., Oleksak, J. E., and Hartsough, E. J. (2021). Cell Adhesion Molecules in Plasticity and Metastasis. Mol. Cancer Res. 19 (1), 25–37. doi:10.1158/1541-7786.MCR-20-0595

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Y.-x., Sun, J.-x., Zhao, J.-h., Yang, Y.-c., Shi, J.-x., Wu, Z.-h., et al. (2017). Non-coding RNAs Participate in the Regulatory Network of CLDN4 via ceRNA Mediated miRNA Evasion. Nat. Commun. 8 (1), 289. doi:10.1038/s41467-017-00304-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Sotillo, R., Schvartzman, J.-M., Socci, N. D., and Benezra, R. (2010). Mad2-induced Chromosome Instability Leads to Lung Tumour Relapse after Oncogene Withdrawal. Nature 464 (7287), 436–440. doi:10.1038/nature08803

PubMed Abstract | CrossRef Full Text | Google Scholar

Stine, Z. E., Walton, Z. E., Altman, B. J., Hsieh, A. L., and Dang, C. V. (2015). MYC, Metabolism, and Cancer. Cancer Discov. 5 (10), 1024–1039. doi:10.1158/2159-8290.CD-15-0507

PubMed Abstract | CrossRef Full Text | Google Scholar

Subramanian, A., Tamayo, P., Mootha, V. K., Mukherjee, S., Ebert, B. L., Gillette, M. A., et al. (2005). Gene Set Enrichment Analysis: a Knowledge-Based Approach for Interpreting Genome-wide Expression Profiles. Proc. Natl. Acad. Sci. 102 (43), 15545–15550. doi:10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

Sumazin, P., Yang, X., Chiu, H.-S., Chung, W.-J., Iyer, A., Llobet-Navas, D., et al. (2011). An Extensive microRNA-Mediated Network of RNA-RNA Interactions Regulates Established Oncogenic Pathways in Glioblastoma. Cell 147 (2), 370–381. doi:10.1016/j.cell.2011.09.041

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, C.-C., Zhu, W., Li, S.-J., Hu, W., Zhang, J., Zhuo, Y., et al. (2020). FOXC1-mediated LINC00301 Facilitates Tumor Progression and Triggers an Immune-Suppressing Microenvironment in Non-small Cell Lung Cancer by Regulating the HIF1α Pathway. Genome Med. 12 (1), 77. doi:10.1186/s13073-020-00773-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomczak, K., Czerwińska, P., and Wiznerowicz, M. (2015). Review the Cancer Genome Atlas (TCGA): an Immeasurable Source of Knowledge. Contemp. Oncol. (Pozn) 19 (1A), 68–77. doi:10.5114/wo.2014.47136

CrossRef Full Text | Google Scholar

Vergoulis, T., Vlachos, I. S., Alexiou, P., Georgakilas, G., Maragkakis, M., Reczko, M., et al. (2012). TarBase 6.0: Capturing the Exponential Growth of miRNA Targets with Experimental Support. Nucleic Acids Res. 40 (Database issue), D222–D229. doi:10.1093/nar/gkr1161

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Xu, D., Huang, H., Cui, Y., Li, C., Zhang, C., et al. (2018). Detection of Dysregulated Competing Endogenous RNA Modules Associated with clear Cell Kidney Carcinoma. Mol. Med. Rep. 18 (2), 1963–1972. doi:10.3892/mmr.2018.9189

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Yu, H., Zhao, Z., and Jia, P. (2015). EW_dmGWAS: Edge-Weighted Dense Module Search for Genome-wide Association Studies and Gene Expression Profiles. Bioinformatics 31 (15), 2591–2594. doi:10.1093/bioinformatics/btv150

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., et al. (2021). clusterProfiler 4.0: A Universal Enrichment Tool for Interpreting Omics Data. The Innovation 2 (3), 100141. doi:10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiao, F., Zuo, Z., Cai, G., Kang, S., Gao, X., and Li, T. (2009). miRecords: an Integrated Resource for microRNA-Target Interactions. Nucleic Acids Res. 37 (Database issue), D105–D110. doi:10.1093/nar/gkn851

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, D., Wang, L., Pang, S., Cao, M., Wang, W., Yu, X., et al. (2021). The Functional Characterization of Epigenetically Related lncRNAs Involved in Dysregulated CeRNA-CeRNA Networks across Eight Cancer Types. Front. Cel Dev. Biol. 9, 649755. doi:10.3389/fcell.2021.649755

CrossRef Full Text | Google Scholar

Xu, J., Li, Y., Lu, J., Pan, T., Ding, N., Wang, Z., et al. (2015). The mRNA Related ceRNA-ceRNA Landscape and Significance across 20 Major Cancer Types. Nucleic Acids Res. 43 (17), 8169–8182. doi:10.1093/nar/gkv853

PubMed Abstract | CrossRef Full Text | Google Scholar

Yoshihara, K., Shahmoradgoli, M., Martínez, E., Vegesna, R., Kim, H., Torres-Garcia, W., et al. (2013). Inferring Tumour Purity and Stromal and Immune Cell Admixture from Expression Data. Nat. Commun. 4, 2612. doi:10.1038/ncomms3612

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, J., Liu, L., Xu, T., Xie, Y., Zhao, C., Li, J., et al. (2019). miRspongeR: an R/Bioconductor Package for the Identification and Analysis of miRNA Sponge Interaction Networks and Modules. BMC Bioinformatics 20 (1), 235. doi:10.1186/s12859-019-2861-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, M., Cheng, L., and Zhang, Y. (2020). Characterization of Dysregulated lncRNA-Associated ceRNA Network Reveals Novel lncRNAs with ceRNA Activity as Epigenetic Diagnostic Biomarkers for Osteoporosis Risk. Front. Cel Dev. Biol. 8, 184. doi:10.3389/fcell.2020.00184

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, W. Y., Zhang, M. M., Liu, C., Kang, Y., Wang, J. O., and Yang, X. H. (2019). Long Noncoding RNA LINC00473 Drives the Progression of Pancreatic Cancer via Upregulating Programmed Death‐ligand 1 by Sponging microRNA‐195‐5p. J. Cel Physiol 234 (12), 23176–23189. doi:10.1002/jcp.28884

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, X., Liu, J., and Wang, W. (2014). Construction and Investigation of Breast‐cancer‐specific ceRNA Network Based on the mRNA and miRNA Expression Data. IET Syst. Biol. 8 (3), 96–103. doi:10.1049/iet-syb.2013.0025

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, Y., Zhou, B., Pache, L., Chang, M., Khodabakhshi, A. H., Tanaseichuk, O., et al. (2019). Metascape Provides a Biologist-Oriented Resource for the Analysis of Systems-Level Datasets. Nat. Commun. 10 (1), 1523. doi:10.1038/s41467-019-09234-6

PubMed Abstract | CrossRef Full Text | Google Scholar