Genome Components

The HiSeq and MiSeq metagenomes were built using 20 sets of bacterial whole-genome shotgun reads. These reads were found either as part of the GAGE-B project [21 (link)] or in the NCBI Sequence Read Archive. Each metagenome contains sequences from ten genomes (Additional file 1: Table S1). For both the 10,000 and 10 million read samples of each of these metagenomes, 10% of their sequences were selected from each of the ten component genome data sets (i.e., each genome had equal sequence abundance). All sequences were trimmed to remove low quality bases and adapter sequences.
The composition of these two metagenomes poses certain challenges to our classifiers. For example, Pelosinus fermentans, found in our HiSeq metagenome, cannot be correctly identified at the genus level by Kraken (or any of the other previously described classifiers), because there are no Pelosinus genomes in the RefSeq complete genomes database; however, there are seven such genomes in Kraken-GB’s database, including six strains of P. fermentans. Similarly, in our MiSeq metagenome, Proteus vulgaris is often classified incorrectly at the genus level because the only Proteus genome in Kraken’s database is a single Proteus mirabilis genome. Five more Proteus genomes are present in Kraken-GB’s database, allowing Kraken-GB to classify reads better from that genus. In addition, the MiSeq metagenome contains five genomes from the Enterobacteriaceae family (Citrobacter, Enterobacter, Klebsiella, Proteus and Salmonella). The high sequence similarity between the genera in this family can make distinguishing between genera difficult for any classifier.
The simBA-5 metagenome was created by simulating reads from the set of complete bacterial and archaeal genomes in RefSeq. Replicons from those genomes were used if they were associated with a taxon that had an entry associated with the genus rank, resulting in a set of replicons from 607 genera. We then used the Mason read simulator [22 ] with its Illumina model to produce 10 million 100-bp reads from these genomes. First we created simulated genomes for each species, using a SNP rate of 0.1% and an indel rate of 0.1% (both default parameters), from which we generated the reads. For the simulated reads, we multiplied the default mismatch and indel rates by five, resulting in an average mismatch rate of 2% (ranging from 1% at the beginning of reads to 6% at the ends) and an indel rate of 1% (0.5% insertion probability and 0.5% deletion probability). For the simBA-5 metagenome, the 10,000 read set was generated from a random sample of the 10 million read set.

Given a pairwise alignment without gaps in genomes and , we compute a pairwise substitution score using a substitution matrix, which defaults to the HOXD matrix [40] (link). The HOXD matrix appears to discriminate well between homologous and unrelated sequence in a variety of organisms, even at high levels of sequence divergence.
The substitution matrix score quantifies the log-odds ratio that a pair of nucleotides share common ancestry, but does not account for the inherent repetitive nature of genomic sequence. Our desire to discriminate between alignment anchors that suggest positional homology and alignments of regions with random similarity or paralogy requires that we somehow consider repetitive genomic sequence in our anchoring score [41] .
We combine the traditional substitution score for a pair of nucleotides with an adjustment for the multiplicity of -mer seeds at the aligned positions:
where is the number of occurrences of the spaced seed pattern that matches the subsequence of at . The product of approximates the number of possible ways that sites in and with the same seeded -mers as and could be combined. For example, consider a repeat element present in both genomes with copy number in genome and copy number in . There are possible pairs of repeats. When a pair of nucleotides in a repeat element have a positive substitution score, the product down-weights the score.
In summary, this scoring scheme assigns high scores to well-conserved regions that are unique in each genome and does not consider gap penalties.

