Estimating the Fraction of Non-Coding RNAs in Mammalian Transcriptomes

Abstract

Recent studies of mammalian transcriptomes have identified numerous RNA transcripts that do not code for proteins; their identity, however, is largely unknown. Here we explore an approach based on sequence randomness patterns to discern different RNA classes. The relative z-score we use helps identify the known ncRNA class from the genome, intergene and intron classes. This leads us to a fractional ncRNA measure of putative ncRNA datasets which we model as a mixture of genuine ncRNAs and other transcripts derived from genomic, intergenic and intronic sequences. We use this model to analyze six representative datasets identified by the FANTOM3 project and two computational approaches based on comparative analysis (RNAz and EvoFold). Our analysis suggests fewer ncRNAs than estimated by DNA sequencing and comparative analysis, but the verity of our approach and its prediction requires more extensive experimental RNA data.

Keywords

randomness test fraction model putative non-coding RNA

Introduction

The rapid progress in large-scale screening of cellular transcriptional output (termed the transcriptome) has dramatically increased the repertoire of transcripts expressed in mammalian cells (Bertone et al. 2004; Carninci et al. 2003; Carninci et al. 2005; Cheng et al. 2005; Kampa et al. 2004; Kapranov et al. 2002; Numata et al. 2003; Ota et al. 2004). Full-length sequencing data (Carninci et al. 2005; Ota et al. 2004) indicate that about 50% of the transcripts are likely to be non-coding RNAs (ncRNAs). Known ncRNAs are involved in a variety of cellular functions, such as gene regulation, rRNA modification, splicing, RNA editing, mRNA degradation. The currently-identified ncRNAs are far less abundant than protein-coding RNAs in mammalian transcriptomes. The large number of putative ncRNAs, however, leads us to speculate that the number of ncRNAs will continue to grow. For those putative ncRNAs, although experimental techniques, such as Northern blot, microarray analysis, and RT-PCR, are providing expression data (Carninci et al. 2005; Cheng et al. 2005; Kampa et al. 2004; Ota et al. 2004; Ravasi et al. 2006), most biological functions remain unknown. Determining the identity of these putative ncRNAs is important: Do these transcripts correspond to genuine ncRNAs with biological functions, or to other RNAs that may be biological or experimental artifacts, non-functional transcripts, or transcriptional noise (Huttenhofer et al. 2005; Johnson et al. 2005; Soares and Valcarcel, 2006).

The FANTOM (Functional Annotation of Mouse) database, which aims to generate the transcriptional landscape of the mouse genome, has identified 34,030 manually annotated putative ncRNAs in the FANTOM3 release, among which 2886 sequences are annotated by the most stringent criteria (Carninci et al. 2005). Only a small population of the putative ncRNAs is known ncRNAs, such as miRNAs and snoRNAs, and others are unknown RNA transcripts.

Computational approaches such as RNAz (Washietl et al. 2005b; Washietl et al. 2005a) and EvoFold (Pedersen et al. 2006) have employed evolutionary conservation of secondary structures, thermodynamic stabilities, and phylogenetic sequence analysis to predict tens of thousands of conserved RNA secondary structural elements in the human genome. The largest dataset predicted by RNAz contains 91,676 RNA structures conserved at least in four mammals (human, mouse, rat, and dog) in non-coding regions (coding exons are removed). The EvoFold program predicts 48,479 RNA structures from various genomic locations (coding, UTR, intronic, and intergenic), of which 517 sequences are classified as ncRNA candidates by the program. However, less than 1% of the predictions by RNAz and EvoFold are known ncRNAs (Pedersen et al. 2006; Washietl et al. 2005a) and others are new RNA folds.

Despite various predictions concerning the biological functions of putative ncRNAs (Huttenhofer et al. 2005; Johnson et al. 2005; Mattick and Makunin I.V. 2006; Soares and Valcarcel, 2006), progress has been slow in identifying associated biological roles in cells. An exhaustive functional characterization of unknown RNA transcripts by experiments is prohibitive. However, determining the fraction of genuine non-coding RNAs in those putative ncRNAs from either experimental or computational data will advance our understanding of the composition of mammalian transcriptomes and the general importance of ncRNAs for cellular function. An approach we explore here is a systematic assessment based on statistical features for the known ncRNA sequences.

Because biological sequences are not purely random, statistical tests might help screen sequences of interest. Previous works on statistical properties of nucleotide sequences have been geared toward characterization of coding/non-coding regions (Almirantis, 1999; Herzel and Grosse, 1997; Kugiumtzis and Provata, 2004; Peng et al. 1992), sequence complexity (Abel and Trevors, 2005; Adami and Cerf, 2000), and evolutionary patterns (Dehnert et al. 2005b; Dehnert et al. 2005a). Moreover, it was found that ribosomal RNAs could be discriminated from random sequences by using statistical measures (Almirantis, 1999). Thus, quantifying sequence characteristics may reveal some aspects of functional features (such as general class type) based on sequence properties.

Here, we use relative z-scores of missing motifs to analyze characteristic global features for specific nucleotide sequence classes. We find that the relative z-score, derived from the monkey test (Marsaglia and Zaman, 1993) for assessing random number generators (RNGs), helps classify six nucleotide sequence classes into three clusters in decreasing degree of randomness: (1) genome/intergene/intron, (2) mRNA/ncRNA, and (3) repeat. We use this characteristic relative z-score of the ncRNA class to predict the ncRNA fraction in putative ncRNAs as determined in the FAN-TOM3 database and by computational programs RNAz and EvoFold.

Materials and Methods

Nucleotide sequence classes

