Genome-wide read mapping

Of 260 million high-quality mRNA-seq reads obtained, 142 million (53%) were matched to a unique site in the rat genome; these were termed uniquely matching reads (UMR) and were used to quantify the expression of known genes (Fig. 1) and to discover new exons (Fig. 2). There were 104 million (39%) nonmatching reads (NMR), which were used for splice unction (SJ) discovery (Fig. 3), and there were 20 million (8%) multiply matching reads (MMR) (i.e., aligning to more than one site in the reference genome), which were useful in completing the annotation of complex genes (Fig. 4).

Open in a separate window

Figure 2.

Mapping novel exons in nonannotated regions of the rat genome. (A) Typical UMR distribution. UMR mapped across nonannotated regions of the genome in nonrandom patterns, often resembling exon–intron structure. Depicted are UMR (each gray line = one read) mapped across a nonannotated region of chromosome 5. For comparison, UMR across an annotated region of chromosome 10 are shown along with the annotation of known exons illustrating correspondence. (B) Aggregation of UMR into read “clusters” resembling exons. A 100-bp “sliding window” was moved across the genome demarcating the beginning (UMR present, i.e., “filled” window) and end (UMR absent, i.e., “empty” window) of “UMR clusters.” Resulting UMR clusters differed in read density. Clusters consisting of high piles of hundreds of reads provide strong evidence for an exon, while clusters with few reads appear indeterminate. (C) Classifier for UMR clusters: exon versus no-exon. The newly defined UMR clusters were dichotomized into “exons” and “no exons” using average read density (read coverage) as a classifier. The ROC curve shown here was obtained by applying the classifier to an exploratory subset of 10,367 (annotated) genes, varying read density from 0.25 to 100. Sensitivity and specificity are plotted for 26 different read densities. The classifier was found to be a very precise test with an area under the curve of 0.97. At a read density of 4 (inflection point of the curve indicated as red circle), the sensitivity of the test was 92% and the specificity was 97%. The favorable test characteristics were confirmed in the validation set (sensitivity 91%, specificity of 97%). Applying the procedure to the full data set, 123,066 UMR exon clusters were found. (D) Location of UMR exon clusters. 10,464 new exons were found (i.e., UMR clusters with density >4 in a nonannotated region); the remaining clusters overlapped known exons belonging either to the 10,367 genes quantified above or to other annotated genes.

Open in a separate window

Figure 3.

Agnostic splice junction (SJ) mapping. (A) Agnostic, i.e., de novo genome-wide SJ discovery from split reads (25 bp + 25 bp). NMR were split into 25-mer halves, which were matched independently to the genome reference sequence. Reads whose halves each matched uniquely to the genome defined SJs de novo; these were termed uniquely mapping splice junction reads (sjUMR). (B) SJ validation on known genes and their application to define connections of novel exons and exon borders. (Green) SJ-spanning uniquely mapping reads (sjUMR; mapped as described in the Methods section). Most SJs were supported by multiple sjUMR (28.6 on average), as depicted here by piles of green lines. (Top panel) A known gene (Pitpnb ENSRNOG00000000665) is shown here, demonstrating how sjUMR mapped precisely to the border of exons as annotated in the reference genome. (Bottom panel) sjUMR connecting novel UMR clusters to groups, thereby defining new genes. In some instances, UMR cluster borders fall into introns (as a result of the inclusiveness of the sliding-window algorithm), as depicted here (right sides of second and third cluster). In those cases, exon borders were defined precisely by sjUMR. Note (top panel, second exon) that very short (<50 bp) exons that were undetectable by UMR could be defined through mapping of SJ reads.

Open in a separate window

Figure 4.

