BACKGROUND
transposon mutagenesis is an experimental method frequently used for surveying bacterial genomes for essential regions, including genes, as well as individual protein domains, regulatory elements, and non-coding rnas that are required for survival.
for example, the himar <dig> transposon inserts randomly into ta nucleotides  <cit> .
those locations that lack a transposon insertion suggest either that the location is essential  or that it was not represented in the library of transposon mutants .
typically, around 15% of the genes in the genome of prokaryotic organisms are essential  <cit> .
knowledge of which genes are essential can be very useful for drug discovery against pathogens .
while the original method used hybridization on dna microarrays  <cit> , deep sequencing has made analysis of transposon insertion libraries much more efficient   <cit> .
short reads are obtained from the genome on either side of each transposon insertion, using amplification with a transposon-specific primer, and then mapped to the genome of the organism, revealing which locations withstood transposon insertions .
the data from a tn-seq experiment can be analyzed in several ways.
first, reads can be used to ascertain the presence or absence of insertions in a gene.
the probability that a gene lacking insertions is essential depends on the diversity of the transposon library , and can be quantified using the binomial  <cit> , negative-binomial distribution  <cit> , or extreme value distribution  <cit> .
alternatively, the number of reads at each site  can be analyzed instead of the mere presence or absence of insertions.
it can be argued that the read count carries additional information because it reflects the abundance of certain clones in the library, and hence the degree to which a region of the genome is essential.
zhang et al.
described a non-parametric test that quantifies the significance of the sum of read counts within a sliding window  along the genome to detect essential regions  <cit> .
both analysis approaches have challenges, depending on the quality of the transposon library and sequencing dataset.
methods that only look at the presence or absence of insertions can be susceptible to spurious reads, such as isolated reads that map to an essential region of the genome only because they have base-call errors.
however, it is difficult to set a threshold for a minimum number of reads, since other sites with a single read might be legitimate.
on the other hand, methods based on read counts are susceptible to several sources of variability, including spikes in the data, where there is a massive over-representation of reads at an isolated site.
the distribution of read counts is usually observed to follow a geometric distribution, but in some datasets, a few sites might have orders-of-magnitude more reads, possibly due to an artifact such as a pcr amplification bias.
this could highly influence statistics based on read counts.
it should also be noted that, even in essential genes, transposon insertions are often observed to be tolerated at the extreme n- and c-terminus of the open-reading frame .
previously, ad hoc methods were used, such as excluding insertions in the n- and c-terminal 5–20% of the orf  <cit> .
however, both the sliding window approach  <cit>  and the extreme value distribution  <cit>  based on the length of the longest sub-sequence of ta sites without insertions are designed to be robust in spite of insertions at the termini of essential genes, and have been used to identify individual essential domains within genes  <cit> .
in this paper, we describe a novel method for analyzing tn-seq data using hidden markov models .
hmms are useful for analyzing sequential datasets, in which a sequence of observed values is explained by an underlying state sequence .
for example, the genome of an organism can be viewed as an alternating sequence of essential and non-essential regions.
we show how an hmm can be designed to incorporate information from read counts at individual ta sites to infer the probability distribution over states, and then use the viterbi algorithm to infer the most likely state sequence .
the sequential-dependence of the model  helps disambiguate the interpretation of each site, thereby coupling neighboring sites together.
the resulting state transition model affords a 'smoothing’ of the read-count data, where, for example, ta sites with no insertions in non-essential regions  are tolerated because neighboring sites have insertions.
however, if a consecutive sequence of ta sites with no insertions is long enough, the most probable state sequence, as determined by the viterbi algorithm, switches locally to essential, providing a different labeling of that region.
the incorporation of read-counts in this hmm requires defining appropriate likelihood functions.
we use the geometric distribution to capture the conditional probability of read-counts in non-essential regions, reflecting the fact that sites with high read counts  are observed with much lower frequency than those with lower read counts.
furthermore, the transition probabilities of the hmm must be carefully defined so that the minimum length of essential regions matches our expectations.
a major contribution of this paper is to show how to calibrate these parameters so that the performance of the hmm will be reasonable and robust across a range of datasets, including those with high or low insertion density , and those with high or low mean read counts .
in addition, we extend the hmm with two extra states, one representing regions with particularly “low” read counts, and one representing regions with higher than average read-counts .
genes belonging to the former class of genes have been characterized before in m. tuberculosis and referred to as “growth-defect” genes  <cit> , as these are genes whose disruption leads to impaired growth of the organism.
we continue this convention here, labeling those genes with depressed read-counts as “growth-defect”  to be consistent with the prior literature.
growth-defect regions are not completely devoid of insertions , but have a lower number of insertions than non-essential regions , suggesting that these clones did not grow as well and had lower abundance due to competition with other clones in the library.
similarly, the latter class of genes  are labeled “growth-advantage” genes.
these could represent genes that have a metabolic cost  and are not necessary for growth in vitro.
the addition of these two states to our hmm allows it to distinguish regions in tn-seq data with suppressed or unusually high read counts in a statistically rigorous way.
methods
the hmm in this application is defined in a straightforward way .
we are given a sequence of observations, c <dig> .cn, which represent read counts at each ta site throughout the genome.
we assume a generative model in which the read count at each site is determined by the local state of each site, which is hidden .
each ta site is assumed to be in one of four states: qes , qgd , qne , qga .
from a given sequence of observations , we want to infer the most probable state sequence q <dig> .qn that could have generated it, based on the joint probability of counts and states:
arg maxqi..qnp
hmms are based on the markov property, i.e.
that observations and successor states only depend on the current state and are conditionally independent of previous history:
p=p
p=p
thus, because of this conditional independence, the total joint probability can be written as:
p=p∏pp
the model we propose depends critically on specifying an appropriate likelihood function for read counts.
in tn-seq experiments, the distribution of read counts can be approximated through a geometric distribution, in that sites with lower counts are more common, and sites with high counts  are much more rare.
an example histogram in shown in figure  <dig> .
thus we model the likelihood function  for qne as geometric:
p∼ciθ
the function is parameterized by θ, which represents the bernoulli probability of insertion for the geometric distribution.
the maximum-likelihood estimate for this parameter is θ=1/c¯, where c¯ is the mean read count at non empty ta sites.
we also use geometric distributions as likelihood functions for the other states.
for qes, we set θ very near to  <dig> , making sites with  <dig> counts highly probable, but also allowing sites with 1– <dig> reads .
for qgd we set θ to be θgd=1/ , reflecting the fact that the growth-defect state must represent approximately ∼100× lower read counts than qne but cannot be less than  <dig> .
for the growth-advantage state, qga, we set θ using five times the mean read count , to capture sites with significantly more insertions  locally than what is observed on average in the genome.
the net effect is that the overlapping densities of the four likelihood functions produce four distinct regions where each one dominates individually, as shown in figure  <dig>
another critical aspect of our model is the definition of the state transition probabilities, as these determine the degree of smoothing of the hmm.
let the transition matrix be defined as tab=p.
the basic assumption is that the probability of self-transition, taa, should be nearly  <dig> for all states, while tab should be nearly  <dig> for a≠b .
this assumption controls the rate at which the hmm transitions from state to state, requiring a significant change in read-counts to justify a transition and smoothing over spurious reads.
for simplicity, we use a fully symmetric matrix, and we allow any state to transition to any other state .
the magnitude of taa determines the tendency of the model to stay in one state for a certain number of steps before being forced into another state that better fits the data.
this depends on several factors, including: a) the expected minimum length of essential regions , and b) the relative magnitudes of the likelihood functions, which are competing to explain the read counts.
to estimate the expected minimum length of essential regions, we utilize the geometric distribution.
the geometric distribution describes the probability of observing a run of successes in a row, which can be used to model the distribution of run lengths.
this depends on the insertion probability in non-essential regions.
because the insertion density of the library will include essential regions with insertion probabilities which are not representative of non-essential regions.
to alleviate this bias, we estimate the insertion probability, pins, empirically by discarding regions with  <dig> or more ta sites in a row lacking insertions, and calculating the insertion density in the remaining areas.
once the insertion probability is estimated, the minimum length of essential regions, r∗, is taken to be the smallest run such that the geometric probability is less than  <dig>  < <dig> ).
typically r∗ is in the range of 5– <dig> ta sites, depending on the dataset.
the self-transition probability is then set as follows:
taa=1-λner∗
where λne represents the likelihood of observing a read-count of zero in a non-essential region.
the rationale for this formula is that the cost of staying in a state such as qne through a region devoid of insertions, must balance the penalty incurred for observing sites with  <dig> read counts ) and the number of such ta sites in a row which are likely to be observed in non-essential regions .
we will show empirically in the results section that this adaptive method for setting the transition probabilities leads to an appropriate assignment of state labels for a variety of types of datasets, and we will examine the resulting length distribution of states produced.
finally, given this definition of the hmm, we use the viterbi algorithm to calculate the most probable state sequence for a given set of read counts  <cit> .
briefly, the viterbi algorithm is a dynamic programming algorithm in which the probability of each state at step i is calculated based on the state-probability distribution from the previous step:
p=maxp×p×p
after computing this incrementally for i= <dig> .n, a back-trace is made from the most probable terminal state qn∗ to extract the sequence of states based on which states were used for updates at each step.
because the viterbi algorithm requires the multiplication of small probabilities, and the state sequence for analyzing transposon insertions is large, an hmm may incur underflow problems.
to overcome this issue, the probabilities are normalized at each iteration, as described by rabiner et al.
<cit> .
RESULTS
the hmm method was applied to a transposon mutant library of m. tuberculosis, constructed by griffin et al.
<cit> .
this library was grown on minimal media and  <dig> % glycerol, and was sequenced on an illumina gaii sequencer with a  <dig> bp read length, resulting in approximately  <dig> million reads.
the reads were mapped to the h37rv genome, and the read counts at each location in the genome were quantified .
the h37rv genome is  <dig> , <dig> bp in length, with a gc-content of  <dig> %.
it contains a total of  <dig>  ta sites, spaced on average  <dig> bp apart.
the overall insertion density, defined as ta sites with at least one insertion , is  <dig> %  of all possible insertion sites.
the average read-count at these locations is c¯= <dig> .
the mean read count was used to calculate the θ parameter for the emission probabilities of the four states as described above.
using these parameters, the most likely sequence of states responsible for the observations was obtained through the viterbi algorithm.
this sequential ordering of states provides an assessment of the essentiality of the entire h37rv genome, regardless of gene boundaries.
table  <dig> contains some statistics for the distribution of states and their observations.
statistics for essentials, growth-defect, non-essential and growth-advantage regions.
mean insertion density is defined the proportion of ta sites containing at least one insertion averaged across the regions belonging to a given state.
mean read counts are defined as the mean value of non-empty read counts , averaged across all regions belonging to the given state.
a total of  <dig> % of the genome is labeled by the essential state .
this is close to the expectations for bacterial organisms, where roughly 10%-15% of the genome is considered to be essential  <cit> .
the majority of sites are labeled non-essential , with a small percentage of sites labeled as growth-defect and growth-advantage .
essential states averaged a very small number of insertions and read counts , demonstrating that the hmm is associating the essential state with stretches devoid of insertions, though these locations can occasionally contain insertions with a very small number of reads so long as as the observations at neighboring sites are consistent with essentiality.
in contrast non-essential regions have a mean insertion density of 70%, and mean read counts of  <dig> in this dataset.
growth-defect regions have some insertions but these are dramatically reduced .
insertion density in growth-advantage regions is almost saturated , and mean read counts are on average >3× larger.
as can be seen in figure  <dig>  both the mean read counts and insertion frequencies among the states increase with the levels of non-essentiality , reflecting the fact that the hmm is successfully separating regions with average read counts and insertions from those with counts significantly lower or significantly higher than average.
analysis of essentiality of individual genes
while the viterbi algorithm does not take into consideration gene boundaries when determining the labeling of states, it is often necessary to determine the essentiality of individual genes in the genome.
to determine individual calls of essentiality, each gene is assigned the essentiality class belonging to the most frequent state found within its boundaries.
however, because genes may contain a mixture of essential and non-essential domains, genes are also classified as essential if they contain sub-sequences of sites belonging to the qes state, which are statistically longer than expected.
thus a gene is also classified as essential if it has at least n sites labeled as qes, where n is 3σ above the expected maximum run length for the gene, based on the extreme value distribution  <cit> .
the essentiality assignments obtained through the hmm method can be validated by comparing to those obtained by sassetti et al with the transposon site hybridization  method  <cit> , which used a completely different experimental methodology for read-out .
this method has been used to assess the essentiality of m. tuberculosis in vivo and in vitro  <cit> , by quantifying hybridization to dna microarrays imprinted with representative oligos for each gene.
table  <dig> shows a comparison between these two methods.
due to the significantly different methodologies, a true comparison between these methods is difficult.
for instance, sassetti et al.
recognized that trash probes for essential genes may actually hybridize to adjacent non-essential regions, particularly if the genes are small.
while the hmm does not depend on hybridization, it may have a difficulty transitioning from one state to another depending on the size of the gene.
in addition, libraries used by these methods were grown on different media and therefore are likely to identify genes that are involved in pathways that correspond to the specific growth media used.
comparison of essentiality predictions with the trash method.
despite these limitations, there is significant agreement in their assessment of essentiality, with  <dig> % of essential and non-essential genes in concordance with the previous results .
approximately half of the genes labeled as 'growth-defect’ by the hmm were previously determined to be essentials, and half as non-essentials, reflecting the borderline nature of these genes and the utility of having an intermediate category.
these are discussed further below.
<dig> genes were called 'growth-advantaged’ due to an excess of transposon insertions, and all of these were previously categorized as non-essentials.
sassetti et al.
<cit>  also defined a set of  <dig> 'growth-defect’ genes.
importantly, these were not characterized by experimentally determining growth rates in individual transposon-insertion mutants.
rather, they were identified as genes that matched the criterion for 'non-essential’ on the first plating of the library , but which had much lower ratios upon re-plating .
the interpretation of these genes is that transposon insertions were not lethal, but that the mutants had a slower growth rate, resulting in gradual depletion in the library due to competition during culturing.
in the experiment from which the dataset we use was derived  <cit> , the dna for sequencing was extracted from the library immediately after selection, thus corresponding to the 'first plating’.
consistent with this, most of these genes  exhibited transposon insertions in our dataset and were categorized by the hmm as non-essential.
we speculate that, if the library had been expanded after selection, clones with insertions in these genes would have gradually decreased in abundance.
although the methods disagree on essentiality of some genes, some of these disagreements may be due to differences in the growth media, as well as the different interpretations of essentiality.
for example glpk, a glycerol kinase, is necessary for glycerol metabolism , but it is not necessary when the library is grown on glucose .
in addition, these differences can also be due to the fact that we identify genes containing essential domains as “essential”, while this distinction was not made in the original trash experiments.
in fact, all of the genes classified as essential by the hmm and as non-essential by the trash method are devoid of insertions in the majority of their ta sites or contain stretches that are significantly longer than expected, suggesting these genes are essential in this library on glycerol.
among these genes are ppm <dig>  and ppp , which independent experiments have shown contain essential domains  <cit> .
in addition to the trash method, we compare our results to those obtained with the reads-based method developed by zhang et al.
<cit> .
this method is capable of assessing the essentiality of the entire genome by looking at the read counts that fall within windows of 400– <dig> bp, and estimating a p-value for each of these windows in the genome to quantify how these regions deviate from expectations.
our results correlate well with the results obtained by window-based method, with a  <dig> % match in the classification of genes .
in addition, the essential and growth-defect states had ta sites with an average p-value of  <dig> , and non-essential and growth-advantage states an average p-value of  <dig>  .
performance on other datasets
to demonstrate that the hmm works on other datasets, we ran it on a tn-seq dataset from h. influenza .
the h. influenza kw <dig> genome is less than half the size of m. tuberculosis  but significantly more at-rich , so there are more ta sites  but they are spaced more closely .
the tn-seq dataset contains  <dig>  reads, hitting only  <dig> % of the ta sites, with a mean read count of  <dig>  .
running the hmm on this lower-density dataset results in  <dig> genes being labeled as essential,  <dig> as non-essential,  <dig> as growth-defect, and  <dig> as growth-advantage.
this distribution is very close to the assignments determined by gawronski et al.
<cit> , who found  <dig> essentials , and  <dig> growth-defect genes .
the overlap  between the essential genes detected by both their method and ours was 94% , and the intersection between their list of growth-defect genes and ours was 60% .
the overlap between essential genes found by the hmm method and those found by gawronski et al.
significantly larger than the overlap between the trash method described above .
this high level of agreement between the two comparisons suggests that the quality of the data used in the analysis  contributes significantly to the quality of the analysis.
in addition, we applied the hmm method to three modified datasets, constructed to represent libraries of different sizes and different volumes of sequencing data.
these datasets were constructed by modifying the original h37rv library analyzed before, to emulate cases where transposon mutant libraries may be sparse or where the amount of sequencing performed on the library is lower .
the first dataset was constructed by setting the read counts at random ta sites to zero , thus lowering the mean insertion density of the dataset while keeping the magnitude of the remaining read-counts the same.
this dataset emulates libraries with significantly less diversity of insertions.
the second dataset was constructed by randomly perturbing approximately one-half of the reads, lowering the magnitude of these reads while keeping the total number of insertions equal.
this dataset represents libraries for which the amount of sequencing performed is significantly less, producing read counts with lower magnitudes.
the final dataset was a combination of these two operations, resulting in a dataset with both lower insertion density and lower mean read count.
statistics about the distribution of reads and insertions in these datasets are shown in table  <dig>
insertion density and mean read count for the glycerol, low density, low reads and, low reads and density datasets.
as can be seen in table  <dig>  the hmm is robust, and capable of adapting to libraries with very different insertion densities and mean read counts, providing results which are generally consistent with each other.
the fraction of the genome labeled as essential is approximately the same in all four datasets .
although the decreased density will result in longer stretches of the genome without a transposon insertion, the hmm is capable of adapting its parameters to become more conservative in designating regions without insertions as essential.
statistics for percentage of sites labeled as essential, growth-defect, non-essential and growth-advantage in the glycerol, low density, low reads and, low reads and density datasets.
growth-defect and growth-advantage genes
one of the principle advantages of our 4-state hmm is that it can distinguish local regions of the genome with significantly depressed or elevated read counts .
the former could represent genes whose disruption is not lethal but could lead to a growth-defect, resulting in a lower representation of clones in the library, and thus a lower abundance of sequencing reads  <cit> .
by analogy, regions with significantly greater than average reads could represent genes whose disruption leads to a growth advantage.
in the h37rv dataset, there were  <dig> genes labeled as qgd , and  <dig> genes labeled as qga .
these are discussed in turn below.
among the genes labeled as growth-defect, there are several notable ones for which a biological explanation can be made .
one of these is pbpa, a penicillin-binding protein in mtb.
mutants have shown decreased growth rates and defective cell septation when pbpa is knocked out m. smegmatis <cit> .
in addition, the wild-type phenotype was restored by complementing in pbpa from m. tuberculosis, suggesting that pbpa plays an important role in cell-division and disruption of this gene might lead to impaired growth in m. tuberculosis.
in fact, this region contains an average insertion density of  <dig> , and an average read-count of  <dig>  significantly below the global insertion density  and read-counts .
notable regions labeled as 'growth-defect’  by the hmm, in which there is a suppression of reads, suggesting an under-representation of transposon insertions in the library.
note that a lower read count is reflected in both a low read count  and a low insertion density .
recent structural and enzymatic studies have shown that bfrb and its ortholog, bfra, are not completely interchangeable.
although they are both ferritin proteins used for iron storage, bfrb has a 20-aa c-terminal extension that enhances its iron oxidation activity  <cit> .
thus growth of bfrb mutants might be hindered because bfra cannot perform this function as efficiently.
in fact, data from the original trash experiments shows that bfrb had a much lower hybridization ratio  compared to bfra , suggesting clones with insertions in bfrb were less competitive.
many genes in the mycobactin biosynthesis cluster  are also labeled as growth-defect genes, suggesting that transposon mutants are viable but grow more slowly than wild-type.
because mtb has only one  iron acquisition system, which is mycobactin-dependent, these biosynthetic genes are essential in iron-depleted environments and non-essential in those environments that are rich in iron.
indeed, it has been shown that mycobactin-deficient mutants of mtb, the growth rate is dependent on the iron concentration  <cit> .
in the original trash experiments , mbtb was specifically shown to be cause a slow-growth phenotype when disrupted, with insertion mutants gradually decreasing in abundance in the library with successive platings  <cit> .
another interesting growth-defect gene is glpx.
glpk , which is the first step in glycerol incorporation, is essential as expected .
glpx is a fructose- <dig> -bisphosphatase, which also should be required when grown on gluconeogenic substrates by circumventing a non-reversible step in glycolysis pathway to generate glucose  <cit> .
in mtb the unexpected non-essentiality of glpx for growth on glycerol has been previously noted  <cit> .
one possible explanation is that rv2131c , an inositol monophosphatase, might also have partial fructose- <dig> -bisphosphatase activity  <cit> .
icl  is also identified as a growth-defect gene in this dataset.
this is one of the two enzymes on the glyoxylate shunt, which has been shown to be critical for infection, based on attenuation of knockouts in mice  <cit> .
as anticipated, icl is essential for growth on fatty-acid substrates like acetate  <cit> .
however, recent evidence suggests that the glyoxylate shunt might play a role even in growth on other carbon sources such as carbohydrates.
for instance, icl knockouts have displayed a growth-defect  on glucose  <cit> .
more recently, it has been shown that inhibitors of malate synthase  are active against cultures whether grown on acetate or glucose  <cit> .
thus, the fact that the hmm labels icl as a growth-defect region in this dataset obtained from growth on glycerol is consistent with these findings and suggests that icl plays an unexpected metabolic role in mtb even when growing on carbon sources other than fatty acids.
another gene identified as belonging to the growth-defect category is tres, which is involved in the trehalose pathway.
trehalose is one of the principle carbohydrates synthesized in mycobacteria.
it is used in producing cell-wall glycolipid components , and is inter-converted with other sugars like glucose and maltose.
the latter are polymerized into intracellular glycogen  and capsular glucan.
several genes in this network have been shown to be essential in vitro, including galu, glga, glgb, pep <dig>  and glge .
however, tres is labeled as a growth-defect gene.
tres is responsible for interconverting trehalose and maltose  <cit> .
it is possible that the organism is sensitive to perturbations of this network .
in fact, it was previously shown that transposon-insertion mutants of tres/rv <dig> display a slow-growth phenotype  <cit> .
as noted before, our 4-state hmm is also capable of detecting regions with unexpectedly high read-counts that might confer growth-advantages to the organism when disrupted.
the  <dig> most notable growth-advantage regions are shown in table  <dig> .
one region of the genome that stands out is the pdim locus, rv2930-rv <dig>  this locus contains genes involved in the biosynthesis of phthiocerol dimycocerosate , including fadd <dig> and ppsabcde.
in addition, other genes outside this locus believed to be involved in pdim biosynthesis, like papa <dig> and mas, are identified as well.
these genes contain read counts well above the global average .
fadd <dig> itself has a mean read count of  <dig>  more than three times the average throughout the genome.
ppsde had a mean read count of  <dig>  and ppsabc a mean read count of  <dig>  pdim is a cell-wall associated glycolipid that modulates the immune response in the host  <cit> .
although it is required for virulence , it is not required for survival in vitro  <cit> .
in fact, biosynthesis of pdim requires resources and imparts a metabolic cost, hence disruption of this pathway is advantageous to cells.
due to the increased metabolic cost, it is widely observed that m. tuberculosis stocks maintained in the lab frequently lose the ability to synthesize pdim via acquisitions of mutations in these genes, often leading larger colony sizes  <cit> .
this growth advantage and consequent selection effect likely explains why clones with transposon insertions in the pdim locus are over-represented in the library.
regions labeled as 'growth-advantage’  by the hmm, in which there is an excess of reads, suggesting an over-representation of transposon insertions.
discussion and
CONCLUSIONS
the hmm described in this work enables the characterization of essentiality throughout an entire bacterial genome from sequencing data of transposon mutagenesis experiments.
although several computational methods have previously been proposed for analyzing tn-seq data, including some based on presence/absence of reads  <cit>  as well as non-parametric models that take quantitative read counts into consideration  <cit> , an hmm provides several advantages over these methods.
for example, an hmm provides a smoothing over adjacent sites that couples them together to help disambiguate the interpretation of read counts at individuals sites.
another advantage of using the hmm is that it is not restricted to annotated gene boundaries, and can identify independent regulatory regions, non-coding rnas, and protein domains that are required for survival.
while methods that depend on a sliding-window  are also capable of assessing essentiality over the entire genome by locally averaging over adjacent ta sites, an hmm formalizes this process in a statistically rigorous way.
in addition, by assessing essentiality among regions in the genome, the hmm can also tolerate insertions in the n- and c- termini of genes, without the need of discarding insertions at these locations in an ad-hoc manner as some methods have done previously  <cit> .
a potential limitation of our method is that it does not take into consideration the doubling-rate or expansion time of the library when estimating the parameters of the model.
this can affect the under- and over-representation of mutants, and therefore the number of reads these genes will contain in the sequence data.
because the hmm depends on the read counts for individual genes, it may be susceptible to libraries that are constructed in different ways.
indeed, the ability to detect essential genomic regions from transposon-insertion sequencing data is highly dependent on the quality of the dataset.
in practice, tn-seq datasets can span a range from hundreds of thousands to millions of reads, but below some point, there are not enough reads to discriminate essential regions confidently.
similarly, well-saturated libraries can have insertions at > 50% of ta sites, but other datasets are more sparse , again increasing the difficulty in distinguishing essential regions.
additionally, some sites may contain reads that are orders-of-magnitude larger due to pcr amplification or the development of “hotspots” due to the interactions of the transposon and the organism’s replication machinery  <cit> .
these problems could be alleviated by comparing cultures before and after the passage in the media or comparing datasets derived using different transposons with specificity to different sites.
while not directly filtered by our hmm, we showed that the parameter estimation equations we propose work on a wide range of real datasets.
in particular we show that they work on dense datasets , as well as sparse ones , and even on artificial datasets down to 25% insertion density.
in all cases, the hmm is stable in that it outputs about the same proportion of essential regions, so as the volume of data decreases, the hmm adapts its predictions and becomes increasingly conservative.
one of the major advantages of using an hmm to analyze transposon mutagenesis data is that additional states may be introduced to capture distinct types of genomic regions .
in this work we add states to capture regions whose disruption leads to a growth defect or a growth advantage.
van opijnen et al.
<cit>  have shown that relative abundance of insertions in a gene can be correlated quantitatively with growth rate , though it depends on the number of generations the culture is allowed to grow.
eventually, clones with insertions in growth-defect genes will be depleted from the library due to competition .
sassetti et al.
found substantial reductions in abundance, even after a second round of plating, where hybridization ratios for certain genes dropped from non-essential  to essential .
while our model only distinguishes one class of growth-defect genes, it could be expanded to more states, discerning finer gradations of growth impairment  <cit> .
a total of  <dig> growth-defect genes, and  <dig> growth-advantage genes in m. tuberculosis h37rv were identified by our 4-state hmm, several of which have been shown to be biologically valid, as knockout mutants have been shown to grow slower  than the parental strain.
identifying these finer distinctions of essentiality  can enrich our understanding of the biological roles of genes.
for example, we found that icl  is labeled as growth-defect .
historically, icl has been viewed as essential in m. tuberculosis specifically for growth on fatty-acid substrates, and non-essential otherwise.
however, that view is too limiting.
icl knock-out strains have in fact been observed to display a growth defect when grown on glucose  <cit> , and the suppression of reads we observed in icl in the transposon mutagenesis data is consistent with this , suggesting that this gene plays an additional role that is not well-appreciated.
availability and requirements
the hmm model utilized in this paper is made available online in the tn-hmm package.
project name: tn-hmm
project home page:http://saclab.tamu.edu/essentiality/hmm/
operating system: linux, windows
programming language: python
license: creative commons attribution-noncommercial
any restrictions to use by non-academics: none
competing interests
the authors declare that they have no competing interests.
authors’ contributions
tri and mad were jointly responsible for the development of the algorithm and the writing of the manuscript.
both authors read and approved the final manuscript.
supplementary material
additional file 1
table s <dig>  excel spreadsheet  containing a list of growth-defect regions identified by the hmm in m. tuberculosis h37rv grown on glycerol .
click here for file
additional file 2
table s <dig>  excel spreadsheet  containing a list of growth-advantage regions identified by the hmm in m. tuberculosis h37rv grown on glycerol .
click here for file
acknowledgements
this work was supported in part by the nih grant u <dig> ai <dig> .

 transposon mutagenesis is an experimental method frequently used for surveying bacterial genomes for essential regions, including genes, as well as individual protein domains, regulatory elements, and non-coding rnas that are required for survival. for example, the himar transposon inserts randomly into ta nucleotides. those locations that lack a transposon insertion suggest either that the location is essential or that it was not represented in the library of transposon mutants. typically, around 15% of the genes in the genome of prokaryotic organisms are essential. knowledge of which genes are essential can be very useful for drug discovery against pathogens. while the original method used hybridization on dna microarrays, deep sequencing has made analysis of transposon insertion libraries much more efficient. short reads are obtained from the genome on either side of each transposon insertion, using amplification with a transposon-specific primer, and then mapped to the genome of the organism, revealing which locations withstood transposon insertions. the data from a tn-seq experiment can be analyzed in several ways. first, reads can be used to ascertain the presence or absence of insertions in a gene. the probability that a gene lacking insertions is essential depends on the diversity of the transposon library, and can be quantified using the binomial, negative-binomial distribution, or extreme value distribution. alternatively, the number of reads at each site can be analyzed instead of the mere presence or absence of insertions. it can be argued that the read count carries additional information because it reflects the abundance of certain clones in the library, and hence the degree to which a region of the genome is essential. zhang et al. described a non-parametric test that quantifies the significance of the sum of read counts within a sliding window along the genome to detect essential regions. both analysis approaches have challenges, depending on the quality of the transposon library and sequencing dataset. methods that only look at the presence or absence of insertions can be susceptible to spurious reads, such as isolated reads that map to an essential region of the genome only because they have base-call errors. however, it is difficult to set a threshold for a minimum number of reads, since other sites with a single read might be legitimate. on the other hand, methods based on read counts are susceptible to several sources of variability, including spikes in the data, where there is a massive over-representation of reads at an isolated site. the distribution of read counts is usually observed to follow a geometric distribution, but in some datasets, a few sites might have orders-of-magnitude more reads, possibly due to an artifact such as a pcr amplification bias. this could highly influence statistics based on read counts. it should also be noted that, even in essential genes, transposon insertions are often observed to be tolerated at the extreme n- and c-terminus of the open-reading frame. previously, ad hoc methods were used, such as excluding insertions in the n- and c-terminal 5–20% of the orf. however, both the sliding window approach and the extreme value distribution based on the length of the longest sub-sequence of ta sites without insertions are designed to be robust in spite of insertions at the termini of essential genes, and have been used to identify individual essential domains within genes. in this paper, we describe a novel method for analyzing tn-seq data using hidden markov models. hmms are useful for analyzing sequential datasets, in which a sequence of observed values is explained by an underlying state sequence. for example, the genome of an organism can be viewed as an alternating sequence of essential and non-essential regions. we show how an hmm can be designed to incorporate information from read counts at individual ta sites to infer the probability distribution over states, and then use the viterbi algorithm to infer the most likely state sequence. the sequential-dependence of the model helps disambiguate the interpretation of each site, thereby coupling neighboring sites together. the resulting state transition model affords a'smoothing’ of the read-count data, where, for example, ta sites with no insertions in non-essential regions are tolerated because neighboring sites have insertions. however, if a consecutive sequence of ta sites with no insertions is long enough, the most probable state sequence, as determined by the viterbi algorithm, switches locally to essential, providing a different labeling of that region. the incorporation of read-counts in this hmm requires defining appropriate likelihood functions. we use the geometric distribution to capture the conditional probability of read-counts in non-essential regions, reflecting the fact that sites with high read counts are observed with much lower frequency than those with lower read counts. furthermore, the transition probabilities of the hmm must be carefully defined so that the minimum length of essential regions matches our expectations. a major contribution of this paper is to show how to calibrate these parameters so that the performance of the hmm will be reasonable and robust across a range of datasets, including those with high or low insertion density, and those with high or low mean read counts. in addition, we extend the hmm with two extra states, one representing regions with particularly “low” read counts, and one representing regions with higher than average read-counts. genes belonging to the former class of genes have been characterized before in m. tuberculosis and referred to as “growth-defect” genes, as these are genes whose disruption leads to impaired growth of the organism. we continue this convention here, labeling those genes with depressed read-counts as “growth-defect” to be consistent with the prior literature. growth-defect regions are not completely devoid of insertions, but have a lower number of insertions than non-essential regions, suggesting that these clones did not grow as well and had lower abundance due to competition with other clones in the library. similarly, the latter class of genes are labeled “growth-advantage” genes. these could represent genes that have a metabolic cost and are not necessary for growth in vitro. the addition of these two states to our hmm allows it to distinguish regions in tn-seq data with suppressed or unusually high read counts in a statistically rigorous way. methods the hmm in this application is defined in a straightforward way. we are given a sequence of observations, c.cn, which represent read counts at each ta site throughout the genome. we assume a generative model in which the read count at each site is determined by the local state of each site, which is hidden. each ta site is assumed to be in one of four states: qes, qgd, qne, qga. from a given sequence of observations, we want to infer the most probable state sequence q.qn that could have generated it, based on the joint probability of counts and states: arg maxqi..qnp hmms are based on the markov property, i.e. that observations and successor states only depend on the current state and are conditionally independent of previous history: p=p p=p thus, because of this conditional independence, the total joint probability can be written as: p=p∏pp the model we propose depends critically on specifying an appropriate likelihood function for read counts. in tn-seq experiments, the distribution of read counts can be approximated through a geometric distribution, in that sites with lower counts are more common, and sites with high counts are much more rare. an example histogram in shown in figure. thus we model the likelihood function for qne as geometric: p∼ciθ the function is parameterized by θ, which represents the bernoulli probability of insertion for the geometric distribution. the maximum-likelihood estimate for this parameter is θ=1/c¯, where c¯ is the mean read count at non empty ta sites. we also use geometric distributions as likelihood functions for the other states. for qes, we set θ very near to, making sites with counts highly probable, but also allowing sites with 1– reads. for qgd we set θ to be θgd=1/, reflecting the fact that the growth-defect state must represent approximately ∼100× lower read counts than qne but cannot be less than. for the growth-advantage state, qga, we set θ using five times the mean read count, to capture sites with significantly more insertions locally than what is observed on average in the genome. the net effect is that the overlapping densities of the four likelihood functions produce four distinct regions where each one dominates individually, as shown in figure another critical aspect of our model is the definition of the state transition probabilities, as these determine the degree of smoothing of the hmm. let the transition matrix be defined as tab=p. the basic assumption is that the probability of self-transition, taa, should be nearly for all states, while tab should be nearly for a≠b. this assumption controls the rate at which the hmm transitions from state to state, requiring a significant change in read-counts to justify a transition and smoothing over spurious reads. for simplicity, we use a fully symmetric matrix, and we allow any state to transition to any other state. the magnitude of taa determines the tendency of the model to stay in one state for a certain number of steps before being forced into another state that better fits the data. this depends on several factors, including: a) the expected minimum length of essential regions, and b) the relative magnitudes of the likelihood functions, which are competing to explain the read counts. to estimate the expected minimum length of essential regions, we utilize the geometric distribution. the geometric distribution describes the probability of observing a run of successes in a row, which can be used to model the distribution of run lengths. this depends on the insertion probability in non-essential regions. because the insertion density of the library will include essential regions with insertion probabilities which are not representative of non-essential regions. to alleviate this bias, we estimate the insertion probability, pins, empirically by discarding regions with or more ta sites in a row lacking insertions, and calculating the insertion density in the remaining areas. once the insertion probability is estimated, the minimum length of essential regions, r∗, is taken to be the smallest run such that the geometric probability is less than < ). typically r∗ is in the range of 5– ta sites, depending on the dataset. the self-transition probability is then set as follows: taa=1-λner∗ where λne represents the likelihood of observing a read-count of zero in a non-essential region. the rationale for this formula is that the cost of staying in a state such as qne through a region devoid of insertions, must balance the penalty incurred for observing sites with read counts ) and the number of such ta sites in a row which are likely to be observed in non-essential regions. we will show empirically in the results section that this adaptive method for setting the transition probabilities leads to an appropriate assignment of state labels for a variety of types of datasets, and we will examine the resulting length distribution of states produced. finally, given this definition of the hmm, we use the viterbi algorithm to calculate the most probable state sequence for a given set of read counts. briefly, the viterbi algorithm is a dynamic programming algorithm in which the probability of each state at step i is calculated based on the state-probability distribution from the previous step: p=maxp×p×p after computing this incrementally for i=.n, a back-trace is made from the most probable terminal state qn∗ to extract the sequence of states based on which states were used for updates at each step. because the viterbi algorithm requires the multiplication of small probabilities, and the state sequence for analyzing transposon insertions is large, an hmm may incur underflow problems. to overcome this issue, the probabilities are normalized at each iteration, as described by rabiner et al..  the hmm method was applied to a transposon mutant library of m. tuberculosis, constructed by griffin et al.. this library was grown on minimal media and % glycerol, and was sequenced on an illumina gaii sequencer with a bp read length, resulting in approximately million reads. the reads were mapped to the h37rv genome, and the read counts at each location in the genome were quantified. the h37rv genome is, bp in length, with a gc-content of %. it contains a total of ta sites, spaced on average bp apart. the overall insertion density, defined as ta sites with at least one insertion, is % of all possible insertion sites. the average read-count at these locations is c¯=. the mean read count was used to calculate the θ parameter for the emission probabilities of the four states as described above. using these parameters, the most likely sequence of states responsible for the observations was obtained through the viterbi algorithm. this sequential ordering of states provides an assessment of the essentiality of the entire h37rv genome, regardless of gene boundaries. table contains some statistics for the distribution of states and their observations. statistics for essentials, growth-defect, non-essential and growth-advantage regions. mean insertion density is defined the proportion of ta sites containing at least one insertion averaged across the regions belonging to a given state. mean read counts are defined as the mean value of non-empty read counts, averaged across all regions belonging to the given state. a total of % of the genome is labeled by the essential state. this is close to the expectations for bacterial organisms, where roughly 10%-15% of the genome is considered to be essential. the majority of sites are labeled non-essential, with a small percentage of sites labeled as growth-defect and growth-advantage. essential states averaged a very small number of insertions and read counts, demonstrating that the hmm is associating the essential state with stretches devoid of insertions, though these locations can occasionally contain insertions with a very small number of reads so long as as the observations at neighboring sites are consistent with essentiality. in contrast non-essential regions have a mean insertion density of 70%, and mean read counts of in this dataset. growth-defect regions have some insertions but these are dramatically reduced. insertion density in growth-advantage regions is almost saturated, and mean read counts are on average >3× larger. as can be seen in figure both the mean read counts and insertion frequencies among the states increase with the levels of non-essentiality, reflecting the fact that the hmm is successfully separating regions with average read counts and insertions from those with counts significantly lower or significantly higher than average. analysis of essentiality of individual genes while the viterbi algorithm does not take into consideration gene boundaries when determining the labeling of states, it is often necessary to determine the essentiality of individual genes in the genome. to determine individual calls of essentiality, each gene is assigned the essentiality class belonging to the most frequent state found within its boundaries. however, because genes may contain a mixture of essential and non-essential domains, genes are also classified as essential if they contain sub-sequences of sites belonging to the qes state, which are statistically longer than expected. thus a gene is also classified as essential if it has at least n sites labeled as qes, where n is 3σ above the expected maximum run length for the gene, based on the extreme value distribution. the essentiality assignments obtained through the hmm method can be validated by comparing to those obtained by sassetti et al with the transposon site hybridization method, which used a completely different experimental methodology for read-out. this method has been used to assess the essentiality of m. tuberculosis in vivo and in vitro, by quantifying hybridization to dna microarrays imprinted with representative oligos for each gene. table shows a comparison between these two methods. due to the significantly different methodologies, a true comparison between these methods is difficult. for instance, sassetti et al. recognized that trash probes for essential genes may actually hybridize to adjacent non-essential regions, particularly if the genes are small. while the hmm does not depend on hybridization, it may have a difficulty transitioning from one state to another depending on the size of the gene. in addition, libraries used by these methods were grown on different media and therefore are likely to identify genes that are involved in pathways that correspond to the specific growth media used. comparison of essentiality predictions with the trash method. despite these limitations, there is significant agreement in their assessment of essentiality, with % of essential and non-essential genes in concordance with the previous results. approximately half of the genes labeled as'growth-defect’ by the hmm were previously determined to be essentials, and half as non-essentials, reflecting the borderline nature of these genes and the utility of having an intermediate category. these are discussed further below. genes were called'growth-advantaged’ due to an excess of transposon insertions, and all of these were previously categorized as non-essentials. sassetti et al. also defined a set of'growth-defect’ genes. importantly, these were not characterized by experimentally determining growth rates in individual transposon-insertion mutants. rather, they were identified as genes that matched the criterion for'non-essential’ on the first plating of the library, but which had much lower ratios upon re-plating. the interpretation of these genes is that transposon insertions were not lethal, but that the mutants had a slower growth rate, resulting in gradual depletion in the library due to competition during culturing. in the experiment from which the dataset we use was derived, the dna for sequencing was extracted from the library immediately after selection, thus corresponding to the'first plating’. consistent with this, most of these genes exhibited transposon insertions in our dataset and were categorized by the hmm as non-essential. we speculate that, if the library had been expanded after selection, clones with insertions in these genes would have gradually decreased in abundance. although the methods disagree on essentiality of some genes, some of these disagreements may be due to differences in the growth media, as well as the different interpretations of essentiality. for example glpk, a glycerol kinase, is necessary for glycerol metabolism, but it is not necessary when the library is grown on glucose. in addition, these differences can also be due to the fact that we identify genes containing essential domains as “essential”, while this distinction was not made in the original trash experiments. in fact, all of the genes classified as essential by the hmm and as non-essential by the trash method are devoid of insertions in the majority of their ta sites or contain stretches that are significantly longer than expected, suggesting these genes are essential in this library on glycerol. among these genes are ppm and ppp, which independent experiments have shown contain essential domains. in addition to the trash method, we compare our results to those obtained with the reads-based method developed by zhang et al.. this method is capable of assessing the essentiality of the entire genome by looking at the read counts that fall within windows of 400– bp, and estimating a p-value for each of these windows in the genome to quantify how these regions deviate from expectations. our results correlate well with the results obtained by window-based method, with a % match in the classification of genes. in addition, the essential and growth-defect states had ta sites with an average p-value of, and non-essential and growth-advantage states an average p-value of. performance on other datasets to demonstrate that the hmm works on other datasets, we ran it on a tn-seq dataset from h. influenza. the h. influenza kw genome is less than half the size of m. tuberculosis but significantly more at-rich, so there are more ta sites but they are spaced more closely. the tn-seq dataset contains reads, hitting only % of the ta sites, with a mean read count of. running the hmm on this lower-density dataset results in genes being labeled as essential, as non-essential, as growth-defect, and as growth-advantage. this distribution is very close to the assignments determined by gawronski et al., who found essentials, and growth-defect genes. the overlap between the essential genes detected by both their method and ours was 94%, and the intersection between their list of growth-defect genes and ours was 60%. the overlap between essential genes found by the hmm method and those found by gawronski et al. significantly larger than the overlap between the trash method described above. this high level of agreement between the two comparisons suggests that the quality of the data used in the analysis contributes significantly to the quality of the analysis. in addition, we applied the hmm method to three modified datasets, constructed to represent libraries of different sizes and different volumes of sequencing data. these datasets were constructed by modifying the original h37rv library analyzed before, to emulate cases where transposon mutant libraries may be sparse or where the amount of sequencing performed on the library is lower. the first dataset was constructed by setting the read counts at random ta sites to zero, thus lowering the mean insertion density of the dataset while keeping the magnitude of the remaining read-counts the same. this dataset emulates libraries with significantly less diversity of insertions. the second dataset was constructed by randomly perturbing approximately one-half of the reads, lowering the magnitude of these reads while keeping the total number of insertions equal. this dataset represents libraries for which the amount of sequencing performed is significantly less, producing read counts with lower magnitudes. the final dataset was a combination of these two operations, resulting in a dataset with both lower insertion density and lower mean read count. statistics about the distribution of reads and insertions in these datasets are shown in table insertion density and mean read count for the glycerol, low density, low reads and, low reads and density datasets. as can be seen in table the hmm is robust, and capable of adapting to libraries with very different insertion densities and mean read counts, providing results which are generally consistent with each other. the fraction of the genome labeled as essential is approximately the same in all four datasets. although the decreased density will result in longer stretches of the genome without a transposon insertion, the hmm is capable of adapting its parameters to become more conservative in designating regions without insertions as essential. statistics for percentage of sites labeled as essential, growth-defect, non-essential and growth-advantage in the glycerol, low density, low reads and, low reads and density datasets. growth-defect and growth-advantage genes one of the principle advantages of our 4-state hmm is that it can distinguish local regions of the genome with significantly depressed or elevated read counts. the former could represent genes whose disruption is not lethal but could lead to a growth-defect, resulting in a lower representation of clones in the library, and thus a lower abundance of sequencing reads. by analogy, regions with significantly greater than average reads could represent genes whose disruption leads to a growth advantage. in the h37rv dataset, there were genes labeled as qgd, and genes labeled as qga. these are discussed in turn below. among the genes labeled as growth-defect, there are several notable ones for which a biological explanation can be made. one of these is pbpa, a penicillin-binding protein in mtb. mutants have shown decreased growth rates and defective cell septation when pbpa is knocked out m. smegmatis. in addition, the wild-type phenotype was restored by complementing in pbpa from m. tuberculosis, suggesting that pbpa plays an important role in cell-division and disruption of this gene might lead to impaired growth in m. tuberculosis. in fact, this region contains an average insertion density of, and an average read-count of significantly below the global insertion density and read-counts. notable regions labeled as'growth-defect’ by the hmm, in which there is a suppression of reads, suggesting an under-representation of transposon insertions in the library. note that a lower read count is reflected in both a low read count and a low insertion density. recent structural and enzymatic studies have shown that bfrb and its ortholog, bfra, are not completely interchangeable. although they are both ferritin proteins used for iron storage, bfrb has a 20-aa c-terminal extension that enhances its iron oxidation activity. thus growth of bfrb mutants might be hindered because bfra cannot perform this function as efficiently. in fact, data from the original trash experiments shows that bfrb had a much lower hybridization ratio compared to bfra, suggesting clones with insertions in bfrb were less competitive. many genes in the mycobactin biosynthesis cluster are also labeled as growth-defect genes, suggesting that transposon mutants are viable but grow more slowly than wild-type. because mtb has only one iron acquisition system, which is mycobactin-dependent, these biosynthetic genes are essential in iron-depleted environments and non-essential in those environments that are rich in iron. indeed, it has been shown that mycobactin-deficient mutants of mtb, the growth rate is dependent on the iron concentration. in the original trash experiments, mbtb was specifically shown to be cause a slow-growth phenotype when disrupted, with insertion mutants gradually decreasing in abundance in the library with successive platings. another interesting growth-defect gene is glpx. glpk, which is the first step in glycerol incorporation, is essential as expected. glpx is a fructose- -bisphosphatase, which also should be required when grown on gluconeogenic substrates by circumventing a non-reversible step in glycolysis pathway to generate glucose. in mtb the unexpected non-essentiality of glpx for growth on glycerol has been previously noted. one possible explanation is that rv2131c, an inositol monophosphatase, might also have partial fructose- -bisphosphatase activity. icl is also identified as a growth-defect gene in this dataset. this is one of the two enzymes on the glyoxylate shunt, which has been shown to be critical for infection, based on attenuation of knockouts in mice. as anticipated, icl is essential for growth on fatty-acid substrates like acetate. however, recent evidence suggests that the glyoxylate shunt might play a role even in growth on other carbon sources such as carbohydrates. for instance, icl knockouts have displayed a growth-defect on glucose. more recently, it has been shown that inhibitors of malate synthase are active against cultures whether grown on acetate or glucose. thus, the fact that the hmm labels icl as a growth-defect region in this dataset obtained from growth on glycerol is consistent with these findings and suggests that icl plays an unexpected metabolic role in mtb even when growing on carbon sources other than fatty acids. another gene identified as belonging to the growth-defect category is tres, which is involved in the trehalose pathway. trehalose is one of the principle carbohydrates synthesized in mycobacteria. it is used in producing cell-wall glycolipid components, and is inter-converted with other sugars like glucose and maltose. the latter are polymerized into intracellular glycogen and capsular glucan. several genes in this network have been shown to be essential in vitro, including galu, glga, glgb, pep and glge. however, tres is labeled as a growth-defect gene. tres is responsible for interconverting trehalose and maltose. it is possible that the organism is sensitive to perturbations of this network. in fact, it was previously shown that transposon-insertion mutants of tres/rv display a slow-growth phenotype. as noted before, our 4-state hmm is also capable of detecting regions with unexpectedly high read-counts that might confer growth-advantages to the organism when disrupted. the most notable growth-advantage regions are shown in table. one region of the genome that stands out is the pdim locus, rv2930-rv this locus contains genes involved in the biosynthesis of phthiocerol dimycocerosate, including fadd and ppsabcde. in addition, other genes outside this locus believed to be involved in pdim biosynthesis, like papa and mas, are identified as well. these genes contain read counts well above the global average. fadd itself has a mean read count of more than three times the average throughout the genome. ppsde had a mean read count of and ppsabc a mean read count of pdim is a cell-wall associated glycolipid that modulates the immune response in the host. although it is required for virulence, it is not required for survival in vitro. in fact, biosynthesis of pdim requires resources and imparts a metabolic cost, hence disruption of this pathway is advantageous to cells. due to the increased metabolic cost, it is widely observed that m. tuberculosis stocks maintained in the lab frequently lose the ability to synthesize pdim via acquisitions of mutations in these genes, often leading larger colony sizes. this growth advantage and consequent selection effect likely explains why clones with transposon insertions in the pdim locus are over-represented in the library. regions labeled as'growth-advantage’ by the hmm, in which there is an excess of reads, suggesting an over-representation of transposon insertions. discussion and  the hmm described in this work enables the characterization of essentiality throughout an entire bacterial genome from sequencing data of transposon mutagenesis experiments. although several computational methods have previously been proposed for analyzing tn-seq data, including some based on presence/absence of reads as well as non-parametric models that take quantitative read counts into consideration, an hmm provides several advantages over these methods. for example, an hmm provides a smoothing over adjacent sites that couples them together to help disambiguate the interpretation of read counts at individuals sites. another advantage of using the hmm is that it is not restricted to annotated gene boundaries, and can identify independent regulatory regions, non-coding rnas, and protein domains that are required for survival. while methods that depend on a sliding-window are also capable of assessing essentiality over the entire genome by locally averaging over adjacent ta sites, an hmm formalizes this process in a statistically rigorous way. in addition, by assessing essentiality among regions in the genome, the hmm can also tolerate insertions in the n- and c- termini of genes, without the need of discarding insertions at these locations in an ad-hoc manner as some methods have done previously. a potential limitation of our method is that it does not take into consideration the doubling-rate or expansion time of the library when estimating the parameters of the model. this can affect the under- and over-representation of mutants, and therefore the number of reads these genes will contain in the sequence data. because the hmm depends on the read counts for individual genes, it may be susceptible to libraries that are constructed in different ways. indeed, the ability to detect essential genomic regions from transposon-insertion sequencing data is highly dependent on the quality of the dataset. in practice, tn-seq datasets can span a range from hundreds of thousands to millions of reads, but below some point, there are not enough reads to discriminate essential regions confidently. similarly, well-saturated libraries can have insertions at > 50% of ta sites, but other datasets are more sparse, again increasing the difficulty in distinguishing essential regions. additionally, some sites may contain reads that are orders-of-magnitude larger due to pcr amplification or the development of “hotspots” due to the interactions of the transposon and the organism’s replication machinery. these problems could be alleviated by comparing cultures before and after the passage in the media or comparing datasets derived using different transposons with specificity to different sites. while not directly filtered by our hmm, we showed that the parameter estimation equations we propose work on a wide range of real datasets. in particular we show that they work on dense datasets, as well as sparse ones, and even on artificial datasets down to 25% insertion density. in all cases, the hmm is stable in that it outputs about the same proportion of essential regions, so as the volume of data decreases, the hmm adapts its predictions and becomes increasingly conservative. one of the major advantages of using an hmm to analyze transposon mutagenesis data is that additional states may be introduced to capture distinct types of genomic regions. in this work we add states to capture regions whose disruption leads to a growth defect or a growth advantage. van opijnen et al. have shown that relative abundance of insertions in a gene can be correlated quantitatively with growth rate, though it depends on the number of generations the culture is allowed to grow. eventually, clones with insertions in growth-defect genes will be depleted from the library due to competition. sassetti et al. found substantial reductions in abundance, even after a second round of plating, where hybridization ratios for certain genes dropped from non-essential to essential. while our model only distinguishes one class of growth-defect genes, it could be expanded to more states, discerning finer gradations of growth impairment. a total of growth-defect genes, and growth-advantage genes in m. tuberculosis h37rv were identified by our 4-state hmm, several of which have been shown to be biologically valid, as knockout mutants have been shown to grow slower than the parental strain. identifying these finer distinctions of essentiality can enrich our understanding of the biological roles of genes. for example, we found that icl is labeled as growth-defect. historically, icl has been viewed as essential in m. tuberculosis specifically for growth on fatty-acid substrates, and non-essential otherwise. however, that view is too limiting. icl knock-out strains have in fact been observed to display a growth defect when grown on glucose, and the suppression of reads we observed in icl in the transposon mutagenesis data is consistent with this, suggesting that this gene plays an additional role that is not well-appreciated. availability and requirements the hmm model utilized in this paper is made available online in the tn-hmm package. project name: tn-hmm project home page:http://saclab.tamu.edu/essentiality/hmm/ operating system: linux, windows programming language: python license: creative commons attribution-noncommercial any restrictions to use by non-academics: none competing interests the authors declare that they have no competing interests. authors’ contributions tri and mad were jointly responsible for the development of the algorithm and the writing of the manuscript. both authors read and approved the final manuscript. supplementary material additional file 1 table s excel spreadsheet containing a list of growth-defect regions identified by the hmm in m. tuberculosis h37rv grown on glycerol. click here for file additional file 2 table s excel spreadsheet containing a list of growth-advantage regions identified by the hmm in m. tuberculosis h37rv grown on glycerol. click here for file acknowledgements this work was supported in part by the nih grant u ai.
