LncRNA2Function: a comprehensive resource for functional investigation of human lncRNAs based on RNA-seq data

Background The GENCODE project has collected over 10,000 human long non-coding RNA (lncRNA) genes. However, the vast majority of them remain to be functionally characterized. Computational investigation of potential functions of human lncRNA genes is helpful to guide further experimental studies on lncRNAs. Results In this study, based on expression correlation between lncRNAs and protein-coding genes across 19 human normal tissues, we used the hypergeometric test to functionally annotate a single lncRNA or a set of lncRNAs with significantly enriched functional terms among the protein-coding genes that are significantly co-expressed with the lncRNA(s). The functional terms include all nodes in the Gene Ontology (GO) and 4,380 human biological pathways collected from 12 pathway databases. We successfully mapped 9,625 human lncRNA genes to GO terms and biological pathways, and then developed the first ontology-driven user-friendly web interface named lncRNA2Function, which enables researchers to browse the lncRNAs associated with a specific functional term, the functional terms associated with a specific lncRNA, or to assign functional terms to a set of human lncRNA genes, such as a cluster of co-expressed lncRNAs. The lncRNA2Function is freely available at http://mlg.hit.edu.cn/lncrna2function. Conclusions The LncRNA2Function is an important resource for further investigating the functions of a single human lncRNA, or functionally annotating a set of human lncRNAs of interest.


Background
Thousands of human long non-coding RNAs (lncRNAs) have been identified and emerging studies have revealed that lncRNAs play important roles in a wide range of biological processes [1,2] and diseases [3,4]. However, functions of most human lncRNAs are still elusive. Functions of a lncRNA may be determined by loss-and gain-offunction biological experiments [5,6]. However, this is not straightforward since it is difficult to knock down a lncRNA expressed as multiple isoforms. Alternatively, computational exploration of human lncRNA functions is helpful to guide further studies on lncRNAs.
Currently, computational investigation of lncRNA functions is still at its early development stage, since it is a considerable challenge due to the characteristics of lncRNAs, e.g., many lncRNA gene sequences are not conserved and do not contain conserved sequence motifs [7], which makes it difficult to infer potential functions of lncRNAs based on their sequences alone. In addition, few available molecular interaction data of new identified lncRNAs also hamper the lncRNA functional annotations [8,9].
Since genes with similar expression patterns across multiple conditions may share similar functions [10] or be involved in related biological pathways [11], identifying protein-coding genes that are co-expressed with lncRNAs may help to assign functions to the lncRNAs. By analyzing lncRNA-mRNA co-expression pattern, Guttman et al. identified several sets of mouse lncRNAs associated with protein-coding gene sets of distinct GO functional categories [12]. In addition, two recent studies separately constructed a mouse co-expressed lncRNA-mRNA network using mouse microarray data and assigned functions to 340 and 1,625 mouse lncRNAs [13,14].
Despite accumulating insights into the mouse lncRNA functions, more than 10,000 human lncRNAs remain to be functionally characterized. Firstly, given a single human lncRNA gene, it needs to be established whether it executes crucial biological functions. Secondly, given a set of human lncRNA genes such as differential lncRNAs between cancer and normal samples, it is an important downstream task to identify significantly enriched function terms. Thirdly, given an important functional term such as a Wnt signalling pathway, how to know which lncRNAs may be involved in the pathway.
Here, based on the expression correlation between lncRNAs and protein-coding genes inferred from RNAseq data of 19 human normal tissues, we functionally annotated 9,625 human lncRNAs with significantly enriched functional terms among the co-expressed protein-coding genes, and developed a user-friendly web interface for the lncRNA community to obtain the lncRNAs associated with a specific functional term, the functional terms associated with a specific lncRNA, or to assign functions to a set of human lncRNAs of interest.

Data sources
We downloaded: (1) genomic coordinates of all human lncRNA genes and protein-coding genes from the GEN-CODE V15 [15], (2) paired-end RNA-Seq data of 19 human normal tissues from the Human Body Map 2 project (ArrayExpress accession no. E-MTAB-513) and another study (GEO accession no. GSE30554), (3) GO assignments for the proteins of the human UniProtKB Complete Proteome from the website of the Gene Ontology Project [16], (4) 4,380 human biological pathways from the ConsensusPathDB database which integrated 12 pathway databases [17].