Integrating information from UMR, sjUMR, MMR, and NMR read classes for complex genome (re-)annotation. (A) UMR- and sjUMR-based definition of new exon groups. A novel gene candidate was identified on chromosome 9 through several UMR clusters (satisfying the classifier criterion of a read density > 4), which were connected through splice junctions defined by sjUMR. Expression of exons was coregulated (suppressed by 70% after SNL). An open reading frame encoding 409 amino acids was noted lacking a start and stop codon, suggesting that an incomplete gene fragment was identified consisting of nine exons. (B) sjUMR validation of a low-read density UMR cluster as a novel exon. Consistent with the limited sensitivity of the UMR cluster classifier (91%), a UMR read cluster with a read density of 3.9 was initially not classified as an exon (because its density was <4); but sjUMR-based SJs defined it unambiguously as a novel exon connecting the above gene fragment to a series of 18 3′-located UMR clusters. Note that the cluster density was low because a region of scattered intron reads was added to the read cluster by the sliding window, an imprecision that was rectified by agnostic splice site mapping, which defined precise exon borders. Through the step depicted in this panel, the candidate gene was extended to 28 exons. (C) NMR-derived contig bridge. A faulty or missing section in the reference genome can be indicated by runs of N (ambiguous bases) as encountered 3′ of the 28 exons assembled above (50 scattered N). Such missing sections in the reference are also an important reason why some mRNA-seq reads do not match. Accordingly, we found that in such cases, contigs of NMR can be assembled bridging faulty sequences, as shown here. As a result of the NMR contig-bridge, the gene candidate could be expanded to 49 exons. (D) MMR candidate exons discriminated by sjUMR. Exon copy numbers >1 in the reference sequence result in clusters of MMR (blue) matching collectively to more than one site. A case of two such MMR clusters in close proximity (within <2 kb; note duplicated region indicated by yellow bar) was observed 5′ of the above gene candidate. As shown here, the ambiguity was resolved by sjUMR; this led to a 5′ extension of the gene candidate by another 15 exons, which resulted in a complete gene with start and stop codon consisting of 64 exons. (E) A summary of the novel 64-exon gene is shown encoding a 3313-amino-acid-long protein (14,084-bp mRNA). While the study was under way, a homologous mouse protein with important CNS relevance was reported (“large previously unknown protein”) (Lu et al. 2009), termed UNC-80, which serves as a substance-P and neurotensin coreceptor. Down-regulation of the rat UNC-80 rat homolog after SNL may contribute to allodynia through alterations in peptide activity in the DRG.

We obtained 2.6 × 10⁸ (274,622,530) 50-bp-long cDNA sequence reads equaling 13 Gb, or ~4.7 whole-genome equivalents, of sequence data. Of these 7,083,741 (2.58%) were discarded because of poor base call quality scores or unknown bases. The remaining 267,538,789 reads were aligned to the full rat reference genome (RGSC 3.4; mitochondrial genome included) using the ELAND aligner software (Illumina) and allowing for ≤2 mismatches per 32 consecutive bases. We took advantage of the extended read length of 50 bp by requiring that the first and the last 32 bp both mapped to the same site (offset by 18 bp) in order to declare a match. We found that 142,572,750 (53.3%) of all processed reads matched a unique genomic site (termed UMR), and were used to quantify the expression of known genes (Fig. 1) and to discover new exons (Fig. 2). Additionally, 20,408,306 (7.6%) reads matched to multiple (i.e., >1) sites in the genome (termed MMR). Finally, 104,557,733 (39%) were NMR. Nonmatching mRNA-seq reads may occur for several reasons, including pre-mRNA splicing, which creates sequences that are not present in the genome, i.e., spanning splice junctions (alignment of NMR to splice sites is discussed in a separate section below); exclusion of the Y chromosome and omission of other sequences from the current rat reference genome assembly RGSC 3.4; sequencing errors; and genetic differences (single nucleotide polymorphisms, insertions, deletions) between the rat strain employed (Sprague Dawley) and the reference genome strain (Brown Norway). Accordingly, the set of NMR was used for SJ discovery (Fig. 3), as described in detail below.

The average read density across known mRNA transcripts was 1.25-fold higher for 3′ exons than 5′ exons. A 3′ bias is to be expected with any mRNA detection method that employs a poly(A) purification step. In our case, the median difference was small, only 25%, because the library construction protocol employed (chemical mRNA fragmentation followed by random priming for cDNA synthesis) was designed to achieve even coverage.

Transcriptome reprogramming

