Mumame is implemented in Perl and consists of two commands: mumame, which performs read alignment to a database of mutations, and mumame_build which builds the database for the former command. The mumame_build command takes a FASTA sequence file and a list of mutations (CSV format) as input. For each entry in the mutation list, it finds the corresponding sequence(s) in the FASTA file, either by sequence identifier or by CARD ARO accessions (Jia et al. 2017). It then excerpts a number of residues upstream and downstream of the mutation position (by default 20 residues for proteins and 55 for nucleotide sequences) and creates one wildtype version and one mutated version of the sequence excerpt with unique sequence IDs. For cases where multiple mutations can occur close to each other on the same sequence, the software attempts to create all possible combinations of mutations (if memory permits; in some situations this is not possible because the number of combinations increases exponentially). The software tool also generates a mapping file between sequence IDs in the database and mutation information from the list.

The main mumame command takes any number of input files containing DNA sequence reads in FASTA or FASTQ format and aligns those against the Mumame database using Usearch (Edgar 2010). For this read alignment, the software runs Usearch in usearch_global mode with target coverage set to 0.55 (by default; any value ≥0.51 should be feasible for target coverage). The output is then matched to the wildtype or mutation information in the Mumame database, and data is collected for each input file and combined into a single output table. The software uses a best-hit strategy, with hits sorted by identity to the reference sequence. As every sequence in the database is present in two versions, one variant with and one without the point mutations, the top hit will always discriminate between the two variants. In addition, as only the immediate region around the mutation site is included in the database, there cannot be any spurious hits adding to wildtype or mutated variants, providing Mumame with the maximum possible degree of precision. The choice of Usearch for read alignment was made because it constitutes the most versatile sequence search software in that it is, unlike most other sequence search tools for large-scale data, able to align nucleotide sequences to both protein and nucleotide sequence databases. However, the software design of Mumame is fairly flexible, allowing implementation of other software for the sequence search process with relatively small effort.

The main output of Mumame is a file with the suffix “.table.txt”. This file contains the reads from each library aligned to the mutation database, with mutation counts in the first set of columns and wildtype counts in the second set of columns. The last line of this file contains the total number of reads in each library, which can be used, e.g., for normalization purposes. The software also saves the output from the Usearch run and, optionally, the read alignments to the database. The output table generated by Mumame can be analyzed using the R script (R Core Team 2016) supplied with the Mumame package. The script reads the read counts for all mutation positions detected, both for wildtype and mutated sequences. The script also takes into account the total library sizes, either explicitly as normalization factors in the count model or implicitly in the proportional model as total library size cancels out (see below). The script assesses if there are significantly different proportions of mutations between different sample groups through a generalized linear model. Alternatively, an overdispersed Poisson generalized linear model (GLM) accounting for the discrete nature of the data and the differences in sequencing depth can be used (Jonsson et al. 2016; Bengtsson-Palme et al. 2017). These two tests were selected after investigating the performance of the two GLM tests, the Student’s t-test on the mutation proportions and the Chi-Square test on the total counts in a simulation study. We simulated a set of different numbers of replicates, effect sizes, sequencing depths and average gene abundances (Suppl. material 1: Table S1), 100 times for each combination of conditions and assumed that counts were Poisson distributed. The simulations showed that the two GLM models overall perform better in terms of detecting significances, particularly when the sequencing depths and effect sizes are small, while generating the expected proportions of false positive detections. The Poisson model is preferable when the number of counts for a targeted gene is low in all sample groups. However, this model performs poorly when estimating effect sizes with small numbers of counts. This is due to situations where one group has zero counts for all replicates. In practice, this means that even if a result is deemed significant by the model, the estimated effect size should only be considered an indication of the directionality of the effect when counts are all zeros in one group.

