BACKGROUND
rifampin, one of the principal drugs used in tuberculosis treatment, is a semi-synthetic antibiotic that inhibits transcription by preventing rna synthesis.
isolates of mycobacterium tuberculosis resistant to rifampin occur at low to moderate frequencies in many regions of the world  <cit> .
mutations in rpob, the gene encoding the β subunit of dna-dependent rna polymerase, are associated with rifampin resistance.
in the laboratory, drug resistance is quantified in terms of minimum inhibitory concentration , which is defined as the minimum concentration of the antibiotic in a given culture medium below which bacterial growth is not inhibited.
several studies have been conducted where mic of rifampin has been measured and partial dna sequences have been determined for rpob in different isolates of m. tuberculosis  <cit> .
however, no model has been constructed to predict rifampin resistance based on sequence information alone.
such a model might provide the basis for quantifying rifampin resistance status based exclusively on dna sequence data and thus eliminate the requirements for time consuming culturing and antibiotic testing of clinical isolates.
tree-based statistical methods  have generated very accurate models relating amino acid sequence of short  peptides to their binding by major histocompatibility complex  class i molecules with higher accuracy than artificial neural networks  <cit> .
both tree-based models and aggregation of such models through random forests  have proven to be quite successful in other problems involving sequence data as covariates such as hiv- <dig> replication capacity  <cit>  and cytidine to uridine rna editing in plant mitochondria  <cit> .
the success of tree-based statistical models and random forests in these problems involving covariates derived from sequence data motivated our application of these models to the problem of rifampin resistance in m. tuberculosis.
the response variable is a set of continuously distributed values for mic, which makes the problem one of regression.
these data are used to answer the following questions: what proportion of the variance in mic is attributable to sequence differences in positions 511– <dig> of the β subunit of rna polymerase of m. tuberculosis?
what particular positions, and what distribution of amino acids at those positions, are associated with most of the variance in mic?
alternatively, the response variable could be cast in discrete terms: resistant or susceptible.
this is possible by assuming a threshold value for mic above which an isolate is considered resistant to rifampin.
among the specific questions we can answer with such a model are the following: what particular positions, and what distribution of amino acids at those positions, allow for distinguishing rifampin-susceptible and rifampin-resistant isolates of m. tuberculosis?
what is the misclassification error rate associated with susceptibility prediction for these data?
we address these questions and evaluate the ability to predict mic from protein sequence data  using tree-based regression and classification methods.
we find that these methods generate highly accurate models of rifampin resistance.
RESULTS
the data set used in the study consists of  <dig> observations with  <dig> distinct genotype-phenotype combinations .
the most frequent combination has  <dig> occurrences, and there are  <dig> unique  combinations.
mic for rifampin varies from  <dig>  μg/ml to >  <dig> μg/ml.
the  <dig> sequences are distributed among  <dig> genotypes,  <dig> of which occur uniquely in the data set.
the plurality genotype is represented by  <dig> samples;  <dig> samples differ from the plurality by one amino acid; and the remaining  <dig> samples differ from the plurality by two amino acids.
some genotypes defined by the partial sequence of rpob are associated with several different phenotypes .
also, some genotypic states are associated with large effects, while some have little or no effect on mic phenotype.
finally, some changes in mic are not associated with changes in the sequence region examined.
these genotype data are short  partial sequences of a single gene, and thus they may not contain all phenotypically relevant genetic information.
indeed, there is evidence that amino acid changes outside of the examined region are associated with changes in mic for rifampin  <cit> .
additionally, the sample size is small , which will decrease power.
nonetheless these data are typical of studies surveying the genetic variation associated with antibiotic resistance and of genotype-phenotype data in general.
thus they make an appropriate subject of investigation.
regression tree analysis
the regression tree for the relationship of rpob amino acid sequence and mic has two splits defining three terminal nodes .
at each node in the tree, the mic prediction given  is the mean of all isolates at that node.
the first split of the topmost node  consists of the entire sample and is based on the amino acid at position  <dig>  with those sequences having serine  going to the left child node, and those having leucine  or tryptophan  going to the right child node.
the best split for each node is that which gives the largest decrease in the error.
here error is measured as the deviance, which for a continuous variable is a constant multiple of the residual sums-of-squares.
reported values were determined using 10-fold cross-validation.
moving down the tree the error decreases, as the sum of the deviance for each pair of child nodes is less than the deviance of the parent node.
given the hierarchical nature of trees and the criterion used to choose splits, the first split, that based on position  <dig>  explains the highest proportion of the overall phenotypic variance.
this bisection of the data results in groups that differ one hundredfold in mean mic .
the subsequent partition based on position  <dig>  the most variable in this region, results in a > 354-fold difference in mic.
the proportion of the variance in mic explained across all splits involving the two amino acid positions  is  <dig> .
all proportions of variances explained by the model as reported here are those estimated through cross-validation and are not based on re-substitution, and thus represent appropriately conservative estimates.
classification tree analysis
from a clinical perspective it may be most relevant to consider the level of drug resistance as a two-state categorical variable  rather than as a continuously distributed variable.
in clinical practice, if an isolate of m. tuberculosis is determined to be rifampin resistant then rifampin is replaced with another antibiotic.
although blood serum concentration of rifampin reaches levels of  <dig> –  <dig> μg/ml about  <dig>  –  <dig> hours after ingestion  <cit> , a clinically relevant mic value for dichotomizing the mic values would be lower than this peak.
we conservatively define mic values ≤  <dig> μg/ml as susceptible and values >  <dig> μg/ml as resistant, a definition consistent with conventional standards  <cit> .
with this dichotomization we can explore the use of tree-based statistical classification to predict rifampin resistance in a way that is more relevant to clinical practice.
the predictor variables are again the unordered categorical designations of amino acids at polymorphic positions.
the classification tree for these data  has two splits based on two of the  <dig> variable amino acid positions.
at each node in the tree, the prediction of rifampin susceptibility status  is given for all isolates at that node.
the first split is based on position 531; those isolates with serine  are predicted to be susceptible, and those with leucine  or tryptophan  are predicted to be resistant.
the class counts for the full data set are given at each node.
for example, the root node  contains all  <dig> cases of which  <dig> are resistant are resistant to rifampin, and the remaining  <dig> isolates are susceptible to rifampin.
the proportion of correctly classified observations across all splits as determined by re-substitution of the observations on the cross-validation pruned subtree is  <dig> .
comparing this tree to the pruned regression subtree  reveals that the two split definitions in each tree are identical.
both the regression and classification tree models are significant  based on permutation tests.
random forest analyses
the random forest analysis, which aggregates results over many tree models, each constructed on subsamples of the data, produced markedly better models as compared to the single tree-based models.
the random forest analysis of the mic data as a continuous variable, a regression problem, produced a model that explained  <dig>  of the variance.
the random forest analysis of the mic data as discrete classes , a classification problem, produced a model that correctly classified  <dig>  of the observations with corresponding sensitivity of  <dig>  and specificity of  <dig> .
although both the regression and classification random forest results are markedly better than the single tree-based models, they do lack the ease of interpretation of a tree model.
however, variable importance can be assessed in random forests by measuring the increase in group purity based individual models containing the variable.
as might be expected, the results for both regression and classification are similar and identify the same amino acid positions as being most important in determining response to rifampin as did the single tree models: primarily  <dig> and  <dig>  and much less so for  <dig> and  <dig> .
discussion
analysis of genotype and phenotype data poses several significant challenges.
data characteristics such as mixture of variable types, high dimensionality, interactions between variables, and preponderance of unordered categorical variables render many candidate analytical methods inappropriate or ineffective.
tree-based statistical models adeptly deal with these all these challenges and do so in a way that produces readily interpretable results.
through the analyses described above, we have learned several things that were not previously apparent.
we have distinguished phenotypically relevant from phenotypically irrelevant changes in genotype by establishing the relative importance of the polymorphic sequence positions, and amino acids at those positions, as they affect susceptibility to rifampin.
for example, although they are polymorphic, changes at positions  <dig>   <dig>   <dig>   <dig> and  <dig> did not significantly affect mic for rifampin.
the hierarchical importance of changes, and their contextual/conditional relationships, are depicted in the resulting tree diagrams in a readily interpretable manner.
inherent in the tree structure is the fact that earlier splits explain more variation in phenotype then subsequent splits.
for example, the first split, at position  <dig>  explains more variation then does the split based on position  <dig>
the models can be used to predict mic for rifampin where genotype is known, as well as provide the basis for hypothesis testing involving future empirical work.
furthermore models can be refined to yield improved predictions by incorporating additional data as they become available.
improved models may be possible with additional data: full length sequence of rpob may include sequence features that are responsible for some variation in mic values for rifampin, and sequence data from additional strains might lead to even more general models.
as demonstrated above, the relationship of genotype to phenotype can be quantified using tree-based statistical models and aggregations thereof.
our approach has been to use types of models in the analysis of genotype-phenotype relationships because they offer distinct advantages compared to other methods and allow for rigorous and ready interpretation of results.
tree-based and random forest analyses are readily applicable to other forms of genotypic information including data that take the general form of visualized fragments  such as microsatellites, restriction fragment length polymorphisms , amplified fragment length polymorphisms , and similar data.
tree-based and random forest analyses can also be applied directly to dna sequence data including single nucleotide polymorphisms .
in general, tree-based statistical and random forest models are applicable to all cases where the goal is to examine the relationship between genotype and phenotype.
CONCLUSIONS
relatively simple models provided accurate predictions of rifampin resistance in m. tuberculosis.
these models demonstrated that only a few variable positions in the β subunit of dna-dependent rna polymerase were responsible for most of the variation in rifampin resistance.
such models might provide the basis for quantifying rifampin resistance status based exclusively on dna sequence data and thus eliminate the requirements for time consuming culturing and antibiotic testing of clinical isolates.
more generally, the results of this study demonstrate the usefulness of tree-based statistical models and random forests in genetic analysis.

 rifampin, one of the principal drugs used in tuberculosis treatment, is a semi-synthetic antibiotic that inhibits transcription by preventing rna synthesis. isolates of mycobacterium tuberculosis resistant to rifampin occur at low to moderate frequencies in many regions of the world. mutations in rpob, the gene encoding the β subunit of dna-dependent rna polymerase, are associated with rifampin resistance. in the laboratory, drug resistance is quantified in terms of minimum inhibitory concentration, which is defined as the minimum concentration of the antibiotic in a given culture medium below which bacterial growth is not inhibited. several studies have been conducted where mic of rifampin has been measured and partial dna sequences have been determined for rpob in different isolates of m. tuberculosis. however, no model has been constructed to predict rifampin resistance based on sequence information alone. such a model might provide the basis for quantifying rifampin resistance status based exclusively on dna sequence data and thus eliminate the requirements for time consuming culturing and antibiotic testing of clinical isolates. tree-based statistical methods have generated very accurate models relating amino acid sequence of short peptides to their binding by major histocompatibility complex class i molecules with higher accuracy than artificial neural networks. both tree-based models and aggregation of such models through random forests have proven to be quite successful in other problems involving sequence data as covariates such as hiv- replication capacity and cytidine to uridine rna editing in plant mitochondria. the success of tree-based statistical models and random forests in these problems involving covariates derived from sequence data motivated our application of these models to the problem of rifampin resistance in m. tuberculosis. the response variable is a set of continuously distributed values for mic, which makes the problem one of regression. these data are used to answer the following questions: what proportion of the variance in mic is attributable to sequence differences in positions 511– of the β subunit of rna polymerase of m. tuberculosis? what particular positions, and what distribution of amino acids at those positions, are associated with most of the variance in mic? alternatively, the response variable could be cast in discrete terms: resistant or susceptible. this is possible by assuming a threshold value for mic above which an isolate is considered resistant to rifampin. among the specific questions we can answer with such a model are the following: what particular positions, and what distribution of amino acids at those positions, allow for distinguishing rifampin-susceptible and rifampin-resistant isolates of m. tuberculosis? what is the misclassification error rate associated with susceptibility prediction for these data? we address these questions and evaluate the ability to predict mic from protein sequence data using tree-based regression and classification methods. we find that these methods generate highly accurate models of rifampin resistance.  the data set used in the study consists of observations with distinct genotype-phenotype combinations. the most frequent combination has occurrences, and there are unique combinations. mic for rifampin varies from μg/ml to > μg/ml. the sequences are distributed among genotypes, of which occur uniquely in the data set. the plurality genotype is represented by samples; samples differ from the plurality by one amino acid; and the remaining samples differ from the plurality by two amino acids. some genotypes defined by the partial sequence of rpob are associated with several different phenotypes. also, some genotypic states are associated with large effects, while some have little or no effect on mic phenotype. finally, some changes in mic are not associated with changes in the sequence region examined. these genotype data are short partial sequences of a single gene, and thus they may not contain all phenotypically relevant genetic information. indeed, there is evidence that amino acid changes outside of the examined region are associated with changes in mic for rifampin. additionally, the sample size is small, which will decrease power. nonetheless these data are typical of studies surveying the genetic variation associated with antibiotic resistance and of genotype-phenotype data in general. thus they make an appropriate subject of investigation. regression tree analysis the regression tree for the relationship of rpob amino acid sequence and mic has two splits defining three terminal nodes. at each node in the tree, the mic prediction given is the mean of all isolates at that node. the first split of the topmost node consists of the entire sample and is based on the amino acid at position with those sequences having serine going to the left child node, and those having leucine or tryptophan going to the right child node. the best split for each node is that which gives the largest decrease in the error. here error is measured as the deviance, which for a continuous variable is a constant multiple of the residual sums-of-squares. reported values were determined using 10-fold cross-validation. moving down the tree the error decreases, as the sum of the deviance for each pair of child nodes is less than the deviance of the parent node. given the hierarchical nature of trees and the criterion used to choose splits, the first split, that based on position explains the highest proportion of the overall phenotypic variance. this bisection of the data results in groups that differ one hundredfold in mean mic. the subsequent partition based on position the most variable in this region, results in a > 354-fold difference in mic. the proportion of the variance in mic explained across all splits involving the two amino acid positions is. all proportions of variances explained by the model as reported here are those estimated through cross-validation and are not based on re-substitution, and thus represent appropriately conservative estimates. classification tree analysis from a clinical perspective it may be most relevant to consider the level of drug resistance as a two-state categorical variable rather than as a continuously distributed variable. in clinical practice, if an isolate of m. tuberculosis is determined to be rifampin resistant then rifampin is replaced with another antibiotic. although blood serum concentration of rifampin reaches levels of – μg/ml about – hours after ingestion, a clinically relevant mic value for dichotomizing the mic values would be lower than this peak. we conservatively define mic values ≤ μg/ml as susceptible and values > μg/ml as resistant, a definition consistent with conventional standards. with this dichotomization we can explore the use of tree-based statistical classification to predict rifampin resistance in a way that is more relevant to clinical practice. the predictor variables are again the unordered categorical designations of amino acids at polymorphic positions. the classification tree for these data has two splits based on two of the variable amino acid positions. at each node in the tree, the prediction of rifampin susceptibility status is given for all isolates at that node. the first split is based on position 531; those isolates with serine are predicted to be susceptible, and those with leucine or tryptophan are predicted to be resistant. the class counts for the full data set are given at each node. for example, the root node contains all cases of which are resistant are resistant to rifampin, and the remaining isolates are susceptible to rifampin. the proportion of correctly classified observations across all splits as determined by re-substitution of the observations on the cross-validation pruned subtree is. comparing this tree to the pruned regression subtree reveals that the two split definitions in each tree are identical. both the regression and classification tree models are significant based on permutation tests. random forest analyses the random forest analysis, which aggregates results over many tree models, each constructed on subsamples of the data, produced markedly better models as compared to the single tree-based models. the random forest analysis of the mic data as a continuous variable, a regression problem, produced a model that explained of the variance. the random forest analysis of the mic data as discrete classes, a classification problem, produced a model that correctly classified of the observations with corresponding sensitivity of and specificity of. although both the regression and classification random forest results are markedly better than the single tree-based models, they do lack the ease of interpretation of a tree model. however, variable importance can be assessed in random forests by measuring the increase in group purity based individual models containing the variable. as might be expected, the results for both regression and classification are similar and identify the same amino acid positions as being most important in determining response to rifampin as did the single tree models: primarily and and much less so for and. discussion analysis of genotype and phenotype data poses several significant challenges. data characteristics such as mixture of variable types, high dimensionality, interactions between variables, and preponderance of unordered categorical variables render many candidate analytical methods inappropriate or ineffective. tree-based statistical models adeptly deal with these all these challenges and do so in a way that produces readily interpretable results. through the analyses described above, we have learned several things that were not previously apparent. we have distinguished phenotypically relevant from phenotypically irrelevant changes in genotype by establishing the relative importance of the polymorphic sequence positions, and amino acids at those positions, as they affect susceptibility to rifampin. for example, although they are polymorphic, changes at positions and did not significantly affect mic for rifampin. the hierarchical importance of changes, and their contextual/conditional relationships, are depicted in the resulting tree diagrams in a readily interpretable manner. inherent in the tree structure is the fact that earlier splits explain more variation in phenotype then subsequent splits. for example, the first split, at position explains more variation then does the split based on position the models can be used to predict mic for rifampin where genotype is known, as well as provide the basis for hypothesis testing involving future empirical work. furthermore models can be refined to yield improved predictions by incorporating additional data as they become available. improved models may be possible with additional data: full length sequence of rpob may include sequence features that are responsible for some variation in mic values for rifampin, and sequence data from additional strains might lead to even more general models. as demonstrated above, the relationship of genotype to phenotype can be quantified using tree-based statistical models and aggregations thereof. our approach has been to use types of models in the analysis of genotype-phenotype relationships because they offer distinct advantages compared to other methods and allow for rigorous and ready interpretation of results. tree-based and random forest analyses are readily applicable to other forms of genotypic information including data that take the general form of visualized fragments such as microsatellites, restriction fragment length polymorphisms, amplified fragment length polymorphisms, and similar data. tree-based and random forest analyses can also be applied directly to dna sequence data including single nucleotide polymorphisms. in general, tree-based statistical and random forest models are applicable to all cases where the goal is to examine the relationship between genotype and phenotype.  relatively simple models provided accurate predictions of rifampin resistance in m. tuberculosis. these models demonstrated that only a few variable positions in the β subunit of dna-dependent rna polymerase were responsible for most of the variation in rifampin resistance. such models might provide the basis for quantifying rifampin resistance status based exclusively on dna sequence data and thus eliminate the requirements for time consuming culturing and antibiotic testing of clinical isolates. more generally, the results of this study demonstrate the usefulness of tree-based statistical models and random forests in genetic analysis. 
