Rapid quantification of sequence repeats to resolve the size, structure and contents of bacterial genomes
© Williams et al.; licensee BioMed Central Ltd. 2013
Received: 26 April 2013
Accepted: 3 August 2013
Published: 8 August 2013
The numerous classes of repeats often impede the assembly of genome sequences from the short reads provided by new sequencing technologies. We demonstrate a simple and rapid means to ascertain the repeat structure and total size of a bacterial or archaeal genome without the need for assembly by directly analyzing the abundances of distinct k-mers among reads.
The sensitivity of this procedure to resolve variation within a bacterial species is demonstrated: genome sizes and repeat structure of five environmental strains of E. coli from short Illumina reads were estimated by this method, and total genome sizes corresponded well with those obtained for the same strains by pulsed-field gel electrophoresis. In addition, this approach was applied to read-sets for completed genomes and shown to be accurate over a wide range of microbial genome sizes.
Application of these procedures, based solely on k-mer abundances in short read data sets, allows aspects of genome structure to be resolved that are not apparent from conventional short read assemblies. This knowledge of the repetitive content of genomes provides insights into genome evolution and diversity.
KeywordsK-mer Genome assembly Repetitive elements Bacterial evolution
New sequencing technologies, by generating hundreds of millions of reads from multiplexed samples, allow the rapid and simultaneous acquisition of many genome sequences, and foster comparative analyses of closely related organisms. Despite achieving high coverage, the short reads generated by many sequencing platforms permit only partial assembly of genomes, due largely to the presence of numerous classes of repetitive sequences. Only regions of unique, single copy sequence and repeat regions shorter than the read length can be accurately assembled, yielding draft genomes that consist of unordered contigs separated by gaps of unspecified size . These technological limitations conceal the size and organization of a genome because the number, contents and distribution of repeat arrays remain unknown. To resolve the total size and repeat organization of a genome typically demands a complete genome assembly, which is ultimately obtained by linking contigs through combinatorial strategies and additional sequencing , which may sometimes be aided by additional finishing approaches, such as optical mapping .
The number, types and configuration of repeated sequences varies greatly within and across species. Even bacterial genomes, which consist mostly of single-copy protein-encoding genes, can contain multiple classes of repetitive sequences that can prevent their assembly and analysis. Bacterial genomes encode as many as 15 ribosomal DNA operons, which can each span several kilobases and are of nearly identical sequence within most organisms, complicating the assembly of their genomes . Even more problematic are the shorter repetitive elements and sequence repeats present in bacterial genomes [5, 6]. For example, most bacteria harbor multiple classes of insertion sequence (IS) elements, which average about a kilobase in length and whose copy numbers are highly variable, sometimes reaching hundreds of copies per genome [7, 8]. In addition, there are several other types of repetitive sequences, including the 21-to-65-nt repetitive extragenic palindromic (REP) sequences , the 127-bp palindromic enterobacterial repetitive intergenic consensus (ERIC) sequences , and the bacterial interspersed mosaic elements (BIME) , each of which can be dispersed throughout the genome or exist in tandem arrays .
Heterogeneity in the numerous classes of repetitive sequences contributes to the genome size variation that has been observed in many bacterial species. For example, isolates of Escherichia coli can differ by up to 30% in genome size, with sequenced strains ranging from 4.5 to 5.9 Mb [13–15]. Some of this size diversity is also attributable to the differential accumulation of plasmids and prophages. Accessory elements can also complicate the repeat structure of genomes because the multiple prophages within a genome can encode the same genes, and plasmids are often maintained in very high copy numbers.
Because high-coverage, short-read data contain information about repeat content that is not apparent in draft assemblies, we developed a procedure to estimate the size and repeat content of genomes from raw sequence data. This method relies on the analysis of the frequencies of overlapping fixed-length sequences (k-mers) , thereby exploiting a computationally efficient and scalable analytical technique that underlies many search and assembly algorithms. To date, k-mer-based algorithms have served as the basis for efficient similarity searching , vector removal , graph-based, short-read assembly [18, 19], short-read error correction [20, 21], estimation of genome repeat structure , comparison of whole genomes , identification of anomalous genome regions [24, 25] and binning of sequence fragments from different sources [26, 27].
We show that direct analysis of the abundances of unique 21-mers (icosihenamers) recovered from sets of whole-genome shotgun short read sequences yields robust estimates of total genome size and of the proportion of the genome represented by repeats of each copy number. Furthermore, sequence coverage can be estimated, sequence quality can be characterized, and sample contamination can often be diagnosed. Therefore, resolution of the repeat architecture of a genome can facilitate comparisons of genome size, structure and complexity, help the assembly of accurate genome sequences, and improve confidence in genome assemblies.
R esults and discussion
Analysis of short-read datasets from novel E. coli strains
Our choice of k-mer length, 21, was determined empirically and reflects the minimum size for these genomes that can reliably distinguish small repeats from single-copy sequences. Abundance analysis using smaller k-mers is sensitive to shorter repeats, resulting in larger proportions of sequence inferred to occur in multiple copies (Additional file 1: Table S1) and k-mers shorter than 15 nucleotides in length are more likely to match by chance alone. The wider range of repeat lengths spanned by longer k-mers causes abundance analysis to be less sensitive to short repeats, such that smaller fractions of the genome will be classified as part of a repeat. Discussions of k-mer length have treated it as a parameter to be optimized [28, 29] and as a sensitivity/specificity tuning parameter for search and assembly [28, 30]. It should be noted that a k-mer length of 21 is biologically relevant: it is short enough to resolve small repeated elements that occur in E. coli, such as REP and ERIC sequences [10, 24] but long enough to distinguish between single copy protein-encoding regions by spanning the sequence that defines such regions as unique.
There are several methods to efficiently enumerate k-mers, including dense representations , hashing , and probabilistic data structures . Technical limitations on computer memory and storage space make dense representations unreasonable for k-mers longer than k = 18; however, the hashing implementation used here (i.e., Jellyfish []) is applicable to k-mers as long as 31, and can in principle be applied to very large datasets (1012 bp). Due to reverse-complement degeneracy, the total number of possible k-mers is approximately (½) 4 k , and for k = 15 or greater, there are many more possible k-mers (>500 million) than basepairs in the largest observed bacterial and archaeal genomes (<13 Mb). The scarcity of long k-mers, together with the fact that most microbial sequence is non-redundant, has the consequence that most k- mers that do occur are present only once, contain sequence from single copy genes, and map to unique locations of the genome.
The multi-modal distributions of 21-mers at abundances above that of the principal peak reflect primarily the several classes of repetitive sequences that are present at different copy numbers in a genome, although unevenness in read coverage of the target genome can cause distortion of the peak shape. The successive peaks occur at multiples of the principal peak’s abundance, each of which reflects the repeat copy number of the sequence contained in those 21-mers, and the area under each peak is determined by the amount of unique sequence at that copy number. In all five E. coli genomes analyzed, the first peak is the largest, consistent with the expectation that complete E. coli genome sequences, and bacteria in general, are relatively gene-rich and consist primarily of single copy sequences.
Size and repeat structure of E. coli genomes estimated by icosihenamer (21-mer) analysis
Genome sequence (bp)
11 × -20×
21 × -79×
Relating abundance histograms to genome repeat structures
The number of 21-mers that reside under each peak in an abundance histogram corresponds closely to the number of basepairs of sequence at a particular repeat copy-number. To estimate the amount of unique sequence under each peak and to gain insight into the repeat structure of each target genome, we applied a maximum likelihood estimator to model the distribution of 21-mers in each histogram spectrum (red lines in Figure 1). This approach is similar to but has advantages over the method used by Li & Waterman  to estimate repeat structure from sequence reads. In our implementation, the inference of sequence repeats is guided by the natural property of 21-mer abundances to occur at integer multiples of the modal k-mer abundance (the principal peaks in Figure 1). Whereas our fitting procedure exploits this property to stabilize the likelihood optimization, the implementation of Li & Waterman  uses the expectation-maximization algorithm. In addition, our approach was evaluated using high-coverage, short-read raw sequencing data and is optimized for current technologies.
The amount of unique sequence, as well as the total amounts of the genome sequence accounted for at each repeat copy-number, vary among the five strains indicating that the strains differ in their repeat structure and numbers of multicopy elements (Table 1). The choice of k limits the shortest repeat that can be resolved. Our choice of k = 21 allows resolution of the small repeated elements that typically occur in E. coli and other bacteria. Resolution of the repeat structure of a genome also provides an estimate of the total genome size (Table 1, bottom row). Despite different values of k providing different interpretations of repeat structure (Additional file 1: Table S1), estimates of genomes size are almost unaffected. Other choices for k from 15 to 29 changed the genome size estimates among our five E. coli datasets by less than 0.25% (Additional file 2: Table S2).
A common characteristic of the repeat structures of these E. coli genomes is a trend towards less unique sequence at higher copy numbers, except around seven times the abundance of the principal peak. The peaks at this copy number correspond to the rrn operons and any other sequence repeated a similar number of times. The rrn operons are present in seven copies in virtually all strains of E. coli and were observed in the 21-mers counted seven times in the completely sequenced E. coli DH1 genome (Table 1).
Despite the similarities in repeat structure, 21-mer abundance analysis using Kmerspectrumanalyzer allows differences in repeat structure between these five novel E. coli genomes to be resolved. The amount of unique sequence at six to eight copies per genome ranges from 4.8 kb in strain A_03_34 to 16.1 kb in strain E_01_37. While E. coli rrn operons are approximately 5.4 kb and expected to be present at seven copies, they are often non-identical at a few positions. Because 21-mer abundance analysis only resolves identical sequence repeats, the two smallest estimates of sequence at six to eight copies were 4.8 kb in strain A_03_34 and 4.9 kb in strain B_04_28 may not be an underestimation but an accurate reflection of the identical sequence among rrn operons. Any differences among the rrn operon sequences would be included in the total genome size estimate but at copy numbers less than seven. The estimates of sequence at six to eight copies in strains C_04_22, D_04_22 and E_01_37 were 9.1, 8.9 and 16.1 kb respectively and are greater than can be explained by the presence of seven rrn operons alone. These larger estimates are evidence of additional repeated sequences at abundances similar to that of the rrn operons.
There is little sequence repeated at high copy numbers within strains A_03_34 and B_04_28: the repeats in nine or more copies are represented by less than 121 bp of unique sequence in both cases. In contrast, strains D_04_22 and E_01_37 possess 1.3 kb and 3.0 kb respectively, of sequence repeated nine or more times. The 1.3 kb of high copy-number sequence from strain D_04_22 contains two protein-coding regions sharing high similarity with sequences annotated as IS911 transposases (GenBank accession number AY555729.1), and some of the 3-kb sequence at high copy number in strain E_01_37 shares high sequence similarity to plasmid pCE10B from E. coli O7:K1 strain CE10 (GenBank accession number: CP003036.1).
Accuracy of 21-mer analysis for determining microbial genome sizes
The accuracy of estimates of total genome size obtained by our fitting procedure was tested in three ways. First, we assessed the sensitivity for discriminating the genome sizes from the whole-genome shotgun short-read datasets of five novel strains of E. coli, whose genome sizes were also estimated by totaling the lengths of I-CeuI restriction fragments resolved by pulsed-field gel electrophoresis (PFGE). Next, we tested the range of microbial genome sizes over which the procedure remains accurate by comparing the size estimates made from publicly available short-read datasets to the lengths of the corresponding complete genome sequences. Finally, we assessed the robustness of these genome size estimates by varying the read depths on target genomes and varying error content by quality-score-based read trimming.
Total genome sizes of five E. coli strains estimated by PFGE and icosihenamer analysis
Genome size (Mb)
21-mer analysis of sequencing reads
To assess the accuracy of the PFGE estimates, we also sized E. coli strain MG1655, for which the genome sequence is available , by PFGE, and compared fragment lengths to those expected based on the I-CeuI sites in the sequenced genome. Except for the largest I-CeuI restriction fragment, fragment length sizes estimated by PFGE were very close to those predicted by in silico digestion of the genome and are listed in Additional file 3: Table S3 along with each fragment size for the novel E. coli strains. Due to the vagaries of resolving and sizing fragments > 1 Mb by PFGE, we estimated the size of the largest I-CeuI in the E. coli MG1655 genome to be 2.793 Mb, whereas it is actually 2.498 Mb according to published sequence .
Although the genome sizes determined by PFGE were similar to those estimated by 21-mer analysis of the short-read data, there remains some discrepancy between the two methods. One source of the difference is that the PFGE only accurately assesses the sizes of linear DNA fragments. Bacterial genomes often contain small circular plasmids that remain intact after I-CeuI digestion and these circular elements are not included when tabulating genome size. The reported lengths of complete genome sequences may differ from size estimates by 21-mer analysis of read-sets because sequences, such as those encoded on multicopy plasmids, are represented in read data in proportion to their actual molecular quantities. In contrast, published genome sizes do not accommodate such differences in copy number leading to disparities in the canonical genome size and actual DNA content of the cell.
Some of the differences, particularly the overestimations of genome size by 21-mer analysis, may be caused by the presence of multicopy plasmids: these elements are always counted only once when reporting a completed genome size but are included in proportion to their actual copy numbers in 21-mer analysis. The over-estimate of the 7.6 Mb Cylindrospermum stagnale PCC 7417 genome (Additional file 4: Table S4) by only 6,600 bp is consistent with single copy plasmids and a small degree of error, and the genome size estimated for Nostoc sp. PCC 7524 could accommodate two copies of its smallest plasmid (6,361 bp) per chromosome. Similarly, the genome size estimated for Thermovirga lienii DSM 17291 accommodates two or three copies of its 31,872 bp plasmid in addition to its chromosome. E. coli KO11FL carries a single plasmid, pRK2, which has been shown experimentally to be maintained in E. coli strains at between 25 to 40 molecules per chromosomal equivalent [35, 36]. We estimate the genome size of E. coli KO11FL to be 181,588 bp greater than its published chromosome size (5,021,812 bp), yielding ≈ 34 copies of the 5,360 bp pRK2 plasmid, in agreement with its known copy numbers. Figure 2 also displays the genome size estimates of the five novel E. coli strains as sized by PFGE demonstrating close agreement between the assessment of accuracy within a single species and across genome sizes of three orders of magnitude.
Our method is most accurate when read coverage over the entire genome is relatively even, since it provides k-mer abundance spectra with mixed Poisson distributions. One limitation in our approach is how well the proposed model can account for systematic errors in contemporary sequencing methods and still retain accuracy. Uneven read coverage broadens the peaks in the k-mer abundance spectra, sometimes causing errors in the assignment of sequences to a particular read coverage. Other factors that bias the k-mer abundance spectra include random sequencing errors and sequences derived from non-target DNA. Sequencing errors replace high-abundance sequences among reads with comparatively rare sequences creating many novel low-abundance k-mers. The fitting procedure implemented in Kmerspectrumanalyzer is designed to exclude these low-abundance k-mers by applying a heuristic low-coverage cutoff. We tested whether the genome sizes of the five novel E. coli, as estimated by 21-mer abundance analysis, were affected by removal of read positions likely to contain errors. We applied a trimming algorithm that retained the longest contiguous read positions with a less than 10% likelihood of being incorrect (i.e., a quality score of 10). This trimming procedure had little effect on the estimates of overall genome size, which changed by a maximum of only 0.19% (Additional file 5: Table S5). Therefore, Kmerspectrumanalyzer is sufficiently robust against sequencing error to negate any need for the pre-processing of read data for error mitigation.
The 21-mer abundance analysis described here would also be appropriate to the whole-genome shotgun sequencing datasets of eukaryotes of known ploidy. Although a k mer length of 21 is sufficient to ensure sparse sampling for microbial genomes, greater values for k may be necessary for the largest eukaryotic genomes. However, genomes with more complex repeat structure and lower sequencing coverage depths, can present challenges to the estimation approach presented here..In principle, this 21-mer abundance analysis can be applied to any sequencing methodology; however, there are some additional limitations. Extremely low coverage datasets (less than ≈ 10×) and datasets with very high error rates (>5%) contain insufficient sampling of the true k-mers to support numerical inferences. While numerical inferences are supported above 10× coverage depth, accuracy is compromised below about 75× coverage depth (Additional file 6: Figure S1). Among the four example datasets we tested, only the estimates for the Owenweeksia hongkongensis DSM 17368 reads (SRA run number: SRR190843) were relatively inaccurate between 75–400× coverage; however, it never differed by more than 4% from the completely assembled genome length. k-mer spectrum approaches can measure sequencing error, characterize heterozygosity, strain variation or mixtures of organisms  within a sample in an annotation-independent, scalable way, and therefore, are applicable to a several sequencing applications, including the interpretation of metagenomic data sets and diagnosing technical aspects of the sequencing procedure.
Sequence repeats severely hamper the assembly of most genomes, and these repeats continue to obscure genome structure even with the high read depths afforded by new sequencing technologies. We provide a simple and rapid means to resolve the repeat structure and total size of a genome by directly analyzing the abundance of distinct k-mers among short reads. By obtaining the genome size and repeat structure of environmental isolates of E. coli from 76-bp Illumina reads, we demonstrate that the sensitivity of this method is great enough to resolve differences among bacterial strains. Total genome size estimates corresponded well with those obtained for the same strains by long-range RFLP mapping on pulsed-field gels. In addition, inferences of genome size from short-read datasets are not limited to strains within species but are accurate across a wide range of genome sizes from 0.005 to 9 Mb. The fitting procedure introduced here depends only on the k-mer spectrum, which summarizes the sequence redundancy of the data set but preserves none of the sequence content.
Previous work on plants  and humans  has exploited k-mer indices of raw sequence data to make inferences about the genome structure, particularly about the frequencies and contents of repetitive elements. More current analyses based on high-throughput sequencing data have proposed a variety of heuristics for estimating genome size by identifying features in the k-mer spectrum  or by using arbitrary cutoffs in the k-mer spectrum . However, such approaches have failed to exploit the fact that repetitive sequences must occur at integer multiples of single copy sequences. By incorporating this information, our method implements a more biologically accurate model of the k-mer abundance distribution, adapting the procedure of  to high-coverage (>30×), short-read sequences.
Accurate estimates of total unique sequence from raw reads, as provided by Kmerspectrumanalyzer, can also aid in the optimization of read assemblies. Whereas the accuracy of an assembly is generally more important than its assembled length, the extent to which the total length of assembled contigs accounts for the actual length of the sequenced genome is a potential metric for the completeness of an assembly. The length of the target sequence is usually unknown, but it can be obtained directly from raw reads by 21-mer analysis, thereby allowing an absolute measure of assembly success. Moreover, abundance distributions in k-mer spectra are also useful in identifying samples that are dominated by PCR artifacts, samples contaminated with sequencing adapters, and samples that contain mostly positive-control calibration genomes (e.g., PhiX174, in the case of Illumina) without performing assembly or similarity searches.
We have observed that the reliability of many widely used statistics for summarizing the k-mer spectrum, including the number of unique k- mers, the fraction of k- mers that are unique, or the fraction of k- mers above an arbitrary or model-fitted threshold, are sensitive to changes in sampling depth. Two summary statistics that appear to be stable descriptions of certain informative properties of the k-mer spectrum are (i) the rank order distribution of k-mers and (ii) the cumulative fraction of the data consumed by rank ordered k-mers. In addition to genome repeat structure resolution, the Kmerspectrumanalyzer package includes tools to visualize the rank-order distribution of k- mers and assist in understanding sequence datasets, even in cases where fitting the spectrum to a mixed over-dispersed Poisson model fails.
Whole-genome shotgun sequencing applies technologies that generate hundreds of millions of short reads and has recently become the most widely used tool in genome analysis . We introduce a straightforward methodology that provides information about the repeat structure of genomes that is ordinarily missing from assemblies of short reads. This additional information offers new insights about genome diversity and evolution that can be gained through the analysis of novel datasets or through the re-analysis of the large volumes of archived short read.
k-mer counting in sequence reads
The abundances of all overlapping 21-bp sequences present in a set of whole-genome shotgun short-read sequences were counted using Jellyfishk-mer counting library (vers. 1.1.5). No read trimming or error-correcting algorithms were applied. The frequencies of different k-mers (in this case 21-mers) at each abundance value contained in a set of sequences are plotted as a k-mer abundance spectrum (sometimes referred to as an abundance or coverage histogram).
A repeated sequence in a sampled genome affects the shape of these k-mer abundance spectra depending on its length and copy number. A sequence of length l will contain (l – k + 1) different k- mers if it does not contain repeats of length greater than k – 1. If the same sequence occurs n times in a genome, shotgun sequencing would sample these k-mers n times more often than those that are single-copy in a genome, resulting in (l – k + 1) k-mers with abundances n times higher than the average read-depth based on the number of sequencing reads. Therefore, repeated sequences in the genome results in higher abundances of the corresponding k-mers. These collections of k-mers at higher-than-normal abundances appear as multiple peaks at different positions along the x-axis of the k-mer abundance spectrum.
Relating k-mer abundance to genome size and repeat abundance
We modeled the abundance distribution of k- mers as a mixture of over-dispersed Poisson (negative binomial) distributions in which the mixture coefficients represent amounts of sequence at each copy number and the overdispersion parameter accommodates uneven read-depth across the target genome. We refer to the peak in the abundance spectrum that consists of the k-mers corresponding to single-copy sequence in the target genome as the principal peak. Peaks in the abundance spectrum that have greater abundance than the principal peak correspond to sequences with different levels of abundance (i.e., repeat copy numbers), which are expected to occur at integer multiples of the principal abundance.
NegBinomial is the mu-alpha parameterization of the negative binomial, where NegBinomial (x; mu, alpha) gives the negative binomial distribution with mean = mu and variance = (1 + alpha)*mu2. When comparing plots of the actual 21-mer spectrum of E. coli K12 to the P N and P NO models calculated from 21-mer counts in the reference genome, we found that the mixed-Poisson model showed multiple peaks that fail to match the shape and general character of the peaks in the actual spectrum, whereas the over-dispersed adequately models the shape of the peaks of the spectrum.
Here, the sum includes the full range of observed k-mer abundances, including values of i, for which z i is zero, and P NO ( ) is the sum of terms defined above.
Fitting and implementation
The first three terms in the equations above, a 1 , c, and s, describe the height, abundance, and width of the principal peak in the k-mer spectrum. Twenty-nine additional terms, fitted sequentially, describe the height of peaks centered at integer multiples of the abundance of the principal peak ranging from 2× to 30×. Although sequencing errors are not explicitly modeled, by excluding points in the abundance spectrum at abundances less than half of the fitted abundance of the main peak, low-abundance sequencing artifacts effectively do not affect model fitting. We tested this assertion by estimating genome sizes with and without the quality score-based trimming of read positions. DynamicTrim from the SolexaQA package version 2.2  was used to trim reads. The fitting procedure also employs several heuristics, including the masking of parts of the spectrum that are out of the range of abundances being fitted, weighted-least-squares optimization to set initial values, non-negativity constraints on the sequence size parameters, and successive fitting of low-order (and lower-dimensional) models before high-order models, all of which serve to stabilize the maximum likelihood optimization on real data.
Short-read datasets were downloaded from the NCBI Sequence Read Archive  and the corresponding complete genome sequences were obtained from the NCBI RefSeq database . For this analysis, the smallest genome size was that of phiX174 (5,386 bp) and the largest was Niastella koreensis GR20-10 (9,033,684 bp). To fit the phiX174 genome, the procedure required a manual low-abundance cutoff of 10,000× and was constrained to fit only 1 term in the mixture model because of the extremely high coverage in that dataset (100,000×); otherwise, all of the microbial genomes were analyzed with the same procedures and parameters. To aid in the interpretation of the k-mer spectra (and the estimation of repeat copy numbers) for the five E. coli strains, we also queried and tabulated the amounts of repeated sequence in the complete genome sequence of E. coli DH1 (NC_017625).
The fitting procedure was implemented in Python 2.7.2 using Numpy, Scipy, and Matplotlib. An open-source implementation of this tool, Kmerspectrumanalyzer, and the scripts used to retrieve, process, and produce the numerical data in these analyses are available at http://github.com/MG-RAST/kmerspectrumanalyzer.git. The tool is available as a module in the KBase sequence analysis framework , allowing end-users to perform these analyses on a third-party, scalable computing infrastructure rather than their own hardware. Motivated by previous work  our procedure is intended to be conveniently reproducible.
Sampling, isolation and strain characterization of environmental Escherichia coli
Strains of Escherichia coli were isolated from liquid samples collected in July 2008 at the Central Contra Costa Sanitary District Treatment Plant in Martinez, California. Samples were titrated to a final concentration of 15% glycerol and stored at −80°C. An aliquot of each sample was diluted 1:100 in LB broth, and 100 μl was plated on MacConkey agar and incubated overnight at 37°C.
Colonies of E. coli were initially selected based on colony morphology and then typed genetically by sequencing fragments of three diagnostic loci, fumC, gyrB and adk, used in the multilocus sequence typing (MLST) analysis of E. coli. MLST proceeded by the colony PCR  using primer pairs: fumC forward 5′-TCA CAG GTC GCC AGC GCT TC-3′ and fumC reverse 5′-GTA CGC AGC GAA AAA GAT TC-3′; gyrB forward 5′-TCG GCG ACA CGG ATG ACG GC-3′ and gyrB reverse 5′-ATC AGG CCT TCA CGC GCA TC-3′, adk forward 5′-ATT CTG CTT GGC GCT CCG GG-3′ and adk reverse 5′-CCG TCA ACT TTC GCG TAT TT-3′. PCR fragments were verified by agarose gel electrophoresis, and prepared for sequencing through the addition of 0.2 μl exonuclease and 0.2 μl calf intestinal phosphatase. Sanger sequencing of purified PCR products was performed at the University of Arizona Genetics Core for Sanger sequencing. Strains that differed in nucleotide sequence from any of the >3000 MLST-typed E. coli at any of the three diagnostic loci were deemed as “unique” and stored in LB broth supplemented with 15% glycerol and stored at −80°C.
Whole-genome shotgun sequencing
Unique strains were streaked onto MacConkey agar and grown overnight at 37°C. Individual colonies were transferred to 2 ml LB broth and grown overnight at 37°C. DNA was isolated by disrupting pelleted cells in 1 ml of TES buffer containing 50mM NaCl, 50 mM Tris–HCl, 50 mM EDTA, 5% SDS (pH 7.6) followed by mechanical lysis with 0.1 mm zirconia/silica beads (BioSpec). The aqueous phase was removed, and treated with an equal volume of phenol/chloroform/isoamyl alcohol pH 7.9 (Ambion), followed by chloroform extraction and the isopropanol precipitation of DNA. Purified DNAs were quantified with Quant-iT PicoGreen (Invitrogen) and submitted to the Yale Center for Genome Analysis for library preparation and sequencing. Whole-genome shotgun sequencing was performed on an Illumina HiSeq 2000 generating 76-bp paired reads from the ends of 155-bp fragments. Short-read data were processed with the CASAVA 1.8.2 package.
Sizing genomes by pulsed-field gel electrophoresis
Estimates of the genome sizes of the five newly sequenced strains of E. coli and of a control strain (E. coli MG1655) were determined by pulsed-field gel electrophoresis following methods described in . In short, cells were grown in 5 ml of LB broth and treated with 180 μg/ml chloramphenicol to align chromosomes. Cells were harvested by centrifugation, washed and resuspended in 0.5 ml of TEN (10 mM Tris [pH 7.5], 100 mM EDTA [pH 8], 250 mM NaCl), and mixed with 0.75 ml of 1.5% low melt temperature agarose (RPI, Mount Prospect, IL) in TEN. Agarose plugs were incubated for 21 hours in lysis solution (0.1% lysozyme, 0.002% RNase, 0.5% Sarkosyl, 10 mM Tris [pH 7.5], 100 mM EDTA [pH 8], 250 mM NaCl), with subsequent overnight incubation at 45°C in 250 mM EDTA containing 0.1% proteinase K and 1% Sarkosyl. To inactivate excess proteases, agarose plugs were incubated in 1 mM PMSF for 1 hr, and washed and stored in 10mM Tris, 100 mM EDTA (pH 8) at 4°C. Agarose plugs were washed five times, each for 20 min, in 50 volumes of distilled H2O, and equilibrated in NEBuffer 4 with 0.1 mg/ml BSA (New England BioLabs, Ipswich, MA). Five units of I-CeuI restriction endonuclease (NEB) was added to initiate digestion, and after incubation overnight at 37°C, EDTA was added to each sample to a final concentration of 0.1 M to terminate digestion.
Electrophoresis was performed on a CHEF-DR II apparatus (Bio-Rad Laboratories, Richmond, CA) in 0.5× Tris-borate-EDTA at 14°C. To resolve DNA fragments in the 20 to 200 kb range, pulse durations were ramped from 5 sec to 12.5 sec; for DNA fragments in the 400 to 800 kb range, pulse durations were ramped from 60 sec to 100 sec. In both cases, electrophoresis proceeded for 24 hr in a 0.9% agarose gel at 6 V/cm. To resolve DNA fragments in the 2000 to 3000 kb range, pulse durations were ramped from 600 sec to 960 sec for 90 hr in a 0.7% agarose gel at 2.4 V/cm. Gels were stained and photographed digitally, and TIFF files of these images were loaded using the tifffile.py Python module (version 2013.01.18, ). Fragment sizes were estimated by interpolation to standards of known size using second order splines implemented in Scipy. Python 2.7.2 source code is included in the Kmerspectrumanalyzer repository in the ‘pfge_analysis’ folder .
Repetitive extragenic palindromic
Enterobacterial repetitive intergenic consensus
Bacterial interspersed mosaic elements
The authors wish to thank Travis Harrison from the KBase development team for developing a wrapper for Jellyfish. This work was supported by grant number 23536 from the Templeton Foundation to HO, and the Offices of Science and of Biological and Environmental Research of the US Department of Energy, as part of the DOE Systems Biology Knowledgebase (under Contract No. DE-AC02-06CH11357 to WT and FM).
- Ricker N, Qian H, Fulthorpe R: The limitations of draft assemblies for understanding prokaryotic adaptation and evolution. Genomics. 2012, 100: 167-175. 10.1016/j.ygeno.2012.06.009.View ArticlePubMedGoogle Scholar
- Eppinger M, Rosovitz MJ, Fricke WF, Rasko DA, Kokorina G, Fayolle C, Lindler LE, Carniel E, Ravel J: The complete genome sequence of Yersinia pseudotuberculosis IP31758, the causative agent of far east scarlet-like fever. PLoS Genet. 2007, 3: e142-10.1371/journal.pgen.0030142.PubMed CentralView ArticlePubMedGoogle Scholar
- Latreille P, Norton S, Goldman B, Henkhaus J, Miller N, Barbazuk B, Bode H, Darby C, Du Z, Forst S, Gaudriault S, Goodner B, Blair HG, Slater S: Optical mapping as a routine tool for bacterial genome sequence finishing. BMC Genomics. 2007, 8: 321-10.1186/1471-2164-8-321.PubMed CentralView ArticlePubMedGoogle Scholar
- Rainey F, Ward-Rainey N, Janssen P, Hippe H, Stackebrandt E: Clostridium paradoxum DSM 7308(T) contains multiple 16S rRNA genes with heterogeneous intervening sequences. Microbiology. 1996, 142: 2087-2095. 10.1099/13500872-142-8-2087.View ArticlePubMedGoogle Scholar
- Lupski J, Weinstock G: Short, interspersed repetitive DNA sequences in prokaryotic genomes. J Bacteriol. 1992, 174: 4525-4529.PubMed CentralPubMedGoogle Scholar
- Bachellier S, Hofnung M: The members of the so-called BEE95 family of dispersed enterobacterial intergenic elements are already known under the name RSA sequences. Mol Microbiol. 1996, 19: 1385-10.1111/j.1365-2958.1996.tb02481.x.View ArticleGoogle Scholar
- Lawrence JG, Ochman H, Hartl DL: The evolution of insertion sequences within enteric bacteria. Genetics. 1992, 131: 9-20.PubMed CentralPubMedGoogle Scholar
- Wagner A: Periodic extinctions of transposable elements in bacterial lineages: Evidence from intragenomic variation in multiple genomes. Mol Biol Evol. 2006, 23: 723-733. 10.1093/molbev/msj085.View ArticlePubMedGoogle Scholar
- Stern M, Ames G, Smith N: Repetitive extragenic palindromic sequences: A major component of the bacterial genome. Cell. 1984, 37: 1015-1026. 10.1016/0092-8674(84)90436-7.View ArticlePubMedGoogle Scholar
- Wilson L, Sharp P: Enterobacterial repetitive intergenic consensus (ERIC) sequences in Escherichia coli: Evolution and implications for ERIC-PCR. Mol Biol Evol. 2006, 23: 1156-1168. 10.1093/molbev/msj125.View ArticlePubMedGoogle Scholar
- Gilson E, Saurin W, Perrin D, Bachellier S, Hofnung M: The BIME family of bacterial highly repetitive sequences. Res Microbiol. 1991, 142: 217-222. 10.1016/0923-2508(91)90033-7.View ArticlePubMedGoogle Scholar
- Aranda-Olmedo I, Tobes R, Manzanera M, Ramos JL, Marqués S: Species-specific repetitive extragenic palindromic (REP) sequences in Pseudomonas putida. Nucleic Acids Res. 2002, 30: 1826-1833. 10.1093/nar/30.8.1826.PubMed CentralView ArticlePubMedGoogle Scholar
- Rasko D, Rosovitz M, Myers G, Mongodin E, Fricke W, Gajer P, Crabtree J, Sebaihia M, Thomson N, Chaudhuri R, Henderson I, Sperandio V, Ravel J: The pangenome structure of Escherichia coli: Comparative genomic analysis of E. coli commensal and pathogenic isolates. J Bacteriol. 2008, 190: 6881-6893. 10.1128/JB.00619-08.PubMed CentralView ArticlePubMedGoogle Scholar
- Ogura Y, Ooka T, Iguchi A, Toh H, Asadulghani M, Oshima K, Kodama T, Abe H, Nakayama K, Kurokawa K, Tobe T, Hattori M, Hayashi T: Comparative genomics reveal the mechanism of the parallel evolution of O157 and non-O157 enterohemorrhagic Escherichia coli. Proc Natl Acad Sci U S A. 2009, 106: 17939-17944. 10.1073/pnas.0903585106.PubMed CentralView ArticlePubMedGoogle Scholar
- Touchon M, Hoede C, Tenaillon O, Barbe V, Baeriswyl S, Bidet P, Bingen E, Bonacorsi S, Bouchier C, Bouvet O, Calteau A, Chiapello H, Clermont O, Cruveiller S, Danchin A, Diard M, Dossat C, El-Karoui M, Frapy E, Garry L, Ghigo J, Gilles A, Johnson J, Le-Bouguénec C, Lescat M, Mangenot S, Martinez-Jéhanne V, Matic I, Nassif X, Oztas S, Petit M, Pichon C, Rouy Z, Ruf C, Schneider D, Tourret J, Vacherie B, Vallenet D, Médigue C, Rocha E, Denamur E: Organised genome dynamics in the Escherichia coli species results in highly diverse adaptive paths. PLoS Genet. 2009, 5: e1000344-10.1371/journal.pgen.1000344.PubMed CentralView ArticlePubMedGoogle Scholar
- Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ: Basic local alignment search tool. J Mol Biol. 1990, 215: 403-410.View ArticlePubMedGoogle Scholar
- White JR, Roberts M, Yorke JA, Pop M: Figaro: a novel statistical method for vector sequence removal. Bioinformatics. 2008, 24: 462-467. 10.1093/bioinformatics/btm632.PubMed CentralView ArticlePubMedGoogle Scholar
- Idury RM, Waterman MS: A new algorithm for DNA sequence assembly. J Comput Biol. 1995, 2: 291-306. 10.1089/cmb.1995.2.291.View ArticlePubMedGoogle Scholar
- Pevzner PA, Tang H: Fragment assembly with double-barreled data. Bioinformatics. 2001, 17: S225-S233. 10.1093/bioinformatics/17.suppl_1.S225.View ArticlePubMedGoogle Scholar
- Chaisson M, Pevzner P, Tang H: Fragment assembly with short reads. Bioinformatics. 2004, 20: 2067-2074. 10.1093/bioinformatics/bth205.View ArticlePubMedGoogle Scholar
- Simpson JT, Durbin R: Efficient de novo assembly of large genomes using compressed data structures. Genome Res. 2012, 22: 549-556. 10.1101/gr.126953.111.PubMed CentralView ArticlePubMedGoogle Scholar
- Li X, Waterman MS: Estimating the repeat structure and length of DNA sequences using L-tuples. Genome Res. 2003, 13: 1916-1922.PubMed CentralPubMedGoogle Scholar
- Sims GE, Jun S-RR WGA, Kim S-HH: Alignment-free genome comparison with feature frequency profiles (FFP) and optimal resolutions. Proc Natl Acad Sci U S A. 2009, 106: 2677-2682. 10.1073/pnas.0813249106.PubMed CentralView ArticlePubMedGoogle Scholar
- Higgins C, Ferro-Luzzi G, Ames W: A novel intercistronic regulatory element of prokaryotic operons. Nature. 1982, 298: 760-762. 10.1038/298760a0.View ArticlePubMedGoogle Scholar
- Monier A, Claverie J-MM, Ogata H: Horizontal gene transfer and nucleotide compositional anomaly in large DNA viruses. BMC Genomics. 2007, 8: 456-10.1186/1471-2164-8-456.PubMed CentralView ArticlePubMedGoogle Scholar
- Brady A, Salzberg SL: Phymm and PhymmBL: metagenomic phylogenetic classification with interpolated Markov models. Nat Methods. 2009, 6: 673-676. 10.1038/nmeth.1358.PubMed CentralView ArticlePubMedGoogle Scholar
- Liu J, Wang H, Yang H, Zhang Y, Wang J, Zhao F, Qi J: Composition-based classification of short metagenomic sequences elucidates the landscapes of taxonomic and functional enrichment of microorganisms. Nucleic Acids Res. 2013, 41: e3-10.1093/nar/gks828.PubMed CentralView ArticlePubMedGoogle Scholar
- Jun S-R, Sims GE, Wu GA, Kim S-H: Whole-proteome phylogeny of prokaryotes by feature frequency profiles: An alignment-free method with optimal feature resolution. Proc Natl Acad Sci U S A. 2010, 107: 133-138. 10.1073/pnas.0913033107.PubMed CentralView ArticlePubMedGoogle Scholar
- Gruenheit N, Deusch O, Esser C, Becker M, Voelckel C, Lockhart P: Cutoffs and k-mers: implications from a transcriptome study in allopolyploid plants. BMC Genomics. 2012, 13: 92-10.1186/1471-2164-13-92.PubMed CentralView ArticlePubMedGoogle Scholar
- Campagna D, Romualdi C, Vitulo N, Del-Favero M, Lexa M, Cannata N, Valle G: RAP: a new computer program for de novo identification of repeated sequences in whole genomes. Bioinformatics. 2005, 21: 582-588. 10.1093/bioinformatics/bti039.View ArticlePubMedGoogle Scholar
- Kelley DR, Schatz MC, Salzberg SL: Quake: quality-aware detection and correction of sequencing errors. Genome Biol. 2010, 11: R116-10.1186/gb-2010-11-11-r116.PubMed CentralView ArticlePubMedGoogle Scholar
- Marçais G, Kingsford C: A fast, lock-free approach for efficient parallel counting of occurrences of k-mers. Bioinformatics. 2011, 27: 764-770. 10.1093/bioinformatics/btr011.PubMed CentralView ArticlePubMedGoogle Scholar
- Pell J, Hintze A, Canino-Koning R, Howe A, Tiedje JM, Brown CT: Scaling metagenome sequence assembly with probabilistic de Bruijn graphs. Proc Natl Acad Sci U S A. 2012, 109: 13272-13277. 10.1073/pnas.1121464109.PubMed CentralView ArticlePubMedGoogle Scholar
- Blattner FR, Plunkett G, Bloch CA, Perna NT, Burland V, Riley M, Collado-Vides J, Glasner JD, Rode CK, Mayhew GF, Gregor J, Davis NW, Kirkpatrick HA, Goeden MA, Rose DJ, Mau B, Shao Y: The complete genome sequence of Escherichia coli K-12. Science. 1997, 277: 1453-1462. 10.1126/science.277.5331.1453.View ArticlePubMedGoogle Scholar
- Turner PC, Yomano LP, Jarboe LR, York SW, Baggett CL, Moritz BE, Zentz EB, Shanmugam KT, Ingram LO: Optical mapping and sequencing of the Escherichia coli KO11 genome reveal extensive chromosomal rearrangements, and multiple tandem copies of the Zymomonas mobilis pdc and adhB genes. J Ind Microbiol Biotechnol. 2012, 39: 629-639. 10.1007/s10295-011-1052-2.View ArticlePubMedGoogle Scholar
- Stepánek V, Valesová R, Kyslík P: Cryptic plasmid pRK2 from Escherichia coli W: sequence analysis and segregational stability. Plasmid. 2005, 54: 86-91. 10.1016/j.plasmid.2004.12.006.View ArticlePubMedGoogle Scholar
- Angly F, Willner D, Prieto-Davó A, Edwards R, Schmieder R, Vega-Thurber R, Antonopoulos D, Barott K, Cottrell M, Desnues C, Dinsdale E, Furlan M, Haynes M, Henn M, Hu Y, Kirchman D, McDole T, McPherson J, Meyer F, Miller R, Mundt E, Naviaux R, Rodriguez-Mueller B, Stevens R, Wegley L, Zhang L, Zhu B, Rohwer F: The GAAS metagenomic tool and its estimations of viral and microbial average genome size in four major biomes. PLoS Comput Biol. 2009, 5: e1000593-10.1371/journal.pcbi.1000593.PubMed CentralView ArticlePubMedGoogle Scholar
- Kurtz S, Narechania A, Stein J, Ware D: A new method to compute K-mer frequencies and its application to annotate large repetitive plant genomes. BMC Genomics. 2008, 9: 517-10.1186/1471-2164-9-517.PubMed CentralView ArticlePubMedGoogle Scholar
- Healy J, Thomas EE, Schwartz JT, Wigler M: Annotating large genomes with exact word matches. Genome Res. 2003, 13: 2306-2315. 10.1101/gr.1350803.PubMed CentralView ArticlePubMedGoogle Scholar
- Wu Y-WW, Ye Y: A novel abundance-based algorithm for binning metagenomic sequences using l-tuples. J Comput Biol. 2011, 18: 523-534. 10.1089/cmb.2010.0245.PubMed CentralView ArticlePubMedGoogle Scholar
- Loman N, Constantinidou C, Chan J, Halachev M, Sergeant M, Penn C, Robinson E, Pallen M: High-throughput bacterial genome sequencing: An embarrassment of choice, a world of opportunity. Nat Rev Microbiol. 2012, 10: 599-606. 10.1038/nrmicro2850.View ArticlePubMedGoogle Scholar
- Cox MP, Peterson DA, Biggs PJ: SolexaQA: At-a-glance quality assessment of Illumina second-generation sequencing data. BMC Bioinforma. 2010, 11: 48-10.1186/1471-2105-11-48.View ArticleGoogle Scholar
- The National Center for Biotechnology Information Sequence Read Archive.http://www.ncbi.nlm.nih.gov/sra,
- The National Center for Biotechnology Information Reference Sequence Database.http://www.ncbi.nlm.nih.gov/refseq/,
- The KMERSPECTRUMANALYZER software repository.http://github.com/MG-RAST/kmerspectrumanalyzer.git,
- The Department of Energy Systems Biology Knowledgebase.http://kbase.us,
- Our approach to replication in computational science, C. Titus Brown Web Log.http://ivory.idyll.org/blog/replication-i.html,
- Wirth T, Falush D, Lan R, Colles F, Mensa P, Wieler LH, Karch H, Reeves PR, Maiden MCJ, Ochman H, Achtman M: Sex and virulence in Escherichia coli: an evolutionary perspective. Mol Microbiol. 2006, 60: 1136-1151. 10.1111/j.1365-2958.2006.05172.x.PubMed CentralView ArticlePubMedGoogle Scholar
- Bisercić M, Ochman H: Natural populations of Escherichia coli and Salmonella typhimurium harbor the same classes of insertion sequences. Genetics. 1993, 133: 449-454.PubMed CentralPubMedGoogle Scholar
- Bergthorsson U, Ochman H: Heterogeneity of genome sizes among natural isolates of Escherichia coli. J Bacteriol. 1995, 177: 5784-5789.PubMed CentralPubMedGoogle Scholar
- Python Modules by Christoph Gohlke, University of California, Irvine.http://www.lfd.uci.edu/~gohlke/,
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.