deepBase v2.0: identification, expression, evolution and function of small RNAs, LncRNAs and circular RNAs from deep-sequencing data

1Eukaryotic genomes encode thousands of small and large non-coding RNAs (ncRNAs), such as microRNAs (miRNAs), long non-coding RNAs (lncRNAs) and circular RNAs (circRNAs). In contrast to the well-known miRNAs, which have been found to be key players in various biological processes, the functions of the majority of lncRNAs and circRNAs have not been fully investigated.

2To better understand the functions of these ncRNAs, a systematic catalog of ncRNA transcripts and their expressions across tissues and evolutional conservation across species is necessary. To date, many databases have been developed to curate computationally predicted and experimentally verified ncRNAs, such as LncRNAdb, GENCODE, lncRNAtor, ChIPBase, NONCODE, miRBase, starBase and circBase. However, these databases either focus on limited species or contain insufficient expression or evolution data.

3Recent advances in high-throughput next-sequencing technology have produced large numbers of short and long RNA sequences, and enable the detection and profiling of known and novel ncRNAs at unprecedented sensitivity and depth. Many studies have identified thousands of sRNAs, lncRNAs and circRNAs from deep-sequencing datasets in various species, including human, mouse, fruitfly, nematode and zebrafish. With the increasing amount of deep-sequencing data available, there is a great need to integrate these large-scale data sets to explore the expression, evolution and function of diverse ncRNAs.

4To meet above-mentioned needs, we have updated deepBase to version 2.0 (deepBase v2.0) (Figure 1). The deepBase v2.0 facilitates the integrative, interactive and versatile display of, as well as the comprehensive annotation and discovery of sRNAs, lncRNAs and circRNAs (Figure 1). In deepBase v2.0, we constructed the most comprehensive expression profiles and evolutional patterns of diverse ncRNAs.

5A total of 558 small RNA deep-sequencing datasets and 478 RNA-seq datasets were compiled in deepBase v2.0 (Table 1). The raw data were downloaded from NCBI GEO and SRA databases. All sequencing datasets and their detailed information are listed in Supplementary Table S1. The genome sequences and transcript sequences were downloaded from NCBI Reference Sequences (RefSeq), UCSC Bioinformatics websites and the following specialized databases: human (UCSC hg19), mouse (UCSC mm10), chicken (Gallus gallus, v4) and Ciona intestinalis (JGI v2.0) genome sequences were downloaded from the UCSC Bioinformatics website. Chimp (Pan troglodytes, panTro4), Gorilla (Gorilla gorilla gorilla, gorGor3), Rhesus (Macaca mulatta, rheMac3), Rat (Rattus norvegicus, rn6), Cow (Bos taurus, bosTau7), Opossum (Monodelphis domestica, monDom5), Platypus (Ornithorhynchus anatinus, ornAna1), X. tropicalis (Xenopus tropicalis, xenTro3), Zebrafish (Danio rerio, danRer7) and C. elegans (Caenorhabditis elegans, ce10) genome sequences were downloaded from the UCSC Bioinformatics website. The Caenorhabditis remanei and Caenorhabditis briggsae genome sequences were downloaded from WormBase, the Drosophila melanogaster genome sequences were download from Flybase. The Bombyx mori genome sequences were download from silkDB. All known miRNAs were downloaded from miRBase (release 21). Other ncRNAs were downloaded from Ensembl Gene Release 76.

6The raw data downloaded from NCBI SRA or GEO databases were classified into different clades, species, tissues and cell-lines according to the description on the website or related literature. For small RNA annotation, after remove the 3′ adapter sequences, the clean reads were then mapped back to their corresponding reference genome using bowtie program. Bedtools intersect tool was used to inspect how many reads could matched to the region of annotated small RNA genes. The reads which were perfectly located in the annotated gene regions were used for expression evaluation.

