GWAMAR: Genome-wide assessment of mutations associated with drug resistance in bacteria

Problem setting

We consider a set $\mathcal{S}$ of closely related bacterial genomes. Typically, this is a set of strains within the same species of bacteria.

Then, we represent the available drug resistance information as a set of drug-resistance profiles $\mathcal{R}$, where each drug resistance profile $r\in \mathcal{R}$ is represented as a vector:

$r:\mathcal{S}\to \left\{S,I,R,?\right\}$.

Here, S, I, R denote that a given strain is known to be drug susceptible, intermediate-resistant, or resistant, respectively. Using question mark ${}^{‵}{?}^{\prime }$ we indicate that the drug resistance status of a strain is unknown. We call a drug resistance profile complete if it does not contain question marks.

Analogously, we represent the genotype data as a set of mutations $\mathcal{M}$, where each mutation $m\in \mathcal{M}$ is represented as a triple (g, p, v), where g,p,v denote the gene family identifier, the position of the mutation in its corresponding multiple alignment, and the mutation profile, respectively. The mutation profile is represented as a vector:

$v:\mathcal{S}\to \sum \cup \left\{?\right\}$.

Here Σ denotes the set of amino acides (for protein-coding genes) or nucleotides (for promoters and rRNA coding genes). We also assume that Σ contains the gap '-' symbol. Using question mark ${}^{‵}{?}^{\prime }$ we indicate that the corresponding amino acid or nucleotide is unknown. Analogously, we call a mutation profile complete if it does not contain question marks.

Analogous to mutation profiles, we call a binary mutation profile complete if it does not contain question marks.

Finally, we define the aim of the tool as: To produce an ordered list of associations between the phenotype and genotype data (represented as drug resistance and mutation profiles) such that the top-scored associations are the most likely to be real.

The pipeline of GWAMAR

Figure 1 illustrates the overall design of the tool. Data preprocessing for GWAMAR consists of several steps which may be performed by our previously developed tool, eCAMBer [14]. These preprocessing steps comprise:

The input drug-resistance profiles, typically, are collected from the articles or databases which provide drug-resistance information for the strains of interest. The set of identified point mutations, the set of drug-resistance profiles and and the phylogenetic tree constitute the input for GWAMAR.

In the next step, GWAMAR computes binary mutation profiles for each mutation profile. This step significantly reduces the number of genetic profiles to be scored. Finally, GWAMAR computes several statistical scores to associate drug-resistance profiles to the mutation profiles, including mutual information, odds ratio, hypergeometric test, weighted support (which is our previously published approach [13]), and the tree-generalized hypergeometric score (our new approach here).

Tree-generalized hypergeometric score

As a part of this work we also introduce a new association score, called tree-generalized hypergeometric score (TGH ). This score is a modification of the CCTSWEEP score introduced in the paper [16]. In this section, we consider a subset of strains S for which a given drug-resistance profile r and a binary mutation profile b are complete, i.e. do not contain question marks. Moreover, we assume that r does not contain any intermediate-resistant strains. In all our computational experiments we transform the intermediate-resistant strains into resistant strains.

In order to present the formal definition of TGH, we first define an auxiliary concept called coloring. For a given tree T , we call a subset c of its nodes a coloring, if it satisfies the following two conditions:

• each path from a leaf to the root contains at most one node from c,

• each internal node in c has a sibling node which does not belong to c.

Here we also introduce a function L which, for each node ω, returns the set of descendants of the node, including the node itself. We say these nodes are visible from ω. Additionally, the function L applied to a coloring c returns the union of all nodes visible from nodes in c.

Let C _T denote the set of all colorings of T . Then, for each complete drug-resistance profile r there exists exactly one coloring ĉ such that the set of leaves visible from ĉ equals the set of drug-resistant nodes in r. We say this coloring is induced by the drug-resistance profile. Analogously, for each complete binary mutation profile b there is exactly one induced coloring $\overline{c}$.

Intuitively, the coloring induced by a given complete drug-resistance profile will contain the set of nodes in which drug-resistance was acquired (assuming a model in which drug-resistance cannot be reversed). Analogously, the coloring induced by a given binary mutation profile will contain the set of nodes in which the mutation was acquired.