Workflow of LncRNA2Function
The schematic workflow of lncRNA2Function is shown in Figure 1. Firstly, RNA-Seq reads sequenced in 19 human normal tissues were firstly mapped to human genome (hg19) using tophat with the default parameters [18], and expression values of all human lncRNA and protein-coding genes in the 19 tissues were computed using cufflinks with the default parameter [19]. Secondly, the Pearson Correlation Coefficients (PCC) of all lncRNA-mRNA gene pairs were computed, and a set of significantly co-expressed protein-coding genes was thus obtained for each human lncRNA (significant: the absolute value of the Pearson correlation coefficient >0.9 and adjusted P-value < 0.05). Thirdly, each lncRNA was functionally annotated with significantly enriched GO terms and biological pathways among the set of co-expressed protein-coding genes. Finally, a web interface was developed to facilitate researchers to browse or search the functions associated with a given lncRNA or lncRNAs associated with a specific function, or to functionally annotate a set of lncRNA genes of interest.

GO and pathway enrichment analysis of human lncRNAs
Given a single human lncRNA gene, we obtained a set of protein-coding genes that were significantly co-expressed with the lncRNA. The lncRNA was then functionally annotated with significantly enriched GO and pathway terms among the set of co-expressed protein-coding genes. The enrichment analysis was separately executed for each term (denoted as T), and a P-value of each term was calculated by the hypergeometric test: Herein, N is the number of all protein-coding genes in human genome, M is the number of protein-coding genes that were annotated in the functional term T, n is the number of protein-coding genes that were significantly co-expressed with the lncRNA, and m is the number of protein-coding genes that were both significantly co-expressed with the lncRNA and annotated in the functional term T.
For each GO term, protein-coding genes directly belong to it as well as those belong to any of its offspring terms are all considered as its annotated genes. Since the statistical analysis is not appropriate to problems with small sample size, those GO and pathway terms with less than 5 annotated protein-coding genes and those lncRNAs with less than 5 co-expressed protein-coding genes were excluded form the enrichment analysis.
Given a set of human lncRNA genes of interest, LncRNA2Function first identify a set of protein-coding genes, each of which are significantly co-expressed with one or more of the given lncRNAs across 19 human normal tissues. Then, the set of lncRNAs are functionally annotated with the enriched GO and pathway terms among the set of co-expressed protein-coding genes. If researchers input a large number of lncRNAs, the LncRNA2Function may obtain thousands of co-expressed protein-coding genes, some of which are co-expressed with only one of the lncRNAs. To improve the accuracy of functional assignments to the set of lncRNAs, users can select the protein-coding genes that are co-expressed with at least K lncRNAs (the K can be assigned based on the size of the set of lncRNAs. The default K is 1).
There are two commonly used methods for controlling false discovery rate (FDR), the Benjamini-Yekutieli (BY) method [20] and the Benjamini-Hochberg (BH) method [21]. The former is suitable for positively related multiple hypothesis tests whereas the later is suitable for independent multiple hypothesis tests [22]. Since the hierarchical GO terms are often dependent, we chose the BY method to correct the P-values from the GO enrichment analysis, and the BH method to correct the P-values from the pathway enrichment analysis. The significant cut-off of corrected P-value was set as 0.05.

Results and discussion
Functional annotations of a single human lncRNA We obtained 5,232,299 significantly co-expressed pairs between 9,625 human lncRNA genes and 10,919 proteincoding genes. Each of the 9,625 lncRNAs was functionally annotated with significantly enriched GO terms and biological pathways among its co-expressed protein-coding genes. Consequently, we obtained 614,174 associations between 5,735 lncRNA genes and 3,890 GO terms, and 240,050 associations between 6,062 lncRNAs and 3,034 biological pathways. To understand the major functions of lncRNAs, we ranked GO biological processes and biological pathways according to the number of lncRNAs associated with each of them. Among the top ranked 200 GO biological processes and pathways, we found that lncRNAs play roles in many important biological processes, including defense response to bacterium, DNA packaging, meiosis, developmental process, metabolic process, cell cycle process, cell adhesion, cell differentiation, Jak-STAT signaling pathway and PI3K-Akt signaling pathway. A part of the enriched functions of lncRNAs have been validated by published studies [23][24][25][26].

Case studies
Due to the lack of a large gold standard dataset of known human lncRNA functions, five well-studied lncRNAs were used as the examples to show the usefulness of LncRNA2Function.