7To accurately annotate lncRNA, we proposed a pipeline to process transcriptome annotation profiles generated by RNA-seq (usually assembled by Tophat and Cufflinks) using several filters described as follow: Transcript length filter: Multi-exonic transcripts whose length exceeded 200nt were kept.Known non-lincRNA annotation filter: Transcripts were excluded from further consideration if they overlap, or reside within 50nt upstream or downstream of, any exon of any transcript from the following annotation sets: coding genes annotated in RefSeq, Ensembl or UCSC knownGene; pseudogenes annotated in Ensembl or Yale Pseudogenes. Those which have any overlap spanning at least 80% of miRNA precursors annotated in miRBase or miRNA/tRNA/rRNA/snRNA/snoRNA annotated in Ensembl were also eliminated. The above rules apply in a strand-specific manner, thus allowing for calling antisense lncRNA.Coding potential filter: Three open-reading-frames-related metrics were used to estimate whether the remaining transcripts exhibit significant protein-coding potential. The getorf program from the EMBOSS suite was used to find ORFs in the candidates and therefore the length and coverage of ORFs and were calculated. ORFScan scores were calculated by the txCdsPredict program from UCSC. Transcripts with ORF <100 aa, ORF coverage less than 30% and txCdsPredict score <800 were classified as lncRNAs. Furthermore, we classified an additional subset of transcripts (ORF coverage between 30% and 90% and txCdsPredict score <800) that may include short ORFs and may serve as either lncRNAs or small peptides as TUCP (transcripts of uncertain coding potential), which was defined and used in several papers. The cutoff for these parameters was determined based on our survey of the known functional lncRNAs collected in the lncRNAdb database.Known protein domains filter: Each remaining transcript was mapped to unitRef90 database with BlastX and transcripts with E-value < 1E–30 were removed. The resulting transcripts were translated in all three forward frames and fed to HMMER v3.0 to examine the occurrence of any PfamA and PfamB known protein family domains collected in the Pfam database. Transcripts with a Pfam hit showing both the full sequence E-value < 1E–5 and the single best domain E-value < 1E–5 were removed.

8Transcript length filter: Multi-exonic transcripts whose length exceeded 200nt were kept.

9Known non-lincRNA annotation filter: Transcripts were excluded from further consideration if they overlap, or reside within 50nt upstream or downstream of, any exon of any transcript from the following annotation sets: coding genes annotated in RefSeq, Ensembl or UCSC knownGene; pseudogenes annotated in Ensembl or Yale Pseudogenes. Those which have any overlap spanning at least 80% of miRNA precursors annotated in miRBase or miRNA/tRNA/rRNA/snRNA/snoRNA annotated in Ensembl were also eliminated. The above rules apply in a strand-specific manner, thus allowing for calling antisense lncRNA.

10Coding potential filter: Three open-reading-frames-related metrics were used to estimate whether the remaining transcripts exhibit significant protein-coding potential. The getorf program from the EMBOSS suite was used to find ORFs in the candidates and therefore the length and coverage of ORFs and were calculated. ORFScan scores were calculated by the txCdsPredict program from UCSC. Transcripts with ORF <100 aa, ORF coverage less than 30% and txCdsPredict score <800 were classified as lncRNAs. Furthermore, we classified an additional subset of transcripts (ORF coverage between 30% and 90% and txCdsPredict score <800) that may include short ORFs and may serve as either lncRNAs or small peptides as TUCP (transcripts of uncertain coding potential), which was defined and used in several papers. The cutoff for these parameters was determined based on our survey of the known functional lncRNAs collected in the lncRNAdb database.

11Known protein domains filter: Each remaining transcript was mapped to unitRef90 database with BlastX and transcripts with E-value < 1E–30 were removed. The resulting transcripts were translated in all three forward frames and fed to HMMER v3.0 to examine the occurrence of any PfamA and PfamB known protein family domains collected in the Pfam database. Transcripts with a Pfam hit showing both the full sequence E-value < 1E–5 and the single best domain E-value < 1E–5 were removed.

12All transcripts that passed all the filters above were classified based on their loci and treated as a set of lncRNAs in the following study. All the lncRNAs in deepBase v2.0 are named according to the following rules. The names of all the lncRNAs in deepBase v2.0 consist of three parts, ‘A-B-C’. Amongst, ‘A’ is composed by three letters, which represents the shorted name of species, ‘B’ is the gene ID and ‘C’ is the transcript ID. For example, hsa-lncRNA10932–8 means the 8th transcripts of gene lncRNA10932 in human.

13Human, mouse and C. elegans circRNA genes were downloaded from circBase v1.0 or obtained from the Supplementary Data of the original articles. In addition, we used the method ‘find_circ’ to identify novel circRNAs from raw RNA-seq data of ENCODE and modENCODE. In brief, raw reads were first aligned to the genome. Then, unmapped sequencing reads were used for back-splice sites identification.

14DESeq was used to evaluate the expression of small RNAs. FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values were calculated by Cufflinks and used to evaluate the expression of lncRNAs. Cufflinks includes a program, Cuffnorm, which can be used to generate tables of expression values that are properly normalized for library size.

