Homology search consists in searching for new members of known ncRNAs families. The important, and increasing, number of publicly available ncRNA annotations makes homology search the first step in any ncRNA annotation process.

Search tools based on sequence alignment, such as BLAST [24], allow for very fast searches of sequences (query sequences) in genomes (target sequences). The sequence similarity between query and target sequences plays a major role in the success rate of BLAST searches. Significant candidates (i.e. those with E-values < 0.1) will present sequence similarities above 84% with a minimum length of 16 nts. Thus, the use of BLAST for homology search becomes more effective for ncRNA families with large conserved sequences, such as rRNA or snRNAs, or when searching within closely related species. When the target and query sequences come from distant species, or present low sequence conservation, with potentially large insertions, pure sequence alignment methods loose efficiency. A way to get around this limitation is to explore the statistical signals imprinted on the sequence by the RNA structural constraints.

ncRNAs present complex three-dimensional structures of packed double-stranded helical regions connected by tertiary interactions that are generally mediated by single-stranded regions [25]. The set of double-stranded and single-stranded regions is called secondary structure. RNA helices consist of stacks of A–U; G–C and G–U base pairs. Those base pairs are structurally equivalent, and the substitution of one base pair by another one has, frequently, minimal or no impact on the molecular structure. The accumulation of base pair substitutions in a RNA sequence can render two homologous sequences very dissimilar. However, when observed from the point of view of a multiple sequence structural alignment, it generates a pattern of covariation between the base-paired positions that can be detected in the respective columns of the sequence alignment [26,27]. This dependency between paired positions in alignments is used to build covariance models, used by ncRNA search tools such as INFERNAL [28], to search for ncRNAs in large sequences. Additionally, if sufficiently diverse sequences are included in the alignments, the known positions of insertions can also be included into the models. Covariance models, for most of the known ncRNAs families, are curated and maintained in the RFam database [29] and could be readily applied in our search.

The results produced by any search tool must be automatically filtered in order to exclude candidates less likely to be real ncRNAs. Candidates conflicting with known annotations or with low score should be discarded from the candidate list. Additionally, all retained candidates should be structurally aligned with known homologues in order to facilitate human validation, required as the last step of the annotation process (Fig. 1A).

The Génolevures database [22] contains 9 budding yeast genomes in which only tRNAs and rRNAs were fully annotated and were not considered in this work. Some other ncRNA families were partially annotated (see Table 1, column “Original Annotation”). The Ashbya Genome Database [23] contains the fully annotated genome of the Ashbya Gossypii (also known as Eremothecium gossypii) yeast.

We performed a BLAST search on those 10 genomes using as queries the S. cerevisiae ncRNAs and all the originally annotated ncRNAs from the original databases. This set of query sequences corresponds to the closest species for which we had reliable and ready to use ncRNA annotations. We believe that using a larger and more distant set of ncRNA sequences would increase the low score candidates with little improvement on the amount of genes found (e.g. of the 86 S. cerevisiae query sequences 37% produced true positive candidates when BLASTed against Candida glabrata, but only 6% when BLASTed against the more distant Debaryomyces hansenii and Yarrowia lypolitica–see Table SI 1). The BLAST search produced 1540 candidates with E-values < 0.1. A first INFERNAL search, using the 83 covariance models corresponding to the ncRNA families present in S. cerevisiae genome, produced 1250 candidates with E-values < 0.5. The candidates were included in the structural alignment of the corresponding family and manually validated according to the criteria described in the section “Candidate Validation”. After validation we retained 554 BLAST candidates and 602 INFERNAL candidates as ncRNAs.

We were interested on how E-values behave as classifiers in this specific data set. Fig. 2 shows the ROC curve [30], after the manual classification of True and False positives (TP and FP), for both tools. INFERNAL E-values were slightly better discriminators than the E-values of BLAST. While both tools will achieve 99% of True Positive Rates (TPR) at similar E-values (0.07 and 0.08 for INFERNAL and BLAST, respectively), the False Discovery Rates (i.e., the proportion of FPs in selected candidates) is significantly lower for INFERNAL (0.35) than for BLAST (0.61).

