Transcriptome profiling and co-expression network analysis of lncRNAs and mRNAs in colorectal cancer by RNA sequencing

Background Long non-coding RNAs (lncRNAs) are widely involved in the pathogenesis of cancers. However, biological roles of lncRNAs in occurrence and progression of colorectal cancer (CRC) remain unclear. The current study aimed to evaluate the expression pattern of lncRNAs and messenger RNAs (mRNAs). Methods RNA sequencing (RNA-Seq) in CRC tissues and adjacent normal tissues from 6 CRC patients was performed and functional lncRNA-mRNA co-expression network was constructed afterwards. Gene enrichment analysis was demonstrated using DAVID 6.8 tool. Reverse transcription quantitative polymerase chain reaction (RT-qPCR) was used to validate the expression pattern of differentially expressed lncRNAs. Pearson correlation analysis was applied to evaluate the relationships between selected lncRNAs and mRNAs. Results One thousand seven hundred and sixteenth differentially expressed mRNAs and 311 differentially expressed lncRNAs were screened out. Among these, 568 mRNAs were up-regulated while 1148 mRNAs down-regulated, similarly 125 lncRNAs were up-regulated and 186 lncRNAs down-regulated. In addition, 1448 lncRNA–mRNA co-expression pairs were screened out from 940,905 candidate lncRNA-mRNA pairs. Gene enrichment analysis revealed that these lncRNA-related mRNAs are associated with cell adhesion, collagen adhesion, cell differentiation, and mainly enriched in ECM-receptor interaction and PI3K-Akt signaling pathways. Finally, RT-qPCR results verified the expression pattern of lncRNAs, as well as the relationships between lncRNAs and mRNAs in 60 pairs of CRC tissues. Conclusions In conclusion, these results of the RNA-seq and bioinformatic analysis strongly suggested that the dysregulation of lncRNA is involved in the complicated process of CRC development, and providing important insight regarding the lncRNAs involved in CRC. Supplementary Information The online version contains supplementary material available at 10.1186/s12885-022-09878-6.


Background
Colorectal cancer (CRC), including colon cancer and rectal cancer, is one of the most common malignant tumors. The progression of CRC is a multi-step process and can be categorized into four stages (Dukes staging system) based on the extent of tumor invasion [1,2]. According to the latest global cancer statistics 2018, CRC has risen to the rank third of malignant tumors and when it comes to the cancer mortality, CRC ranks second, ahead of the stomach cancer and liver cancer [3]. An upward trend in morbidity rate was observed in China, rank fourth in men and third in women [4]. In previous studies, several molecular mechanisms such as the oncogene p53, APC [5], gene methylation [6,7] Open Access † Mingjie Li and Dandan Guo contributed equally in this study.
*Correspondence: wyaslyy@126.com and non-coding RNA regulation [8][9][10] were shown to contribute to the occurrence and development of CRC. Additionally, high-throughput screening of the expression changes between CRC tumor tissues vs. adjacent normal tissues revealed a lot of diagnostic and prognostic biomarkers [11][12][13]. However, the comprehensive understanding of the progression and prognosis of CRC patients remains a formidable challenge due to the genetic heterogeneity and complex genomic alterations found in this cancer [14,15].

Sample information
Twelve samples (harboring 6 CRC tissues and 6 paired adjacent normal tissues) used in RNA-Sequencing (RNA-Seq) were collected from six Chinese patients who were diagnosed with stage II b or IIIb CRC. The raw sequencing data is secondary analyzed, and the 6 pairs of CRC tissues were divided into two groups (group 1 and group 2, corresponding to clinical stage II and III, Table  S1) based on their clinical stages. 60 pairs of CRC tissues used in expanded validation cohort were collected at Fujian Provincial Hospital from June 2015 to August 2017. We received the written informed consents from patients, and this study was reviewed and approved by the ethics committee of Fujian Provincial Hospital (No. K2012-009-01).