Case study 1: HOTAIR
The HOTAIR is a well-studied lncRNA. Rinn et al. found that the HOTAIR interacts with the Polycomb repressive complex 2 (PRC2) to modify chromatin and repress transcription of the HOX genes, which regulate development [27]. Niinuma et al. revealed that overexpression of HOTAIR was strongly associated with highrisk grade and metastasis among gastrointestinal stromal tumors (GIST) specimens, and knockdown of HOTAIR suppressed GIST cell invasiveness [28]. In addition, Gupta et al. demonstrated that the lncRNA HOTAIR is increased in expression in primary breast tumors and metastases, and enforced expression of HOTAIR in epithelial cancer cells leaded to altered histone H3 lysine 27 methylation, gene expression, and increased cancer invasiveness and metastasis in a manner dependent on PRC2. Conversely, loss of HOTAIR can inhibit breast cancer invasiveness [26].
To examine whether our LncRNA2Function can functionally annotate the lncRNA HOTAIR with development and metastasis-related functional terms, we applied the LncRNA2Function to the HOTAIR, and found that it was annotated with 99 GO biological processes and 33 pathways (The significant Corrected P-value cutoff is 0.05). Of the 99 GO biological processes, 77.8% (77/99) are involved in the development and morphogenesis as expected (The top 20 GO development-related biological processes are shown in Table 1), and 9.1% (9/99) are involved in the cell invasion and metastasis, such as cell migration (GO:0016477), cell adhesion (GO:0007155), biological adhesion (GO:0022610) and cell motility (GO:0048870). In addition, Of the 33 biological pathways, 72.7% (24/33) are involved in the cell invasion and metastasis (Table 2), such as focal adhesion, beta1 integrin cell surface interactions, NCAM1 interactions, Syndecan-1mediated signaling events, PI3K-Akt signaling pathway and cell surface interactions at the vascular wall. Taken together, these results demonstrated that our LncRNA2-Function can successfully recall the known functions of a well-studied lncRNA HOTAIR and suggested that it is applicable to infer potential functions of new identified lncRNAs.

Case study 2: HCP5
The lncRNA HCP5 was found to be associated with AIDS [29][30][31]. Rodriguez-Novoa et al. analyzed a total of 245 HIV patients and found a good correlation between HLA-B*5701 and HCP5 (negative and positive predictive values of 100% and 93%, respectively). Colombo et al. analyzed that 1,103 singles infected with human immunodeficiency virus (HIV) and concluded that HCP5 genotyping could serve as a simple screening tool for ABC-HSR, particularly in settings where sequence-based HLA typing is not available.
To assess whether the HCP5 can be correctly predicted to have immune-related functions, we applied our LncRNA2Function to it and found that HCP5 was annotated with 549 GO biological processes terms and 270 biological pathways. As expected, most of them are indeed immune system and response functional terms, which are strongly associated with the development of AIDS. The top 20 GO biological terms assigned to the HCP5 are shown in Table 3 while the top 20 biological pathways assigned to the HCP5 are shown in Table 4.

Case study 3: HULC
The lncRNA HULC is highly upregulated in liver cancer and plays an important role in tumorigenesis [32]. Depletion of HULC resulted in a significant deregulation of several genes involved in liver cancer [33], and colorectal carcinomas that metastasize to the livers but not to lymph nodes experience an up-regulation of HULC in all the samples tested (n = 8), with a strong-to-moderate expression in six out of eight [34].
To examine whether the HULC was predicted to have liver-related functions, we analyzed it using our lncRNA2-Function. Expectedly, LncRNA2Function also works well to functionally annotate the HULC. The results showed that it was annotated with 373 GO biological processes and 383 biological pathways (the significant P-value cutoff is 0.05). Of the 373 GO biological processes and 383 pathways, over 80% are involved in the known liver-related biological functions, such as metabolic function, bile secretion, lipid transport and homeostasis, cholesterol homeostasis, regulation of blood coagulation, protein-lipid complex subunit organization, detoxification, Immune defense and complement activation. The Figure 2 shows the top 25 GO functional terms assigned to the HULC, and the Table 5 shows the top 20 pathways enriched in protein-coding genes that are co-expressed with the liverrelated lncRNA HULC.