15The sequences of all identified long non-coding RNAs were extracted from genomes of corresponding species using Bedtools v2.17.0. As described by Ulitsky et al., we performed similarity searches with each lncRNA sequence in one species against that in the rest species using NCBI BLASTN with parameters ‘-task blastn -word_size 6 -evalue 1e–5 -strand plus’ and the resulting hits with the lowest e-value was kept as the most conserved lncRNA.

16Expression correlation (Pearson correlation coefficient) between lncRNAs and protein-coding genes was estimated using coexpression program (the program is available from the authors upon request), which was written in C language and ALGLIB library, and p-value was adjusted with the False Discovery Rate (FDR) correction. GO ontology data for the Ensembl Genes was downloaded from the Ensembl website. There are thirteen functional categories for human genome. Enrichment analysis of these pathways in the dataset was determined using a hypergeometric test with Bonferroni and FDR correction.

17deepBase provides genome-wide identification of small RNAs in multiple types, including microRNAs, snoRNAs, rRNA-derived small RNAs. In total, there are 18 508 small RNAs are identified in 51 studies of 10 species (Table 1). Users who want to see the results of these small RNAs can use the Browse page. For example, deepBase v2.0 integrates a study which sequenced small RNAs in 60 samples of tumour or normal human cervical tissues. After our systematic re-annotation of all types of small RNAs in these samples, the ‘Browse’ page shows 3,176 small RNAs in the result table. By ordering the ‘Sample’ column of the table, users can find 14 small RNAs that could be detected in all of the 60 samples, including 12 miRNAs, one snoRNA (SNORD85) and one snRNA (RNU2–59P). By sorting the ‘Expression’ column of the table, the highest expressed miRNA hsa-mir-21 ranks on the top of the table. The ‘View’ button on the last column can link to the detail information of this miRNA, which displays the distribution of all the sequencing reads on the miRNA precursor in each sample. Most of the reads correspond to the mature arm of the miRNA precursor, and the highest one has 10 005 sequencing copies.

18In the Browse page, we also provide integration of lncRNAs identified by RNA-seq, organized by different species and corresponding studies (Table 1). About half of the lncRNAs were not annotated in GENCODE (Supplementary Figure S1). For instance, we applied our lncRNA identification pipeline (detailed in Methods) to accurately annotate lncRNAs in mouse from RNA-seq data generated by a previous study. 14 462 mouse lncRNAs were identified, among which 11 220 were not annotated in GENCODE. We compared the result of lncRNAs in deepBase v2.0 with other three studies. Supplementary Table S2 shows the number of overlapped lncRNAs in deepBase v2.0 and novel lncRNAs that is not reported in the other studies. According to our results, although varied in amounts, the lncRNAs in different species shared common features (Supplementary Figures S2–S5). For example, the transcript length of lncRNAs is generally shorter than that of protein-coding genes but longer than other ncRNAs (Supplementary Figure S3). Users can view basic information of each lncRNA in the query result, and pressing the ‘Predict’ button of each entry will call the web-based lncFuntion program (introduced in the later session) to predict possible functions of the lncRNA.

19One of the most important features of miRNAs is their special expressions in different stages, tissues or cells. deepBase v2.0 provides the expression profiles of small RNAs in each sample. Users can find their interested miRNA name in the small RNA ‘Expression’ page. The heatmap will display the normalized expression number of miRNA in each sample. Click the name of the miRNA will link to the detail expression result page. Users can find that hsa-mir-122 is especially highly expressed in normal liver tissues compared to other samples.

20Spatial and temporal specific expression is one of the well-characterized properties of lncRNAs. In the Expression page, expression profiles of lncRNAs are visualized through heatmaps, providing an easy way for users to spot tissue-specific lncRNAs. For instance, H19, a maternally imprinted lncRNA, was reported to be highly expressed in placenta and down-regulated in all tissues except skeletal muscle immediately after birth. This phenomenon can be observed through our expression heatmap of H19 generated from data of the Human Body 2.0 project.

21We conducted evolution conservation analysis of lncRNAs in 14 species. Many studies show that most lncRNAs are lineage-specific and in fast phylogenetic evolution. Out of 18,964 lncRNAs in human we identified, 16 926 have at least one ortholog in primates. However, only 754 lncRNAs present to be conserved in Platypus. 13 800 human lncRNAs have not previously been annotated, among which hsa-lncRNA12238 from antisense is highly but only conserved in primates, indicative of weak evolution constraints. hsa-lncRNA12238 presents intensively high expression in human testis. To verify that, we then predict the function of hsa-lncRNA12238 in Function entry, and GO items show strong correlation with such reproductive process as spermatogenesis, sperm-egg recognition and so forth.