Library preparation and sequencing
Total RNA was extracted from tissues with TRIzol as per the manufacturer's protocol (Invitrogen, USA). A total of 3 μg RNA per sample was used as initial material for the RNA sample preparations. Ribosomal RNA was removed and the sequencing library was generated using Hieff NGS ® MaxUp rRNA Depletion Kit (Yeasen, China) following manufacturer's recommendations. Libraries from CRC tissue and adjacent normal tissues were analyzed on a single Genome Analyzer IIx lane (Illumina, USA) using 115 bp sequencing. Raw RNA-seq data were filtered by fastx_toolkit-0.0.14 (http:// hanno nlab. cshl. edu/ fastx_ toolk it/) according to the following criteria: 1) reads containing sequencing adaptors were removed; 2) nucleotides with a quality score lower than 20 were trimmed from the end of the sequence; 3) reads shorter than 50 were discarded; and 4) artificial reads were removed.

Reads mapping and transcript abundance estimation
The H. sapiens reference genome (GRCh37) was downloaded in Ensemble database (Human-download DNA sequence). The original transcriptome reads sequenced were aligned against the reference genome using TopHat v1.3.1, and bam (binary SAM) file alignment results were output. The pre-built GRCh37 index was downloaded from the TopHat homepage and used as the reference genome. The aligned read files were processed by Cufflinks v1.0.3, which uses the normalized RNA-seq fragment counts to measure the relative abundances of transcripts. The unit of measurement is Fragments Per kilo-base of exon per million fragments mapped (FPKM). Confidence intervals (CI) for FPKM estimated were calculated using a Bayesian inference method.

Differentially expressed gene testing
The downloaded Ensemble GTF file (GRCh37) was submitted to Cufflinks v2.2.1 along with the original alignment (SAM) files produced by TopHat. Cufflinks re-estimates the abundance of the transcripts listed in the GTF file using alignments from the SAM file and concurrently tests for differential expression with the default parameters. Only the comparisons with q_value less than 0.05, |log 2 FC| ≥ 1, Max FPKM (N, T) ≥1 and test status marked as "OK" in the Cufflinks output were regarded as differential expression. Meanwhile, since we hope to study the overall gene expression in colorectal cancer tissues, genes expressed separately in stage II or III respectively were excluded, which may better reflect the commonality of this sequencing.

Functional enrichment analysis and lncRNA-mRNA co-expression network
DAVID v 6.8 is a web-based functional annotation tool. The unique lists of differentially expressed genes and all the expressed genes (FPKM> 0) were submitted as the gene list and background list, respectively. The cut-off value of the False Discovery Rate (FDR) was 0.05, and only the results from the Gene ontology analysis (GO) and Kyoto Encyclopedia of Genes and Genomes pathway analysis (KEGG) were selected as functional annotation categories. Pearson correlation analysis was used to estimate co-expression relationships between lncRNAs and mRNA. A set of co-expressed lncRNA-related genes were filtered with a Pearson coefficient threshold of 0.95 and p < 0.01. Cytoscape 3.2.1 tool was applied to construct the lncRNA-mRNA network.

Validations of differentially expressed lncRNAs
The differentially expressed lncRNAs were verified by Reverse transcription quantitative polymerase chain reaction (RT-qPCR) using SYBR ® Premix Ex Taq ™ reagent (TAKARA, Japan) on ABI ViiA ™ 7 (Applied Biosystems, USA) per the manufacturer's instructions. The selection criteria for validation included, 1) The gene expression level was relatively high for detection; 2) The gene expression pattern was consistent in the 6 tumor tissues (all higher than/all lower than the matched normal tissues); and 3) Higher differential expression ratio in cancer/normal tissues. Primer sequences were listed in Table S2. In addition, the correlationship between MIR4435-1HG (an up-reguated lncRNA) and COL4A1, SATB2-AS1 (a down-regulated lncRNA) and SGK2, were confirmed using Pearson correlation analysis in 60 samples collected. Gene expression levels were normalized to glyceraldehyde-3-phosphate dehydrogenase (GAPDH). All the RT-qPCR reactions were performed in triplicate. Expression data was expressed as mean ± SD and P < 0.05 was considered statistically significant.

Characterization of sequencing and mapping
All 12 samples were subjected to massively parallel paired-end cDNA sequencing. On average, 16 Gb (14.2-19.6Gb) datum were obtained from CRC tissues and adjacent normal tissues. We used TopHat tool to align the reads to the Ensemble reference human genome GRCh37. The proportion of reads that mapped to the Ensemble reference genes ranged from 82.7 to 90.9% for the twelve samples. Correlation coefficients of expression levels between different samples are shown in Fig. 1. After grouping the samples, the scatter relationship between tumor tissues and normal tissues was shown in Fig. 2. The average coverage of our sequencing depth was approximately 108(94-137) times of human transcriptome and the details of the mapping results were listed in Table 1. This sequencing received 18,489 mRNAs and 9753 lncRNAs, accounting for 89 and 70% of annotated genes (mRNA:20730, lncRNA:13869). The mRNA and lncRNA expression level of FPKM≥1 were 12,773 and 1669, accounting for 62 and 12% respectively ( Table 2).

Differentially expressed lncRNAs and mRNAs in CRC tissues
FPKMs were calculated for normalization of the expression level of lncRNAs and mRNAs. 1716 differentially expressed mRNAs and 311 differentially expressed lncRNAs were found in 6 pairs of CRC tissues vs. adjacent normal tissues. Among these, 568 mRNAs were up-regulated while 1148 mRNAs down-regulated, similarly 125 lncRNAs were up-regulated while 186 lncRNAs down-regulated. In group I, 903 differentially expressed mRNAs and 153 differentially expressed lncRNAs were screened out. Among them, 296 mRNA were up-regulated and 607 mRNAs down-regulated while 56 lncRNAs were up-regulated and 97 lncRNAs down-regulated. In group II, 566 differentially expressed mRNAs and 126 differentially expressed lncRNAs were found. Among them, 174 mRNAs were up-regulated and 392 mRNAs    down-regulated while 37 lncRNAs were up-regulated and 89 lncRNAs down-regulated (Fig. 3).

Functional enrichment analysis and mRNA-lncRNA co-expression network
We constructed a co-expression network of the dysregulated lncRNAs and mRNAs. 1448 lncRNA-mRNA co-expression pairs were screened out from 940,905 candidate lncRNAs and mRNAs (Fig. 4). GO analysis and KEGG revealed that these co-expression mRNAs were closely correlated with cell adhesion, collagen adhesion, cell differentiation and formation of extracellular matrix organization, and mainly enriched in fatty acid degradation, butanoate metabolism and PI3K-Akt signaling pathway (Table S3 and S4). It is public knowledge that PI3K-Akt signaling pathway had a profound effect on CRC progress. Naturally, as depicted at Fig. 5, we performed the mapping analysis for PI3K-Akt signaling pathway. According to co-expression analysis, many lncRNAs were enriched on important nodes of the PI3K/ Akt signaling pathway (Fig. 5, FDR < 0.05).

Discussion
As one of the most malignant tumors, CRC is becoming a great social burden in the world. It was reported that there would be 18.1 million new cancer cases and 9.6 million new cancer deaths worldwide in 2018, among which CRC ranked the 4th in incidence and the 2nd in mortality, seriously endangering people's healthy and property safety [3]. Improvement of this severe situation mainly depends on identification of biomarkers for early diagnosis and development of therapies for CRC treatment. Here, the differentially expressed mRNAs and lncRNAs were screened out by using RNA-seq for 6 pair of CRC tissues. Based on the sequencing results, differential lncRNA-mRNA co-expression network and gene list enrichment analysis revealed the potential regulatory roles of lncRNAs in the development of CRC. Finally, the expression patterns of 10 lncRNAs, as well as correlativity between selected lncRNAs and mRNAs, were detected in an expanded tissues sample set to verify the reliability of RNA-seq. Protein-coding genes make up only 1.5-2% of the human genome, while the non-coding genes consist of almost 98%. LncRNA, a class of RNA with length more than 200 bp, is now attracting wide attention. It was once considered sort of transcriptional noises due to deletion of protein-coding regions. But now, accumulating evidences showed that lncRNAs were generally involved in many human cancers, such as glioma, gastric cancer, breast cancer, liver cancer, endometrial cancer and so on [16]. However, the underlying functional roles and mechanisms of most lncRNAs remain elusive. In last decade, a lot of lncRNAs were identified for early diagnosis and prognosis monitoring of CRC. Through the bioinformatics database and large-scale verification, Xu et al., identified the differentially expressed lncRNA-SNHG11 as an appropriate candidate for early diagnosis of CRC patients [17]. A prognostic risk formula including three lncRNAs (LINC01602, AP003555.2 and AP006284.1) was successfully established to evaluate the prognosis of CRC patients, these three-lncRNAs signature presented a great potential of being the independent biomarker for the prognosis of CRC patients [18]. LINC01133 was detected down-regulated in CRC tissues and Kaplan-Meier survival analysis revealed patients with high-LINC01133 had a better survival outcome [19]. Encouragingly, several lncRNAs mentioned above were included in our differentially expressed genes set, which also confirmed the effectiveness of the current sequencing. Based on those studies, we also hope to further analyze the impact of these dysregulated lncRNAs on early diagnosis and prognosis of CRC patients in the future. BBOX1-AS1, an aberrant expressed anti-sense lncRNA depicted in this study, presented increasing status in CRC cell lines. Knockdown of BBOX1-AS1 inhibited the progression of CRC cell, including cell proliferation, migration, invasion and conversely promoted apoptosis of tumor cells by sponging miR-361-3p/SH2B1 regulatory axis [20]. Consistent with our study, lncRNA DPP10-AS1 was shown to be significantly decreased in CRC tumor tissues, along with changes in colon cancer stem cell properties. In vitro and in vivo studies uncovered that DPP10-AS1, worked as a tumor suppressor, inhabited proliferation, migration and invasion but facilitated apoptosis of Fig. 4 LncRNA-mRNA co-expression network. The red nodes in the network represented lncRNAs while the blue nodes were co-expressed mRNAs. LncRNAs and mRNAs with correlation coefficients greater than or equal to 0.95 were selected, and then a network was constructed using Cytoscape 3.3.1 tool CRC cells through the potential miR-127-3p/ADCY1 axis [21]. Another lncRNA MIR503HG in the validation set of RT-qPCR was widely known for its tumor suppressor-like role in CRC. Rescue test uncovered that overexpression of miR-107 reversed the antitumor effect of MIR503HG on CRC cells by potential mechanism of epithelial-mesenchymal transformation [22]. It was worth mentioning that MIR503HG was decreased in tumor tissues and cells in their study, which was contrary to the current study (Fig. 5B). On one hand, the sample set of this study might be insufficient. As was well-known, with the increase of sample size, the average expression level of the gene in the population tended to its true level. On the other hand, as mentioned above, there existed large differences in tumor heterogeneity of CRC patients, and even different parts of the same piece of tissue are expressed differently due to cell composition and genetic heterogeneity.
Drug resistance was one of the main obstacles in the therapy of CRC, and understanding of chemoresistance will greatly improve the treatment and prognosis of patients. Accumulating evidences suggested that lncR-NAs might play significant roles in the chemoresistance. In vivo and vitro studies validated that lncRNA-HAND2-AS1 inhabited the proliferation and 5-FU resistance in 5-FU-resistant CRC tumor cells [23]. Targeted lncRNA therapy has a profound prospect and may be an alternative option for CRC patients accompanied by chemotherapy resistance.
Recently, RNA-seq can be used to distinguish differences in gene expression between different time points and different groups, especially transcriptome differences between normal and tumor tissues. RNA-seq is characterized by high throughput and high repeatability and has been applied to the transcriptome study of animals and plants [24,25]. Meanwhile, RNA-seq also provides a new means for accurate and early diagnosis of clinical diseases, and provides a feasible prospect for improving clinical treatment [26]. In the current study, RNA-seq was performed in 6 pairs of CRC tissues and differentially expressed lncRNAs were identified, which laid a foundation for our subsequent research on diagnostic biomarkers and expanded the understanding of the pathogenesis of CRC. However, the current RNA-seq technology also faced some technical difficulties. Ribosomal RNA and short-reads RNA were removed prior to sequencing in this study, and thus limiting the number of reads and the accuracy of these RNA expression levels, which might introduce potential errors and ultimately affected experimental results. Expression data could not accurately reflect gene expression levels because many genes possessed various isforms and some are rarely detected. Therefore, there still existed uncertainty about genes whose expression was significantly altered when using the cuffdiff tool. Beyond all doubt, the emergence of new technologies is inevitably in possession of disadvantages as well as advantages, and RNA-Seq does the same. But certainly, it is worth expecting that there will be more strategies to reduce these shortcomings in the future.
Based on the co-expressed mRNAs, differentially expressed lncRNAs were enriched in PI3K-Akt signaling pathway, which was known as a complicated pathway in the progress of CRC (Fig. 5). The high-frequency activation of PI3K signaling pathway has been validated sharing a close relationship with occurrence and development of CRC, contributing important values for the diagnosis and therapy of CRC [27]. Studies have confirmed that abnormal activation of PI3K occurred in the early, advanced and metastatic stages of CRC [28][29][30], and experiments in vivo and in vitro clearly prompted that the regulation of signal transduction system by targeted-PI3K signaling pathway had a certain improvement effect on CRC therapies [31]. In the current study, according to the construction of Fig. 6 The expression values of 10 lncRNAs in an expanded 60 samples cohort. a, b, c, d and e The expression levels of 5 up-regulated lncRNAs. f, g, h, i and j The expression levels of 5 down-regulated lncRNAs. All the RT-qPCR tests were performed in triplicate and data were shown as means ± SD. Lower Δct values indicated higher expression. ***P < 0.001; ****P < 0.0001. P < 0.05 was considered statistically significant lncRNA-mRNA co-expression network, these lncRNAs were found to be related to cell proliferation, cycle, DNA repair event and even pathological epithelialmesenchymal transformation, which played a crucial role in malignant metastasis and drug resistance of tumor cells. Analysis of the co-expression of lncRNA and mRNA was helpful to predict the role of lncRNA in the development of various diseases, especially in cancer, and laid a foundation for revealing their mechanism of carcinogenesis and progression. Besides, CASC15 and RP11-334E6.12 were co-expressed with GF and RTK, pivotal upstream molecule of PI3K, might make an influence on the PI3K/Akt signal and ultimately changed the biological behaviors by regulating the metabolism and proliferation of CRC cells (Fig. 5). Beyond that, as conspicuously demonstrated at Fig. 5, more than 10 aberrantly expressed lncRNAs were enriched in extracellular matrix (ECM) pathway, accompanied by PI3K/Akt pathway, might have significant impacts on tumor migration and metastasis. We firmly believed that aforementioned results would shed a light on subsequent works. CRC is a common malignancy worldwide, with a particularly high incidence in China. The etiology, pathophysiology, and underlying molecular mechanisms of CRC remains largely unknown, and further study on roles of candidate lncRNAs needs to be fully carried out. LncRNA is becoming a research hotspot in CRC, and we hope that this study will not only expand the library of lncRNA markers in CRC, but further revealing the mystery of pathogenic mechanisms. Dysregulated lncRNAs are expected to be advisable biomarkers for early/advanced diagnosis and prognosis. We are confident that lncRNA will become the new favorite of targeted biotherapy.

Conclusions
In summary, we comprehensively investigated the RNA profiles of CRC tumor and adjacent normal tissues through the RNA-seq technology and identified differentially expressed mRNA and lncRNAs. The abnormal expression of lncRNAs were verified and furthermore, related lncR-NAs were fortunately enriched in the PI3K-Akt signaling pathway, which undoubtedly laid a solid foundation for the future work on the roles of lncRNAs in CRC. Fig. 7 The expression correlation between selected lncRNAs and relative mRNAs. a and (b) The correlationships between MIR4435-1HG and COL4A1, SATB2-AS1 and SGK2 based on RNA-sequencing. c and (d) The correlationships between MIR4435-1HG and COL4A1, SATB2-AS1 and SGK2, based on RT-qPCR in 120 samples. Y axes in Fig. 7A and B were the FPKM while C and D were Δct