As a first step in our analysis, the expression of 10,367 Ensembl known protein-coding genes (out of 17,738 annotated) was quantified. Raw UMR counts for each gene (provided in Supplemental Table 1) correlated highly when controls were compared (biological replicates). For these pairs Pearson correlation coefficients were always r = 0.99 (Supplemental Table 2). The correlation among SNL-control pairs was lower (average 0.82), as was to be expected if gene expression was altered after SNL. For subsequent analyses, read counts were normalized to the total number of reads obtained for each sample (multiples of 10⁶). Results are shown in Figure 1B–F (the remaining 7371 Ensembl known protein-coding genes were expressed at low levels). As shown in Figure 1, when SNL animals were compared with controls at the 2-wk time point, 1268 to 1415 genes were induced, and 772 to 775 genes were suppressed. A false discovery rate (FDR) of 0.5% (52 of 10,367 genes) in each group (i.e., “induced” and “suppressed”), resulted in a true-positive estimate of 1289 (12.4%) induced and 721 (7.0%) suppressed genes (statistical procedures and computations are detailed further in Methods). In other words, nearly 20% of the L4 DRG transcriptome was reprogrammed to a new expression level in response to SNL of the neighboring (L5) nerve root.

Alterations persist long term

Examining the same set of 10,367 genes, we found that 1432 to 1631 genes were induced and 1041 to 1160 were suppressed 2 mo after SNL. Fewer total mRNA sequences were obtained for this time point (one lane compared with four for the early time point), resulting in a slight increase in the measurement variance and an observed FDR of 0.7% (51 to 87 false-positive genes vs. 52 expected). Accounting for the observed (i.e., higher) FDR, 1458 genes (14.1%) were induced and 1027 (7.5%) were suppressed at 2 mo.

Of all genes induced ≥2-fold at 2 mo, 96% were induced ≥1.4-fold at 2 wk. Similarly, of all genes suppressed by a factor of ≥2 (i.e., to less or equal than half of controls) at 2 mo, 97% were suppressed by at least 1.4-fold at 2 wk. No genes were found to be induced or suppressed by a factor of ≥2 at one of the two time points but regulated in the opposite direction at the alternate time point. While the study was not designed as a time-series experiment and not powered to detect single genes (or very few) displaying dynamic regulation, the results strongly suggest that the large majority of transcriptome alterations persisted stably from 2 wk to 2 mo. Overall, >20% of measured genes were altered at 2 mo, demonstrating that transcriptional reprogramming persisted in neuropathic pain matching the chronicity of the condition (Fig. 1E,F).

Extent of DRG gene induction and suppression detectable by mRNA-seq

Genes were induced up to >100-fold and were suppressed as much as 10-fold in SNL animals (compared to sham controls). Ten genes were induced 32- to 121-fold, 42 genes were induced ≥10-fold, and 212 were induced ≥4-fold; 165 genes were suppressed ≥4-fold (at 2 wk and/or 2 mo). The magnitude and breadth of changes documented appear generally more extensive than alterations detectable in comparable microarray studies of the DRG in pain (Bonilla et al. 2002; Costigan et al. 2002; Valder et al. 2003; Davis-Taber and Scott 2006; Rodriguez Parkitna et al. 2006). The principle of quantification by mRNA-seq is digital counting, which is free of background hybridization problems (or similar artifacts limiting sensitivity) and does not suffer from saturation of probes (known to diminish precision when quantifying high-expressing genes in microarrays) (Wold and Myers 2008). The resulting superior dynamic range of mRNA-seq established by the original studies (Cloonan et al. 2008; Marioni et al. 2008; Mortazavi et al. 2008; Nagalakshmi et al. 2008; Sultan et al. 2008; Wilhelm et al. 2008; Wold and Myers 2008; Wang et al. 2009) was found here to be immensely informative when applied to the nervous system.

Mapping nonannotated portions of the DRG transcriptome