We compare the degree of randomness of the following nucleotide sequence classes: genome, intergene, intron, mRNA (or coding sequence), and ncRNA. The genome, intergene, intron, and mRNA classes are generated from RefSeqs (Pruitt et al. 2005) in Supplementary Table S1. The ncRNA class contains 7,698 representative ncRNA genes from Noncode (Liu et al. 2005), RNAdb (Pang et al. 2005), Rfam (Griffiths-Jones et al. 2005), and European ribosomal RNA database (Wuyts et al. 2004) (Table 1). Noncode and RNAdb have low redundancy, while the remaining three sources have multiple copies. In order to determine the effect of multiple copies in randomness analysis, we create two other versions that have fewer or no rRNAs, tRNAs, or spliceosomal RNAs. We also create an additional class, the repeat sequence class, as a control: a long concatenated “sequence” is generated by repeating a 1,024-nt random sequence 2,048 times to form a 2,097,152-nt repetitive sequence.

Table 1.

The composition of the training ncRNA class and other two constructs.

Group	Training set	Version 2	Version 3	Source
Noncode ncRNA	4251	4251	4251	Noncode
RNAdb ncRNA	2204	2204	2204	RNAdb
rRNA	129	0	10	EID^a
tRNA	78	0	10	Rfam
spliceosomal RNA	28	0	10	Rfam
Total sequence number	6690	6455	6485
Relative z-score (mean)	6.93	6.77	7.33
Relative z-score (std)	0.33	0.11	0.17

EID: European ribosomal RNA database.

To generate the fraction model, we use the ncRNA class and mouse genomic, intergenic, and intronic sequences (Supplementary Table S2).

Mouse genomic and intergenic sequences have lengths of 10⁸ nt or more, so we cut them into 2000-nt segments which can conveniently be manipulated and shuffled.

Putative ncRNA datasets

We analyze six putative ncRNA datasets identified by experimental and computational methods (Table 2). The two experimental datasets come from the FANTOM3 database of full-length cDNAs that do not code for proteins (Carninci et al. 2005). The computational datasets are formed by predictions from programs RNAz (Washietl et al. 2005b) and EvoFold (Pedersen et al. 2006). The sequences predicted by RNAz are grouped into three datasets based on different number of conserved organisms and p-values: (1) the set1.P0.5 dataset is predicted with sequences conserved at least in human, mouse, rat, and dog at p > 0.5; (2) the set1.P0.9 dataset is predicted at p > 0.9; (3) the set2.P0.5 dataset is predicted with sequences conserved at least in human, mouse, rat, dog, and chicken at p > 0.5. The sequences predicted by EvoFold are conserved in eight species (human, chimpanzee, mouse, rat, dog, chicken, pufferfish, and zebrafish).

Table 2.

The six putative ncRNA datasets.

Dataset	Sequences	Total length (nt)	Ave length (nt)
FANTOM3 putative ncRNA	34,030	67,856,244	1,994
FANTOM3 stringent putative ncRNA	2,886	4,535,792	1,572
RNAz set1.P0.5	91,676	12,474,689	136
RNAz set1.P0.9	35,985	5,475,570	152
RNAz set2.P0.5	20,391	2,798,941	137
EvoFold	48,479	1,869,205	39

The randomness (monkey) test

The monkey test assesses sequence randomness based on a χ² distribution of goodness-of-fit measure (Σ(OBS-EXP)²/EXP) for the overlapping k- and k-1-letter words (Marsaglia, 2005; Marsaglia and Zaman, 1993). This test can assess both uniformity and independence of random sequences. We employ an applied version of the monkey test, the DNA test (Marsaglia and Zaman, 1993), in our study. In the DNA test, the number of missing k-letter words in a long sequence approximately follows a normal distribution whose mean and variance depend on the word and alphabet sizes.

Thus, the z-score can quantify the degree of randomness:

z = \frac{m i s s i n g - m e a n}{s t d}

(1)

where the number of missing words is an observed value, mean is the theoretical average number of missing words, and std is the standard deviation.

In our analysis, we use the default parameters of the DNA test: it counts 10-letter words for a 2,097,152 letter sequence; the expected number of missing words behaves like a normal distribution with mean 141,909 and standard deviation 339. Figure 1 illustrates the conversion from nucleotide sequences to numeric sequences. Different word size and sequence length can also be applied to the DNA test.

Figure 1

Application of the DNA monkey test to biological sequence analysis. The DNA test counts 10-letter missing words from a 4-letter alphabet (A, C, G, T/U). Three overlapping 10-letter words are shown in blue, green, and red, respectively.

Applying the Monkey Test to Dna Sequences

The DNA test requires long input sequences on the order of 10⁶ nt (2,097,152 nt) (Marsaglia and Zaman, 1993). We adopt the following procedure to analyze sequences whose lengths are shorter than 2,097,152 nt (Fig. 2): (1) we randomly shuffle original nucleotide sequences in a group by the Mersenne Twister RNG (Matsumoto and Nishimura, 1998); (2) we concatenate shuffled sequences into one sequence, our “concatenated sequence;” (3) we cut that concatenated sequence into segments in the length of 2,097,152 nt; (4) we generate a random sequence with the same dinucleotide composition as the corresponding concatenated biological sequence; (5) we submit both biological and random sequences to the DNA test. Following this procedure, we generate at least 100 concatenated sequences and corresponding random sequences for a given sequence group and submit them to the DNA test. These pairs of z-scores result in relative z-scores which will be described in detail in Results. The distribution of the relative z-score of a given sequence group is summarized by mean, standard deviation, and range.

Figure 2

