Genome sequencing of bacteria: sequencing, de novoassembly and rapid analysis using open source tools

Sequencing output

Genomes were sequenced on two runs using four multiplex identifier (MID) adaptors for each run. Rapid Library (RL) sequencing yielded poor results because of the relatively low read lengths, so the sequencing was repeated using General Library Preparation (GLP). The sequencing outputs are presented in Additional file 1: Table S1. The average read length was shorter using the RL kit than the GLP kit, with fewer total bases and less total coverage. Pooled data with ~ 30 fold coverage from both sequencing runs were used in all downstream data analysis.

Reference mapping of re-sequencing and annotation/re-annotation

The automated annotation pipelines RAST, IMG, and PGAAP were evaluated with reference-mapped contigs found using the MIRA3 and Mosaik assemblers. These annotations were compared to the original annotation of the A. bokrumensis SK2 (NC_008260 genome from GenBank [gi:19683]); in addition, the raw sequence of SK2 was re-annotated with these three annotation pipelines (Additional file 1: Table S2). Re-annotation of the finished genome sequence increased the number of Coding Sequences (CDSs) by about 5–6%, and the number of tRNA and rRNA operons did not change. The length of the CDSs did not change remarkably and the average length of CDSs in the finished genome was 993 ± 666 (mean ± standard deviation, the largest CDS being 10,803 bp), and in the re-annotated genome 934 ± 677 (the largest being 10,677 bp). In contrast, when using truly re-sequenced data, both the MIRA3 and Mosaik assemblers predicted considerably more CDSs. The CDS size clearly decreased using PGAAP annotation, e.g. 11 contigs from MIRA3 mapping resulted in an average CDS of 490 ± 337 (the largest comprised of only 2,925 bp, which is statistically different from the finished genome at p < 0.01 by Tukey pairwise comparison of one-way ANOVA). This occurred because in the reference mapping mode 100% of the consensus contigs contained many artificial indels, and many more artificial stop codons were introduced (data not shown). These artificial indels originated mostly from homopolymeric regions or other sequencing errors of different pyrosequence reads from the same fragment. Different independent reads contained different numbers of homopolymeric nucleotides of high quality, which led to a consensus of less than 100%, because re-mapping algorithms are optimized to discover SNPs. For example, one aconitase A CDS was fragmented into five pieces using both the IMG and RAST annotation pipelines (Additional file 1: Figure S1). PGAAP behaved even worse as it tended to annotate only the longest CDS with the same name/function, so the remaining sequence regions were left unannotated, and therefore the proportion of non-coding parts increased substantially.

A full comparison of the matching annotation results (CDSs) obtained from the genome assemblers and annotation pipelines used is presented in Figure 1a; CDSs were compared on the basis of sequence match at 100% similarity – to compare CDSs with exact start and stop matches. All pipelines, including the re-annotation pipeline, predicted more CDSs than the original sequencing (NC_008260 in GenBank). Some of these CDSs were unique to each pipeline and were not identified by the others. The number of matches identified increased by 52.8% when re-sequenced and de novo assembled data was used compared with the re-annotated sequence data from the finished genome (Figure 1b, 769 and 1507 CDS). This was corroborated by comparing raw sequences; While both the reference mapping and the de novo assembly contained sequencing errors (Figure 2), these errors influenced the consensus of reference-mapped reads more severely.

Further manual comparison revealed differences in gene calling that could affect the accuracy and annotation details of genome analyses. The tRNA and rRNA operons were used as examples (Additional file 1: Table S2). Annotation of MIRA3 mapping by RAST was in good agreement with the NC_008260 re-annotation; all three rRNA operons were annotated accurately, and the number of tRNAs was 43 (originally 42). IMG-ER annotation of the same mapping was very similar to RAST: three rRNAs and 43 tRNAs. In contrast, PGAAP failed to find the rRNA operons and only three 5S rRNAs were annotated; the rest of the operon was found in the assembly but was not recognized by PGAAP. Annotation results for tRNAs were similar to RAST and IMG-ER. The greater number of contigs (18) that resulted from Mosaik assembly induced suboptimal results, e.g. RAST failed to annotate any rRNA operons and predicted one fewer tRNA. An exact sequence-based comparison of all three rRNA operons in the re-annotation of NC_008260 and re-sequencing/annotation of our A. borumensis SK2 are provided in Additional file 1: Figure S1.