22We provide a web-based tool named IncSeeker to implement our filtering pipeline of finding lncRNAs transcriptome assembled from RNA-seq data. Users can upload their transcriptome data in BED12 or GTF format, and our web-server will process the file by a series of tunable criterion and filter steps, including transcript length, exon number, coding ability and potential coding domain. A strict set of lncRNAs will then be displayed on the result page. To evaluate the performance of lncSeeker, we used lncSeeker to predict lncRNAs in GENCODE gene sets. There are 11 324 high-confidence lincRNA transcripts and 95 309 protein-coding transcripts in GENCODE Release 19 (GRCh37.p13). The result showed that the accuracy of lncSeeker was 96.82%, the sensitivity and specificity were 70.08% and 100% (Supplementary Table S3), respectively. We further compared lncSeeker with other four tools related with lncRNAs identification, getorf, txCdsPredict, CPAT and CNCI. The result showed lncSeeker had the highest accuracy and specificity.

23A web-based tool, lncFunction, was also developed to predict lncRNA functions from coding and non-coding co-expression networks. There are six options on the query page, and the option ‘Input the Target LncRNA’ is required. The running parameters, expression values from selected study and every outcome in the three categories of function prediction are available for users to download.

24Although a few dozen of lncRNAs and circRNAs have been characterized to some extent and reported to function in important cellular processes, such as differentiation, development, proliferation, self-renewal, pluripotency, carcinogenesis and progression, the functions of most annotated ncRNAs are unknown. In this study, by analyzing a large set of sRNAs, lncRNAs and circRNAs identified from deep-sequencing datasets, we have annotated ncRNAs with a broad range of structural, expression, and evolutionary features.

25Currently, there is an increasing amount of databases developed for ncRNAs, including NONCODE, lncRNAdisease, ncFANs, lncRNA2Function, FAME, miRo, miRGator and lncRNAwiki. NONCODE is a comprehensive database which contains all kinds of ncRNAs, but no evolution information for these ncRNAs. lncRNAdisease is a specialized database which collects disease-associated lncRNAs and their interaction entries. On comparison with ncFANs, which is used to predict lncRNAs from ‘Affymetrix microarray datasets’, our tool is designed to identify lncRNAs from ‘RNA-Seq’ based datasets. Compared with other pipeline/tools, our pipeline integrates more tools as well as built-in scripts and annotates and identifies lncRNAs more reliably. Compared with our previous deepBase v1.0 which was designed to focus on small RNAs, deepBase v2.0 was designed to identify and annotate lncRNAs and circRNAs, and to predict the functions of lncRNAs. deepBase v2.0 also provides more comprehensive expression and evolution profiles of lncRNAs, circRNAs and small RNAs. The distinctive features of deepBase v2.0 include the following: (i) deepBase v2.0 has provided the most comprehensive expression analysis of sRNAs and lncRNAs from 1036 RNA-Seq datasets from 19 species. The constructed gene expression profiles of both ncRNAs and protein-coding genes from our analyses are a valuable resource for understanding the similarities and differences of transcriptional regulation of protein-coding genes and ncRNAs across different tissue/cell-line types. (ii) To the best of our knowledge, this is the first attempt to construct evolutional patterns of lncRNAs and circRNAs across several evolutional clades. Conservation patterns provided by our database may help biologists to select important ncRNAs for further functional validation. (iii) Two web-based tools, lncSeeker and lncFunction, can be used to identify high-confidence lncRNAs, and to predict lncRNA functions. We expect that access to these tools will enable more researchers to search for functions of novel lncRNAs in the ever-increasing amounts of deep-sequencing data.

26Next-generation sequencing technologies play vital roles in improving our understanding of functional genomics. As sRNA-Seq and RNA-Seq technology is applied to a broader set of species, cell lines, tissues and conditions, we will continually maintain and update the database. The integration of transcriptome datasets from the deepBase database, and the cancer genomics data and clinical information from the Gene Expression Omnibus (GEO), The Cancer Genome Atlas (TCGA) and International Cancer Genome Consortium (ICGC), will improve our understanding of expression and function of ncRNAs in diseases.

27deepBase v2.0 is freely available at http://deepbase.sysu.edu.cn/ or http://biocenter.sysu.edu.cn/deepBase/. The deepBase data files can be downloaded and used in accordance with the GNU Public License and the license of primary data sources.