Sequence manipulation scheme for randomness analysis.

Thermodynamic Analysis

In addition to randomness test, we use a complementary thermodynamic analysis (the melting temperature and energy landscape analyses) to screen short putative ncRNAs in the FANTOM3 dataset. The melting temperature, defined by the peak of the heat capacity, is predicted by RNAheat (Hofacker et al. 1994; Mccaskill, 1990) from the Vienna RNA folding Package (Hofacker et al. 1994) (Version 1.6). The energy landscape, defined by the base pair dissimilarity “distance” between the minimum energy structure and each suboptimal structure, is measured by the “Valley Index” (Kitagawa et al. 2003). We use a 90% confidence ellipse in 2D plots of free energy vs. melting temperature and of free energy vs. Valley Index to determine stability of test structures as described in our previous work (Laserson et al. 2005). This method implies an error rate of 10%, i.e. random sequences pass the test in 10% of cases. The thermodynamic tests are applied to ten ncRNA families, namely tRNA, 5S rRNA, 5.8S rRNA, 6S RNA, SRP RNA, SL1 RNA, U6 RNA, UnaL2, snoRNA, and His3 in Rfam (Griffiths-Jones et al. 2005) and 151 putative ncRNAs (<400 nt) in FANTOM3 (Carninci et al. 2005).

Results

The relative z-score measures the degree of randomness

We employ the DNA test to generate a randomness measure, the z-score (Eq. 1) (Marsaglia and Zaman, 1993). Biological sequences often fail the DNA test because the number of missing motifs for biological sequences is often greater than the theoretical mean value (Eq. 1). Random sequences generated with nonuniform nucleotide compositions also fail the DNA test. Since dinucleotide composition has become standard to measure the background effect (Goni et al. 2004; Workman and Krogh, 1999), we use the random dinucleotide sequence as a control to the corresponding biological sequence. Thus, we define the relative z-score for the degree of randomness to be:

z_{r e l a t i v e} = \frac{z_{b i o s e q}}{z_{d i - r a n d o m}}

(2)

where z_bioseq is the z-score for a biological sequence and z_di-random is the z-score for the random dinucleotide control sequence. Higher-order background information (e.g. tri-nucleotide, tetra-nucleotide) can be considered in the future for model improvement. We also name the z-score (e.g. z_bioseq and z_di-random) the absolute z-score to differentiate it from the relative z-score. A relative z-score close to 1 indicates that the dinucleotide composition is the major contributor to non-randomness; a value larger than 1 indicates that other sequence factors affect the non-randomness; a value near zero means that single-nucleotide composition rather than dinucleotide more likely causes the non-randomness. Below, we show that the relative z-score can distinguish various classes of biological sequences.

The relative z-score classifies nucleotide sequence classes into three clusters

We then use the relative z-score to assess the degree of randomness of six nucleotide sequence classes: genome, intergene, intron, mRNA, ncRNA, and artificial repeat. The relative z-scores are obtained by running the DNA test on one hundred 2,097,152-nt concatenated sequences for each class from the three phylogenetic domains.

In the three-domain collection (Fig. 3a), the genome class displays the broadest distribution in sequence randomness (Table 3): 86% of relative z-scores fall into the region of 0.9–5.0. The intergene and intron classes have narrow relative z-score distributions overlapping with the peak of the genome distribution. The mRNA and ncRNA classes overlap with one another and have a lower degree of randomness than the previous three classes. The artificial repeat class spans a wide range and has the lowest degree of randomness among the six nucleotide sequence classes (Fig. 3e). As expected, regular motif patterns of the repeat class lead to a low degree of randomness. Thus, we find that the relative z-score partitions the six nucleotide sequences from the three domains into three clusters: (1) genome, intergene, and intron; (2) mRNA and ncRNA; and (3) repeat sequences.

Table 3.

Relative z-scores of the six nucleotide sequence classes in the three domains.

Class	Collection		Archaea		Bacteria		Eukarya
	Range	Mean (std)	Range	Mean (std)	Range	Mean (std)	Range	Mean (std)
Genome	0.9–11.7	2.9 (2.2)	1.2–11.5	3.6 (2.5)	1.2–15.2	3.0 (2.3)	1.0–5.8	2.2 (1.0)
Intergene	1.3–2.2	1.6 (0.2)	1.2–1.3	1.2 (0.02)	4.1–5.1	4.6 (0.2)	1.4–2.2	1.7 (0.2)
Intron	1.6–2.0	1.8 (0.1)	/	/	/	/	1.6–2.0	1.8 (0.1)
mRNA	5.7–7.3	6.3 (0.2)	2.8–3.5	3.1 (0.2)	4.0–5.1	4.6 (0.3)	3.7–5.1	4.4 (0.3)
ncRNA	6.3–7.7	6.9 (0.3)	/	/	8.7–9.5	9.1 (0.1)	6.1–6.5	6.2 (0.1)
Repeat	24.5–230.8	58.0 (32.9)	24.5–230.8	58.0 (32.9)	24.5–230.8	58.0 (32.9)	24.5–230.8	58.0 (32.9)

Figure 3

The degree of randomness (relative z-scores) for different sequence classes in the three phylogenetic domains by the DNA test: (a) the three-domain collection; (The intron class corresponds only to eukaryotes). (b) Archaeal dataset; (c) Bacterial dataset; (d) Eukaryotic dataset; (e) the repeat (control) dataset. The distribution of the relative z-score of each sequence class is estimated by the density function in R (R Development Core Team 2006). The color legend in the inset of (a) applies to (b), (c), and (d).