Previous “genome-wide” nervous system transcriptome studies have been based on microarrays and, therefore, have been restricted to known mRNA. mRNA-seq is not bound to such limitations. We found 6,224,094 RNA reads (4.7% of UMR) mapping to nonannotated regions of the rat genome. Most of these occurred in dense aggregates of many reads (dozens to thousands) that were separated by long nontranscribed regions (Fig. 2A). To define which read-aggregates resembled novel exons, we developed a classifier. First, mRNA-seq reads were joined into “clusters” by a 100-bp sliding window (Fig. 2B). Next, newly defined clusters were compared with exons of well-expressed protein-coding genes (a subset of the 10,367 genes described above) containing 97,646 exons. One-tenth of the exons were used for test development. Test conditions were systematically varied requiring different average read densities (from <1 to >100 UMR per 50 bp) to declare that a cluster was an exon. Cluster classification results were compared with the actual exon annotation, recording the frequency with which an annotated exon was detected (true-positives) and how often an exon was erroneously predicted in an intron region (false-positives). A receiver operating characteristic (ROC) curve (Fig. 2C) demonstrated that this was an effective classifier (area under the curve 0.97). Optimal precision was achieved at a density cut-off of >4 with a sensitivity of 92% and a specificity of 97% in the exploratory exon set (10% of all gold standard exons). The sensitivity was 91%, and the specificity was 97% when the classifier was applied to the validation set (the remaining 90% of exons). Clusters meeting these criteria were termed “exon clusters.” Applying this procedure to the entire genome, we found 92,628 exon clusters in the 10,367 previously considered known protein-coding genes, 19,974 exon clusters in other annotated genes (e.g., weakly expressed protein-coding genes, RNA genes, or genes predicted on the basis of orthologs), and 10,464 new exon clusters in nonannotated regions of the genome (Fig. 2D). The new exon clusters had a total length of 7,290,658 bp, i.e., taken together they covered 0.27% of the genomic sequence.

We evaluated the possibility that contrast could be used as an additional criterion for classifying read clusters. Within the gold-standard regions, true-positive clusters had on average a steeper flank than false-positives. But the overlap was considerable (as seen in Supplemental Fig. 1), suggesting that only about two out of three clusters would be classified correctly. Improving a strong classifier by combining it with a weak one usually fails to improve the overall result—as one of us previously found in another area of informatics research, speech recognition (P. Beyerlein, unpubl.).

“Classic” genes—i.e., regions with well-defined start, end, and exon boundaries—may not be the only source of RNA. “Dark matter” (Johnson et al. 2005) regions of the genome may be transcriptionally active; for example, it has been reported that up to 94% of the human genome can be transcribed (The ENCODE Project Consortium 2007). We found 968,726 UMR that mapped neither to annotated areas of the genome nor to newly defined clusters. These were labeled “dark matter reads.” Dark matter reads were rare, comprising only 0.71% of UMR, and were dispersed sparsely, covering 1.4% of the genome sequence at an average density of 1.2×. Given this sparsity, statistical inferences cannot be derived from this portion of the data set, leaving DRG dark matter transcription an uncertain possibility that will require even larger sequencing efforts before we can arrive at a critical assessment.

Sequencing depth in the present study was ample for discovery of new exons (resembling those of annotated genes) as demonstrated by the test characteristics provided above (i.e., sensitivity and specificity in the high-90s percentage range) and by the high average read density, 43×, for the newly discovered 10,464 exon clusters.

Agnostic splice junction discovery

Pre-mRNA (cis-)splicing is fundamental to defining transcriptomes. Splicing determines how many different proteins (or regulatory RNA) are derived from a set of neighboring exons, and which of them should be grouped together in a single semantic unit, i.e., a gene. A remaining stated goal of mRNA-seq, for which no published solution yet exists, is to determine pre-mRNA SJs de novo, i.e., solely from the experimental data (mRNA-seq reads) without use of existing annotations or information from other sources or experiments.

In previous mRNA-seq studies, splicing was detected either by mapping reads to known SJs or to limited sets of hypothetical junctions derived by alternative (in silico) assembly of known exons (Pan et al. 2008; Wang et al. 2009). These approaches were the only practicable methods when mRNA-seq reads were very short, i.e., 25–36 bp. Such methods are termed biased because only the limited set of SJs that are predicted to exist can be experimentally detected. Accordingly, previous studies could not discover new SJs in nonannotated sections of the genome nor could they discover unsuspected SJs involving known genes (e.g., splicing of newly discovered exons to known genes). In our study, the length of sequences was longer (50 bp), thanks to improved technology, which allowed us to implement a computationally novel approach for unbiased SJ discovery that overcame the shortcomings of previous approaches. The new approach is termed “agnostic” because it operates free of preconceived notions about exon borders, mapping SJ de novo based solely on the information represented by the mRNA-seq data.