we also show that it detects growth-defect and growth-advantage genes previously shown to impair or enhance growth when disrupted. a 4-state hmm provides an improved way of analyzing tn-seq data and assessing different levels of essentiality that enables not only the characterization of essential and non-essential genes, but also genes whose disruption leads to impairment of growth. we evaluate the performance of a 4-state hmm on a sequence dataset of m. tuberculosis transposon mutants. we propose a novel method for analyzing sequence data from transposon mutant libraries using a hidden markov model, along with formulas to adapt the parameters of the model to different datasets for robustness. this approach allows for the clustering of insertion sites into distinct regions of essentiality across the entire genome in a statistically rigorous manner, while also allowing for the detection of growth-defect and growth-advantage regions. we show that the hmm produces results that are highly correlated with previous assignments of essentiality for this organism.
 knowledge of which genes are essential to the survival of an organism is critical to understanding the function of genes, and for the identification of potential drug targets for antimicrobial treatment. previous statistical methods for assessing essentiality based on sequencing of tranposon libraries have usually limited their assessment to strict'essential’ or'non-essential’ categories. however, this binary view of essentiality does not accurately represent the more nuanced ways the growth of an organism might be affected by the disruption of its genes. in addition, these methods often limit their analysis to open-reading frames. we propose a novel method for analyzing sequence data from transposon mutant libraries using a hidden markov model, along with formulas to adapt the parameters of the model to different datasets for robustness. this approach allows for the clustering of insertion sites into distinct regions of essentiality across the entire genome in a statistically rigorous manner, while also allowing for the detection of growth-defect and growth-advantage regions.  we evaluate the performance of a 4-state hmm on a sequence dataset of m. tuberculosis transposon mutants. we also test the hmm on several synthetic datasets representing different levels of transposon insertion density and sequence coverage. we show that the hmm produces results that are highly correlated with previous assignments of essentiality for this organism. we also show that it detects growth-defect and growth-advantage genes previously shown to impair or enhance growth when disrupted.  a 4-state hmm provides an improved way of analyzing tn-seq data and assessing different levels of essentiality that enables not only the characterization of essential and non-essential genes, but also genes whose disruption leads to impairment of growth. next-generation sequencingsequence analysishidden markov modelsessentiality 