Massively parallel cDNA sequencing (RNA-Seq) provides an unbiased way to study a transcriptome, including both coding and noncoding genes. Until now, most RNA-Seq studies have depended crucially on existing annotations and thus focused on expression levels and variation in known transcripts. Here, we present Scripture, a method to reconstruct the transcriptome of a mammalian cell using only RNA-Seq reads and the genome sequence. We applied it to mouse embryonic stem cells, neuronal precursor cells and lung fibroblasts to accurately reconstruct the full-length gene structures for most known expressed genes. We identified substantial variation in protein coding genes, including thousands of novel 5′ start sites, 3′ ends and internal coding exons. We then determined the gene structures of more than a thousand large intergenic noncoding RNA (lincRNA) and antisense loci. Our results open the way to direct experimental manipulation of thousands of noncoding RNAs and demonstrate the power of ab initio reconstruction to render a comprehensive picture of mammalian transcriptomes. A r t i c l e s cells (NPC) and mouse lung fibroblasts (MLF) and correctly identified the complete annotated full-length gene structures for most expressed, known, protein coding genes. The reconstruction of the three transcriptomes revealed substantial variation in protein coding genes between cell types, including thousands of novel 5′ start sites, 3′ ends or additional coding exons. Many of these variant structures are supported by independent data. We also discovered the gene structure and expression level of over 2,000 noncoding transcripts, including hundreds of transcripts from previously identified lincRNA loci, over a thousand more lincRNAs with similar properties and hundreds of multi-exonic antisense ncRNAs. We show that lincRNAs have no significant coding potential and that they are evolutionary conserved. Our results open the way to direct experimental manipulation of this new class of genes and highlight the power of RNA-Seq along with an ab initio reconstruction to provide a comprehensive picture of cell-specific transcriptomes. RESULTS We used massively parallel (Illumina) sequencing to sequence cDNA libraries from poly(A) + mRNA from ESC, NPC and MLF cells, with 76-base paired-end reads. For the ESC library, we generated a total of 152 million paired-end reads. Using a gapped aligner 13 , 93 million of these were alignable (497 Mb aligned bases, 262-fold average coverage of known protein coding genes expressed in ESC). We obtained similar numbers for the NPC and MLF libraries (Online Methods). In ESC, 76% of these reads mapped within the exonic regions of known protein coding genes, 9% were in introns of known protein coding genes, and 15% mapped in intergenic regions. We found a strong correlation between expression levels of protein coding genes as measured by RNA-Seq and Affymetrix expression arrays (r = 0.88 for all genes; Supplementary Fig. 2 ). We next developed Scripture, a genome-guided method to reconstruct the transcriptome using only an RNA-Seq data set and an (unannotated) reference genome sequence. Scripture consists of five steps (Fig. 1 A r t i c l e s junctions within the transcript, and ~30% of 76-base reads are expected on average to span an exon-exon junction. From the aligned spliced reads, we construct a 'connectivity graph' (Fig. 1d) , where two bases in the genome are connected if they are immediate neighbors either in the genomic sequence itself or within a spliced read. We use agreement with splicing motifs at each putative junction to orient the connection (edge) in the connectivity graph 9,13 (Fig. 1d). (ii) To infer transcripts, we use a statistical segmentation approach 4 and both spliced and unspliced reads to identify paths in the connectivity graph with mapped read enrichment compared to the genomic background (Fig. 1e) . This is done by scoring a sliding window using a test statistic for each region, computing a threshold for genome-wide significance, and using the significant windows to define intervals. (iii) From the paths, we construct a 'transcript graph' connecting each exon in the transcript (Fig. 1f) . Each path through the graph is directed and represents one oriented (strand-specific) isoform of the gene. Alternative spliced isoforms are identified by considering all possible paths in the transcript graph. (iv) We augment the transcript graph with connections based on pairedend reads and their distance constraints, allowing us to join transcripts or remove unlikely isoforms (Fig. 1g, below) . (v) We generate a catalog of transcripts defined by the paths through the transcript graph. Paired-end information, consisting of reads that came from the two ends of the sequenced RNA fragment, provides valuable additional information in the reconstruction. First, the presence of paired ends linking two regions shows that they appear in the same transcript; such a connection might not otherwise be apparent because low expression levels or unalignable sequence might prevent a continuous chain of overlapping sequence reads (spliced or unspliced) across the transcript. We thus augment the transcript graphs with paired-end information, where available, to (indirectly) link nodes in the graph. We use these indirect links (Fig. 1g) to add edges between disconnected graphs, add internal nodes (exons) that might have been missed within a path (transcript) and add extra support for existing edges. This refines the structure of our transcripts and increases our confidence in them, especially in weakly expressed transcripts, which are more likely to have coverage gaps. Second, the distribution of library insert sizes constrains the distance between the paired-end reads; these distance constraints can be used to infer the relative likelihood of some potential transcripts (for example, those in which the paired ends would be much closer or much further than expected). We infer the distribution of insert sizes for a given library from the position of read pairs on transcripts from those genes for which there is only a single transcript model (that is, no detectable alternative splicing; Online Methods). For example, in the ESC library, this distribution matches well with the experimentally determined sizes. Using this distribution, we assign likelihoods to each connection, filtering unlikely ones (Online Methods). We applied Scripture to our mouse ESC RNA-Seq data set and compared our reconstructions to protein coding gene annotations 15 . Scripture identified 16,389 nonoverlapping, multi-exonic transcript graphs that correspond to 15,352 known multi-exonic genes (Online Methods). Of reconstructed genes, 88.4% are covered by a single graph (no fragmentation of the reconstructed transcript) and 8.0% are covered by two transcript graphs (fragmentation of the transcript to two separate pieces in the reconstruction). Focusing on the 13,362 genes with significant expression (P < 0.05 compared with © 2010 Nature America, Inc. All rights reserved. 5 0 6 VOLUME 28 NUMBER 5 MAY 2010 nature biotechnology A r t i c l e s background coverage; see Online Methods), Scripture reconstructed the full-length structure of the longest known splice isoform (from 5′ to 3′ end, including all exons and splice junctions; Fig. 2a ) for 10,355 of them (~78%). All of our reconstructed transcripts for known multi-exonic transcripts also had the correct orientation (strand), allowing us to reconstruct genes that overlap one another on opposite strands (Fig. 2a) . Complete transcript structures were recovered across a very broad range of expression levels (Fig. 2b,c) for both single and multi-exonic genes. For example, Scripture accurately reconstructed the full-length transcript of ~73% of the known protein coding genes at the second quintile of expression, and ~94% of the genes from the top quintile. Furthermore, the average proportion of bases reconstructed for each transcript was high (Fig. 2c) . Even for the bottom 5% of expressed genes, we recovered on average 62% of each of these transcripts' bases (Fig. 2c) . For single-exon genes, we recovered on average 80% of the transcribed bases. We obtained similar results in the other two cell types (19,835 and 20,407 transcript graphs for 14,212 and 13,351 known genes in NPC and MLF, respectively). Most of the genes that were not fully reconstructed are those with low expression; it should be possible to reconstruct most of these by generating more RNA-Seq data. The few highly expressed genes that were not fully reconstructed are either the result of alignment artifacts caused by recent processed pseudogenes or stem from novel transcriptome variations, missing from the current annotation (explored in detail below). Given that most of the Scripture reconstructions of protein coding genes were accurate, we next investigated the differences between the reconstructed transcriptome and the known gene annotations (Supplementary Table 1 ). We focused on transcripts with (i) novel 5′ start sites; (ii) novel 3′ ends; and (iii) previously unidentified exons within the transcriptional units of known protein coding genes. In each category, we first discuss below the reconstructed transcripts in ESC and then consider the results for the NPC and MLF. 1. Alternative 5′ start sites are supported by H3K4me3 marks. We found 1,804 transcripts in ESC that match the annotated 3′ end but have an alternative 5′ start site, derived from an extra exon not overlapping the annotated first exon. We distinguish between internal alternative 5′ start sites (1,397 cases; Fig. 3a) , which occur downstream of the annotated start, and external alternative 5′ start sites (407 cases; Fig. 3b ), which occur upstream of the annotated start. Ninety percent of the internal 5′ start sites and 75% of the external 5′ start sites contained an H3K4me3 modification, a mark of the promoter region of genes 16 (Supplementary Fig. 3 ). These alternative start sites are on average 21 kb upstream of the annotated site, substantially revising the annotated promoters. Notably, ~60% of the transcripts with an alternative start site (internal or external) had no reconstructed isoform starting at the annotated 5′ start site. We obtained similar results from NPC and MLF (Fig. 3a, A r t i c l e s 5′ start sites (2,302 supported by H3K4me3 in their respective tissues), and 807 external 5′ start sites in at least one cell type. In particular, 33% of these novel 5′ ends are likely active in ESCs but not in NPCs or MLFs. 2. Alternative 3′ untranslated regions are supported by polyadenylation motifs. There are 551 (~4%) ESC-reconstructed transcripts with an alternative 3′ end downstream of the annotated 3′ end (mean distance 30 kb downstream, Fig. 3c ). Of these, 275 (~50%) showed evidence of a polyadenylation motif within the novel 3′ exon, which is only slightly lower than for annotated 3′ ends (60%) and much higher than for randomly chosen size-matched exons (6%). The frequency of the polyadenylation motif supports the accuracy of the reconstruction. To conservatively distinguish between upstream (early) termination and incomplete reconstruction, we designated novel 3′ ends only in those cases that did not overlap any of the known exons in the annotated transcript and that contained complete 5′ start sites. We identified 759 transcripts with upstream 3′ ends in ESC (Fig. 3d) ; 44% of them contained a polyadenylation motif, supporting their biological relevance. For most (90%) of these transcripts, Scripture also reconstructed an isoform that contained the annotated 3′ end. We obtained similar results for NPC and MLF (Fig. 3c,d , right; Supplementary Table 1) . Altogether, we identified 940 downstream 3′ ends and 1,850 upstream 3′ ends in at least one cell type. 3. Additional coding exons are highly conserved and preserve ORFs. We found 534 transcripts in ESC with at least one extra, previously unannotated internal coding exon spliced into annotated protein coding transcripts (Fig. 3e) . These transcripts contained 588 novel internal exons, ranging in length from 6 bp to 3.5 kb (median, 111 bp; 20-80% quantiles, 60-224 bp). Of these extra exons, 322 (54.5%) were present in all versions of the reconstructed transcript in ESC. Most (83%) of these novel exons maintain the reading frame of the transcript and are as highly conserved as known coding exons (Supplementary Fig. 4) , consistent with their coding capacity. We validated the presence of the novel exons within five of five tested transcripts, using reverse transcription followed by PCR (RT-PCR) followed by Sanger sequencing (Online Methods). We obtained similar results in MLF (124 transcripts, 144 exons) and NPC (325 transcripts, 363 exons) (Fig. 3e, right) . A majority of exons (~70%) were present in all versions of the reconstructed transcript within a cell type. Altogether, we identified 960 novel internal exons in at least one cell type (Fig. 3e, right) . We next turned to identifying the gene structures of transcripts expressed from known lincRNAs loci. We had previously identified 317 lincRNA loci on the basis of K4-K36 domains in ESC cells 4 . When applied to ESC RNA-Seq data, Scripture reconstructed multi-exonic gene structures for 250 (78.8%) of them (Fig. 4) . This is comparable to the proportion (78.5%) reconstructed for protein coding genes with K4-K36 domains in ESC. Scripture reconstructed 87% (160 of 183) of ESC lincRNAs for which we previously identified an RNA hybridization signal from tiling microarrays. We discuss possible reasons for the few remaining discrepancies in Supplementary Note 2. The reconstructed lincRNA transcripts in ESC have on average 3.7 exons, an average exon size of 350 bp and an average mature spliced size of 3.2 kb (in comparison, protein coding genes have on average 9.7 exons, exon length of 291 bp and length of 2.9 kb). The Scripture-identified strand information for each lincRNA is consistent with that inferred from the location of H3K4me3 modification and with the orientation determined from a strand-specific RNA-Seq library which we generated independently (Online Methods). Most lincRNAs likely represent 5′ complete transcripts based on overlap with H3K4me3 (82%) and 3′ complete transcripts based on presence of a polyadenylation motif (~50%, comparable to 60% for protein coding genes and far above background of 6%). Similarly, Scripture successfully reconstructed lincRNA gene structures for K4-K36 lincRNA loci in MLF and NPC (232 of 289 in MLF and 224 of 270 in NPC). Most are likely 5′ complete (69% in MLF and 81% in NPC based on overlap with H3K4me3) and many may be 3′ complete based on detectable 3′ polyadenylation sites (18% in MLF and 37% in NPC). In addition, we successfully reconstructed another 116 lincRNAs previously identified only in mouse embryonic fibroblasts but which were now reconstructed in at least one of the other three cell types. Altogether, we identified gene structures for 609 previously defined lincRNA loci in at least one of the three cell types. In addition to the previously identified lincRNAs, we found another 1,140 multi-exonic transcripts that map to intergenic regions (591 in ESC, 318 in MLF, and 528 in NPC; Fig. 5 ). Most of these transcripts do not seem to encode proteins, and are designated as noncoding, on the basis of their codon substitution frequency (CSF) scores 17, 18 (Online Methods) across the mature (spliced) RNA transcript (88%; Fig. 5a ) and on the lack of an open reading frame (ORF) larger than 100 amino acids (80%; Fig. 5b ). Careful review of the © 2010 Nature America, Inc. All rights reserved. 5 0 8 VOLUME 28 NUMBER 5 MAY 2010 nature biotechnology A r t i c l e s remaining ~12% revealed 66 loci that are likely to be novel protein coding genes (high CSF score, ORF >200 amino acids and very high evolutionary conservation; Supple mentary Fig. 5) . Most of the novel lincRNA loci were not identified in our previous study owing to the stringent criteria we imposed when using chromatin maps to identify lincRNAs. Specifically, we required that a K4-K36 domain extend over at least 5 kb and be well separated from the nearest known gene locus 4 . Indeed, most novel intergenic transcripts (76%) were enriched for a K4-K36 domain (a comparable proportion as that for expressed protein coding genes) but failed to meet the size and distance criteria or could not be identified at a genome-wide significance level (without knowing their locus a priori). On average, the genomic loci of the novel lincRNAs are closer to neighboring genes and have smaller sizes (~3.5 kb average), and the transcripts are shorter (859 bp). Of the lincRNAs that did not have a chromatin signature that reached genomewide significance, ~40% showed chromatin modifications enriched at a nominal significance level (compared to 57% for protein coding genes). On average, the lincRNAs are expressed at levels that are readily detectable, albeit somewhat lower than those of protein coding genes. The median expression level of the reconstructed lincRNAs, as estimated by reads per kilobase of exonic sequence per million aligned reads (RPKM; see Online Methods) was approximately one-third of the expression of protein coding genes (Fig. 5d) , with ~25% of lincRNAs having expression levels higher than the median level for protein coding genes (Fig. 5d) . The novel lincRNAs identified in this study are expressed at somewhat lower levels than those from chromatin identified loci, consistent with the fact that chromatin enrichment is positively correlated with expression levels (Fig. 5d) . We compared the novel lincRNA genes to a collection of ~35,000 mouse cDNA and found evidence that ~43% of our lincRNAs were present in this collection 1 . This is comparable to the reported fraction (40%) of known transcripts covered by the same cDNA catalog 1 . The remaining lincRNAs are found in our study but not in the comparison catalog. These were likely previously missed owing to the different cell types and limited coverage of the previous study 1 . The reconstructed full-length gene structures of lincRNAs allow us to accurately assess their evolutionary sequence conservation in each exon and in small windows. To this end, we identified the orthologous sequences for each lincRNA across 29 mammals and estimated conservation by a metric (ω; Online Methods) reflecting the total contraction of the branch length of the evolutionary tree connecting them 19 . We calculated ω over the entire lincRNA transcript, as well as over individual exons. On the basis of our high-resolution gene structures, the lincRNA sequences show greater conservation than random genomic regions or introns (Fig. 5c) , comparable to eight known functional lincRNAs [20] [21] [22] , and lower than protein coding exons. The results are consistent with our previous estimates of conservation 4 . Interestingly, conservation levels are indistinguishable between the chromatin-defined lincRNAs 4 and the novel ones identified only in this study (Fig. 5c) , consistent with membership in the same class of functional large ncRNA genes. These conservation levels are considerably higher than those reported for a previous catalog of large noncoding RNAs 1 . We also determined the specific regions within each lincRNA that are under purifying selection and thus likely to be functional, by computing ω within short windows (Online Methods). On average, 22% of the bases within the lincRNAs lie within conserved patches (comparable to the value of 25% for the eight known functional lincRNAs, much higher than the 7% for intronic bases and lower than the 77% for protein coding bases, Supplementary Fig. 6 ). These conserved patches provide a critical starting point for functional studies 23 . A substantial fraction (~41%) of the novel lincRNAs reconstructed in at least one cell type show evidence for expression in at least two of the A r t i c l e s three cell types. This is comparable to the 45% of the previously identified lincRNAs present in at least two out of the three cell types. In contrast, 80% of expressed protein coding genes are expressed across two of the three cell types. This is not merely a result of the lower overall expression of lincRNAs, as the fraction of cell type-specific lincRNAs is higher than that of tissue specific protein coding genes in every expression quantile (Supplementary Fig. 7) . Thus, lincRNAs are likely to be more tissue specific than protein coding genes. A substantial portion of lincRNA loci also produce alternative spliced isoforms. For example, within ESC we identified two or more alternative spliced isoforms for 25% of lincRNA genes, comparable to the 30% for protein coding genes (15% of lincRNAs in MLF have alternative spliced isoforms, and 14.7% in NPC). Altogether, 28.8% of the 1,749 lincRNA loci had evidence for alternative isoforms in any of the three cell types. Scripture reconstructed hundreds of transcripts that overlap known protein coding gene loci but are transcribed in the opposite orientation and likely represent antisense transcripts. To determine orientation, we required that any identified antisense transcript be multi-exonic (Online Methods). Using these criteria, we identified 201 antisense multi-exonic transcripts in ESC (Fig. 4b) ; these transcripts had an average five exons per transcript and an average transcript size of 1.7 kb. On average, the antisense transcripts overlapped the genomic locus of the sense protein coding gene by 1,023 bp (83% of the locus), and most (64%) overlapped at least one sense exon, but this overlap was substantially lower (766 bp, 48% of the transcript exons). Some of these antisense transcripts (79, ~40%) were identified by a previous cDNA sequencing study 1, 24 , but most (122, ~60%) were previously unidentified. Most (~85%) antisense transcripts were non-protein coding by both ORF analysis (Fig. 5b) and CSF scores (Fig. 5a) . Four of the newly identified antisense transcripts had a large, conserved open reading frame and are likely novel, previously unannotated protein coding genes. We validated the reconstructed ESC antisense transcripts by three independent sets of experimental data. (i) Most of the antisense loci carried an H3K4me3 mark at their 5′ end (Fig. 4b) , consistent with their independent and antisense transcription (for example, 64% of the 164 transcripts where it was possible to detect an independent H3K4me3 mark because the 5′ end of the antisense transcript did not overlap the 5′ ends of the sense gene). (ii) We generated and sequenced a strand-specific library in ESC (17.5 million Illumina reads; Online Methods), and found a significant (P < 0.05) number of reads on the antisense strand in >90% of cases (the remaining are likely missed in this limited sequencing owing to lower expression). (iii) We confirmed five of five tested antisense transcripts using RT-PCR to unique exons of the antisense transcript (Online Methods) followed by Sanger sequencing. We obtained similar results for antisense transcripts in MLF and NPC (112 and 202 multi-exonic antisense transcripts, respectively). Altogether, we identified 469 antisense transcripts expressed in at least one cell type, only 125 of which (27%) were previously identified in large-scale sequencing of mouse cDNAs 24 . The remaining 344 (73%) were unidentified by the previous study, likely reflecting the distinct cell types used in that study and the limited coverage of previous catalogs. The 469 antisense transcripts are expressed at levels comparable to those of the novel lincRNAs (Fig. 5d) but show substantially lower sequence conservation. Indeed, the antisense ncRNAs showed very little evolutionary conservation as estimated by the ω metric for the portions that do not overlap protein coding exons on the sense strand, suggesting that the antisense ncRNAs are a distinct class from the lincRNAs (Fig. 5c) . Despite the availability of the genome sequence of many mammals, a comprehensive understanding of the mammalian transcriptome has been an elusive goal. In particular, the computational tools needed to reconstruct all full-length transcripts from the wealth of short read data were largely missing. A recent study proposed to overcome this limitation experimentally by using very long reads (for example, 454 sequencing) as a scaffold for short read reconstruction 25 . This is applicable, albeit at a substantial cost, for highly expressed genes but would require extraordinary depth to cover more weakly expressed ones. Here we present Scripture, a new computational method to reconstruct a mammalian transcriptome with no prior knowledge of gene annotations. Scripture relies on longer reads that span splice junctions to connect discontiguous (spliced) segments and resolve multiple splice isoforms, and uses paired-end information to refine these transcripts. Scripture can identify short but strongly expressed transcripts as well as transcripts with much lower expression for which there is aggregate evidence along the entire transcript length. Although Scripture does rely on a reference genome sequence, many of its components can also be used in the development of methods for assembly of transcripts from read data only. We applied Scripture to RNA-Seq data from pluripotent ESCs and differentiated lineages and showed that we can accurately reconstruct most expressed, annotated protein coding genes, at a broad range of expression levels, as well as uncover many new isoforms in the protein coding transcriptome. This variation may have key regulatory roles, defining new cell type-specific promoters, untranslated regions and protein coding exons. We used Scripture's sensitivity and resolution to reconstruct the gene structures and strand information of hundreds of lincRNAs and multi-exonic antisense transcripts, many of which are only moderately expressed. Scripture identified over a thousand lincRNAs across the three cell types studied. Most of the lincRNAs identified were not previously found by classical large-scale cDNA sequencing 1 . Many of these lincRNAs could not be reliably identified solely on the basis of chromatin structure owing to their proximity to protein coding genes or their short genomic lengths. Overall, we found that the ratio of expressed protein coding to noncoding genes in these cell types was ~10:1 but that the total number of RNA molecules was more heavily biased toward the protein coding fraction (~30:1), results similar to previous observations 26 . Scripture identifies precise gene structures for most previously found lincRNA loci (as well as for the newly discovered ones), a prerequisite for further studies. For example, we used these to identify the specific regions within each lincRNA that are under purifying selection (conservation), a starting point for experimental and computational investigation. Taken together, our results highlight the power of ab initio reconstructions to annotate a genome, to discover transcriptional variation within known protein coding genes and to provide a rich catalog of precise gene structures for noncoding RNAs. The next step is clearly to apply this approach to a wide range of mammalian cell types, to obtain a comprehensive picture of the mammalian transcriptome. Methods and any associated references are available in the online version of the paper at http://www.nature.com/naturebiotechnology/.
Scripture REF construct a splice graph from the mapped reads and reconstructs transcripts corresponding to all possible paths in this graph.
Ab initio reconstruction of cell type–specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs
Nature Biotechnology
Nature Biotechnology
Biology Medicine