Figure 2 (A) presents an example of colorings induced by a given drug-resistance profile (large red nodes) and a given binary mutation profile (small orange nodes) for a flat tree. In that situation the colorings may be interpreted as independent drawing of balls as in the standard hypergeometric distribution model. Knowing this property of TGH we proposed its name as it generalizes the standard hypergeometric test in the case of a flat tree. Figure 2 (B) presents another example of colorings induced by the same pair of profiles, but for a tree which is not flat. In this model the dependencies between different strains are captured by the topology of the tree.

Let us now assume we want to compute the TGH score for a pair of complete drug-resistance profile r and complete binary mutation profile b. Let us additionally assume that the size of coloring $\overline{c}$ induced by b equals n. Morover, let the number of nodes in coloring $\overline{c}$ visible from the coloring ĉ equals k. This value can be interpreted as the number of times the considered mutation was acquired not earlier than the resistance was acquired.

Now, let V _T (n) denote the number of colorings of size n:

$V_T\left(n\right)=\#\left\{c\in C_T:\left|c\right|=n\right\}$

V _T (n) may be interpreted as the total number of binary mutation profiles for which the induced coloring is of the same size as for $\overline{c}$.

Then, let B _{T, ĉ} (k, n) denote the number of colorings of size n, such that exactly k nodes of that coloring are visible from nodes of coloring ĉ.

$B_{T,ĉ}\left(k,n\right)=\#\left\{c\in C_T:|L\left(ĉ\right)\cap c|=k\mathsf{\text{and}}\left|c\right|=n\right\}$

Here, the value B _{T, ĉ} (k, n) may be interpreted as the number of binary mutation profiles such that their induced coloring has n elements, out of which k is visible from the nodes in ĉ.

Finally, for the complete drug-resistance profile r and complete binary mutation profile b, which induce colorings ĉ and $\overline{c}$, respectively, we define the TGH score as follows:

$H_T\left(r,b\right)=-\mathsf{\text{log}}\left(\frac{{\sum }_{i=k}^n{B}_{T,ĉ}\left(i,n\right)}{V_T\left(n\right)}\right)$

Here, we take the negative logarithm to have consistent property for all considered scoring methods, such that the higher the score the more likely drug-resistance profile r is associated with binary mutation profile b.

Time complexity

Let D denote the number of drug-resistance profiles considered. Additionally, let N denote the number of considered strains and M denote the number of binary mutation profiles. Finally, let K denote the maximal number of children of an internal node in the tree. Then, the time complexity of the algorithms we implemented to compute the hypergeometic score, the mutual information, odds ratio, and weighted support is O(D N M ). In order to compute TGH, we implement a dynamic programing algorithm which computes the values B _ω,ĉ (k, n) for each internal node ω, k and n. This strategy gives an algorithm with complexity O(D·N ^{K− 1}·N² +D·N·M ). For the brevity of the presentation we skip details of the algorithm.

Mtu173 case study

The first case study is based on the set of 173 fully sequenced strains of M. tuberculosis with publicly available data.

The preprocessing steps of preparing the genotype data were performed using eCAMBer, our previously published tool to support comparative analysis of multiple bacterial strains [13].

In particular, first, we used eCAMBer to download the genome sequences and annotations from the PATRIC database [15]. Next, we applied eCAMBer to unify the genome annotations of protein-coding genes and to identify the clusters of genes with high sequence similarity. Then, for the subset of 4379 such identified gene clusters with genes present in at least 90% of the strains, we computed multiple alignments using MUSCLE [17]. The multiple alignments were computed for amino-acid sequences of protein-coding genes, as well as nucleotide sequences of their promoter regions, and rRNA genes. Finally, based on the computed multiple alignments, we identified 118913 non-synonymous point mutations. The set of identified point mutations constituted the input genotype data for GWAMAR.