After examining the randomness trends in all three domains, we analyze that in the three domains of life separately. In Eukarya, the six nucleotide sequence classes still form three clusters as they do in the three-domain collection (Fig. 3d). However, the genome/intergene/intron cluster has a narrower range than the one in the three-domain collection and the mRNA/ncRNA cluster is not as compact as that one in the three-domain collection. In Archaea, the genome, intergene, and mRNA classes form a cluster, with the mRNA class being less random than the intergene class (the intron and ncRNA classes are not available) (Fig. 3b). In Bacteria, the genome, intergene, and mRNA classes form a cluster, and the ncRNA class forms another one (Fig. 3c). The limited ncRNA dataset for bacteria, containing only a small number of available ncRNAs except for tRNAs and rRNAs, produces a relative z-score distribution mostly reflecting the randomness features of tRNA and rRNA families. Therefore, the clustering pattern among the six nucleotide sequence classes is maintained in Eukarya but changed in Archaea and Bacteria due to sequence bias in datasets. It also shows that the ncRNA class is less random than the other biological sequence classes in the three-domain collection and separate domains.

Overall, the relative z-score can partition two sequence clusters—genome/intergene/intron (more random) and mRNA/ncRNA—in the three-domain collection and Eukarya. For the ncRNA class, it is characterized by a distinct relative z-score in the three-domain collection, Bacteria and Eukarya. Its lower randomness than the genome/intergene/intron cluster may be explained by low motif diversity of the ncRNA class likely caused by the sequence conservation within some ncRNA families and RNA's preference for specific recurrent motifs.

Putative ncRNAs are not all functional

We now assess the six putative ncRNA datasets listed in Table 2 using the relative z-score. The total length of the EvoFold dataset is 1,869,205 nt which is shorter than the required length (2,097,152 nt) of the DNA test (Marsaglia and Zaman, 1993), so randomly selected ncRNAs are added to reach the length requirement. The relative z-score (1.437) of this “pseudo” EvoFold dataset is an estimate of the true value. Another “pseudo” EvoFold dataset created with additional genomic sequences has almost the same relative z-score (1.436). The DNA test result shows that none of the six datasets have a relative z-score close to the ncRNA class (Fig. 4). Instead, the six datasets form non-overlapping relative z-score distributions which fall in the genome/intergene/intron cluster. In order of decreasing degree of randomness, we have EvoFold, RNAz set2.P0.5, FANTOM3 putative, RNAz set1. P0.5, FANTOM3 stringent, and RNAz set1.P0.9.

Figure 4

The degree of randomness of the six putative ncRNA datasets measured by the DNA test. The relative z-score distribution of the six datasets is denoted as follows: (a) EvoFold, (b) RNAz set2.P0.5, (c) FANTOM3 putative, (d) RNAz set1.P0.5, (e) FANTOM3 stringent and (f) RNAz set1.P0.9.

Three explanations to the much smaller relative z-scores of these sequences compared to the ncRNA class are possible. First, these putative ncRNAs may consist of a mixture of real ncRNAs and other types of RNA transcripts (i.e. “RNA noise”). Second, while most of these putative ncRNAs may indeed have biological functions like known ncRNAs, the relative z-score is a poor indicator (e.g. ncRNAs are group I intron-like with a low relative z-score). Third, most of these putative ncRNAs may not have biological functions. Below, we develop a model to probe each possibility in turn.

The Fraction Model for Estimating the Proportion of Ncrnas in Test Datasets

Whatever the fraction of true ncRNAs in these datasets, we attempt to formulate an approximate fraction model based on the assumption that the Fantom3 and computational datasets are mixtures of true ncRNAs and “RNA noise”. We thus construct a fraction model with various ncRNA to RNA noise ratios to match the relative z-scores of the test datasets.

In our fraction model, the ncRNAs are sampled from the training ncRNA dataset defined in Methods which represents diverse RNA families in Noncode, RNAdb and several prominent families from Rfam and European ribosomal RNA database (Table 1). This choice of training dataset is reasonable because the Fantom3 dataset is highly heterogenous (broad length distribution, multiple experimental data sources (Carninci et al. 2005)) and the computational datasets represent whole genome scans. We thus assume the ncRNAs in the test datasets to be as diverse as those in the training dataset. Further analysis shows that the sequence properties of the training and test datasets are quite similar (Supplementary Figs. S1–S4). On the other hand, we may miss a part of new ncRNAs whose relative z-scores are different from the representative score of the training dataset.

To model RNA noise or background sequences, we use mouse intronic, intergenic, and genomic sequences. RNA noise may include biological and experimental artifacts, non-functional transcription, and transcriptional noise (Huttenhofer et al. 2005; Johnson et al. 2005; Soares and Valcarcel, 2006). Because our analysis shows that the putative ncRNAs fall in the cluster of genome/intergene/intron (Fig. 4), such RNA noise may come from genomic, intergenic, and/or intronic sources. Though RNA noise could also originate from mRNAs, this possibility is less likely because existing protein-coding sequences have been filtered (Carninci et al. 2005) or removed (Washietl et al. 2005a) a priori. Thus, we use mouse genomic, intergenic, and intronic sequences (Supplementary Table S2) to model RNA noise in mammalian transcriptomes under various ratios.

