Copy number variations primed lncRNAs deregulation contribute to poor prognosis in colorectal cancer

Multi-omics integration analysis

The expression profiles of PCGs, CNVs and 450k methylation were combined with prognosis status based on univariate cox proportional hazards model. Eventually, 2118 genes, 5015 CNV regions and 7083 CpG sites were obtained with the significant threshold p<0.05. Subsequently, 289 CRC patients sharing in three omics were classified into 5 molecular subtypes using iCluster. Five kinds of subtypes C1, C2, C3, C4, C5 consisted of 23, 44, 23, 60, 139 CRC patients respectively (Table 1). C3 group was harboring the worst survival probabilities (overall survival), and C2 subtype owned the optimal prognosis status (p=0.0056) (Figure 1A). Progression free survival (PFS) among five subtypes was almost consistent with overall survival, similarly, C3 group got hold of the worst survival probabilities (PFS) (Supplementary Figure 1). The distribution of TNM stage among five subtypes clearly showed that the proportion of advanced patients (stageIII+stageIV) in C3 group was much larger than other four subtypes accounting for about 60% (Supplementary Figure 2). Simultaneously, we described the spectrum of gene mutation status among five subtypes, and selected the top 10 genes with the highest mutation rate in each subtype. There are total 20 mutation related genes acquired and this phenomenon highlighted that the most common mutation related genes in five subtypes had higher coincidence. We made further observations revealing that these 20 genes were differently distributed in five subtypes and the mutation frequency of the same gene in different subtypes is also distinguished (Figure 1B). Obviously, the frequency of gene mutations in C5 subtype was significantly higher than other subtypes, and even in some C5 subtypes samples, these 20 gene mutations almost existed. In addition, TTN, APC, TP53 and KRAS mutations were especially more common than other genes, which was highly shining upon that these four genes mutations may imply pivotal role in carcinogenesis. These findings suggested that our strategy for subtypes classification based on the expression profiles of PCGs, 450k methylation and CNVs data could predict distinct prognostic situation of CRC and feature certain regulatory relationship among genomic, epigenomic and transcription level.

Differentially expressed lncRNAs and mRNAs across subtypes

DE-lncRNAs in five subtypes and carcinoma and normal tissues were 2668, 2890, 2838, 2634, 2265, 2718, respectively (Table 2). DE-PCGs in five subtypes and carcinoma and normal tissues were 3702, 3989, 3934, 3647, 3013, 3642, respectively (Table 2 and Figure 2G). It was evidently exhibiting that C5 subtype was holding less differential mRNAs and lncRNAs, on the contrary, C2 and C3 subtypes were holding more differential mRNAs and lncRNAs than other subtypes. A total of 4253 DE-lncRNAs and 5808 differential expression of PCGs were obtained. The upregulated lncRNAs were remarkably more than down-regulated lncRNAs (Figure 2A–2F). Whereafter, we downloaded lncRNAs closely related to disease from the LncRNA Disease and Lnc2Cancer database and 611 lncRNAs were obtained, which were compared to 4253 subtypes specific lncRNAs (Figure 2H). There were 120 lncRNAs contained in previous DE-lncRNAs and the significance of these lncRNAs was tested by hypergeometric test (p=0.003929) (Supplementary Table 1). In order to understand the clustering situation of these DE-lncRNAs, gene set enrichment analysis (GSEA) was employed based on fold change of each lncRNA in five molecular subtypes and carcinoma and normal tissues (Figure 3A–3F). The results indicated that these DE-lncRNAs were always clustered in the gene set with larger difference multiples. There was almost no intersection of lncRNAs between all five molecular subtypes and carcinoma and normal tissues. Most lncRNAs were merely enriched in one subtype, such as, subtype 1 having 237 lncRNAs; subtype 2 having 202 lncRNAs; subtype 3 having 226 lncRNAs; subtype 4 having 102 lncRNAs; subtype 5 having 53 lncRNAs and carcinoma and normal tissues having 61 lncRNAs, respectively. The remaining lncRNAs were focused on two kinds of subtypes or three kinds of subtypes. It was obviously seen carcinoma and normal tissues shared a great deal of lncRNAs with other subtypes (Figure 3G). This phenomenon highly suggested that these identified lncRNAs are distinctively enriched in specific molecular subtypes.