The input phenotype data was collected from over 20 publications issued together with the fully sequenced genomes. Additional file 1 shows the drug-resistance status of the strains retrieved from literature together with references. Based on the drug-resistance information collected for Ciprofloxacin and Ofloxacin, we introduced a new drug-resistance profile for Fluoroquinolones. A strain is categorized as susceptible to Fluoroquinolones if it was susceptible to at least one of the drugs, but not resistant to any of them. Similarly, a strain was categorized as resistant to Fluoroquinolones if it was resistant to at least one of the drugs, but not susceptible to any of them. If none of the cases was satisfied for a strain, then the drug-resistance status of the strain was categorized as unknown. We restrict further analysis to the set of six drugs or drug families: Fluoroquinolones, Ethambutol, Isoniazid, Pyrazinamide, Rifampicin and Streptomycin.

The input phylogenetic tree was reconstructed using the maximum-likelihood approach implemented in the PHYLIP package. As an input for the tool we used the set of all the identified mutations concatenated into one multiple alignment file. Additional file 2 presents the reconstructed phylogenetic tree.

Literature search did not provide us any additional support for the remaining three mutations (A505T in embB, D869G in embB and G1108C in rrs), which haven't been reported in TBDReaMDB. These mutations may potentially be false positives or real drug-resistance-associated mutations.

Mtu_broad case study

Similar to the previous case study, based on the drug-resistance information available in the database for Ciprofloxacin, Ofloxacin, Levofloxacin and Moxifloxacin, we introduced a new drug-resistance profile for Fluoroquinolones. A strain was categorized as susceptible to Fluoroquinolones if it was susceptible to at least one of the drugs, but not resistant to any of them. Similarly, a strain was categorized as resistant to Fluoroquinolones if it was resistant to at least one of the drugs, but not susceptible to any of them. If none of the cases was satisfied for a strain, then the drug-resistance status of the strain was categorized as unknown. We restrict further analysis to the set of six drugs or drug families: Fluoroquinolones, Ethambutol, Isoniazid, Pyrazinamide, Rifampicin and Streptomycin.

Similarly, as in the previous case study, the phylogenetic tree was reconstructed using the maximum-likelihood approach implemented in the PHYLIP package. As an input for the tool we used the set of all available mutations concatenated into one multiple alignment file.

A closer look at the mutation D89G/N, which is not present in TBDReaMDB, reveals that the mutation has recently been associated with resistance to Fluoroquinolones [20].

The set of associations provides some additional support for the four mutations which were present in TBDReaMDB, but not categorized as high-confidence.

Assessment of accuracy

Here we use the two datasets described above to assess the accuracy of the various association scores, viz: mutual information, odds ratio, hypergeometric score, weighted support and TGH. The CCTSWEEP software might contain some bugs and did not produce correct output. Its authors had not responded to our queries in time. So we omitted it from our experiments.

As our gold standard we use the set of 88 drug-resistance associations classified as high-confidence in the Tuberculosis Drug Resistance Mutation Database (TBDReaMDB) [5].

In both case studies, as the set of positives, we assume the subset of high-confidence mutations present in TBDReaMDB, which are also present in the genotype data. This is the set of mutations which may be potentially detected using the available datasets. There are 39 and 75 of such associations for the mtu173 and mtu_broad datasets, respectively. The set of negatives is constructed by random sampling from the whole set of identified putative associations except for the associations which are classified as positives. The number of sampled negatives equals the total number of mutations present in genes which has at least one high-confidence mutation. We use this approach in order to significantly reduce the probability of classifying as negatives associations which are real but not present in the database.

Figure 3 presents precision and recall curves for different association scores. Consistently, results on both datasets suggest that tree-aware scores outperform tree-ignorant scores.

Rifampicin resistance and compensatory mutations

It is commonly known that Rifampicin resistance in M. tuberculosis gets acquired by point mutations within the RRDR region which corresponds to the Rifampicin binding spot [21]. However, these mutations are often associated with deleterious effect on bacteria fitness [22]. This effect may be potentially reversed by compensatory mutations. Recently, three new papers have been released, focusing on identifying putative compensatory mutations within rpoA, rpoB and rpoC genes [23–25].

Figure 4 presents the distribution of the putative compensatory mutations identified in these recent studies, put together with mutations identified based on our two case studies.

Interestingly, several mutations identified by our approach have also been reported in at least one of the papers [23–25].

Additional file 3 contains the detailed table with the complete list of putative compensatory mutations.