We consider three submodels based on the ncRNA/noise partitioning: (1) ncRNA/intron; (2) ncRNA/intergene; (3) ncRNA/genomic sequence. In each submodel, eleven ncRNA/noise ratios (0:10, 1:9, 2:8,…, 10:0) are applied to create concatenated sequences whose relative z-scores correspond to specific ncRNA fractions (from 0% to 100% with a resolution of 10%). The data points for the ncRNA fraction vs. the relative z-score defines a function f (z) where z is the mean of relative z-scores for a given ncRNA fraction. As expected, the three functions generated from the three submodels are monotonically increasing and converging to 1 as the fraction of ncRNAs increases (Fig. 5). The two functions f₂(z) and f₃(z) almost overlap with one another because the intergenic sequences cover about 68% of the genomic sequences (Supplementary Table S2). The ncRNA/intron submodel predicts lower ncRNA fractions than the other two submodels because the intron class has a larger average relative z-score than the genome and intergene classes.

Among the three submodels considered above, the third (ncRNA/genome) appears to simulate the putative ncRNA datasets most accurately. This expectation is supported by the following evidences: transcription of the mammalian genomes is wide (e.g. over 60% of the mouse genome is transcribed (Carninci et al. 2005)); and putative ncRNAs are transcribed from diverse genomic locations including coding regions, 5’- and 3’-UTRs, introns, and intergenic regions (Carninci et al. 2005; Cheng et al. 2005; Pedersen et al. 2006; Washietl et al. 2005a). Therefore, based on given relative z-scores of putative ncRNA datasets, we can use the fraction model (the ncRNA/genome submodel preferred) to predict the ncRNA fractions in mammalian transcriptomes.

In all three submodels, each ncRNA fraction has a relative z-score distribution represented by error bars which quantify the uncertainty of our prediction of ncRNA fractions (Fig. 5). At each data point, the error of f (z) is estimated by the standard deviation of z_mean. Based on the eleven data points for each submodel, the average error of ncRNA fraction is about ±5%. Given this resolution of our model, we predict that less than 5% ncRNAs exist in a dataset with the mean relative z-scores below 2.

Figure 5

The model used to assess the ncRNA fraction in the FANTOM3, RNAz and EvoFold datasets. The mean relative z-score of the six datasets is shown in dashed lines in the same order as Fig. 4. Error bars show standard deviations of relative z-scores. The four predictions for the FANTOM3, the FANTOM3 stringent, the RNAz set1.P0.5 and the RNAz set1.P0.9 datasets by f₃(z) are highlighted in bullets.

Predicting the ncRNA Fraction in the Six Representative Datasets

We collect six putative ncRNA datasets from the three sources: the FANTOM3 database, and the predictions of RNAz and EvoFold. The overlap between one another is quite small: the RNAz program only found 781 conserved RNA folds out of more than 15,000 putative ncRNAs in the FAN-TOM2 database (Washietl et al. 2005a); about 6000 folds have at least partial overlaps between the predictions of RNAz and EvoFold.

We use the mean relative z-score of the two FANTOM3 putative ncRNA datasets to interpolate the ncRNA fraction with functions f_1–3(z) (see Fig. 5 and Table 4). The relative z-score of the FANTOM3 dataset (dashed line c in Fig. 5) is out of the predicting region of f₁(z), which suggests that its ncRNA fraction may be too small to be detected by the ncRNA/intron submodel, so the ncRNA fraction is estimated to be less than 5%. The preferred model suggests that the FANTOM3 stringent dataset with higher confidence for ncRNA annotation may contain a much higher fraction of genuine ncRNAs (47%) than the FANTOM3 dataset (18%).

Table 4.

The ncRNA fractions predicted by the model. Relative z-scores are shown in mean values and standard deviations.

Dataset	Relative z	f₁(z)			f₂(z)		f₃(z)
		%	Sequence #	%	Sequence #	%	Sequence #
Fantom putative	2.63 (0.05)	<5%	<1,701	22%	7,487	18%	6,125
Fantom stringent	3.32 (0.06)	9.7%	280	50%	1,443	47%	1,356
RNAz set1. P0.5	3.08 (0.04)	<5%	<4,584	42%	38,503	39%	35,754
RNAz set1. P0.9	3.49 (0.04)	21%	7,557	53%	19,072	52%	18,712
RNAz set2. P0.5	2.21 (0.01)	<5%	<1,020	<5%	<1,020	<5%	<1,020
EvoFold^a	~1.44 (0.02)	<5%	<2,424	<5%	<2,424	<5%	<2,424

The relative z-score of the EvoFold dataset is estimated by concatenated sequences mixed with EvoFold predictions and known ncRNAs because the total length of the EvoFold dataset is shorter than the required length of the DNA test.

We also predict the ncRNA fraction of the four computationally-identified datasets predicted by RNAz (the set1.P0.5, set1.P0.9, and set2.P0.5 datasets) and EvoFold (the EvoFold dataset) (Table 4). Our preferred model predicts that set1.P0.5 and set1.P0.9 contain 39% and 52% genuine ncRNAs, respectively. The other two datasets (RNAz set2. P0.5 and EvoFold) have high sequence randomness (i.e. small relative z-scores); they locate out of the predicting region of all three submodels. We thus propose that these two datasets contain less than 5% genuine ncRNAs. The four computationally-predicted datasets, although all based on comparative analysis, contain different ncRNA fractions.

The corresponding number of genuine ncRNAs in the six datasets also varies greatly: 6,125 for FANTOM3, 1,356 for FANTOM3 stringent, 35,754 for RNAz set1.P0.5, 18,712 for RNAz set1. P0.9, < 1,020 for RNAz set2.P0.5, and <2,424 for EvoFold. Since some conserved RNA folds form clusters and share common RNA transcripts (Pedersen et al. 2006), the actual number of genuine ncRNAs in the computational predictions may be smaller than the above numbers.

Thermodynamic Analysis of the FANTOM3 Putative Ncrnas