We have expanded the Reference Gene Catalog^{8 (link)} to include genetic elements related to stress response and virulence genes; these expansions can be visualized in the Reference Gene Catalog Browser (https://www.ncbi.nlm.nih.gov/pathogens/refgene/). One reason we expanded AMRFinderPlus is to understand the linkages between AMR genes and stress response and virulence genes in food-borne pathogens; thus, the stress response and virulence genes included in the Reference Gene Catalog are composed primarily of E. coli-related genes derived primarily from González-Escalona et al.^{23 (link)} as well as BacMet^{24 (link)}, but also have been supplemented by manual curation efforts for other taxa. Stx gene nomenclature adopts the system of Scheutz et al.^{25 (link)} and the intimin (eae) gene nomenclature uses existing designations in the literature^{26 (link),27 (link)}. Genes are incorporated only if there is literature supporting the function of that protein or closely related sequences that meet the identification criteria. As a major focus of our work is to improve NCBI’s Pathogen Detection system^{16 (link)}, we excluded genes that belonged to organisms not deemed clinically relevant. To remove ‘housekeeping’ proteins that were universally found in one or more taxa in the Pathogen Detection system, sequences were not included if they were found at a frequency of greater than 95% in a survey of 58,531 RefSeq bacterial assemblies belonging to any of the following species: Acinetobacter, Campylobacter, Citrobacter, Enterococcus, Enterobacter, Escherichia/Shigella, Klebsiella, Listeria, Salmonella, Staphylococcus, Pseudomonas, and Vibrio. If genes of particular interest in foodborne pathogens exceeded this threshold, they were excluded in the taxa where they appear to be nearly universal (see “Identifying genomic elements” below). In addition, genes with misidentified functions, such as copper-binding proteins that use copper as a co-factor yet do not confer resistance to copper, also were excluded. As we continue to expand the database, we use similar criteria when adding genes.

In building the training set, we used an advanced method based on Markovian Jensen–Shannon divergence (MJSD) to obtain the core (native) components of all available prokaryotic genomes to ensure the most balanced representation was used in our regression. We were able to significantly reduce the runtime of genome segmentation and clustering algorithm, as implemented in IslandCafe [27 (link)], by introducing a reverse-calculation step during recursive segmentation. MJSD, entropy, and statistical significance were calculated as described in [27 (link)]. Specifically, the information content of a genome sequence, quantified by the entropy function for probability distribution p_i, is obtained as, $H^m\left(p_i\right)=-\sum _{w}P\left(w\right)\sum _{x\in \mathcal{A}}P(x|w){log}_{2}P(x|w)$ , where P(x|w) is the probability of nucleotide x given the preceding oligonucleotide w of length m (m defines the model order, is set to 2 in IslandCafe) and P(w) is the probability of oligonucleotide w. A genome is initially segmented by iterating the computation of entropy and thus MJSD at each position along the genome and identifying the location of highest MJSD of (user-defined) significance in the genome. This process is then iterated for the resulting genomic segments.

To compare the editability of different genomic elements, including the protein-coding gene related elements (5′-UTR, CDS, intron and 3′-UTR) and the repeat-associated elements (SINE, LINE, LTR, DNA transposon, Helitron, tandem repeat and other unclassified repeat loci), we calculated the A-to-I editing density for each type of genomic element by counting the number of A-to-I editing sites located in this element type, out of the total number of transcribed adenosines (RNA depth ≥ 2X) from this element type. The editing density of each element type was first calculated for each sample of a species separately, then the mean editing density across samples was calculated as the representative value for a species (Figure 2C).
We also calculated the editing-level-weighted editing densities for each element type (Figure 1, Figure 2, Figure 3, Figure 4, Figure 5S3C and S3D). To do so, an editing site with for example an editing level of 0.1, would be regarded as 0.1 editing site instead of 1 editing site, when counting the number of editing sites for an element type. Only editing sites and transcribed adenosines with RNA depth ≥10X were used in the weighted analysis.

The Repeat Masker (v4.1.4) (33 ) with -pa 16 -qq options was used to quantify repeat elements from reference genomes of various species. RMBlast (v2.11.0) was used as the repetitive sequence search algorithm, and the search was based on the Dfam (v3.6) database (34 (link)). In addition, TRF (v4.09) (35 (link)) was used to find tandem repeat sequences.