A second INFERNAL search was performed including all the remaining RFam families with exception of viral and miRNA families (719 Rfam families). According to the sensitivity analysis performed previously (Fig. 2), we reduced the expected value cutoff to < 0.1 obtaining 1360 candidates (applying the less restrictive cutoff, E-values < 0.5, we would have obtained 5578 candidates). From this search, only 41 candidates were selected. This surprisingly low number can be explained by the fact that most of the 719 Rfam families do not have homologues in yeasts, thus the resulting candidates are mostly FPs with high E-values (see Supplementary material SI 1). Curiously, all except two of the selected candidates correspond to homologous ncRNAs not found in the first search. A representative example is the box C/D snoRNA snR47 that was identified in Y. lypolitica, D. hansenii, K. thermotolerans and Z. rouxii with the covariance model of the snoRNA SNORD36 that is the mammalian ortholog of snR47. The remaining two candidates were, until now, unknown in the Saccharomycotina clade: a Box C/D snoRNA found in Y. lypolitica (Fig. 3) and a TPP riboswitch found in the 5’ UTR of protein coding genes in three different species (Fig. 4). The snoRNA candidate was found by INFERNAL using a covariance model of an archaeal snoRNA from the Pyroccocus family. The conservation of the characteristics Box C/D, and the complementarity of the putative guide sequences with the ribosome prevent a rapid exclusion of this candidate. On the other hand, the fact that the homologues of the putative ribosomal target sites are not modified in S. cerevisiae raises doubts about the nature of this candidate that only experimental validation could confirm. The TPP riboswitches are the only riboswitches found in eukaryotes [31,32] and the structural alignment with the crystal structures of the bacterial (Escherichia coli) [33] and plant (A. thaliana) [34] TPP riboswitches reveals a striking similarity of key structural nucleotides. The distribution of this regulatory domain across the Saccharomycotina phylogeny (Fig. 4B), totally absent in the S. cerevisiae branch, while present in D. hansenii, A. adeninovorans. and Y. lypolitica, raises interesting questions about the evolution of this type of mechanisms, and show the utility of extending homology search beyond the close group of species.

Table 1 presents the results of the homology pipeline distributed by species and ncRNAs families. The total retained candidates correspond to 79% of all ncRNAs that were expected assuming the S. cerevisiae database as the reference for the ncRNA families present on yeast. At least 60% of the expected ncRNAs were found in all species. Unsurprisingly the homology search was much more effective in the species from the upper branch (from S. cerevisiae to A. gossypii), as 90% of the sequences used as queries came from that branch. Comparing the proportion of ncRNAs found by INFERNAL and BLAST for each species (Fig. 5B), we can observe that BLAST was as sensitive as INFERNAL as long as the searched genomes were close enough. In more distant species, while both tools loose performance, INFERNAL is much more sensitive.

Finally, it was not possible to find the RNA component of the Telomerase complex with any of the used tools. This ncRNA presents a challenge to automatic search programs due to minimal sequence and secondary structure conservation, extensive insertions/deletions and variable size between species [9,35,36]. The telomerase ncRNAs contains some structural features that are common to most known yeasts [36] such as the template region complementary to the Telomeric Repeat Sequence (TRS), a characteristic uridine-rich pseudoknot, two helical regions known to be the binding sites of the yKu and EST1 complexes, and a uridine-rich Sm binding domain. The template region can be detected with a simple BLAST using the TRS of the organism as a BLAST query if this TRS is long enough to produce meaningful hits [9]. If the TRS is too short or unknown, there is no simple way to find the telomerase ncRNA with bioinformatics analysis alone. In this case, the combined search for all structural features occurring in the correct order in the same region of the genome could, eventually, be an alternative approach.

Contrary to the homology search pipeline, a de novo search involves looking for what we do not know, i.e., ncRNAs for which no homologous are available a priori. One of the limitations of the de novo searches is that even the most promising candidates will require experimental evidence to be validated as bona fide ncRNAs. In general, de novo ncRNA search tools [37,38] rely on the same set of assumptions: (i) the homologous sequences of a ncRNA share the same overall secondary structure; (ii) alignments of ncRNAs reveal the covariation patterns resulting from compensatory mutations and, consequently, allow the inference of the secondary structure; (iii) when applying standard folding algorithms to the alignment, the resulting minimum free energy (MFE) will be lower than the average MFE obtained by folding random sequences with the same nucleotide composition. Although none of these assumptions is always true–in particular the third one [39,40], together they can be considered as good indicators for the acceptance of a predicted candidate as a ncRNA.