In addition to developing the fraction model to assess putative ncRNAs, we also employ thermodynamic tests (involving free energy, melting temperature and energy landscape) to analyze 151 short (<400 nt) putative ncRNAs in FANTOM3; only short RNAs are used for this analysis to reduce errors in secondary folding algorithms. Among these short putative ncRNAs, 23% are thermodynamically stable (Fig. 6). We also analyze known ncRNA families: tRNA, 5S rRNA, 5.8S rRNA, 6S RNA, SRP RNA, SL1 RNA, U6 RNA, UnaL2, snoRNA, and His3. Except for the non-structural snoRNA family, the other nine ncRNA families have higher passing rates than the FANTOM3 subset. The average passing rate of the ten ncRNA families is 60%. This analysis shows that ncRNAs are more stable than both the random sequences and the FANTOM3 subset. It also suggests that 23% of short putative ncRNAs in FANTOM3 may have biological function, in agreement with the relative z-score analysis (18%).

Figure 6

Thermodynamic analysis of selected sequences of the FANTOM3 putative ncRNA dataset (<400 nt) and ten known ncRNA families. The passing rate, tested sequence number and dataset name are shown above the passing rate bar.

Systematic errors in the fraction model

Since the FANTOM3 data are assembled from diverse sources (tissues, experimental conditions, various laboratories, etc), the simplest and reasonable model for the FANTOM3 data is a collection of diverse RNA families similar to the composition of our training dataset. Still, systematic errors of the fraction model can arise due to possible differences in the composition of ncRNAs in training and test datasets. For example, the test datasets could be enriched with specifically low or high relative z-score ncRNAs. To simulate such datasets, we generate four biased fraction models containing 10% and 25% of group I intron (low relative z-score, data not shown) and rRNA (high relative z-score, data not shown) sequences in the ncRNA partition of the ncRNA/genome model, labeled f₄, f₅, f₆ and f₇, respectively. Figure 7 shows that the four biased models have different error ranges. For example, the predicted ncRNA fraction in the FANTOM3 stringent dataset increases from 47% to 56% and 72% for f₄ and f₅ (10% and 25% group I intron), but decreases to 42% and 34% for f₆ and f₇ (10% and 25% rRNA) (Supplementary Table S4). Thus, for datasets with enriched 10% random-like (e.g. group I intron) or structural ncRNAs (e.g. rRNA), the error would be ~10%. The two enriched 25% models are less likely due to the high diversity of our test datasets (Supplementary Figs. S2–S4). These four biased models can partially simulate scenarios where the test datasets have compositions different from that of the training ncRNA dataset. Expected increase of experimentally characterized ncRNAs in the near future will help improve these fraction models.

Figure 7

The six submodels used to simulate systematic errors in our fraction model. The f₄ and f₅ are constructed by 10% and 25% group I intron sequences in the ncRNA partition of the ncRNA/genome model, respectively; and the f₆ and f₇ are constructed by 10% and 25% rRNA sequences in the ncRNA partition, respectively. The submodels f₈ and f₉ are ncRNA/genome models using version 2 and version3 ncRNA reference dataset (Table 1), respectively. The mean relative z-score of the six test datasets is shown in dashed lines in the same order as Fig. 4.

We have also estimated the uncertainties arising from changing the composition of training dataset when some ncRNA families are removed. Table 1 shows two additional training datasets (versions 2 and 3) with low and no representation of rRNA, tRNA and spliceosome families. These datasets lead to only 0.5–1.5 standard deviations from the mean relative z-score of the original training dataset, implying errors of less than 5% in predicted ncRNA fraction.

Discussion

We have proposed here a fraction model for assessing ncRNA content of sequence datasets based on a sensitive relative z-score for measuring the degree of sequence randomness. Our fraction model relies on the relative z-score to help distinguish the genome/intergene/intron cluster from the mRNA/ncRNA cluster. It assumes that the relative z-score can discriminate the ncRNA class and the genome/intergene/intron cluster (Fig. 3a), and that the training and test datasets have similar sequence conservation patterns (Figs. S1–S4). These are reasonable first-order approximations but only with much more data can these assumptions be validated.

Clearly, sequence randomness reflects some functional features. For example, ncRNAs likely contain recurrent motifs (e.g. GNRA, UNCG) underlying RNA's modular architecture (Hendrix et al. 2005; Leontis et al. 2006). Of course, there are numerous other aspects that must be considered. Moreover, all statistical features analyzed here rely on currently available datasets. Though many more RNAs await discovery, the modularity of RNA and reliance on a few recurring motifs suggest that our approach is worth considering. It is interesting to recall that when the ribosome structure was solved, only a few new tertiary motifs emerged despite expectations to the contrary.

Other computational studies of ncRNA, most of which are based on comparative genomic analysis, such as QRNA (Rivas and Eddy, 2001), RNAz (Washietl et al. 2005b) and EvoFold (Pedersen et al. 2006), are limited by the requirement of high sequence conservation across species. However, many ncRNAs exhibit low sequence conservation (Pang et al. 2006). In contrast, the relative z-score assesses the randomness degree of any sequences whether conserved or not. Nevertheless, this approach is not applicable to single ncRNA sequences, because those sequences are three to four orders of magnitude shorter than the required length (2,097,152 nt) for reliable statistical analysis by the current approach. Though it may be possible to reduce the sequence length by changing the word size and then recalculating the mean value and standard deviation of Eq. 1, longer sequence lengths are more reliable for the monkey test application.

