EMIRGE overview: iterative mapping and correction of reference SSU sequences

As full-length rRNA genes could not be recovered from the de novo assemblies, we developed a strategy based on mapping of all reads to a large reference database of SSU sequences and iterative determination of the most probable full gene sequences. We chose as a reference database a filtered subset of the SSU sequences contained in Silva [15]. An ideal mapping strategy would not depend on the completeness or correctness of the reference database, and evaluate the probabilities of errors in the mapping and the sequence of the reads. Therefore, we developed a method that models reads as being generated by unknown reference SSU sequences. Each iteration consists of mapping reads to the current reference sequences, computing the probability that each read was generated by each reference sequence, computing estimates of reference sequence abundance, and then correcting the nucleotides of the reference sequences before the subsequent iteration and mapping begins anew. Gradually, the correct reference sequences and the estimates of organism abundance adjust and then stabilize, at which point the iterations stop (Figure ​(Figure2).2). At each iteration, if two reference sequences have evolved to be close to identical, we merge them. Conversely, if the evidence from mapping indicates separate strains mapping onto the same reference, the reference is split into two sequences for future iterations.

Open in a separate window

Figure 2

Convergence of reconstructed SSU sequences and abundance estimates. (a-d) Algorithm convergence for both the simulated simple microbial community (a, b) and natural community (c, d) is shown. In (a, c), the number of nucleotide (nt) changes made in all reconstructed SSU sequences is plotted for each iteration. In (c, d), each line represents a different reconstructed SSU sequence: the prior probability (abundance estimate) of each SSU sequence is plotted for each iteration. Only SSU sequences with ≥ 1% prior probability at convergence are shown.

Central to our approach to determining both the correct SSU gene sequence and computing abundances of each sequence type is a probabilistic approach that acknowledges that we do not know which reference sequence each read should be mapped to. Especially for communities with many closely related strains, and because of sequencing errors, many assignments of reads to specific reference sequences are uncertain. There is even more uncertainty if the read was generated by a sequence not represented in the reference database. The expectation maximization algorithm [22] provides a computational framework for explicitly dealing with this uncertainty. The actual sequences present in the community are unlikely to be present in the SSU database. Thus, the SSU database sequences serve as candidate initial sequences and the probability that each base is correct changes with each iteration. For each mapping cycle and final sequence determination, the base with the highest probability is chosen at each position.

We use the EM algorithm to alternate between an expectation step (E-step), in which the probability for each read being generated by each reference sequence is computed, and a maximization step (M-step), in which we calculate both (i) the probability values for each base in each reference and (ii) a prior probability that each candidate reference generated any read.

The E-step: computing the probability that a SSU sequence generates a read

In the E-step we construct a distribution representing the probability that each reference sequence generated each read of interest. To compute the probability that a specific reference sequence s generated the observed read r, Pr(s|r), we use Bayes' theorem:

The E-step likelihood

To calculate the likelihood, Pr(s|r), for a given read-reference pair, we make the false but simplifying assumption that each base position in the read is independent. The probability that the reference SSU sequence generated the read is thus a product over all mapped positions (k) of the probabilities of observing each base b_k, given the reference SSU sequence:

To compute Pr(b_k|s), we consider the possibility that the reference base n can be one of A, C, T, or G:

The first term in the summation is the probability of observing read base b_k, given that the reference sequence base at position k is n, and can be calculated with the following:

where p_k is the error probability of the called base b_k in the read, as reported by the base-calling software. We make the simplifying assumption that a base mismatch has an equal probability of being one of the three non-match bases. The indicator variable I is one if and only if the base called in the read matches the mapped base in the reference sequence. The second term in the summation above, Pr(n), is calculated based on the current alignment of reads to the reference sequence from the previous M-step in the EM algorithm (see below).

The M-step: computing SSU sequence and abundance probabilities

In the M-step, the model parameters representing the candidate reference sequences and their abundances are updated, based on the current best estimate of Pr(s|r) from the previous E-step.

Correction of the reference sequence

We can calculate Pr(n) for each base position in each reference sequence. Ignoring sequence quality scores and with each read mapped to only one sequence, a maximum likelihood estimator for the probability of, for example, base A at reference sequence position k would simply be a summation of all A bases in reads mapped at that position divided by the total number of reads mapped at that position. However, because the generating reference sequence for each read is unknown, we instead compute the probability of that base at position k based on the current calculated Pr(s|r) from the previous E-step, as well as quality scores of mapped bases:

where b_{k, j} is the aligned base at position k in read j, I is an indicator variable indicating a match of the aligned base with the base under consideration, n, and P_{k, J} is the error probability of the aligned base k in read j. The consensus sequence chosen for the next round of mapping is simply the sequence of bases with the highest probability at each position (Figure ​(Figure1c1c).

Adjusting reference sequence abundances

In each M-step, we also calculate the prior probabilities (abundances) of each reference sequence, based on the current calculated Pr(s|r) from the previous E-step. Again, if the reference sequence generating each read was known, an estimate of these prior probabilities could be obtained by observing the fraction of reads generated by each reference. However, each read is essentially split among several possible 'read-generating' reference sequences probabilistically from the previous E-step. Thus, we compute:

where J is the total number of reads with mappings.

Algorithm initialization and termination

The EM algorithm performs best and avoids local maxima when initialized with reasonable parameters. Therefore, we initialize EMIRGE by first choosing the single best reference sequence for each read, as reported by the mapping software. If a read maps equally to two or more sequences, a reference sequence is chosen at random. When paired reads are available, the reads are mapped together and the probabilities are computed as if a single longer read was used. We then assume for each read-reference pair and begin an M-step. We terminate the algorithm when no further changes are made to the nucleotide sequence of the reference SSU genes, at which time the sequence abundances have also stabilized (Figure 2a, c).

Assessing algorithm performance on simulated microbial communities

To test the ability of EMIRGE to reconstruct correct SSU gene sequences and abundances from metagenomic data, we first simulated realistic error-containing Illumina reads from the simple microbial community constructed in silico (see Materials and methods). After combining SSU genes from the same organism (generating composite sequences in cases where multicopy genes were not identical), and combining two Lactococcus lactis subspecies, which share near identical SSU sequences (> 99% identity), the nine composite SSU genes range in abundance from 2.3% to 26.1%. To further challenge our algorithm, we also introduced random mutations in 10% of the bases in each SSU sequence in the starting reference database. We processed the simulated reads with EMIRGE and examined the nine most abundant reconstructed SSU sequences that emerged from the starting database of approximately 125,000 sequences (Additional file 1). These nine sequences, which together represented 96% of the summed prior probabilities, accurately recover each of the expected SSU sequences (Figure ​(Figure3a).3a). In addition, the observed abundances match the expected abundances excellently (Pearson ρ = 0.998, P-value = 8.5e-10; Figure ​Figure4).4). Thus, on simulated but realistic data, our method accurately reconstructs both community SSU gene sequences and their abundances.

Open in a separate window

Figure 3

Community composition is captured by correctly reconstructed full-length SSU sequences. (a, b) Phylogenetic trees showing algorithm-reconstructed sequences (black diamonds) and their best blast hits, for both the simulated simple (a) and natural (b) microbial communities. Reconstructed sequences are presented with their (arbitrary) algorithm-assigned identifier and their prior probability, which serves as an abundance estimate, after the final round. All reconstructed sequences match to the expected organism in the simulated community (a), and all but two sequences match to metagenomic contigs assembled from traditional Sanger sequencing in the natural community (b). The two novel sulfobacilli in the natural community are presented with their closest blast hit in GenBank. Units are base substitutions per site, and bootstrap values ≥ 50 are shown at the branches.

Open in a separate window

Figure 4

SSU abundance estimates are accurate. For the nine most abundant reconstructed sequences in the simulated simple community, the final prior probability estimated by EMIRGE is plotted against the expected SSU abundances from the associated community members. The algorithm recovers the expected abundances excellently (Pearson ρ = 0.998, P-value = 8.5e-10).

