The tag length used in 5′SAGE technology (19 bp) renders it difficult to identify precisely the genome position of some of the sequence tags. We therefore developed a method to improve 5′SAGE by the generation of 27 mer sequence tags combined with massively parallel sequencing. Figure 1a shows the process used for 5′-end library construction for SOLiD sequencing. cDNA was synthesized after ligating an RNA linker including an EcoP15I site, which cleaves 27 bp downstream of the recognition site leaving a 5′ overhang of two bases,[15] to the 5′ end of mRNAs using the oligo-cap method[16]. The full-length, double-strand cDNA with a biotinylated 5′end was then digested with EcoP15I. A 27 bp fragment from the 5′end of the cDNA was then purified using avidin-bound magnetic beads. DNA fragments were ligated with sequencing adaptors and amplified by PCR, followed by sequencing with the SOLiD system.

We conducted 5′-end sequencing of mRNAs isolated from human colon adenocarcinoma HT-29 cells treated the epigenetic agent, 5Aza, a potent inhibitor of DNA methylation. Using the SOLiD platform, we generated approximately 62 million 25 bp reads (Table 1). Erroneous sequence tags were eliminated using quality values before alignment. The quality value (QV) is a well-known indicator for evaluation on error probability of a base. It is assigned by the SOLiD base (color) caller and estimates the probability each base is called correctly. The probability of a base being wrong p is given by p = 1/10^QV/10 [17]. We removed any read which had 5 or more of the first 22 bases with a QV <9. Although massively parallel sequencing technologies achieve dramatically higher throughputs than capillary sequencing, the ability of the short read lengths to map uniquely to a reference sequence is more sensitive to base calling errors. There is the flexibility to choose a stricter or looser QV criterion depending on how much the user wants to reduce the risk mentioned above. Filtering low quality reads yielded ∼35.9 million reads, and ∼23.0 million (64%) of these were aligned to the genome with at most two mismatches in order to cope with sequencing errors and SNPs. Among mapped reads, ∼16.4 million (72%) were uniquely anchored on the genome. The ratios of unique reads aligned with zero, one, and two mismatches were 36%, 32%, and 32%, respectively. 73% of unique reads were located within 500 bases of the public representative transcription start sites (TSSs) of well-annotated protein-coding genes in the RefSeq database. Among the remaining 27% of unique reads, 11% were localized to RefSeq gene-coding regions within the genome, while 16% may represent small RNAs and annotated genes. In this study, we further analyzed the RefSeq associated reads. The 30 most highly expressed transcripts identified here in HT-29 cells are shown in Table S1. The most highly expressed genes in the untreated (control) cells encoded ribosomal proteins and ornithine decarboxylase antizyme 1. The frequency of these transcripts in the library prepared from 5 Aza treated-cells was similar to the level in the library prepared from the control cells.