Our fraction model, if valid, predicts that less than 52% of putative ncRNAs predicted by FAN-TOM3 and computational approaches are functional. This is not consistent with the speculation that most of the putative ncRNAs are functional (Mattick and Makunin IV. 2006) but agrees with other computational studies. For example, the EvoFold program predicted that 517 out of 48,479 conserved RNA folds are ncRNA candidates (Pedersen et al. 2006), which agrees with our prediction that less than 5% (<2,424 folds) are genuine ncRNAs. The RNAz program screened the dataset of FANTOM2 putative ncRNAs and only identified 781 out of more than 15,000 putative ncRNAs having conserved RNA secondary structures (Washietl et al. 2005a). This number is much less than our predicted number (6,125) in the FANTOM3 dataset partly because the relative z-score assesses both conserved and non-conserved sequences. In addition, false positives from computational predictions can contribute to over-counting of genuine ncRNAs. For example, the high false positive rate, 28.9% (p = 0.5), for the RNAz program suggests that only a part of predictions may be real ncRNAs.

Moreover, our fraction model assumes that putative ncRNAs contain genuine ncRNAs and background noise. At least three potential errors may be introduced into the model: (1) limited amount of training data for ncRNAs; (2) limited source of background noise; and (3) “contamination” of mRNAs in a tested dataset. For the first type of error, as the number and diversity of ncRNA families increase, the verity and precision of our fraction model can be assured and improved. The second type of error arises from the limited knowledge of transcriptome noise. Available experimental data indicate that over 60% of the mammalian genomes are transcribed (Carninci et al. 2005), but annotation is an ongoing process. Finally, as shown in our randomness analysis, the mRNA and ncRNA classes share a same region of randomness in the three-domain collection and Eukarya.

In conclusion, based on a first-level approximation, we suggest that fewer putative ncRNAs in datasets identified by experiments or computational approaches may exist. Clearly, our understanding of the general importance of ncRNAs in mammalian transcriptomes will advance as the number of genuine ncRNAs is better estimated rather than speculated. Our fraction model, if validated, might also be used to investigate ncRNA fraction of other putative ncRNA databases. The relative z-score may also help guide the detection of novel aptamers and ribozymes through design of sequence pools, an area of current work (Gevertz et al. 2005; Kim et al. 2007).

Footnotes

Acknowledgements

The work was supported by Human Frontier Science Program (HFSP) and by a Joint NSF/NIGMS initiative in Mathematical Biology (DMS-0201160). We thank Stefan Washietl and Jakob Skou Pedersen for kindly providing us sequences predicted by the RNAz and EvoFold programs, respectively.

Supplementary Material

Table S4.

The systematic errors caused by biased ncRNA training datasets. Four submodels f₄–f₇ are created to simulated biased training data using 10% and 25% group I intron and rRNA sequences, repectively. The submodel labels are same as Fig. 7.

Dataset	f₅(z)	f₄(z)	f₃(z)	f₆(z)	f₇(z)
Fantom putative	31%	22%	18%	16%	12%
Fantom stringent	72%	56%	47%	42%	34%
RNAz set1 P0.5	60%	47%	39%	35%	28%
RNAz set1 P0.9	78%	61%	52%	47%	38%
RNAz set2.P0.5	<5%	<5%	<5%	<5%	<5%
EvoFold^*	<5%	<5%	<5%	<5%	<5%

Same dataset as described in Table 4.

References

Abel

D.L.

and Trevors

J.T.

2005. Three subsets of sequence complexity and their relevance to biopolymeric information. Theor. Biol. Med. Model, 2.

Adami

and Cerf

N.J.

2000. Physical complexity of symbolic sequences. Physica. D-Nonlinear Phenomena., 137: 62–9.

Almirantis

1999. A standard deviation based quantification differentiates coding from non-coding DNA sequences and gives insight to their evolutionary history. J. Theor. Biol., 196: 297–308.

Bertone

, Stolc

, Royce

T.E.

2004. Global identification of human transcribed sequences with genome tiling arrays. Science, 306: 2242–6.

Carninci

, Kasukawa

, Katayama

2005. The transcriptional landscape of the mammalian genome. Science, 309: 1559–63.

Carninci

, Waki

, Shiraki

2003. Targeting a complex transcriptome: The construction of the mouse full-length cDNA encyclopedia. Genome Res., 13: 1273–89.

Cheng

, Kapranov

, Drenkow

2005. Transcriptional maps of 10 human chromosomes at 5-nucleotide resolution. Science, 308: 1149–54.

Dehnert

, Helm

W.E.

and Hutt

M.T.

2005a. Information theory reveals large-scale synchronisation of statistical correlations in eukaryote genomes. Gene, 345: 81–90.

Dehnert

, Plaumann

, Helm

W.E.

2005b. Genome phylogeny based on short-range correlations in DNA sequences. J. Comput. Biol., 12: 545–53.

10.

Gevertz

, Gan

H.H.

and Schlick

2005. In vitro RNA random pools are not structurally diverse: A computational analysis. RNA, 11: 853–63.

11.

Goni

J.R.

, de la Cruz

and Orozco

2004. Triplex-forming oligonucleotide target sequences in the human genome. Nucleic Acids Res., 32: 354–60.

12.

Griffiths-Jones

, Moxon

, Marshall

2005. Rfam: annotating non-coding RNAs in complete genomes. Nucleic Acids Res., 33: D121–D124.

13.

Hendrix

D.K.

, Brenner

S.E.

and Holbrook

S.R.

2005. RNA structural motifs: building blocks of a modular biomolecule. Q. Rev. Biophys., 38: 221–43.

14.

Herzel

and Grosse

1997. Correlations in DNA sequences: The role of protein coding segments. Physical Review E, 55: 800–10.

15.