Co-expressing modules between lncRNAs and PCGs based on WGCNA

Co-expression between lncRNAs and PCGs was analyzed by hierarchical clustering. We defined the samples with distance more than 150000 as outlier samples for screening, and finally got 492 CRC cases for co-expressing analysis (Figure 4A). We accepted 3 as the final soft threshold in order to ensure the network belonging to non-scale property (Figure 4B, 4C), which can effectively strengthen strong correlation and weaken weak correlation or negative correlation. Next, the expression matrix was transformed into adjacency matrix, and then the adjacency matrix was transformed into topological overlap matrix (TOM) for reduced noise and false correlation. According to the dynamic tree cut criterion linked with TOM structure, the threshold for module partitions was constructed at 30. After determining the initial modules, the next was to establish new modules presenting as different colors, including calculating specific eigengenes of each module; clustering analysis of each module; grouping close modules into one new module with height = 0.25, deepsplit = 2, minModuleSize = 30. A total of 27 modules (Figure 4D) were eventually obtained and the module of grey was signatured with genes which were unable to be classified into other modules. P value represents the significant aggregation degree and fold change represents the aggregation multiple. There was no significant difference between lncRNAs and PCGs among 27 modules (Table 3).

For the interest of finding biologically associated modules, we need to correlate modules with external information, for instance clinical data of array information and functional enrichment analysis of gene information. Therefore, we executed correlation analysis by combining patients’ gender, age, height, weight, BMI, T, N, M, stage with each module (Figure 5A). What could be clearly indicated in the correlation analysis that all modules except white, pink, black, darked were related to at least one phenotype. Besides, four modules: tan, blue, yellow, magenta were finally selected with more than three kinds of phenotypic correlations for functional enrichment analysis and the significant enrichment pathways p<0.05 were picked up. These four kinds of modules were clustered in distinguished pathways and the cross-talk between these four modules was missing (Figure 5B), suggesting each module may drive in different biological function. Tan module was associated with 18 pathways (Figure 6A, Supplementary Table 2) and neuroactive ligand-receptor interaction and calcium signaling pathway were in the highest flight in this module. MAPK signaling pathway also occupied superior gene ratio in this module and had been proved holding significant position in CRC-promoting process. Only two pathways were enriched in yellow module, of which microRNAs in cancer may possess outstanding role (Figure 6B, Supplementary Table 3). Besides, blue module enriched in 47 pathways (Figure 6C, Supplementary Table 4), and the remarkable pathways included neuroactive ligand-receptor interaction, PI3K-AKT signaling pathway and MAPK signaling pathway. Genes in magenta module mainly enriched in Wnt signaling pathway, TGF-beta signaling pathway and DNA replication which were implicated with classical undertakings in carcinogenesis, especially in CRC (Figure 6D, Supplementary Table 5). Aberrant activation of WNT signaling pathway markedly increases the expression of β-catenin, thus promoting normal colon epithelial cells to infinitely proliferate and cancerate [25–27]. Moreover, TGF-beta/smad signaling pathway can advance invasion and metastasis of CRC by means of immune suppression, angiogenesis, augmenting interaction between tumor cells and extracellular matrix [28–30]. Overall, DE-lncRNAs and PCGs co-expression modules constructed by WGCNA were enriched in various pathways implicating important functions in regulating CRC progression.

Identification of CNVs associated lncRNAs in CRC