De novoassembly

Three assemblers (MIRA3, Newbler, and CABOG) were compared on the basis of output statistics of the de novo assembly using their default parameters in the assembly process (Additional file 1: Table S3). All de novo assemblies of the four strains were highly fragmented into relatively small contigs when assembled using CABOG. They contained more contigs in total, and had lower N50 values when compared with the other assemblies. MIRA3 assemblies resulted in longer contigs and higher N50 values when compared with the other assemblies (Additional file 1: Table S3). De novo assembly of the re-sequenced genome of A. borumensis SK2 using MIRA3 covered 98.7% (3,079,251 of 3,120,143 bps) of the reference genome, which is lower than the reference mapping of the same number of reads. In addition, the pairwise similarity of the reference genome and de novo assembly was 95.3% (3,010,624 identical sites). In contrast, de novo assembly using CABOG covered only 81.6% of the total genome (2,545,645 of 3,120,143 bps). Newbler covered slightly fewer base pairs than MIRA de novo assembly - 98.1% (3,061,666 of 3,120,143 bps).

The genome sizes of novel strains subjected to true de novo assembly varied: Flavobacterium strain GOBB3-209 had ~2.4 Mbps (GC content: 33.5%) and GOBB3-C103-3 ~4.2 Mbps (34.6 CG%), and Marinomonas GOBB3-320 had ~4.8 Mbps (44.6% GC%). The number of assembled contigs for the small genome bacterium (genome ~3 Mbps) was comparable to A. borkumensis SK2, which has roughly the same genome size. Therefore, we assume that our de novo assembly might cover ~99% of the total genome of GOBB3-209. Larger genomes (> 4 Mbps) were more fragmented and thus it is difficult to estimate how much of the total genome was covered and annotated using our approach.

Comparison of de novoassembled genomes annotated by automated pipelines

Assembly algorithms

Two algorithms are widely employed in de novo genome assembly. The first is the overlap layout consensus (OLC) or overlap contig consensus approach, and the other is the de Bruijn graph (DBG) or Eulerian path [17]. The latter is more useful for shorter reads (<150 bp) numbering hundreds millions compared to a few million 454 Titanium pyrosequencing reads (>400 bp). We used our data to conduct preliminary tests on the effectiveness of Velvet, which implements a DBG-based approach [18], however, the assemblies were highly fragmented (data not shown). For a comprehensive review of de novo assembly algorithms we refer the reader to Miller et al. [19].

MIRA3 is a hybrid combination of the OLC and greedy algorithms [20], and is an iterative assembler that learns from past mistakes with OLC and greedy components (B. Chevreux, personal communication). The MIRA3 project is being actively developed and has a growing group of users. MIRA3 performed better in all aspects of bacterial genome assembly using Roche 454 reads and currently seems to have the greatest potential. In addition, it can be combined with high quality reads from Illumina technology into hybrid assembly, which potentially evades bottlenecks in de novo genome assembly of relatively simple bacterial genomes and hopefully will allow this step to be fully automated in the near future.

Newbler is a commercial product, and probably uses the OLC approach, developed by Roche Diagnostics, but it is usually freely available to laboratories running Roche 454 sequencing. It is not a free and open source software package, and release descriptions indicate that the originally published algorithm may differ from the current one. In our hands, the available version of Newbler performed almost as well as MIRA3, however, commercial tools cannot be evaluated in full detail. Celera Assembler (CA) is another open source OLC based tool that evolved from a Sanger-era assembler; the revised pipeline for NGS reads including Roche 454 data is generally called CABOG [19].

It should be mentioned that we did not apply paired end (PE) library sequencing and it is clear that such mate pairs would allow more contigs to be closed into scaffolds or at least a list of ordered and orientated longer contigs. We disregarded this kind of approach in the search for the most simple and cost efficient method for assembling bacterial genomes, and in theory a coverage of 30 times should allow an average bacterial genome to be assembled. The bacterial WGS in the contigs we obtained should not lose the correct prediction of too many genes, and from this point of view proper scaffolding is more important for eukaryotic genomes with larger size.