Hofacker

I.L.

, Fontana

, Stadler

P.F.

1994. Fast folding and comparison of RNA secondary structures. Monatshefte fur Chemie, 125: 167–88.

16.

Huttenhofer

, Schattner

and Polacek

2005. Non-coding RNAs: hope or hype? Trends Genet., 21: 289–97.

17.

Johnson

J.M.

, Edwards

, Shoemaker

2005. Dark matter in the genome: evidence of widespread transcription detected by microarray tiling experiments. Trends Genet., 21: 93–102.

18.

Kampa

, Cheng

, Kapranov

2004. Novel RNAs identified from an in-depth analysis of the transcriptome of human chromosomes 21 and 22. Genome Res., 14: 331–42.

19.

Kapranov

, Cawley

S.E.

, Drenkow

2002. Large-scale transcriptional activity in chromosomes 21 and 22. Science, 296: 916–9.

20.

Kim

, GAN

H.H.

and Schlick

2007. A computational proposal for designing structured RNA pools for in vitro selection of RNAs. RNA, 13: 478–92.

21.

Kitagawa

, Futamura

and Yamamoto

2003. Analysis of the conformational energy landscape of human snRNA with a metric based on tree representation of RNA structures. Nucleic Acids Res., 31: 2006–13.

22.

Kugiumtzis

and Provata

2004. Statistical analysis of gene and intergenic DNA sequences. Physica A-Statistical Mechanics and Its Applications, 342: 623–38.

23.

Laserson

, Gan

H.H.

and Schlick

2005. Predicting candidate genomic sequences that correspond to synthetic functional RNA motifs. Nucleic Acids Res., 33: 6057–69.

24.

Leontis

N.B.

, Lescoute

and Westhof

2006. The building blocks and motifs of RNA architecture. Curr. Opin. Struct. Biol., 16: 279–87.

25.

Liu

C.N.

, Bai

B.Y.

, Skogerbo

2005. NONCODE: an integrated knowledge database of non-coding RNAs. Nucleic Acids Res., 33: D112–D115.

26.

Marsaglia

2005. Monkeying with the goodness-of-fit test. Journal of Statistical Software, 14:

27.

Marsaglia

and Zaman

1993. Monkey tests for random number generators. Computers and Mathematics with Applications, 26: 1–10.

28.

Matsumoto

and Nishimura

1998. Mersenne Twister: A 623-dimensionally equidistributed uniform pseudorandom number generator. ACM Trans on Modeling and Computer Simulation, 8: 3–30.

29.

Mattick

J.S.

and Makunin

I.V.

2006. Non-coding RNA. Hum. Mol. Genet., 15: R17–R29.

30.

Mccaskill

J.S.

1990. The equilibrium partition-function and base pair binding probabilities for RNA secondary structure. Biopolymers, 29: 1105–19.

31.

Numata

, Kanai

, Saito

2003. Identification of putative noncoding RNAs among the RIKEN mouse full-length cDNA collection. Genome Res., 13: 1301–6.

32.

Ota

, Suzuki

, Nishikawa

2004. Complete sequencing and characterization of 21,243 full-length human cDNAs. Nat. Genet., 36: 40–5.

33.

Pang

K.C.

, Frith

M.C.

and Mattick

J.S.

2006. Rapid evolution of noncoding RNAs: lack of conservation does not mean lack of function. Trends Genet., 22: 1–5.

34.

Pang

K.C.

, Stephen

, Engstrom

P.G.

2005. RNAdb—a comprehensive mammalian noncoding RNA database. Nucleic Acids Res., 33: D125–D130.

35.

Pedersen

J.S.

, Bejerano

, Siepel

2006. Identification and classification of conserved RNA secondary structures in the human genome. PLoS Computational Biology, 2: 251–62.

36.

Peng

C.K.

, Buldyrev

S.V.

, Goldberger

A.L.

1992. Long-Range Correlations in Nucleotide-Sequences. Nature, 356: 168–70.

37.

Pruitt

K.D.

, Tatusova

and Maglott

D.R.

2005. NCBI Reference Sequence (RefSeq): a curated non-redundant sequence database of genomes, transcripts and proteins. Nucleic Acids Res., 33: D501–D504.

38.

R. Development Core Team, 2006. R: A language and environment for statistical computing. R. Foundation for Statistical Computing.

39.

Ravasi

, Suzuki

, Pang

K.C.

2006. Experimental validation of the regulated expression of large numbers of non-coding RNAs from the mouse genome. Genome Res., 16: 11–9.

40.

Rivas

and Eddy

S.R.

2001. Noncoding RNA gene detection using comparative sequence analysis. BMC Bioinformatics, 2.

41.

Soares

L.M.M.

and Valcarcel

2006. The expanding transcriptome: the genome as the “Book of Sand”. EMBO J., 25: 923–31.

42.

Washietl

, Hofacker

I.L.

, Lukasser

2005a. Mapping of conserved RNA secondary structures predicts thousands of functional noncoding RNAs in the human genome. Nat. Biotechnol., 23: 1383–90.

43.

Washietl

, Hofacker

I.L.

and Stadler

P.F.

2005b. Fast and reliable prediction of noncoding RNAs. Proc. Natl. Acad. Sci. U.S.A., 102: 2454–9.

44.

Workman

and Krogh

1999. No evidence that mRNAs have lower folding free energies than random sequences with the same dinucleotide distribution. Nucleic Acids Res., 27: 4816–22.

45.

Wuyts

, Perriere

, and de Peer

Y.V.

2004. The European ribosomal RNA database. Nucleic Acids Res., 32: D101–D103.