CNV is a frequent form of genomic structural changes, which is closely related to occurrence and deterioration of tumors. We profiled the spectrum of CNV-related lncRNAs on the whole genome in CRC patients The frequency of copy number deletions were largely more than amplifications, suggesting copy number deletions may predominantly implicate in CRC, and majority deletions were found on chromosome 8, in contrast, the amplifications were mainly concentrated on chromosome 20 (Figure 7A). The correlation between the expression profiles of lncRNAs and CNVs was demonstrating a positive relationship and the distribution was significantly higher than the random, p<1e-16 (Figure 7B). Moreover, the frequently changed regions in CRC genome were identified on the basis of GISTIC algorithm. A genome-wide view of the CNVs was shown in Figure 7C and the overall frequent change regions were complex. The focal amplification events (indicated by dashes), surpassing the significance threshold (green line), mainly distributes among chromosomes 8, 11, 12p13, 16q12, 20p11, and the focal deletions events were concentrated on 1p33, 3q26, 4p16, 5q11, 5q22, 20p12. Similarly to previous condition, frequent copy number deletions of lncRNAs were largely more than the copy number amplifications in the whole genome. This finding emphasized the impact of copy number deletions of lncRNAs in CRC, and it may be significantly related with prognosis of CRC patients.

In order to deeply explore the effect of CNV on transcription of lncRNAs, we screened out 17 lncRNAs with more than 7% of the frequency of CNV in 452 CRC samples (Supplementary Table 6). In the light of these 17 lncRNAs, we analyzed the expression difference of each lncRNA between the samples with copy deletion or copy amplification and normal group. A total of 12 lncRNAs were obtained according to the criterion in which the expression in each group was greater than 0. There are four lncRNAs: CASC11, HM13-AS1, ABALON, NKILA highly expressed in amplification group than diploid group (p<0.05), and FAM87A, LOC101927752, KBTBD11-OT1, LOC100287015, LOC101929066, these five lncRNAs were low-expressed in copy number deletion group than normal copy number group (p<0.05) (Figure 8). All these findings indicated that deletions or amplifications of lncRNAs may closely conduct the expression level of lncRNAs. However, these differentially expressed lncRNAs are deserved to explore for the potential clinical significance based on large-scale samples supporting.

LncRNAs-based prognostic biomarkers in CRC patients

As we all know, lncRNAs are holding diving or inhibiting role in tumors and increasing numbers of lncRNAs have been conducted as prognosis associated biomarkers in CRC. In our research, we wondered whether the expression profiles of lncRNAs regulated by copy number deletions or amplifications can impact survival status of CRC patients. Based on the previous five molecular subtypes, we systematically analyzed the copy number alternations of these DE-lncRNAs among five subtypes, and then, 104 lncRNAs conforming to the following criterions were obtained (Supplementary Table 7): CNV change frequency among each sample is more than 0.1%; holding difference at least five subtypes; positively related with CNVs. We then wiped off the samples of whom the expression of each lncRNA was under zero, and subsequently correlated lncRNAs with disease-free survival. We finally acquired 10 lncRNAs with significant prognosis (p<0.05) based on the quartiles method (Figure 9, Supplementary Table 8). LOC101927604, CASC15, CLDN10-AS1, C14orf132, LMF1 and LINC00675 were differentially expressed among all six subtypes (including normal group), and the expression levels of LOC105377267, LINC-PINT, CCDC144NL-AS1 and LOC284454 were differential in five subtypes (Table 4).