The Mumame software is freely available (http://microbiology.se/software/mumame or https://github.com/bengtssonpalme/mumame) and can also be installed via Conda, using the command “conda install -c bengtssonpalme mumame”.

To quantify the abundances of fluoroquinolone resistance mutations in the gyrA and parC genes (Johnning et al. 2015b), we downloaded the CARD database on 2018–05–24 (Jia et al. 2017). We extracted all mutation information regarding the gyrA and parC genes from the “snps.txt” file and created a new file with that information. We then created a new Mumame database using mumame_build with default options. That database was used to align all the reads from the samples generated by Kraupner et al. (2018; data provided by courtesy of Stefan Ebmeyer and Joakim Larsson) to the database using Mumame in the Usearch mode (Edgar 2010) and the following options “-d gyrA_parC -c 0.95”. This study exposed aquatic bacterial communities in an aquarium system to different concentrations of ciprofloxacin (0, 0.1, 1 and 10 µg/L) in triplicates. The communities were then subjected to both amplicon sequencing of the target genes for ciprofloxacin (56,394–290,391 reads per sample after preprocessing) and shotgun metagenomics (109–190 million reads per sample). We analyzed both the shotgun metagenomics data as well as the amplicon sequences derived specifically from Enterobacteriaceae gyrA and parC genes. Prior to this sequence similarity search raw reads were quality filtered using Trim Galore! (Babraham Bioinformatics 2012) with the settings “-e 0.1 -q 28 -O 1”. We then used the R script provided with the Mumame software to compare the numbers of matches to mutated and wildtype sequences in the database. The same database and method combination was used to quantify fluoroquinolone resistance mutations in sequence data from an Indian lake exposed to ciprofloxacin pollution (Bengtsson-Palme et al. 2014; ENA project PRJEB6102; https://www.ebi.ac.uk/ena/data/view/PRJEB6102), as well as in an Indian river upstream and downstream of a wastewater treatment plant processing pharmaceutical waste (Kristiansson et al. 2011; Pal et al. 2016; MG-RAST project 18323; https://www.mg-rast.org/linkin.cgi?project=mgp18323). These samples were preprocessed in the same way as in the Indian lake study (Bengtsson-Palme et al. 2014). These samples were taken as part of a survey of environments close to a pharmaceutical production waste treatment facility and consist of two samples of river sediment from upstream of the treatment plant, one by the outlet, three taken downstream, and one from a nearby lake. These were subjected to shotgun metagenomic sequencing, generating 14.5–37.3 million reads per sample.

Finally, we investigated data from the experiment by Lundström et al. (2016; ENA project PRJEB11402; https://www.ebi.ac.uk/ena/data/view/PRJEB11402), in which aquatic bacterial communities were investigated in an aquarium system exposed to different concentrations of tetracycline (0, 0.1, 1, 10, 100 and 1000 µg/L). The samples were sequenced using shotgun metagenomics to a depth of 122–273 million reads per sample. The experiment used two replicates per concentration, meaning that we can only test for statistically significant trends in the data rather than differences between specific concentrations. To quantify resistance mutations in the sequence data, we created a Mumame database for tetracycline resistance mutations in the 23S rRNA gene targeted by tetracycline. We extracted the mutational information related to tetracycline from the CARD “snps.txt” file and then built the database using mumame_build with the additional option “-n”. We then aligned all reads from the Lundström et al. (2016) data to the Mumame database using the options “-d Tet -c 0.95 -n”. Reads were quality filtered and statistical differences were assessed as above.

To assess the limitations of the method in terms of sequencing depth, the samples from the highest and lowest ciprofloxacin concentrations generated by Kraupner et al. (2018; 10 µg/L and 0 µg/L, respectively) were downsampled to 1, 5, 10, 20, 30, 40, and 50 million reads. Thereafter, the reads from the downsampled libraries were aligned to the fluoroquinolone resistance mutation database using Mumame as above. Statistical differences were assessed at all simulated sequencing depths and average effect sizes calculated for the significantly altered genes.

As a proof-of-concept that our method to identify point mutations in metagenomic sequence data is functional, we used Mumame to quantify the mutations in amplicon data from the gyrA and parC genes. These genes are targets of fluoroquinolone antibiotics, and often acquire resistance mutations attaining high levels of resistance. We quantified such mutations in an amplicon data set specifically targeting these two genes in Escherichia coli. This data set derives from an exposure study with increasing ciprofloxacin concentrations, and enrichments of mutations in the classical fluoroquinolone resistance determining positions S83 and D87 (gyrA) and S80 and E84 (parC) have previously been verified using other bioinformatic methods (Kraupner et al. 2018). This data set, therefore, serves as an ideal positive control for our novel method. We found that Mumame was able to identify the difference between the highest concentration (10 µg/L) and the lower ones reported in the original study (Fig. 1). However, Mumame only reported an average frequency of mutations of around 11–12% for gyrA mutations (Fig. 1A), while the original paper found frequencies of 60–85% (S83) and 30–40% (D87). The A67 position was not quantified in the original paper. The exact reason for the discrepancies is unknown, but it is likely caused by a taxonomic filtration step that selects for E. coli reads used in the Kraupner et al. (2018) study, while Mumame does not perform prior filtering. The decision to exclude filtering was made in order to mimic a situation with true metagenomic data where several target species may co-exist. For parC, Mumame only quantified the S80 position (Fig. 1C), because the E84 mutations were not included in the version of the CARD database used for this study. For position S80, Mumame identified around 35% mutated sequences at the highest concentration of ciprofloxacin, while the original study reported around 50%. When a similar E. coli filtering step was introduced in our analysis, the proportions were much closer to those in the original paper (data not shown).

We next evaluated the performance of Mumame on the real shotgun data that was generated from the same samples as the amplicon libraries. Ideally, this analysis should generate virtually the same result as the amplicon analysis. Indeed, we found similar results for the A67 and S83 gyrA mutations (Fig. 1B). For the D87 mutation, the frequencies were much lower than for the other two mutations, albeit still significantly larger than at the lower concentrations (p < 0.01). For the parC gene, the shotgun metagenomic analysis had large variability within the sample groups, which prevented any statistically significant results (Fig. 1D). This is surprising, given that the total number of reads detected for both the mutant and wildtype variants were higher for the parC gene than for gyrA. Thus, the large degree of variability could potentially be due to features of the parC gene. For example, it is possible that this gene has higher similarity to closely related species, leading to that species replacement could influence the results. Taken together, these results indicate the high noise levels present for individual gene variants even in deeply sequenced shotgun metagenomes from controlled exposure studies.

Figure 1.

Fluoroquinolone resistance mutations in ciprofloxacin-exposed bacterial communities. Total mutation frequencies quantified using Mumame for three known mutations conferring resistance to fluoroquinolone in the E. coli gyrA gene based on amplicon sequencing (A) and shotgun metagenomic data (B) from the same samples. Corresponding data for the S80 mutation in parC is shown in (C) for amplicon data, and (D) for shotgun data.

Noting the much more instable levels of mutations in the shotgun metagenomes, we next investigated the effects of sequencing depth on the ability of our method to detect significantly altered mutation frequencies. For this analysis, we used downsampled data from the shotgun metagenomic library of the ciprofloxacin exposure study (Fig. 2). As expected, we found that the number of significantly altered mutation frequencies detected increased with larger sequencing depth (Fig. 2A). In addition, the average effect size of the significant mutations became gradually lower with larger sequence depth, also in accordance with expectations (Fig. 2B). Importantly, the average effect size of detectable mutation frequency differences seems to decrease linearly with sequencing depth. This means that we can calculate an expected detection limit for the method given the characteristics of the data and experimental setup. At 10 million reads, we expect that the proportion of reads with mutation must be 30–40% higher in the exposed sample in order for it to be detected as significant. This proportion decreases, on average, to 10% for 50 million reads (Fig. 2B). These numbers also depend on other factors, such as the number of replicates per treatment, but nevertheless they can be used as an approximation to aid the design of metagenomic studies or to interpret non-significant results derived from Mumame analyses.

Figure 2.

Influence of sequencing depth on detected mutations and their effect sizes. Relationship between the number of investigated reads and number of mutations with significantly altered frequencies (A) and the average effect size for those mutations (B); as assessed using Mumame on shotgun metagenomic data from a ciprofloxacin exposure experiment.

After validating the method and testing the limit of detection, we used Mumame to quantify resistance mutations in a similar controlled aquarium setup under exposure to the antibiotic tetracycline (Lundström et al. 2016). In this study, no amplicon sequencing of the target gene for tetracycline (23S rRNA) was performed, and thus there was no a priori true result to which we could compare. While Mumame was able to successfully detect tetracycline resistance mutations in the data, we somewhat surprisingly found no enrichment of such mutations in this data (Fig. 3). Notably, this result was obtained despite a very high sequencing depth (on average 181,595,072 paired-end sequences per library). Obtaining a negative result at this sequencing depth suggests that there actually is no enrichment of known E. coli resistance mutations in the samples.

Figure 3.

Resistance mutations in tetracycline-exposed bacterial communities. Frequencies of E. coli tetracycline resistance mutations at exposure to different concentrations of tetracycline, based on shotgun metagenomic data.

As a final investigation of the performance of the method, we also let Mumame quantify the fluoroquinolone resistance mutations in river and lake sediments polluted by antibiotic manufacturing waste, primarily ciprofloxacin (Kristiansson et al. 2011; Bengtsson-Palme et al. 2014; Pal et al. 2016). These libraries are fairly old and were not as deeply sequenced as the other data sets we investigated. While the experimental setup of these studies in terms of number of replicate samples per geographical location does not allow for proper statistical testing, we found overall more fluoroquinolone resistance mutations downstream of the pollution source, at least for the E. coli gyrA and parC genes (Fig. 4). We also detected a few such mutations in other species, but the counts of those were low and the results largely non-informative due to the small number of detections per mutation (Suppl. material 2: Fig. S1). This serves as an example of that mutations can be detected also in shallowly sequenced metagenomic data, but that without proper experimental design, interpretation of the results is difficult or even impossible.

Figure 4.

Resistance mutations in antibiotic-polluted sediments. Relative frequency of gyrA and parC sequences with resistance mutations in samples taken downstream, at, or upstream of the pharmaceutical production wastewater treatment plant, as well as in a lake polluted by dumping of pharmaceutical production waste. The numbers at the top of the bars show the total number of gyrA/parC sequences (wildtype or mutated) identified in each sample.