In the present de novo search pipeline, we performed a whole genome Multiple Sequence Alignment (MSA) between the ten budding yeast genomes, which resulted in a set of local MSAs. The MSAs were filtered to discard too small or known protein coding regions. We then applied RNAz [37] to select the MSAs with higher probability of belonging to a ncRNA. Each selected MSA was then evaluated within its genomic context. Unannotated sequences, belonging to selected MSAs for which at least one of the sequences fall into an annotated region, are automatically given the same annotation. MSAs with no known annotation (those falling in intergenic regions for example) must be manually validated and will, eventually, represent new ncRNAs (Fig. 1B).

The automatic steps of the de novo pipeline, applied to the 10 genomes, produced 630 candidates distributed in the following way: 376 (60%) previously annotated ncRNA genes (273 of which correspond to rRNAs or tRNAs), 13 (2%) new ncRNA annotations of expected genes not found with the homology pipeline and 210 (33%) repetitive elements wrongly classified as ncRNA. The remaining 19 (3%) candidates were considered putative new ncRNA genes. Those candidates occur in intergenic regions with no previous annotations, they have a pairwise similarity higher than 50% between all sequences and display potential secondary structures supported by covariation or, at least, some compensatory mutations. The genomic location of some candidates suggests a potential regulatory role, 14 of the 19 genes occur less than 200 from the 5’ or 3’ UTR of known genes. Notice that all, except one of the candidates, were identified in only two species, a fact which prevents a more detailed analysis of the sequence variation (see Supplementary File SI 1).

Confirmation of the candidates requires experimental validation. However, the fact that 64% of the candidates could be confirmed as real ncRNA supports the assumption that, at least, some of the putative candidates correspond to real ncRNAs genes or regulatory elements. Curiously, the immediate utility of the de novo pipeline was the discovery of genes of already known families, functioning as a complement to the homology pipeline.

The A. gossypii genome was obtained from the Ashbya Genome Database (agd.vital-it.ch) [23]. The A. adeninovorans genome was obtained from (C. Neuvéglise, personal communication) and the P. sorbitophila from (V. Leh, personal communication). All other genomes and the original annotations corresponding to 297 ncRNA sequences that were used as BLAST queries are from the Génolevures Database (www.genolevures.org) [22]. From the yeast genome database (www.yeastgenome.org) [47] we obtained 86 S. cerevisiae ncRNAs also used as BLAST queries. From the RFam–ncRNA families database (rfam.sanger.ac.uk) [29] (version 9.1) we downloaded 802 covariance models for INFERNAL search.

The 383 query sequences were BLASTed against each one of the 10 genomes using the “blastall -p blastn” command (version 2.2.21) with default parameters. No query sequence was BLASTed against its own genome. The obtained candidates with E-Values < 0.1 were retained.

The homology search using INFERNAL (version 1.0) was performed in two steps. First, 83 RFam covariance models, corresponding to the ncRNAs families already known in yeasts were searched and the obtained candidates with E-Values < 0.5 were retained. Second, 719 of the remaining RFam families (corresponding to all the remaining families with exception of the viral and miRNAs families) were searched and the obtained candidates with E-Values < 0.1 were retained. All INFERNAL homology searches were performed using the “cmsearch” command with default parameters.

The retained candidates were automatically aligned with the known homologous fungal sequences using the “cmbuild” and “cmalign” commands from the INFERNAL package. Each alignment was manually validated according to the criteria described in the “Candidate Validation” section above.

For the de novo search, a whole genome MSA was performed using the TBA tool [48] according to the protocol described in “A Practical Guide to Using TBA” (www.bx.psu.edu/miller_lab) with the following tree: “((((((sace cagl) zyro) (sakl klth)) (klla ergo)) (deha piso)) (yali arad))”. The MSAs smaller than 50 nts or overlapping coding regions were discarded. The remaining MSAs were split using a sliding window of 120 nts with a step of 40 nts. We searched the resulting MSAs with RNAz and retained all candidates with reported probability higher than 0.5. Retained candidates were evaluated according to sequence conservation, secondary structure prediction, possible secondary structure signals such as covariation and compensatory mutations and genomic location with respect to neighboring genes.