To gauge how EMIRGE performs on more complex communities, we next simulated reads from a mock community of 67 human gut microbes [28]. The community reconstructed from the full-length SSU sequences reported by EMIRGE was highly similar to the expected community (Figure ​(Figure5;5; Additional file 2). We attempted to quantify this similarity; using the weighted UniFrac statistic [30], the EMIRGE and expected communities were not distinguishable (weighted UniFrac distance = 0.0124, P-value < 0.01). The EMIRGE community was more similar to the expected community when the input reads were longer or sequencing effort was higher. However, even with few paired SSU reads of full length, or with enough reads as short as 36 bp, the community structure could be recovered correctly (Figure ​(Figure6).6). Insert size of the sequencing library had little effect on the ability of EMIRGE to recover the expected community.

Open in a separate window

Figure 5

Accurate SSU sequences and abundance estimates are recovered by EMIRGE for a complex microbial community. Using reads from the complex simulated community, full-length SSU genes reconstructed by EMIRGE with estimated abundances of > 0.5% were aligned and placed in a phylogenetic tree with the expected community members. Estimated EMIRGE sequences and relative abundances (blue) correspond in most cases to expected sequences and expected abundances (red). Grey circles on branches indicate bootstrap values > 80.

Open in a separate window

Figure 6

Effect of sequencing library characteristics on EMIRGE performance. The effects of sequencing effort (x axis), read length, and insert size were evaluated by running EMIRGE on the complex community with varying input. Reconstructed communities were compared to the expected community with the weighted UniFrac distance metric [30]. For the varying insert size experiment, a single sequencing effort was chosen (76-bp read length; 80,000 genomic reads; see Materials and methods).

We also attempted to evaluate how well the short SSU fragments produced by de novo assembly could reconstruct the structure of a microbial community. Assembly produced many more SSU gene fragments than expected for each of the simulated communities, though many of these could be classified to a specific genus (Table ​(Table1).1). We used BLAST [31] to search for the closest known full-length SSU sequence to these short assembled fragments. The best hits were clustered at 97% identity, and each cluster was given an abundance relative to the average k-mer coverage of the matching SSU fragments. The resulting phylogenies produced weighted UniFrac distances to the expected communities that were roughly an order of magnitude larger than the distances produced by EMIRGE (Table ​(Table11).

Simulating short read sequencing reads for simulated microbial communities

For the simple simulated community, we reconstructed reads from the ten evenly distributed genomes used in the in vitro community of Morgan et al. [27]. For the complex community, we first downloaded the SSU sequences for organisms in the 'uneven 1' mock community from Turnbaugh et al. [28]. Because not all of these organisms have available genome sequences, we padded each SSU sequence with 1,000 random bases before simulating reads. This likely makes the de novo assembly problem (though not EMIRGE's method) easier than can be expected in a real community, as there are unlikely to be shared k-mers between taxa in these random padding sequences. We used the wgsim program [43] to simulate 60e6 paired-end, 76-bp, error-free reads from the genomes for the simple community, with an insert length mean of 200 bp and a standard deviation of 25 bp. We simulated reads with varying coverage depth, insert size, and read length for the complex community. Unless otherwise specified, we used a data set with 80,000 paired, 76-bp reads with an insert size mean and standard deviation of 300 ± 30 bp. When simulating other insert sizes, a standard deviation of 10% was used. Many of these simulated reads fell in the padded genome sequence outside of the SSU genes. We assigned quality score vectors to each read by sampling at random from 1 million real quality value vectors from the real microbial community reads, and made mutations in all simulated reads with the probabilities specified by the assigned quality scores. To calculate the expected abundances of SSU rRNA genes in the simple simulated community, we divided the copy number of the SSU gene for any given genome by the sum of copy numbers for all SSU genes in the community. For the complex community, we used the expected abundances given by Turnbaugh et al. [28]. Both simulated datasets are available [44].

Evaluation of SSU gene fragments produced by de novo assembly

Contigs from de novo assembly were searched with BLAST [31] against the Silva SSU database, and contigs with an e-value ≤ 1e-10 were identified as SSU-fragment-containing. These fragments were classified at the genus level using the RDP classifier [32] if the classification had a bootstrap value > 50%. Best BLAST hits to the fragments were clustered at 97% identity using UCLUST [45], and relative abundances for each cluster representative were calculated based on the average k-mer coverage, as reported by Velvet, for each of the fragments hitting a sequence in that cluster. We note that this strategy will not allow the discovery of novel sequences as EMIRGE will; however, for the two simulated datasets all the expected sequences were in the search database.

Implementation details

The reference SSU database was built by first removing sequences < 1,200 bp or > 1,900 bp from the Silva [15] SSU reference 102, and then clustering the remaining sequences at 97% identity with UCLUST [45]. In order to evaluate the ability of EMIRGE to recover novel sequences, for each database sequence, 10% of the sites were chosen at random and mutated to a different base, also chosen at random. For read mapping, we used the short read mapper bowtie [19] in paired-end mode. For initiation, bowtie reported a single best alignment, allowing up to three mismatches in a 20-bp seed, and a maximum sum of quality values in mismatched bases of 300. The minimum and maximum insert size allowed was set to ± 3 times the expected standard deviation around the expected median insert size. For subsequent iterations, bowtie used the same parameters, but reported all mappings in the top strata. Prior to running EMIRGE, reads were trimmed from the 3' end to remove bases with quality scores of 2 or lower, and paired reads were kept if both reads were at least 60 bp long after trimming. For data sets with shorter read lengths, no trimming was done.

At each iteration, SSU sequences were merged into one sequence if the identity of non-gapped positions in a global alignment was greater than 97%. A single SSU sequence (and its prior probability) was divided into two sequences if the second most probable base in more than 4% of all positions had a probability greater than 10%. In this way, sequences that evolved over iterations to be the same were merged, and sequences with evidence from the reads for multiple OTUs were duplicated and allowed to evolve as separate OTUs in future iterations.

EMIRGE was implemented in python. Forty iterations of EMIRGE for the natural community took 6.8 hours on eight Intel Xeon 2.0 GHz cores, with a maximum memory footprint of 900 MB. The read mapping steps accounted for 91% of the run time and is parallelizable. Thus, run time decreases approximately linearly with the number of processor cores available.

Fluorescent in situ hybridization microscopy

Cell fixation and epifluorescent microscopy were performed as described previously [46], except that cells were fixed after thawing from storage at -80°C. We developed a probe for a broader specificity of Sulfobacillus spp. from Bond et al. [47], SUL230 (5'-GGRGCUCGCGGCCCAUUA-3'), and all cells were counterstained with DAPI. Probe stringency was maintained by hybridization at 46°C with 30% formamide.