Reference mapping is a different field of genome assembly related to the exercise of alignment; it should be much easier to map NGS reads to very close reference sequences. We used reference mapping only with our test strain and our initial approach was to use only Mosaik [21] for this purpose. Mosaik is a specific reference-guided assembler using the Smith-Waterman algorithm to align a hashed table of short reads to the reference genome. The purpose of this analysis was to compare the finished genome with re-sequenced 454 reads mapped to the reference and to the de novo assembly; we did not expect surprises in this process. However, Mosaik was initially unable to close all gaps between the contigs (18 contigs) resulting from re-sequenced data, leaving ~ 20,000 bps out (Table 1). Because the preliminary results were somewhat suboptimal, including several regions that were not covered, we tested the reference mapping performances of Newbler and MIRA3 for comparison. Newbler performed better, but as with the Mosaik mapping it left some regions without mapped reads (overall statistics in Table 1). MIRA3 reference mapping revealed that there were still a few weakly-covered regions although the total coverage (30 fold) should be acceptable. Fine-tuning of the MIRA3 reference mapping would allow all gaps to be closed with coverage of very few reads that were not accepted and thus not mapped by the software.

Gene prediction and annotation - tools

In most so-called automated online tools, a truly automatic process with no manual intervention remains a design goal of annotation in terms of predicting tRNA, rRNA and protein coding genes – coding sequences (CDSs). The usual pipelines for finding a gene in a raw DNA sequence involve detection of an ORF, finding the gene and predicting its function by comparing it with genes in existing databases (see below). Automated annotators normally use several HMMs (Hidden Markov Models) and BLAST-based gene prediction methods, e.g. tRNAscan-SE [22], and BLASTp and BLASTn for protein coding and RNA genes, respectively.