We reasoned that some reads would span SJ symmetrically, i.e., the splice site would be exactly in the middle of the read. In these instances the 25-bp sequence on each side of the junction would map to a different location on the same chromosome indicative of the positions and precise borders of two exons belonging to the same gene (Fig. 3A). Accordingly, the two 25-bp ends (each of the two halves) of every 50-bp read obtained in the study were mapped independently against the entire genomic reference sequence, allowing us to identify 4,539,891 reads containing an SJ exactly in the middle (1.7% of all sequences obtained in the study); such reads were termed sjUMR (algorithmic definition provided in Methods).

We found that 97% of 25-bp halves of all sjUMR mapped within the 123,066 UMR clusters defined above, providing a first validation of the approach (examples are shown in Fig. 3B,C). In other words, we confirmed that the overwhelming majority of sjUMR connected bona fide exons and not random positions. Next, sjUMR were grouped into clusters, termed SJ [read] clusters (SJC), each of which consisted of all the sjUMR spanning the same splice connection. From this analysis we identified 99,291 SJC connecting two of the previously defined exon clusters. To validate the predicted exon–exon splice connections, we determined if exons connected by SJC belonged to the same gene. We found that SJC that were located on one side in an exon cluster of an annotated gene connected to another exon cluster of the same gene in 94% of instances. In 3% of cases, they connected to a previously nonannotated exon cluster, and in another 3%, they connected to an exon cluster belonging to an apparently unrelated gene. While this indicates an error rate of 3%–6% “wrong” connections, at least 94% of the connections described by SJC were confirmed by the annotation as true connections, and up to 97% may, indeed, be correct considering that a frequent type of newly discovered exons are those belonging to (and thereby extending) a previously discovered gene.

The majority of spliced introns in mammals (as well as most lower eukaryotes) are flanked by a 5′-GT…AG-3′ consensus motif, i.e., 5′-GT at the beginning of the intron (directly following the last base of the splice donor exon) and AG-3′ at the end (directly preceding the first base of the splice acceptor exon). As an additional means of external validation of the agnostic SJ discovery method, we examined the frequency of the consensus motif. Among 65,108 unselected SJC consisting of 10 or more sjUMR, 97% flanked the 5′-GT…AG-3′ motif; this is virtually indistinguishable from the frequency observed in annotated known protein-coding genes (also 97%), further emphasizing the high degree of validity of our agnostic splice site discovery algorithm.

SJC can consist of sjUMR with different offsets (i.e., different start and end positions) appearing as “jitter” of sjUMR in a graphical representation. While the observation of a jitter may at first appear counterintuitive, it has two discernible causes, which were taken into account for the final splice site prediction, as follows: (1) Mismatches at the end of the sjUMR r₁,…,r₂₅. We allowed for up to two mismatches when aligning reads to the reference genome; these mismatches were permitted to occur at the end of r₁,…,r₂₅ or at the beginning of r₂₆,…,r₅₀. In other words, an sjUMR could be a perfect match if cut into two unequal fragments r₁,…,r₂₃ and r₂₄,…,r₅₀ (or r₁,…,r₂₄ and r₂₅,…,r₅₀; r₁,…,r₂₆ and r₂₇,…,r₅₀; r₁,…,r₂₇ and r₂₈,…,r₅₀); i.e., it is in fact an NMR aligning in a bona fide fashion to the identified SJ; yet, given the 25 bp + 25 bp cut method used, it aligns with mismatches. We resolved this by simply trimming mismatched bases of such sjUMR, which resulted in unambiguous SJs, which, in most cases were supported by additional reads with the correct offset (not requiring trimming). (2) A microhomology may exist between the 5′-end of an intron and the 5′-end of the subsequent exon. In this case, the position of an SJ based on any form of mRNA sequencing is ambiguous, albeit this has no consequence for the predicted mRNA, which is the same regardless of the SJ used. The 5′-GT…AG-3′ motif flanking 97% of introns could be used to resolve the overwhelming majority of such cases with single-base-pair precision.

