Student's t-test, which takes into account both types of variation, to identify genes differentially expressed between the two groups. In doing this, however, we found some important flaws of classic statistics: a t-test can only be applied in a meaningful sense after normalization for the total number of tags in the library and proper variance stabilization. We did this by linear scaling based on the total number of counts and subsequent square root transformation. The square root transformation (approximately) stabilizes the variance of raw counts but not of the scaled data. Hence, we cannot stabilize variance while normalizing for library size at the same time (22). This problem is particularly prominent
in our experiment where one wild-typeand one transgenic sample had, respectively, 3 and 10 times lower numbers of counts than the other samples.
Vencio et al.(10) proposed a Bayesian method for the analysis of
replicated SAGE data that takes into account stochastic effects for the low abundant genes as well as differences in library size. It reports the Bayesian error rate which can be interpreted as the chance that a gene
is found differentially expressed underthe null hypothesis. With a Bayesian error rate of 0.05, we detected 1559 up- and 1620 downregulated canonical tags in the comparison of transgenic with wild-type mice. The distribution of the detected fold-changes can be inferred from the Volcano plot in Figure 2 and ranged between 71-fold (2700089E24Rik, found only once in all wild-type samples, but 19 times in transgenic samples) and 1.13-fold. Differentially expressed tags were found in the entire range of expression levels (Supplementary Figure 5). A list of the 20 most significant tags is given in Table 2. Vencio's test does not consider multiple testing. To estimate the number of false positives obtained, we calculated Bayesian error rates when permuting the samples (Supplementary Figure 6). The number of tags found differentially expressed with an error rate of 0.05 in the permutated situations was 270 ± 103. Thus, the false discovery rate in our list of 3179 differentially expressed genes is estimated to be 8.5%.
View this Table 2. List of the 20 most significantly table: differentially expressed tags [in this window] [in a new window]
Figure 2. Volcano plot of canonical tags. For every tag, the ratio in expression levels of transgenic over wild-type mice (2log scale, x-axis) is plotted against the Bayesian error rate (10log scale, y-axis). The horizontal line indicates
View larger version (10K): the significance threshold applied, the
[in this window] 3179 differentially expressed tags being [in a new window] above that line. The plot shows that the [Download PowerPoint tags with highest average differences
slide] between trasngenic and wild-type mice (far left and right part of the plot) are
not all significant (due to large intragroup variation). The most significant tags (top of the plot) generally display small differences in expression between transgenic and
wild-type but are, due to relatively high expression levels, very accurately measured and therefore display low intragroup variation.
In addition to differentially expressed canonical tags, we detected
differential expression of 2479 noncanonical and 15 mitochondrial tags. Biological implications
By alternative splicing the DCLK gene produces numerous proteins. Recent functional studies from a.o. knock-out mice strongly indicate involvement of the DCLK gene in several molecular pathways. Some are microtubule-associated proteins (23) that may regulate
microtubule-guided transport of SNARE-protein containing synaptic vesicles (24), while the DCLK-short variant has Ca++/calmodulin-dependent protein kinase (CaMK) properties (8,25). In the current study, we
evaluated which biological pathways were affected in the hippocampus by the expression of the DCLK-short isoform. The global test (11) was applied to the DGE data to identify the differential regulation of gene sets, as defined by the Gene Ontology consortium. Unlike commonly used
overrepresentation tests or gene set enrichment analysis, this method uses the gene expression measurements of a particular set of genes, giving optimal power for small sample-size experiments and detection of gene sets where many genes display a small effect. The most significantly affected
pathways are reported in Table 3. Strikingly, the CaMK pathway was the second most significant pathway. Disturbances in the expression of genes in the CaMK pathway are possibly a consequence of transcriptional feedback
mechanisms. Also inline with the function of the DCLK gene, we find indications for disturbed synaptic vesicle transport along microtubules as a result of alterations in gene expression of vesicle SNARE proteins (rank 19) and microtubule plus-end binding proteins (rank 1), potentially affecting neurotransmitter release and axonal outgrowth.
View this table: [in this window] [in a new window] Table 3. Significantly deregulated pathways in DCLK transgenic mice
The effect of sequencing depth on detection of differentially expressed genes
Before the development of deep sequencing technology, constructionof a large-scale SAGE library containing up to 100 000 canonical tags would typically take 1 year and a considerable financial investment. The number of tags in such a library is 60 times smaller than the number of tags we obtained for each group of samples in a 3-day experiment. To illustrate the effect of the increased sequencing depth, we have compared our results
tothe results from a simulated SAGE experiment, which includesonly 1/60 of the original number of DGE reads, randomly taken. The number of detected
differentially expressed genes decreased15-fold, from 3179 with the original number of reads to 200 in the simulated SAGE experiment (Bayesian error rate <0.05). The lowest abundance of a significantly detected differentially expressed transcript was 0.8 t.p.m. in our deep sequencing experiment versus 91 t.p.m. in the simulated SAGE experiment. As noted before (26), many of the genes with most significant changes in expression
are low-abundant genes and would not have beenidentified in a typical SAGE experiment.
Comparison with microarrays and qPCR
The exact same RNA samples had been analyzed previously by five different whole-genome expression microarray platforms: Applied Biosystems, Affymetrix, Agilent, Illumina and home-spotted oligonucleotide arrays
(9). We compared the results from DGE and the microarray experiments after mapping all canonical tags and microarray probes to the ENSEMBL transcript
database. With DGE, we detected15 189 ENSEMBL transcripts with
abundances >2 t.p.m. With most microarray platforms, a lower number of transcripts gave signal above background, except for Agilent, where there may have been considerable background signal caused by
cross-hybridization (Table 4). Affymetrix had the highest percentage of transcripts in common with DGE. In general, less abundant transcripts were more difficult to detect with microarrays. The median expression of 538 transcripts detected by DGE but not by any of the microarray platforms had a median abundance of only 4 t.p.m., while the transcripts that were detected by all platforms had a median abundance of 106 t.p.m.
View this table: [in this window] [in a new window] Table 4. Overlap between DGE and microarrays in detectable transcripts
Figure 3 shows the correlation between absolute transcript abundance and microarray probe intensity. In line with other reports (27–29), we observed a reasonable correlation between the intensity of the microarray hybridization signal and the number of tags sequenced. The correlation was highest for Affymetrix chips (Pearson correlation: 0.63). For the Affymetrix data, intensities of the 11 perfectly matched probes were summarized into a single value. Indeed, the use of 11 different probes per transcript, in contrast to a single probe per transcript for the other platforms, should even out probe-specific hybridization characteristics. The correlation in detected transcripts was higher than previously found for SAGE or MPSS versus Affymetrix (30,31), mainly due to the higher number of tags sequenced with DGE.
Figure 3. Correlation between absolute expression level (DGE) and microarrays signal intensity. Correlation of the tag
View larger version (26K): abundance (square root transformed; [in this window] x-axis) and intensities [normalized as [in a new window] described in (9)] on the five microarray [Download PowerPoint platforms (y-axis) for matching ENSEMBL slide] transcripts, for wild-type sample 1. Pearson correlations are indicated in the graphs. ABI: Applied Biosystems; AFF: Affymetrix; ILL: Illumina; AGL: Agilent; LGTC: home-spotted long oligonucleotide arrays.
Technical replicate measurements were used to compare the precision of DGE with that of microarrays. As a measure for precision we determined the distribution of the differences between independent replicate
measurements of log ratios between wild-type and transgenic samples, as proposed by Irizarry et al. (1). Figure 4A shows the distribution of these differences for DGE and the two microarray platforms with highest and lowest precision (Agilent and home-spotted oligonucleotide arrays, respectively). The distribution of DGE is narrower (interquartile range (IQR): 0.51) than that of Agilent (IQR: 0.61) and home-spotted arrays (IQR: 0.75), indicating that DGE has a higher precision than microarrays.
Figure 4. Assessment of precision and accuracy of DGE. (A) Samples from the wild-type and transgenic pools were sequenced in three different lanes. We calculated the three possible independent log ratios between
transgenic and wild-type samples (technical replicates). As a measure of precision, we determined the pair-wise differences between these technical replicates. The distribution of these differences is plotted as a density
View larger version function (black line). This is also done for
(17K): three technical replicates of wild-type over [in this window] transgenic ratios determined on Agilent (red) [in a new window] and home-spotted (blue) microarrays. We [Download PowerPoint balanced the number of observations per
slide] platform through random selection of 21 886 features. (B) As a measure of accuracy, we
correlated logged ratios of the expression in transgenic versus wild-type mice as obtained