DNA sequence data
Table lists the 12 H. influenzae
clinical strains and the reference strain Rd, a largely non-pathogenic strain, used in the comparative genomic studies described herein, their NCBI locus tags, the location where the sequencing was performed, and their clinical origins. Nine of the clinical strains were sequenced using 454 LifeSciences novel pyrosequencing technology [25
]. The number of sequencing runs, the extent of genomic coverage, and the number of contigs resulting from first and in some cases second pass assemblies are tabulated (Table ).
Bacterial strains and sources used for whole genome sequencing, comparative genomics, and computation of the NTHi core and supragenomes
Sequencing data for the 9 Nthi strains sequenced with 454-technology
Determination of gene clustering parameters
Gene clustering parameters for the grouping of homologs were empirically determined by minimizing the change in the number of clusters per change in the parameters (Figure ). We hypothesize that this minimum point coincides with the best estimate threshold for distinguishing true orthologs from functionally distinct homologs. Some homologs will be more similar than 70%, while some orthologs will be more divergent than 70%, but as a uniform criterion, the threshold is optimized. Visual inspection of the clusters reveals that most clusters are reasonable. Mosaic genes were particularly difficult to cluster due to high levels of rearrangement. In the remainder of the paper, genes in the same cluster are considered to be the same gene.
Figure 1 A plot of the total number of clusters as a function of clustering parameters shows an inflection point near 0.65 identity and 0.70 match length. The inflection, which minimizes the rate of change in the number of clusters per change in parameters, suggests (more ...)
Enumeration of gene clusters and genic relationships among NTHi strains
We identified 2,786 gene clusters among the 13 strains (Table ). Of these, 52% were found in every strain (core genes) and 19% were found in only a single strain (unique genes). The remaining 29% of genes were found in some combination of two or more strains, but not all (distributed genes; Figure ). The number of clusters found per strain varied from 1,686 in PittEE to 1,878 in PittII (Table ). All strains possessed some unique genes not seen in any of the other strains. A pair-wise comparison was performed among all possible strain pairs, which determined the mean number of genic differences between any two strains was 395 with a standard deviation of 94 (Figure ). This analysis also identified minimal and maximal genic differences of 81 and 577, respectively, for the strain pairs 2866:PittII and 2866:PittAA. The number of coding sequences identified per genome by AMIgene did not correlate strongly with genome size. This is likely due to the presence of split open reading frames (ORFs) in the 454 sequenced genomes as an analysis of the 4 completed genomes showed a linear relationship between gene number and genome size with an R2 = 0.910. In contrast, the correlation between total gene clusters and genome size is 0.86, implying that the number of distinct genes found on the genome is linearly related to the genome size.
Figure 2 A histogram of gene clusters observed in exactly N of 13 H. influenzae strains compared to the expected number of genes estimated by the supragenome model (trained on all 13 strains). Over 1,400 genes were observed in all 13 strains, indicating that there (more ...)
Gene identification and clustering results
Figure 3 A pairwise genic comparison of 12 NTHi strains of H. influenzae and the reference strain Rd KW20. The comparison of two strains is found at the intersection of the row and column corresponding to the respective strains. Strains are compared based on the (more ...)
A dendrogram based on non-core genic differences (Figure ) demonstrates the diversity in the NTHi population. A typical strain differs from its nearest neighbor by more than 200 genes. The strains collected from otitis media with effusion (OME) patients at Children's Hospital in Pittsburgh (designated as Pitt strains) show that a genetically diverse population can be isolated contemporaneously from a single geographic location from patients with similar indications. In contrast, two pairs of strains, PittEE/R2846 and PittII/R2866 are relatively similar despite geographically distinct points of isolation. Interestingly, the laboratory strain Rd KW20 is not an outlier among the clinical strains. For comparison, a maximum likelihood tree was generated using sequence from seven multi-locus sequence typing (MLST) housekeeping genes for the same set of 13 strains (Figure ). The topology of the trees is significantly different, both in terms of pairwise groupings and overall structure.
Figure 4 Plotting of relationships among the sequenced NTHi strains by gene sharing and multi-locus sequence typing. (a) A dendrogram based on genic differences among the 13 strains of H. influenzae. While several pairs of strains appear to be closely related, (more ...)
The identified number of new genes and core genes found per addition of each genome (as determined by incremental clustering of the 13 strains) shows an exponentially decaying trend in both cases (Figures and ). Qualitative inspection suggests a diminishing return on new genes found in future sequences, though it is expected that approximately 40 new gene clusters will be found in each of the next few genomes that are sequenced. The number of core genes appears to trend towards a horizontal asymptote near 1,450 genes. A quantitative analysis of these results is developed below in the section 'Mathematical development of a finite supragenome model'.
Figure 5 The expected number of total gene clusters and core gene clusters identified at the addition of each genome to the clustering dataset. Modeling predictions are based on the eight strain training set (see 'Mathematical development of a finite supragenome (more ...)
The observed and expected number of new gene clusters found at the addition of each genome to the clustering dataset. Modeling predictions are based on the eight strain training set (see 'Mathematical development of a finite supragenome model').
Whole genome alignments reinforce the great diversity observed among gene clusters
Whole genome alignments were generated between Rd and each of the 12 clinical strains to quantify genomic insertions and deletions independently of gene identification (Table ). On average, each of the clinical strains had 127 genomic insertions (>90 base-pairs (bp) in length) that did not correspond to any Rd KW20 sequence. Similarly, each clinical strain contained, on average, 147 genomic deletions (>90 bp) when compared to the Rd KW20 strain. The average total length of non-matching sequences between the 12 clinical strains and Rd was 321 kb, approximately 18% of the genome. The quantity of non-matching sequences reasonably accounts for the average of 390 genic differences between strain pairs. Figure shows a genomic region in which two different forms of an insert, homologous to the plasmid ICEhin, have integrated into the same site of two different genomes, but which is wholly absent from the other strains in the alignment. Similarly, a 40 kb contiguous region in Rd shows extensive deletional diversity among seven of the clinical strains, with only two of the clinical strains demonstrating the same local genomic organization (Figure ). Interestingly, the two strains, PittAA and PittEE, that are similar in this region are highly divergent overall (Figure ). Genic diversity also exists on a smaller scale. Figure displays a 20 kb region from 7 clinical strains that shows 5 different combinations of possession and loss of the lic2C gene, the NTHI0683 gene, and the UreABCEFGH operon.
Analysis of inserted and deleted Sequence in 12 strains with respect to Rd KW20
Figure 7 A multi-sequence alignment using 86-028NP as a reference shows varying degrees of homology among 6 strains to a 50 kb region homologous to the plasmid ICEhin1056. The plasmid is integrated in 86-028NP and is partially present in R2866, but absent from (more ...)
Figure 8 A 40 kb region present in Rd KW20 shows two blocks of genomic variation among other strains. The upstream block is bounded on the right by a frame-shifted insertion sequence (IS) element (HI1018). The downstream block (HI1024-HI1032) includes genes with (more ...)
Figure 9 A 20 kb region that demonstrates strain diversity at the level of an individual gene (lic2C), a pair of genes (NTHi0683/4), and a group of seven functionally related genes (urease system). 86-028NP is used as a reference for the alignment, and sequence (more ...)
Global genomic alignments of PittEE against R2846 and R2866 were performed (Figures and ). PittEE and R2846 are very similar at the global level and this is reinforced by the gene cluster analysis, which revealed only 96 genic differences. In contrast, R2866 has a large inversion and several large insertions and deletions with respect to PittEE. This diversity at the global level corresponds to the 377 genic differences identified between these two strains by cluster analysis (Figure ). Global alignments were not visualized for most strains since the ordering of the contigs had not been determined.
Figure 10 A global alignment of R2846 and PittEE as visualized by Mummerplot. A point is placed at the (x,y) coordinate if the x-coordinate of R2846 matches the y-coordinate of PittEE. Green matches indicate a reverse complement match. It can be seen that PittEE (more ...)
Global alignment of R2866 and PittEE shows a large inversion and several regions unique to each strain. The strains are similar across the majority of the genome; however, there is one large inversion as well as several regions unique to each strain.
Codon usage analysis
The codon usage of each gene cluster was compared to the typical H. influenzae
codon usage pattern by the epsilon-score calculated by CodeSquare [26
]. A low epsilon score indicates that a gene's codon usage is similar to typical patterns of the organism, while a high score indicates atypical codon usage. Since the epsilon score is partially dependent on the length of a coding sequence, all scores were normalized by length. The average normalized score is 0 and low values continue to indicate typical codon usage. Figure is a scatter plot of the normalized epsilon scores versus the number of strains in which the gene was found. The range of normalized epsilon values is similar for core, distributed, and unique genes, though the median values are slightly higher for distributed and unique genes (Tables and ). The Mann Whitney U-test was employed to determine the significance of this difference. To eliminate any remaining length bias, only genes with lengths of 200-300 amino acids were analyzed. The median normalized-epsilon value of core genes is significantly smaller than the medians of distributed and unique genes, and as a consequence, these non-core genes are more likely to have foreign origins. Interestingly, there is no significant difference between distributed and unique genes and most of these non-core genes display typical H. influenzae
Figure 12 Codon usage of genes is quantified by a normalized epsilon score . Low epsilon scores indicate that a gene's codon usage is similar to the typical H. influenzae codon usage pattern. The range of epsilon scores is similar for all three classes of genes: (more ...)
Codon usage comparisons of core, contingency and unique genes
Codon usage comparison of core, contingency and unique genes
Phage homology analysis
Phage insertion is a common origin of genomic diversity. The influence of phage was quantified by a homology search between all gene clusters and the NCBI NT database. A gene cluster was said to be 'phage associated' if one of the top ten significant matches was annotated as a sequence of phage origin. Overall, 9.3% of gene clusters were phage associated. The distribution of these genes is not uniform among core and non-core genes. Only 0.3% of core genes were phage associated, while 14.6% and 25.8% of distributed and unique genes, respectively, were phage associated (Table ).
Percentage of genes with probable phage origin per category
Development of a finite supragenome model
The comparative genomic data presented above are supportive of the DGH and reinforces the concept that, at the species level, there is an H. influenzae
supragenome that is much larger than the genome of any single individual strain, and hence many strains must be sequenced to generate an accurate picture of the species supragenome. Among the questions we may ask about the supragenome, the most obvious is, how many strains must be sequenced to observe the entire (or nearly all) of the supragenome?. The problem is similar to determining the read coverage necessary to sequence an entire individual genome using a random shotgun library approach. Lander-Waterman statistics provide an answer in the latter case by using the assumption that reads are independently and randomly sampled from the genome with equal probability. Previously, Tettelin et al
] developed a supragenome model for S. agalactiae
that, like Lander-Waterman statistics, is based on the assumption that contingency genes are independently sampled from the supragenome with equal probability, except in the case of rare genes, which are modeled as unique events that appear only once in the entire global population. The model requires four parameters: the number of core genes, the number of contingency genes, the probability of finding a contingency gene, and the expected number of 'unique' genes found per strain. This model predicted that the supragenome of S. agalactiae
is infinite in size (that is, the expected number of unique genes found in each strain is non-zero). While the model is an insightful attack on the problem, we question the assumption that contingency genes are sampled in the population with equal probability. It is important to compare the existing model against a new model that does not rely on this assumption.
The Supragenome is represented here by a generative model that emits genomes according to a set of probabilistic rules. The supragenome contains N genes that are modeled as Bernoulli random variables with 'success' probabilities that correspond to the population frequency of each gene. A genome is generated by observing the Bernoulli variables: a gene is present if the corresponding trial is a success and otherwise absent. Each gene variable is assumed to be independent of all other genes. This assumption is sometimes violated in real H. influenzae genomes. For example, genomic islands are sets of genes that are not independent. However, we proceed with this assumption since it significantly reduces the complexity of the model and is reasonable in many cases.
The true population frequencies are, in general, unknown. Therefore, population frequencies are also treated in a probabilistic fashion. It is assumed that there are K discrete classes of genes. Each class k has an associated population frequency, μk. All genes in class k will have population frequency μk. Each of the N genes is assigned to a class according to a probability distribution given by the vector π, where πk is the probability that a gene is assigned to class k. Conceptually, πk is the percentage of genes in the supragenome that have population frequency μk. The assignment of a gene to a class is independent of all other gene assignments.
The complete model is depicted in plate notation in Figure . 'Z' is the hidden class variable in which zn corresponds to the class of gene n. 'X' is the observed gene variable, where xn,s corresponds to the presence or absence of gene n in strain s. The outer plate represents the supragenome, while the inner plate represents instances of specific genomes. The model requires 2 × K + 2 parameters: N, K, a mixture coefficient πk for each class, and a Bernoulli probability μk for each class. The number of gene classes, K, and their associated Bernoulli probabilities, μk, are fixed in advance. Care must be taken to choose classes that represent low and high population frequencies. Seven classes were selected for this study (K = 7) with associated probabilities μ = <0.01, 0.1, 0.3, 0.5, 0.7, 0.9, 1.0>. The class with probability 1.00 represents 'core' genes that appear in all strains.
Figure 13 A plate diagram of the H. influenzae supragenome model. Each node in the diagram represents a random variable, and the arrows indicate dependence between the variables. Independent, identically distributed (IID) nodes appear in boxes with an index listed (more ...)
The remaining parameters, N and πk, are selected under a maximum likelihood scheme. Suppose that |S| genomes have been sequenced and a particular gene from class k was observed in n of the |S| strains. The probability of this observation is given by a binomial probability since this result is the sum of independent Bernoulli variables. As a function of πk and N, the probability is given by:
However, we do not know the true gene class, so we must consider a mixture of binomial probabilities:
Now consider the complete set of genes. Let c = <c0, c1, ..., cS>, where cn is the number of genes observed that appear in exactly n of |S| strains. The probability of the total observation is given by a multinomial distribution:
The parameters N and π can be determined by maximizing the log-likelihood of the observation c:
The log-likelihood function was maximized by fixing N and maximizing with respect to π. The maximization was performed using the MATLAB function fmincon with the constraint:
and requiring that the coefficients are between 0 and 1. The maximization was performed for values of N starting at the minimum possible value (the number of genes actually observed) to 6,000. The combination of N and π that maximized the overall log-likelihood was selected as the best parameter estimate.
Supragenome modeling validation and results
The model was validated by training the supragenome parameters using only the first 8 sequenced genomes and comparing the predictions with the observed results for 13 strains. The maximum likelihood number of genes was 3,078. Of these genes, 1,423 are core genes, 417 are contingency genes with population frequency >0.1, and 1,238 are contingency genes with 0.1 population frequency. No genes were predicted in the 0.01 population frequency class. Predictions for the 0.01 class may be inaccurate due to the small sample of 8 genomes. The 1/100 maximum likelihood confidence interval for total genes ranged from 2,975 to 3,681. Figure shows the distribution of the genes among the seven classes.
Figure 14 The distribution of genes among gene classes in the supragenome model trained on 8 or 13 strains. The only significant difference occurs in the rare gene categories with frequency 0.01 and 0.10. A small sample of eight strains is not expected to generate (more ...)
Figure compares model predictions based on 8 strains to actual observations of core genes (shared among the first N strains) and total genes found after sequencing the 9th through 13th strains. In both cases the model predictions follow the observed trends. Figure compares predictions to observations of the number of new genes found in the Nth sequenced strain. Again the model predictions follow the observed trend. Figure compares the best-fit gene distribution (based on 8 strain models) to the observed distribution of genes found in exactly N of 13 strains. Overall, the predicted trends follow the observed distribution; however, the predictions were too low for genes seen in 1 of 13 strains, and too high for genes seen in 2 of 13 strains. This bias may be due to the small sample size (eight strains) used to train the supragenome model. Predictions for genes seen in four to seven strains were also somewhat lower than observed.
The supragenome model predicted an average of 1,776 genes per strain with a standard deviation of 14 genes. Of the 13 strains, the average number of genes was 1,793 with a standard deviation of 62 genes. The model predicted an average of 373 different genes when comparing any two strains with a standard deviation of 17 genes. Among the 13 sequenced strains, the average was 395 with a standard deviation of 91 genes. In both cases the model predication for average was reasonable, while the standard deviation was underestimated by about four-fold. This suggests that the assumptions used for the supragenome model may omit important sources of variation. Genomic islands and other genes that appear together in the genome likely contribute to the total variance.
Altogether, the above results show that the supragenome model generates reasonable predictions for the average properties of the supragenome. To obtain improved predictions, the model was re-trained on all 13 strains. The supragenome class distribution for the extended model is shown in Figure . The results are similar to the model trained on 8 strains, except that the class with population frequency 0.01 is now predicted to contain 2,609 genes, while the 0.10 frequency class was reduced in size to 590 genes. This large change is due to improved resolution of rare genes in the 13 strain training set. The model now predicts 5,230 genes, with a 1/100 likelihood interval ranging from 4,425 to 6,052 (Table ). Nearly all of the increase over the eight strain model is due to the class of rarest genes. Of these genes, 1,437 are core genes, 594 are contingency genes with population frequency >0.1, and 3,199 genes are rare contingency genes with population frequency <0.1. Figures and show the prediction trends for total, core, and new genes observed after sequencing N strains (up to 30 strains).
Maximum likelyhood estimate for size of supragenome and 1/100 likelihood intervals based on 8 and 13 strain training sets
Figure 15 A theoretical plot of the number of new genes expected to be found in the Nth genome for future H. influenzae sequencing projects. The plot was generated using strains isolated in North America, and the extrapolation may not hold for isolates from other (more ...)
Figure 16 A theoretical plot of the number of total genes and core genes expected among N sequenced H. influenzae genomes for future sequencing projects. The extrapolation may not hold for strains isolated outside of North America since the plot was constructed (more ...)