In addition, we verified the prognostic value of identified 10 lncRNAs in GEO (Gene Expression Omnibus) dataset (GSE39582) based on GPL570 platform (https://www.ncbi.nlm.nih.gov/geo/). There were about 8 lncRNAs successfully annotated to probes existing in GEO platform and only 5 lncRNAs, CASC15, LINC-PINT, C14orf132, LMF1 and CCDC144NL-AS1 holding significant prognosis difference (p<0.05) (Figure 10).

Data preparation and processing

Involved original data were downloaded from official TCGA data portal (https://tcga-data.nci.nih.gov), including RNA sequencing data, 450k methylation data, CNV data, DNA mutation data and clinical information. Subsequently, series of managements were performed, for instance, RNA sequence data: primary data were downloaded in “counts” form, and the processed FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) data were transformed to TPM (TranscriptsPerKilobase of exonmodel per Million mapped reads). The gene expression profiles of 458 CRC patients and 41 normal cases were acquired. Then, we defined lincRNA, sense-intronic, sense-overlapping, antisense, processed-transcript, 3primer-overlapping as lncRNAs based on the genecode file. Meanwhile, gene type belonging to protein-coding was assigned to PCGs; Next, administrations of 337 CRC patients holding 450k methylation data were involved with removing of NA probes, cross-reactive CpG sites [18] and unstable CpG sites existing in the sex chromosomes and single nucleotide sites; Then, CNV data and single nucleotide mutation data were pictured, with the removed germline difference and processed by “mutect software”, respectively. Eventually, we also obtained corresponding clinical information of 458 CRC patients. The patients with following time less than 30 days were removed. Underlying these prepared data, we employed R package-iCluster to integrate PCGs, methylation and CNVs for subtypes clustering. Number of clusters was set to 5 and 50 iterations were set to identify stable samples clusters with default parameters [19]. The R code can be downloaded at http://www.mskcc.org/mskcc/html/85130.cfm.

Differentially expressed lncRNAs and mRNAs across five subtypes

The R package-DEseq2 was claimed for determining these differentially expressed lncRNAs and mRNAs between five subtypes and normal tissues specimens, with the criteria of fold change>2, FDR<0.05. Genes with count<1 were rejected among the whole genome profiles. The basic data structure consisted of two tables, including countData and colData. Gene set enrichment analysis (GSEA) was furtherly applied to distinguish the distribution of DE-lncRNAs and the difference multiple was used for ranking.

Weighted gene co-expression network analysis (WGCNA) to identify phenotypes-related modules

WGCNA is often regarded as efficient methodology applying in picking up co-expressed gene modules and identifying potential therapeutic targets associated with clinical phenotypes [20]. In present research, based on the differential expression profiles of lncRNAs and PCGs, co-expression network was constructed using WGCNA. In order to make network conform to non-scale signature, the soft threshold was defined as 3. The non-scale signature means that log(k) is negatively correlated with log(P(k)), and the correlation coefficient is greater than 0.8. The possible modules were identified with the average-linkage clustering and dynamic tree cut on the ground of height=0.25, deepSplit = 2, minModuleSize = 30. Considering the clinical information downloading from TCGA, the correlation between each module and corresponding clinical information was calculated. Immediately, annotations of potential pathways were targeted in order to screen out biologically interested modules based on the userListEnrichment function included in package WGCNA [21].

Copy number expression profiles of whole genomic lncRNAs

GISTIC (Genomic Identification of Significant Targets in Cancer) algorithm is responsible for identifying variant regions that are more prone to drive cancer pathogenesis [22]. It can visualize regions in the genome manifesting amplifications or deletions across thousands of samples. G-score is allocated to each alternational region to evaluate the amplitude of aberration and the frequency of occurrence [23]. Here, GISTIC2.0 software was applied to define the prepared CNVs profiles of all genes underlying 452 CRC samples, and lncRNAs-related CNVs expression profiles were extracted. We defined copy number>1 or <-1 as copy number amplifications or deletions respectively. False Discovery Rate (FDR) q-values were assigned to each alternational region. “Peak regions”, also known as significantly aberrant regions, point to the greatest frequency and amplitude of aberrations. The “peak regions” aim at determining whether the signal is primarily due to broad events (longer than half a chromosome arm), focal events, or significant levels of both [23, 24].

CNV-related lncRNAs acting as prognostic biomarkers in CRC patients

CNVs profiles of DE-lncRNAs among five subtypes were described. Associated lncRNAs were prepared with the criterions as following: 1) the identified lncRNAs must be positively correlated with CNVs (correlation coefficient>0); 2) the identified lncRNAs should have significant difference across at least 5 subtypes (number of DEsubtype≧5); 3) the identified lncRNAs hold copy number change ratio more than 0.1% in all involved cases (CNV frequency>0.001); 4) prognostic-related lncRNAs were marked with p<0.05 (samples were divided into four groups based on quartile method or two groups based on median level). The selected patients must manifest expression value of previous marked lncRNAs greater than 0.

Statistical analysis

All statistical strategies in our research were based on R 3.4.3 (https://cran.r-project.org/bin/windows/base/old/3.4.3/). All workings are dependent on default parameters except to special notes.