Thereafter, the CDSs are assigned to annotations based on various functional resources such as COG clusters [23], Pfam [24], TIGRfam, Gene Ontology etc. Functional annotations may be further “grouped” into metabolically relevant “pathways” such as COG functional categories, Entrez Protein Clusters (ProtClustDB [25], FIGfams-subsystems [10], KEGG [26] and MetaCyc [27] pathway collections, etc. Thereafter, annotated genomes might be maintained by integrated network systems such as RAST-SEED, IMG and others. Such networks allow further comparative analysis to be performed easily with no need for the researcher to have deep expertise in bioinformatic algorithms and tools.

All of the annotation pipelines we tested are second generation tools, which try to combine multiple gene-calling algorithms and knowledge databases for comparison with related species and training sets. Therefore, these tools should perform much better than the early gene-prediction methods such as Glimmer [28] and GeneMark [29].

Draft WGS assembly and preliminary annotation using pipelines proved able to describe the core genomes and different metabolic features of novel environmental miroorganisms fairly well: genome size, basic metabolic pathways, the number of tRNAs, but not rRNAs. For example, two Flavobacterium strains isolated from the same environment and at the same season showed large differences in genome size; strain GOBB3-209 is 45% smaller than GOBB3-C103-3. Therefore, it is not surprising that strain GOBB3-209 lacks several features found in GOBB3-C103-3, e.g. the capsule and extracellular polysaccharide pathway, the denitrification pathway, etc. (data not shown). Re-annotation of these finished genomes after a re-sequencing project (e.g. A. borkumensis SK2 and Roche 454 pyrosequencing in our study) using automated pipelines may be less useful; in our hands, no improvement was observed after re-annotation. For example, even the prediction of rRNA operons was suboptimal; only one copy in all de novo sequenced and assembled genomes was located, although most bacteria described so far have more than one rRNA operon [30].

However, it is worrying that most genomes currently deposited in the public database rely on automated methods, because it has been reported that their performance seems not to have improved over several years [16]. The deposited data sets may contain genomes and gene annotations that differ in their degree of precision and resolution owing to the use of different sequencing methods and annotations. The most severe problem is that erroneous and incomplete annotations are often carried over into the public resources and are difficult to trace and correct afterwards. For example, several hundred CDSs might be removed or added and the start sites corrected (e.g. re-annotation of the uropathogenic Escherichia coli strain CFT073 [31]), totalling more than 1000 changes when such data are re-evaluated carefully. Even when many biochemical, physiological, and genetic data support broad genome similarity, the current automated annotation tools can fail to predict certain metabolic pathways. Therefore, more detailed studies that combine all types of available data are needed [32].

Bottlenecks to be considered

The computational efficiency of de novo assembly algorithms implemented in free and open source software (e.g. MIRA3) no longer seems to be a bottleneck. At least with longer reads with ~30-fold coverage, a reasonable draft genome can be produced within hours and without manual intervention on direct shotgun sequencing.

On the other hand, the tools supplied by bioinformatic service providers such as RAST, PGAAP, and IMG cannot yet be fully automated and involve manual intervention. It is clear that the performance of these tools will improve as more carefully curated and finished genomes become available to aid in the automation process. However, considering the number of genomes from unique species available today (in the range of several thousand [33]), and because the potential abundance of microbes in nature is huge [34], progress in this area will take time. Nevertheless, the bacterial genomes used in our study are relatively well covered by existing knowledge of phylogenies and databases of fully sequenced and finished genomes. In the GOLD database [35] there are closely-related bacteria from the genus Marinomonas: M. mediterranea MMB-1, ATCC 700492 (unpublished), M. posidonica IVIA-Po-181 (unpublished) Marinomonas sp. MWYL1 (unpublished) and Alcanivorax borkumensis SK (finished) [36]; and from the genus Flavobacterium: F. columnare ATCC 49512 (unpublished), F. branchiophilum FL-15[37], F. johnsoniae UW101 [38], and F. psychrophilum JIP02/86 [39]. Therefore, after true finishing of the drafts, the genomes of these organisms can be described with no particular effort, after manual curation, at the same level as phylogenetic relatives of finished genomes. Further justification and re-annotation would be based directly on new biological discoveries.

The usefulness of detecting SNPs has been discussed previously [40], and because the error rate of 454 pyrosequencing is higher than Sanger sequencing and probably more than Illumina (Solexa/Genome Analyser), it might be a suboptimal choice for re-sequencing and reference mapping projects. Pop & Salzberg [41] reasoned that fragmented draft genomes would produce fewer annotated genes because of false stop codons; our observations differ but lead to a similar result with regards to reference mapping. The flaws in reference mapping caused by errors introduced by 454 pyrosequencing led to frame shifts and thus to a greater number of CDSs.

We did not use paired-end (PE) libraries because we wanted to keep the preparation cost of sequencing as low as possible. Our results indicate that the PE approach is useful for scaffolding truly de novo-sequenced data to a sufficient degree of coverage (~30 fold) even in assembly of relatively small and simple bacterial genomes. Complex genomes containing repetitive elements may need more attention, although the same is true for annotation of genomes containing genomic and/or pathogenic islands. Isolation of genomic DNA by standard methods often fragments the chromosomal DNA into smaller pieces with a size limit of a few tens of kbs, so extended PE libraries with very large fragments are not easy to construct and need specific treatment of the genomic DNA prior to sequencing. Therefore, it is not possible to assemble larger repetitive elements correctly compared to the fragments in DNA extracts.

Origins of the strains

Three environmental bacterial strains, Flavobacterium sp. GOBB3-C103-3, Flavobacterium sp. GOBB3-209 and Marinomonas sp. GOBB3-320, were obtained from estuarine enrichment cultures based on samples from a watershed in the northern Baltic Sea [42, 43]. The GenBank accession numbers of the 16S rRNA gene sequences were AF321019, AF321038, and AF321017, respectively. One bacterium, type strain Alcanivorax borkumensis (strain SK2/ATCC 700651/DSM 11573, TaxID: 393595) [NEWT/NCBI]), was chosen as a reference strain because its genome was fully sequenced and annotated in the pre-NGS era (ref PubMed = 16878126).

Growth conditions of strains