Case study 4: H19
H19 is an important lncRNA that play roles in the infertility [35] and multiple cancers such as breast cancer [36,37], cervical cancer [38], liver cancer [39,40] and bladder cancer [41]. For example, Korucuoglu et al. revealed that H19 expression was lower in the infertility group as compared to the control group (4-fold change,   We applied the LncRNA2Function to the lncRNA H19 and found that it was annotated with 6 GO biological processes and 31 biological pathways. The GO terms includes female pregnancy (GO: 0007565), estrogen biosynthetic process (GO:0006703), growth hormone receptor signaling pathway (GO:0060396), cellular response to growth hormone stimulus (GO:0071378) and JAK-STAT cascade involved in growth hormone signaling pathway (GO:0060397), which suggest that H19 may play roles in infertility or breast cancer by participating in these biological processes. In addition, the cancerrelated lncRNA H19 was correctly annotated with many important caner pathways, such as PI3K-Akt signaling pathway, GPCR signaling-G alpha s Epac and ERK pathway, Nuclear signaling by ERBB4 pathway, Akt signaling pathway and JAK-STAT-Core cancer pathway. These results suggest that our LncRNA2Function correctly recall the known functions of H19.

Case study 5: PCA3
The lncRNA prostate cancer antigen 3 (PCA3) is a highly specific biomarker upregulated and plays crucial roles in prostate cancer (PCa) [42][43][44][45]. Clarke et al. found that up-regulation of two new PCA3 isoforms in PCa tissues improves discrimination between PCa and benign prostatic hyperplasia (BPH). In 2012, the US Food and Drug Administration approved the use of the lncRNA PCA3 for the detection of prostate cancer.
To test whether our LncRNA2Function can annotate the PCA3 with prostate-related functions, we applied the LncRNA2Function to the PCA3. LncRNA2Function first identified 77 protein-coding genes that are coexpressed with the PCA3 and then annotated it with only one pathway named 'Regulation of Androgen receptor activity' (corrected P-value: 0.020385). This pathway has 62 genes, which includes 4 protein-coding genes that are co-expressed with the PCA3. These four genes are HOXB13, KLK3, KLK2 and SPDEF that have been validated to be useful in the diagnosis and monitoring of prostatic carcinoma and be suitable target for developing specific cancer therapies. Consequently, lncRNA2Function can correctly predict the functions of PCA3 by its co-expressed protein-coding genes.

Functional annotation for a set of human lncRNAs
High-throughput genomic technologies like lncRNA microarray and RNA-Seq usually generate hundreds of candidate lncRNA genes of interest, such as a cluster of co-expressed lncRNA genes across multiple conditions or a set of differentially expressed lncRNAs between cancer and normal samples. To manually map each lncRNA to functional terms is by far a simple task. Therefore, how to identify significantly enriched functions among the set of lncRNAs is an important downstream task for interpreting high-throughput experimental data.
As a proof-of-concept, a set of liver-specific lncRNAs and a set of heart-specific lncRNAs inferred from RNA-Seq data of 19 human normal tissues were used as examples to show the functionality of our lncRNA2Function system in annotating a set of lncRNAs of interest, respectively. As expected, lncRNA2Function correctly assigned the functional terms to the two distinct sets of lncRNAs. Users can test these two sets or their own lncRNA sets at our 'LncRNA set analyzer' web interface http://mlg.hit. edu.cn/lncrna2function/lncrna_enrich.jsp.

Web interface of LncRNA2Function
To facilitate researchers to access the functional annotations of lncRNA genes, we developed a web interface named 'LncRNA annotation browser', which is a userfriendly interface to browse or search lncRNAs associated with a specific functional term, or functional terms associated with a given lncRNA. To enable researchers to analyze a set of lncRNA genes of their interest, we implemented a web interface titled 'LncRNA set analyzer', which can help investigators to annotate a set of lncRNAs with Gene Ontology and 4,380 biological pathways curated from 12 pathway databases. In addition, we developed a web interface titled 'LncRNA expression viewer' to facilitate investigators to graphically view the expression dynamics of genes across multiple human normal tissues. Users can not only view expression value of a single lncRNA or protein-coding gene across 19 human normal tissues, but also simultaneously view the expression index of both lncRNA and protein-coding genes to learn about whether they are co-expressed across the 19 tissues. Furthermore, we provide a submission page that allows other researchers to submit known functional annotations of lncRNAs that are not documented in our LncRNA2Function system ( Figure 3). They do not have to be an author on the original study to submit a record. Once approved by the submission review committee, the submitted records will be made available to the public in the coming release. LncRNA2-Function is freely accessible at http://mlg.hit.edu.cn/ lncrna2function.

Conclusions
Thousands of human lncRNAs have been identified in recent several years, while the vast majority of the lncRNAs remain to be functionally characterized. In this study, we functionally annotate 9,625 human lncRNAs with the enriched functions among the protein-coding genes that are co-expressed with each lncRNA. Furthermore, we developed a web interface, which facilitates researchers to search the functions of a specific lncRNA or the lncRNAs associated with a given functional term, or annotate functionally a set of human lncRNAs of interest. The lncRNA2Function will become an important tool for investigating functions of human lncRNAs.