Stable patterns of pre-mRNA splicing in chronic pain

Alteration in pre-mRNA splicing is a mechanism of cellular regulation. We considered the possibility that splicing might be altered in the DRG of animals with chronic pain as compared to controls. To test this possibility, we compared the usage of alternative splice sites with the same statistical approach used above for comparison of gene expression levels. Based on a false discovery rate analysis of the 573 best-supported alternative splice events, >99% of alternative splicing choices were changed <2.79-fold. When SNL and control samples were compared, only four to seven alternative splice events fell outside the 99% control range compared with five to six (1% of 573) expected false-positive events. Reads crossing the splice junctions at other locations (22 + 28, 23 + 27, 24 + 26, 26 + 24, 27 + 23, 28 + 22) were then added to the analysis to improve the degree of support by using all data. As a result the 99% confidence interval was tightened (−1.37 to 1.37 on a log₂ scale), but the outcome remained principally unaltered—four to six alternative splice events now fell outside the 99% control range compared with an unchanged expected false-positive number of events (five to six). Thus, no differences were found to suggest that alterations in splicing were a common mechanism of gene regulation in pain.

De novo transcriptome annotation: Integration of all mRNA-seq read classes and agnostic SJ discovery

Applying the newly defined SJC to annotating previously nonannotated portions of the rat genome, we found that 3420 of the 10,464 nonannotated exon clusters (newly discovered as described above) were connected to annotated genes by SJC. We found 421 groups of exons consisting solely of newly discovered nonannotated exon clusters connected by SJC, i.e., transcriptional units resembling new genes. Of these, 90 were homologous to nonannotated regions of the mouse genome and 320 to known mouse genes. No homology was found for the remaining 11 (BLASTN E-value ≤ 0.01). Comparison with the human genome led to a smaller fraction of known homologs, 269, which was to be expected because of the greater phylogenetic distance. Integrating information from UMR, sjUMR, MMR, and NMR read classes provided complex genome annotation/reannotation including discovery of large novel genes. For example, we mapped a 64-exon-long homolog of UNC-80, which was previously not annotated in the rat genome. While our study was under way, a mouse form of UNC-80 was prominently reported as a novel substance-P and neurotensin coreceptor (Lu et al. 2009). Substance P has long been recognized as a critical peptide neurotransmitter at the spinal level in chronic pain. Down-regulation of the rat UNC-80 homolog after SNL may contribute to allodynia through alterations in peptide activity in the DRG. The case of UNC-80 serves as an example of how de novo transcriptome annotation may be highly informative in identification of novel candidate regulators that would not have been identifiable by previous technologies such as microarrays or qPCR.

Induction and suppression of 10,367 known protein-coding genes

For a given sample of mRNA-seq reads S (e.g., all reads obtained for a single DRG), the read count for an individual gene c(g, s) is the sum of read counts for the set E of all exons e₁…e_n annotated by Ensembl (RGSC 3.4, release 50) as belonging to the gene g, i.e.,

A read was counted as belonging to an exon if its start position fell within annotated exon boundaries. “Known protein-coding” (kpc) Ensembl genes were considered in this part of the analysis, requiring a c(g, s) ≥ 100 for at least two samples at the early time point (2 wk), thereby excluding low-expressing genes that contribute only noise. As a result 10,367 genes were included in the analysis. c(g, s) were normalized to the total number of UMR in each respective sample s:

i.e., the count of reads for a gene per 1 million reads in the respective sample. A count of 1 was added to c(g, s) as a discounting step precluding the possibility of division by 0 in subsequent operations.

To validate mRNA-seq for the quantification of gene expression changes, a comparison with quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) was performed. Fold-changes of gene expression measured by Taqman qRT-PCR were highly correlated (r = 0.956, slope 0.984) with mRNA-seq, as shown in Supplemental Figure 3, indicating the validity of mRNA-seq as a method for quantifying relative changes of transcript levels.

