Long non-coding RNAs as pan-cancer master gene regulators of associated protein-coding genes: a systems biology approach

Data acquisition

Gene expression data was obtained from TCGA data portal (https://cancergenome.nih.gov/ downloaded in the year 2016. We selected UNC_IlluminaHiseq_RNAseqV2 platform only to avoid heterogeneity of data and normalized level-3 sample data sets. A total number of samples was 5,601 from 15 different cancer types (primary cancers) and 601 corresponding normal samples. Each sample had associated expression values for 20,531 protein-coding genes.

Significant differential expression analysis of each gene

We analyzed significant differential expression by using unpaired t-test with unequal group variances (Welch approximation) in Multiple Experiment Viewer (MeV) version 4.9.0 (http://mev.tm4.org/) and Linear Model for Microarray (LIMMA) as a package in RStudio program (https://www.r-project.org/). For assessing significant differential expression in the t-test, the default cut-off was a critical p-value of 0.01 and in LIMMA, double filtration approach with the alpha value 0.01 and fold change >2 cut-offs was used. Combining these two results in order to increase accuracy and eliminate false positives, significantly differentially expressed common genes were enlisted using Venn diagram.

Identification of up- and down-regulated genes

Using our significantly differentially expressed common genes list, we applied Comparative Marker Selection version 10.1 (http://software.broadinstitute.org/cancer/software/genepattern/modules/docs/ComparativeMarkerSelection/10) to identify up- and down-regulated genes between the two classes of samples (cancer and normal) by t-test (median centered) at 10,000 permutations and a p-value cutoff of 0.01 with Bonferroni correction.

Correlation of protein and mRNA expression levels

We compared protein expression level (intensity obtained from immunohistochemistry results) published in “Human Protein Atlas (HPA)” (available from https://www.proteinatlas.org/, data accessed in 2016–2017) with mRNA expression level of each of the 70 significantly differentially expressed coding genes. As per Cancer Atlas in HPA, the staining intensity for each protein expression level from immunohistochemistry data is given in the range: high, medium, low and not detected, from the expert manual analysis. Side comparisons with normal samples are also presented.

Functional annotation

For functional annotation, we used GeneCoDis3 web server (http://genecodis.cnb.csic.es/). This web server takes a list of genes as inputs and annotates a gene set using Gene Ontology, KEGG Pathways, and other terms. Significance of annotations to gene set is determined by p-values computed using hypergeometric distribution or chi-square test. For our analyses, we examined only those annotations for which p-values obtained through hypergeometric analysis and corrected by the FDR method of Benjamini and Hochberg were provided for significance. Hypergeometric p-value of <or = 0.05 (default) was used. The results can be seen as the tags cloud where most significant terms that are enriched in the list of genes are presented as fonts of variable sizes. The font size of each tag depends on the number of genes enriched from our gene list that are present in a particular annotation.

Data acquisition

LncRNAtor database (http://lncrnator.ewha.ac.kr/expression.htm) is an online database which integrates expression profile, interacting protein, integrated sequence curation, evolutionary scores, and coding potential. TCGA, GEO, ENCODE, and modENCODE databases provide the relevant datasets. We selected differentially expressed lncRNAs from TCGA RNAseq platform, level-3 count data. Searching for the same cancer types as present in TCGA database, we found nine cancer types which were the same as the cancer types used in coding genes’ work and for which data from TCGA for comparing normal vs tumor was present. In this way, a total of 2,900 samples with 403 normal and 2,497 tumor samples were investigated. Total number of lncRNAs in lncRNAtor was 24,154. Significantly differentially expressed lncRNAs identified by the database authors using DESeq2 package for TCGA level 3 data with adj p-value 0.01 were enlisted from the database. Bioconductor software specifies about DESeq2 as follows: “It estimates variance-mean dependence in count data from high-throughput sequencing assays and tests for differential expression based on a model using the negative binomial distribution” ( Love, Huber & Anders, 2014). We then used Comparative Marker Selection version 10.1 to discriminate up- and down-regulated lncRNA genes between two classes of samples by t-test (median centered) at 10,000 permutations and a p-value cutoff of 0.01. We used lncRNAtor and GeneCodis3 to predict co-expressed mRNAs and their role in biological processes and molecular function in order to assign putative functions to uncharacterized lncRNAs by “guilt-by-association”.

Prediction of lncRNA-mRNA interaction and function

We used the human lncRNA-mRNA interaction database (http://rtools.cbrc.jp/cgi-bin/RNARNA/index.pl) (Terai et al., 2016) which is validated. It contains lncRNA-mRNA and lncRNA-lncRNA interactions between 23,898 lncRNAs and 20,185 mRNA (human) sequences obtained from the GENCODE Project (http://www.gencodegenes.org/releases/19.html). The first step in this pipeline; RACCESS program is used for extracting accessible regions in each RNA sequence, Step 2: tandem repeats are removed using the TANTAN program, Step 3: reverse complementary ‘seed matches’ are detected using LAST, Step 4: the binding energies between query and target RNAs from the seed matches are evaluated using INTARNA. In Step 5, the candidate interaction pairs are ranked according to the interaction energies (SUMENERGY and MINENERGY). Finally, the interaction site secondary structure with minimum interaction energy is predicted in each pair of RNA sequence using RACTIP. Using this database, we analyzed interactions of common significantly differentially expressed mRNAs and lncRNAs enlisted from our studies. Interaction energy (ranking was done using SUMENERGY score which is adequate for long RNA sequences such as lncRNAs and mRNAs), location of the interaction site on mRNA and joint secondary structures of the binding site were predicted. We used Excel for generating plots of SUMENERGY.

Integrated Network analysis

For experimentally verified transcriptional regulatory network (TFs-gene interaction) analysis, we retrieved data from Open-Access Repository of Transcriptional Interactions (ORTI, URL: http://orti.sydney.edu.au/about.html) database, which apart from other databases on mammalian transcription factors (TF) and their regulated target gene (TG) interactions, harbors Human Transcriptional Regulatory Interaction database (HTRIdb). LncRNAs are also included in the network interactions. In another network, we annotated differentially expressed genes as co-expressed, physically interacting and pathway genes using GeneMania for further analysis (Fig. S1). These two types of interaction data were imported into Cytoscape 3.5.1 (http://www.cytoscape.org/download.php). The network analyzer tool—using the directed network with TF as the source and gene as the target (directed TF to gene network)—was used to produce quantitative data for node degree, betweenness centrality, clustering coefficient among others.

Identification of up- and down- regulated genes

We further analyzed our data to identify which of these genes are up-regulated and which ones are down-regulated in cancers as compared to normal tissues. For this purpose, comparative marker selection in GenePattern suite of programs was used. We observed that out of these 70 commonly dysregulated genes, 48 were up-regulated and 4 genes down-regulated in all 15 types of tumors while remaining 18 were significantly differentially expressed genes that differed in their expression pattern in one or two tumor types only and therefore, still taken as showing a uniform pattern of expression (Fig. 3). These have been tabulated in Table 3.

Functional annotation

We annotated the functioning of these 70 common significantly differentially expressed genes using GeneCodis3 web-server implementing Gene Ontology (GO) resources for biological processes. Except two genes, these genes (Table 4) were found to be involved in prominent biological processes such as gene expression, DNA-dependent regulation of transcription, transcription from RNA polymerase-II, DNA replication, mitotic cell cycle and DNA repair in decreasing order of the number of genes involved in each biological process. These are all major processes widely implicated in tumor development and progression.

Correlation with protein expression levels

It is widely known that expression levels of mRNA and protein may not correlate at times due to the stability issues of mRNA and corresponding protein and further post-transcriptional modifications and regulatory processes. Further, if proteins are to be used as targets, it is important that their levels correlate with their mRNA levels too. In our attempt to understand further if mRNA and protein levels correlate in our list of genes, we used Human Protein Atlas (HPA) database. In HPA, protein expression levels are scored based on staining intensity score from immunohistochemistry experiments in cancer tissues. When mRNAs were found to be up-regulated or down-regulated in tumors and corresponding proteins had a score of high and low immunoreactivity, respectively, it was deciphered to be a positive correlation. A negative correlation was understood to be present when protein intensity score was high while mRNA is downregulated and vice versa. In some cases, we were unable to establish a clear correlation and categorized these genes as ‘uncorrelated or not able to correlate’. Thus, out of 70 genes, 34 genes were positively correlated, 9 genes negatively correlated, and remaining 16 genes were found to be uncorrelated (Table 5 and Fig. S2). Data was not available for 11 of these. Of these 34 positively correlated genes, using GeneCodis3, we found that processes/annotations over-represented for majority of genes were: gene expression, transcription from RNA polymerase-II promoter and DNA repair. Mitotic cell cycle and DNA replication were also among the biological processes involved (Table 6).

Involvement of these genes/proteins in cancer development and progression and their mechanistic actions have been detailed in a separate Supplemental File.

The genes MSL3L2 or MSL3P1 and YDJC appear to be novel and their functions in cancer are unknown, to the best of our knowledge.

To further understand their regulation in cancers, we turned our focus on all of these 34 key essential genes. Apart from regular gene regulators, lncRNAs are among the emerging regulators involved in biological processes predominantly through epigenetic regulation and found to be powerful regulators in cancer as well. Further, epigenetic regulation is at the heart of almost all the biological processes over-represented in our gene set. Therefore, in order to gain further regulatory and functional insights, we proceeded for lncRNA studies.

Identification of up- and down-regulated lncRNAs

Analyses of up- and down-regulated status of these 24 common significantly differentially expressed lncRNAs, showed that lncRNAs PVT1, UHRF1, AC005154.5, SNHG11, RP11-1149O23.3, and RP11-368I7.2 were up-regulated in tumors; MAGI2-AS3, MIR22HG down-regulated in tumors and the remaining were tumor-specific differentially expressed lncRNAs (Fig. 5).

Functional analyses of dysregulated lncRNAs

In light of the above studies and our own observations, we proceeded to study the function of lncRNA molecules in more detail. Function of a few lncRNAs such as Xist and NEAT1, both involved in cancer progression, have been experimentally validated (Patil et al., 2016; Qian et al., 2017). While Xist is involved in X-chromosome inactivation (Xiong et al., 2017; Zhang et al., 2017; Song et al., 2016a), NEAT1 plays a regulatory role in initiation and progression of cancers (Fang et al., 2017; Yang et al., 2017a; Li et al., 2016). PVT1 lncRNA functions as an oncogene by acting as a sponge for miRNAs (Li et al., 2018; Liu et al., 2016b; Li, Meng & Yang, 2017) as well as several other mechanisms as has been detailed above. Function of many lncRNAs remains to be discovered. One of the several approaches used to know the likely functions of lncRNAs is “guilt-by-association” with co-expressed protein-coding mRNAs, as these protein-coding genes are likely to be co-regulated by lncRNAs synergistically (Sun et al., 2015). In this manner, role of lncRNAs in specific biological processes and molecular function can be predicted.

We wanted to find out the function of lncRNAs in our narrowed down list. As no de-novo function prediction tool for lncRNAs is available as yet, we used lncRNAtor which predicts lncRNA function from co-expression of protein-coding mRNAs. We characterized each lncRNA for co-expressed mRNAs from lncRNAtor for each of the 9 cancer types. These co-expressed mRNAs were analyzed for their enrichment in biologicals processes and molecular function using GeneCodis3. Most prominent biological processes are common in which these lncRNAs may be involved through “guilt-by-association” are: signal transduction, gene expression, mitotic cell cycle, nerve growth factor receptor signaling pathway, apoptotic process and blood coagulation, in that order. DNA repair and DNA replication processes were also predicted. The predicted molecular functions are in protein binding, nucleotide-binding and ATP binding (Fig. 6 and Fig. S3 ). It is highly interesting to note from our observations that the major biological processes predicted for lncRNAs also correspond to the biological processes predicted for our significantly differentially expressed protein-coding gene list. This shows that these lncRNAs and protein-coding genes from our list may be interacting or regulating each other in one way or another.