In order to assess the reproducibility of the system, we compared the unique reads mapping to RefSeq genes isolated from two independent sequence runs from the same library. These data showed that data generated from technical replicates of the 5′SOLiD library were highly reproducible (Pearson's correlation coefficient r>0.99) (Figure S1a, b). A high degree of reproducibility of the sequence data generated from independent analyses of the same library is essential for accurate detection of subtle changes in gene expression.

We next calculated the number of genes expressed in a single cell type. As human cells contain approximately 300,000 mRNA molecules[1], we evaluated RefSeq associated genes among the 5′SOLiD tags. The level of expression of these genes ranged from ∼0.02 to ∼4,704 transcripts per cell. The unique transcript tags identified in the control and 5Aza library at more than 0.1 copies/cell corresponded to approximately 14,900 and 14,500 RefSeq associated genes, respectively (Table 2). Only transcripts present at greater than 3∼5 copies/cell can be detected by regular SAGE sequencing or DNA array methods [18]. These RNAs represent 95% of the total expressed mRNAs (Figure 1b). However, genes expressed at less than 5 copies /cell accounted for about 60% of the total expressed genes in this study (Figure 1c). Similar numbers of expressed genes and a similar range of copy number were observed in control and 5Aza treated cells.

To evaluate the dynamic range and sensitivity of 5′SOLiD, we compared the 5′SOLiD library to a 5′SAGE (19 base) library prepared from HT-29 cells treated with 5Aza. The 5′SAGE libraries comprised approximately 70,000 tags (Table S2). The sensitivity of the 5′SAGE method appeared to be equivalent to the original SAGE method. Figure 2 shows scatter plots comparing RNA levels determined using the 5′SOLiD and 5′SAGE systems. The X and Y coordinates indicate the number of 5′-end tag-associated genes identified within the control and 5Aza treated cells. A simple comparison of the scatter plots of the data generated from these two different methods, indicates that the dynamic range of the number of unique transcripts identified by 5′SOLiD per single run is 100–1,000 fold greater than that observed for 5′SAGE (Figure 2, Table 3). However, the number of 5′SOLiD tags identified per library sample is one hundred-fold greater than the number of 5′SAGE tags detected per library sample. Therefore to subtract out the differential sensitivity caused by the significantly greater number of SOLiD tag reads relative to the number of capillary tag reads, we selected a random one percent of the 5′SOLiD tag reads so that the selected 5′SOLiD tags were roughly equivalent in number to the 5′SAGE tags. Figure 2b displays the scatter plot of this selected subset of 5′SOLiD tags. The Pearson coefficients of the genes from the selected 5′SOLiD tags were greater than those observed for the 5′SAGE tags, indicating that the 5′SOLiD tags are more likely to fit to the diagonal. In addition, the MA plots generated from each scatter plot are shown in Figures 2d–f. In order to display the wide dynamic range of a method, an MA plot is preferred to a scatter plot because points are depicted along the horizontal line rather than the diagonal line. As indicated, the larger variance relative to 5′SAGE is more easily observed than when using scatter plot analysis for both the 5′SOLiD and 5′SAGE. It shows graphically how the dynamic range of the 5′SOLiD is approximately 100 fold greater than that of the 5′SAGE. More specifically, most of data plotted in Figures 2d (5′SOLiD) and 2f (5′SAGE) fall within the red regions, an index to compare the dynamic range. Observe that the width of this red region in the 5′SOLiD data plot is 2⁽¹⁶⁻⁴⁾ = 1024, while that of the 5′SAGE plot is only 2⁽¹⁰⁻⁴⁾ = 32. It is also remarkable that the distribution of data points within the red regions are less broad when using the 5′SOLiD compared to the distributions of the data points when using the 5′SAGE.

We determined that the expression of 277 genes differed significantly between the control and 5Aza-treated cells, using the criteria of a fold change ≥2 and a Bonferroni adjusted p-value <0.001. The number of these differentially expressed genes corresponds to 1.8% (in 5Aza) of all expressed genes in a single cell type. In contrast, 5′SAGE identified only a few genes exhibiting significant changes in expression levels. Only 41 genes (Bonferroni adjusted p-value <0.001 and greater than 2-fold difference) exhibited a significant difference in expression between control and 5Aza-treated cells.

Table 3 compares the expression of cell cycle-related genes by gene profiling using the 5′SAGE and 5′SOLiD methods. As predicted by the scatter plots of the data generated by the two methods, most cell growth-related genes were present below the limits of detection of 5′SAGE. In contrast, most of these genes were detected with a sensitivity of between 100 to 1,000 tags on the 5′SOLiD platform. For example, whereas the myc proto-oncogene was not detected by 5′SAGE analysis (tag count: cont, 0; 5Aza, 0), it was easily detected using the 5′SOLiD method (cont, 337; 5Aza, 340) (Table 3).

Genes whose expression was altered by greater than 20 fold in response to exposure of cells to 5Aza relative to the control cells are listed in Table S3. A greater number of these genes exhibited an increase rather than a decrease in expression in response to 5Aza.

We next used the Gene Ontology Consortium database to categorize the function of those highly expressed genes whose expression was altered by greater than a 10 fold change in response to 5Aza treatment (Table S4). We observed that genes involved in the immune response, defense response and chemotaxis were differentially expressed in 5Aza-treated cells.

To confirm the results of the 5′SOLiD analysis, we performed quantitative real-time PCR (qPCR) on forty transcripts identified in the previous section. We observed a good correlation between the changes in gene expression in response to 5Aza treatment detected by 5′SOLiD data and qPCR. T and regression analyses between the two methods revealed a high level of consistency between the two sets of data (Pearson's correlation coefficient r = 0.985) (Figure 1d). Figure S2 shows nine representative examples of the correlation between the 5′SOLiD tags and qPCR. It is clear from this analysis that the results of 5′SOLiD method correspond very well to those of qPCR, even for genes expressed at only low copy numbers.

5′SOLiD tags identify the 5′-end of transcripts. Therefore, we analyzed the distribution of 5′-end tags annotated in exons/introns of known transcripts. Although the highest tag density occurs in the canonical position around the 5′-end of well-characterized genes or in the upstream regions of defined TSSs in the databases, some tags were localized to introns, inner exons or proximal to the most 3′ exon (Figure 3a).

To evaluate whether the 5′-end tags within a given gene region are expressed independently, we performed regression analysis of the fold differences between the numbers of tags localized to the promoter+1st exon relative to the number localized to inner exons/introns in the Control and 5Aza -treated libraries. Each linear regression analysis indicated a good fold-change correlation (Figure 3b). This result indicates that the expression of 5′-end tags in the same gene region correlates with those in the other exons and introns. These results suggests that the transcripts within a given gene region are coordinately perturbed by 5Aza.

We next evaluated the number of 5′SOLiD tags that were represented among the known TSS in the DBTSS database, which contains the largest collection of human TSSs, and conversely we evaluated what fraction of the TSS contained with the DBTSS were represented within the 5′SOLiD TSSs data. Figure 3c shows that 80% of the 5′-end tags present within the 5′SOLiD dataset mapped to within 1,000 bp upstream or downstream of the 1,562,911 known TSSs within the DBTSS, an area which comprises only 4.3% of the genome. The rest are potentially novel TSSs revealed by the enormous collection of 5′end tags. Figure 3d shows that only 80% of the known TSSs in the DBTSS are represented within the 5′SOLiD TSSs tags. The remaining 20% of known TSSs in the DBTSS may be derived from many different tissues and samples, which may explain why we did not detect them within the particular cell type used in this study. In addition, the diversity of the transcription start sites was not changed by treatment of cells with 5Aza.

HT-29 cells were cultured in RPMI1640 or McCoy's 5A medium supplemented with 10% fetal bovine serum. Cells were incubated with 5Aza (Sigma) (500 nM) for 72 h; the culture medium was replaced with fresh medium with 5Aza every 24 h. Total RNA was extracted using RNA-Bee prior to analysis by 5′SAGE and RT-PCR. The quality of total RNA was analyzed using an Agilent 2100 Bioanalyzer. Prior to 5′SAGE analysis and 5′SOLiD, we examined the influence of 5Aza on HT-29 cells using bisulfite-PCR to determine the methylation status of CpG islands in the MAGE-1 promoter. We found decreased methylation in the MAGE1 promoter after 5Aza-treatment in HT-29 cells. This preliminary experiment confirmed that epigenetic changes were induced by 5Aza (data not shown).

5′SAGE libraries were created as previously described [4]. Five to ten micrograms of poly(A)+RNA was treated with bacterial alkaline phosphatase (BAP; TaKaRa). The poly(A)+RNA was extracted twice with phenol∶chloroform (1∶1), ethanol precipitated, and then treated with tobacco acid pyrophosphatase (TAP) . Two to four micrograms of the BAP- TAP treated poly(A)+RNA were divided into two aliquots and an RNA linker containing recognition sites for EcoRI/MmeI was ligated using RNA ligase (TaKaRa): one aliquot was ligated to a 5′-oligo 1 (5′-GGA UUU GCU GGU GCA GUA CAA CGA AUU CCG AC -3′) linker, and the other aliquot was ligated to a 5′-oligo 2 (5′-CUG CUC GAA UGC AAG CUU CUG AAU UCC GAC -3′) linker. After removing unligated 5′-oligo, cDNA was synthesized using RNaseH free reverse-transcriptase (Superscript II, Invitrogen) at 12°C for 1 h and 42°C for the next hour, using 10 pmol of dT adapter-primer (5′-GCG GCT GAA GAC GGC CTA TGT GGC CTT TTT TTT TTT TTT TTT-3′).

After first-strand synthesis, RNA was degraded in 15 mM NaOH at 65°C for 1 h. cDNA was amplified in a volume of 100 ul by PCR with 16 pmol of 5′ (5′ [biotin]- GGA TTT GCT GGT GCA GTA CAA –3′) or (5′ [biotin]- CTG CTC GAA TGC AAG CTT CTG-3′ ) and 3′ (5′-GCG GCT GAA GAC GGC CTA TGT-3′) PCR primers. The cDNA was amplified using 10 cycles at 94°C for 1 min, 58°C for 1 min, and 72°C for 2 min. PCR products were digested with the MmeI type IIS restriction endonuclease (University of Gdansk Center for Technology Transfer). The digested 5′-terminal cDNA fragments were bound to streptavidin-coated magnetic beads (Dynal, Oslo, Norway). cDNA fragments that bound to the beads were directly ligated together in a reaction mixture containing T4 DNA ligase in a supplied buffer for 2.5 h at 16°C. The ditags were amplified by PCR using the following primers: 5′ GGA TTT GCT GGT GCA GTA CA 3′ and 5′ CTG CTC GAA TGC AAG CTT CT 3′. The PCR products were analyzed by polyacrylamide gel electrophoresis (PAGE) and digested with EcoRI. The region of the gel containing the ditags was excised and the fragments were self-ligated to produce long concatamers that were then cloned into the EcoRI site of pZero 1.0 (Invitrogen). Colonies were screened with PCR using the M13 forward and reverse primers. PCR products containing inserts of more than 600 bp were sequenced with the Big Dye terminator ver.3 and analyzed using a 3730 ABI automated DNA sequencer (Applied Biosystems, Foster City, CA). All electrophoretograms were reanalyzed by visual inspection to check for ambiguous bases and to correct misreads. In this study, we obtained 19–20 bp tag information.

Five to ten micrograms of poly(A)+RNA was treated with bacterial alkaline phosphatase (BAP;TaKaRa). The poly(A)+RNA was extracted twice with phenol∶chloroform (1∶1), ethanol precipitated, and treated with tobacco acid pyrophosphatase (TAP). Then an RNA linker containing recognition sites for EcoP15I was ligated using RNA ligase (TaKaRa): a 5′-oligo 1 (5′- CUG CCC CGG GUU CCU CAU UCU CU CAG CAG -3′) linker. After removing unligated 5′-oligo, cDNA was synthesized using RNaseH free reverse-transcriptase (Superscript II, Invitrogen) at 12°C for 1 h and 42°C for the next hour. A 5′-end- cDNA library was produced using 10 pmol of dT EcoP adapter-primer (5′- GCG GCT GAA GAC GGC CTA TGT GCA GCA G(T)₁₇ -3′).

After first-strand synthesis, RNA was degraded in 15 mM NaOH at 65°C for 1 h. cDNA was amplified in a volume of 100 ul by PCR using 16 pmol of 5′ (5′ [biotin]- CTG CCC CGG GTT CCT CAT TCT –3′) and 3′ (5′- GCG GCT GAA GAC GGC CTA TGT -3′) PCR primers. The cDNA was produced using 10 cycles at 94°C for 1 min, 58°C for 1 min, and 72°C for 10 min.

PCR products were digested with the EcoP15I type IIS restriction endonuclease (NEB). The digested 5′-terminal cDNA fragments were bound to streptavidin-coated magnetic beads (Dynal, Oslo, Norway). cDNA fragments that bound to the beads were blunted and then directly ligated to the SOLiD linker in a reaction mixture containing T4 DNA ligase in the supplied buffer for 2.5 h at 16°C. The samples were amplified by PCR using the primers: 5′-CCACTACGCCTCCGCTTTCCTCTCTATG-3′ and 5′-CTGCCCCGG GTTCCTCATTCT-3′. The PCR products were purified by polyacrylamide gel electrophoresis (PAGE). The purified libraries were sequenced with the SOLiD system according to the manufacturers' protocol (Applied Biosystems, Foster City, CA). We confirmed the integrity of the cDNA using an Agilent 2100 Bioanalyser prior to construction of the 5′SOLiD libraries.

In the SOLiD Sequencing System, after each cycle of ligation, the color that is detected represents 4 potential two base combinations and the conversion into nucleotide base space is usually done after the sequence is aligned to the reference genome translated into in the colorspace coding. This two base encoding strategy provides higher sequencing accuracy and inherent error checking capability. The details of principles of SOLiD sequencing and advantages of colorspace are summarized at the link below. http://marketing.appliedbiosystems.com/images/Product_Microsites/Solid_Knowledge_MS/pdf/SOLiD_Dibase_Sequencing_and_Color_Space_Analysis.pdf. The algorithm “SOLiD™ System Color Space Mapping Tool (mapreads)” uses spaced seeds for high speed and sensitive global alignment. It transforms the reference to colorspaces then maps the reads. The opposite transformation that changes colorspace to dnaspace (or basespace) should be avoided because a single mismatch in colorspace will affect all the following base calls. In our report, a mismatch means a colorspace differs between a read and the reference. The software can be downloaded from http://solidsoftwaretools.com/gf/project/mapreads/. SOLiD tags were mapped to the human genome (hg17) using the mapping tool developed at ABI. The first 25 color spaces of tags were used for alignment allowing up to two mismatches. “Unique tags”, whose best sequence matches were located on unique positions in the genome, were used in subsequent tag assignment analysis for high accuracy. Unique tags that were aligned with the same start position on the same strand of the genome were treated as identical even if they contained at most two mismatches relative to the genome sequence.

We associated a unique tag with a known transcription start site (TSS) registered in the DBTSS (ver5.2, http://dbtss.hgc.jp/) if the start position of the tag was located within the ±n bp of the TSS in the genome. Similarly, we assigned a unique tag to a RefSeq gene (UCSC, http / hgdownload.cse.ucsc.edu /goldenPath/hg17/database/) if the start position of the tag fell within a range from 500 bp upstream of the TSS to the end of the RefSeq gene. For further analysis, a unique tag was associated with a transcript annotation based upon being mapped within the promoter, an exon, or an intron of a RefSeq gene.

cDNA was prepared from total RNA using the Applied Biosystems cDNA archive Kit and random primers. Multiple reactions containing 10 ug total RNA per 100 ul reaction volume were performed for each sample following the manufacturer's recommendations. TaqMan gene expression assays were performed using a FAM-labeled MGB probe. Each TaqMan assay was performed in triplicate for each RNA sample. 7.5 ng total cDNA(as total input RNA) in a 15 ml final volume was used for each assay. Assays were performed with 2xUniversal Master on Applied Biosystems 7500 (Foster City, CA) using universal cycling conditions (2 min at 50°C, 10 min at 95°C, followed by 40 rounds of 15 s at 95°C and 1 min at 60°C). G3PDH was chosen as the reference gene.

5′SOLiD tags have been deposited at NCBI Short Read Archive under the project accessions SRA002659.