To estimate the number of induced and suppressed genes in SNL animals compared to controls, the FDR concept was applied. The relative change r of the normalized expression of a given gene n(g, s) between two samples r(g, s₁, s₂) was calculated as:

i.e., the fold-change on a log₂ scale, where +1 corresponds to doubling of gene expression and −1 to a reduction by half.

Next, the ratio r(g, s₁, s₂) of each of the 10,367 genes was determined, with s₁ and s₂ corresponding to two biological replicates, i.e., either two control or two SNL animals. Under ideal conditions, all ratios on the log₂ scale should be 0 (i.e., a change of 1×, which is “no change”). Examining the actual distribution of these control ratio distributions, we found that the respective medians were very close to 0, their shape was symmetric, and the average 99% confidence interval was ±0.8096, i.e., only 0.5% of control expression ratios r(g, s₁, s₂) were higher or lower than the confidence interval cut-off, corresponding to a FDR of 52 genes induced and suppressed among the total of 10,367. Next, we determined the number of genes induced and suppressed after SNL by calculating r(g, s₁, s₂) with s₁ and s₂ corresponding to SNL and control samples, repeating the procedure for independent biological replicates. Results are shown in Figure 1E.

Validation of mRNA-seq with qPCR

mRNA-seq was performed using two quality-control standards, the universal human reference RNA (Stratagene) and human brain reference RNA (Ambion). Fold-change differences among 755 transcripts quantified by mRNA-seq were calculated as described above and compared with established Taqman qRT-PCR data (Shi et al. 2006) as shown in Supplemental Figure 3.

Pathway analysis

Data were analyzed through the use of Ingenuity Pathways Analysis (Ingenuity Systems; http://www.ingenuity.com). A network is a graphical representation of the molecular relationships between genes and gene products. Genes or gene products are represented as nodes, and the biological relationship between two nodes is represented as an edge (line). All edges are supported by at least one reference from the literature, from a textbook, or from canonical information stored in the Ingenuity Pathways Knowledge Base. Human, mouse, and rat orthologs of a gene are stored as separate objects in the Ingenuity Pathways Knowledge Base, but are represented as a single node in the network. The intensity of the node color (Supplemental Fig. 2) indicates the degree of up- (red) or down- (green) regulation. Nodes are displayed using various shapes that represent the functional class of the gene product. Annotation of genes is provided in Supplemental Table 12.

Novel gene discovery

A read histogram was constructed for each chromosome c by counting the number of reads covering a given base position x, denoted as h_c(x), resulting in a group of chromosome-spanning histograms h_c(·)(c = 1, … ,21) with single base position resolution covering the genome. For this analysis, data from all samples were pooled:

A sliding window with length l was shifted over the genome integrating the area under the curve of the read histogram h_c(x) in the interval [p, p + l], i.e.,

When w_c(p, p + l) is zero, the consecutive interval of length l has no reads. An interval [a, b] on the chromosome c, for which all values w_c(a, a + l), w_c(a + 1, a + 1 + l), w_c(a + 2, a + 2 + l), …, w_c(b − l, b) are non-zero, is called a read cluster. We denote this read cluster with K_a,b,c. Setting the window length l = 100, we defined 882,913 K_a,b,c. In other words, we identified 882,913 read clusters consisting of one or several neighboring reads; clusters were separated by sections devoid of reads (no read within 100 bp).

Because an individual read cluster by this definition could consist of as little as a single read and as many as thousands, we reasoned that only some, i.e., those with many reads, would be strong predictors of exons. Therefore, we introduced read density as a measure to characterize clusters further. For each read cluster K_a,b,c, we calculated the read density as follows:

Next, we developed a classifier dichotomizing clusters into a set of exon-like clusters F and its complement based on read density, whereby

To establish the test characteristics of the procedure, we tested it for D = 0, …, 100 on a training set as described in the Results section, constructing the ROC curve shown in Figure 2C, which suggested D = 4 as a cut-off for optimal precision. Sensitivity and 1-minus-specificity values for each tested read density are provided in Supplemental Table 7.

Support was provided by the National Institute of Neurological Disorders and Stroke (NINDS) of the NIH (R01NS063022 and R21NS062271) (A.S.B.) and by the Schulze Family Foundation (A.S.B.).