the subsequent partition based on position the most variable in this region, results in a > 354-fold difference in mic. mutations in rpob, the gene encoding the β subunit of dna-dependent rna polymerase, are associated with rifampin resistance in mycobacterium tuberculosis. such a model might provide the basis for quantifying rifampin resistance status based exclusively on dna sequence data and thus eliminate the requirements for time consuming culturing and antibiotic testing of clinical isolates. sequence data for amino acid positions 511– of rpob and associated mic of rifampin for different isolates of m. tuberculosis were taken from studies examining rifampin resistance in clinical samples from new york city and throughout japan. models may be better with improved measurements of mic and more sequence data. the random forest analysis of the mic data as discrete classes produced a model that correctly classified of the observations with sensitivity of and specificity of. the first partition in the data, based on position results in groups that differ one hundredfold in mean mic. we used tree-based statistical methods and random forests to generate models of the relationships between rpob amino acid sequence and rifampin resistance. the proportion of variance explained by a relatively simple tree-based cross-validated regression model involving two amino acid positions is. when considered as a classification problem, a cross-validated tree-based model correctly classified most of the observations and was very similar to the regression model. random forest analysis of the mic data as a continuous variable, a regression problem, produced a model that explained of the variance. several studies have been conducted where minimum inhibitory concentration of rifampin has been measured and partial dna sequences have been determined for rpob in different isolates of m. tuberculosis. however, no model has been constructed to predict rifampin resistance based on sequence information alone. highly accurate regression and classification models of rifampin resistance can be made based on this short sequence region.
 mutations in rpob, the gene encoding the β subunit of dna-dependent rna polymerase, are associated with rifampin resistance in mycobacterium tuberculosis. several studies have been conducted where minimum inhibitory concentration of rifampin has been measured and partial dna sequences have been determined for rpob in different isolates of m. tuberculosis. however, no model has been constructed to predict rifampin resistance based on sequence information alone. such a model might provide the basis for quantifying rifampin resistance status based exclusively on dna sequence data and thus eliminate the requirements for time consuming culturing and antibiotic testing of clinical isolates.  sequence data for amino acid positions 511– of rpob and associated mic of rifampin for different isolates of m. tuberculosis were taken from studies examining rifampin resistance in clinical samples from new york city and throughout japan. we used tree-based statistical methods and random forests to generate models of the relationships between rpob amino acid sequence and rifampin resistance. the proportion of variance explained by a relatively simple tree-based cross-validated regression model involving two amino acid positions is. the first partition in the data, based on position results in groups that differ one hundredfold in mean mic. the subsequent partition based on position the most variable in this region, results in a > 354-fold difference in mic. when considered as a classification problem, a cross-validated tree-based model correctly classified most of the observations and was very similar to the regression model. random forest analysis of the mic data as a continuous variable, a regression problem, produced a model that explained of the variance. the random forest analysis of the mic data as discrete classes produced a model that correctly classified of the observations with sensitivity of and specificity of.  highly accurate regression and classification models of rifampin resistance can be made based on this short sequence region. models may be better with improved measurements of mic and more sequence data. 