Strains were isolated from enrichment cultures when the total bacterial abundance reached ~1–2 × 10⁶ cells ml^-1 (early stationary phase, after approximately two weeks from inoculation). Cultivable bacteria were isolated by spreading 100 μl samples on triplicate ZoBell medium or more nutrient-poor AC medium agar plates with 60 mmol l^-1 riverine dissolved organic carbon. ZoBell medium [40] was made from 800 ml, 0.2 μm filtered (Durapore, Millipore) seawater and 200 ml of MilliQ ultrapure water (Millipore). This was supplemented with 5 g casein pancreatic digest (Trypton, Difco), 1 g yeast extract (Difco) and 0.01 g FePO₄.2H₂O (Aldrich). Plates were incubated in the dark at 15°C until no more colonies appeared (about 5–15 days). Isolates were purified by serially streaking a single cell colony three times on to new agar plates. The purified isolate was grown in ZoBell medium (10 ml) in the dark at 15°C (1–2 days). After the purity of the isolate was confirmed by one more plating it was re-grown in liquid Zobell medium and 0.8 ml of this cell suspension was mixed with 50% glycerol (0.2 ml) and stored at -80°C. Type strain Alcanivorax borkumensis SK2 was routinely maintained on Marine Broth 2216 at 30°C, for plating with agar (15 g l^-1).

Genomic DNA preparation and quantification, 16S rRNA gene amplification, and sequencing

Genomic DNA was extracted from pelleted pure isolates using a Qiagen DNeasy Blood and Tissue kit (cat. 69504). The amount of genomic DNA was measured using a Qubit Fluorometer (Invitrogen). To check the isolate identity, a partial 16S rDNA gene from the isolated bacteria was amplified and directly sequenced using the Applied Biosystem 3730XL according to the manufacturer’s directions. The bacterial universal primers 27F (3’-AGAGTTTGATCATGGCTCAG-5’) and 1492R (3’-TACGGYTACCTTGTTACGACTT-5’) were used for amplification. The PCR product was purified with PCR Kleen Spin Columns (Bio-Rad Inc) and the nucleotide sequences determined from three partial fragments of the 16S rRNA gene covering almost the full length of the gene. Nucleotide sequences were determined from the purified 16S rRNA gene with the primers 27F (3’-AGAGTTTGATCATGGCTCAG-5’), 800R (3’-CCAGGGTATCTAATCC-5’), 1492R (3’-ACGGGCGGTGTGTRC-5’), and 347F (3’-TACGGGAGGCAGCAG-5’).

Whole genome sequencing

Purified DNA from each isolate was directly sequenced using Roche standard technology (Roche 454 GS-FLX system, GS Titanium chemistry by Zürich Functional Genomics Centre). The sequencing libraries were prepared with a GS Rapid Library Kit (Roche, cat. no. 05 608 228 001) together with a GS Rapid Library MID Adaptors Kit (Roche, 05 619 211 001) using 500 ng of DNA. In addition, a General Library Preparation Kit (Roche, 05 233 747 001) together with a Titanium Library MID Adaptors kit (Microsynth) was employed using 3.7 μg DNA from Alcanivorax borkumensis SK2, 4.5 μg from Flavobacterium sp. GOBB3-C103-3 and 5 μg from both Marinomonas sp. GOBB3-320 and Flavobacterium GOBB3-209, according to the manufacturer's protocol. The sequencing reactions were performed using a Roche 454 Genome Sequencer FLX with the GS Titanium Sequencing Kit XLR70 (Roche, cat. no. 05 233 526 001) using the 2-region gasket of the GS Titanium PicoTiterPlate Kit (70x75) (Roche, cat. no. 05 233 682 001), according to the manufacturer's instructions. Imaging and signal processing were done using GS FLX SW v2.3, gsRunProcessor fullProcessing.

Analysis of sequencing data – reference and de novoassembly and rapid annotations

For reference mapping, Mosaik release 1.0.1388 [18], Newbler Version 2.3 (Roche) and MIRA3 version 3.2.0 [44] were used; raw reads from the Roche sequencing of Alcanivorax borkumensis ATCC 700651 DNA isolated in the laboratory and the backbone of Alcanivorax borkumensis SK2 (NC_008260 [gi:19683]) were downloaded from NCBI. De novo assembly was performed using the Roche reads obtained from all four strains; the assemblers were MIRA3, Newbler v. 2.3 (Roche), and modified Celera Assembler (CABOG Version 6.1, [45]). Multiple contigs obtained from the reference mapping and de novo assembly were subjected to three automated annotation tools: RAST [46], NCBI-PGAAP [8], and IMG-ER [7].