The NRF Module in Genomic Island Identification

Novel regions have the potential to affect virulence and niche specificity of pathogenic microorganisms. In 2000, Perna et al. [34] published the complete genomic sequence of the pathogenic E. coli O157:H7 strain EDL933 and compared it to the previously sequenced, non-pathogenic laboratory E. coli K12 strain MG1655. Genomic regions found in E. coli O157:H7 EDL933 but not in K12 were called "O-Islands" and genomic regions present in K12 but not EDL933 were termed "K-Islands". Perna et al. found 177 O-Islands greater than 50 bp in E. coli O157:H7 EDL933, constituting 1.34 Mb and 234 K-Islands greater than 50 bp in the K12 strain MG1655, representing 0.53 Mb. They determined the presence of these islands using a custom modification of the MUMmer program [12] and found that many of the genomic regions in strain EDL933 were bacteriophage-related and suggested that they may play a role in the virulence of the organism.

Using Panseq, we re-analyzed the genome sequence data used in these experiments with the NRF module, checking "Unique among the sequences selected" to generate both sets of genomic islands in a single step. Perna et al. [34] defined islands as regions with less than 90% sequence identity over 90% of the sequence length. We found that genomic islands identified in this previous analyses can extend into conserved regions, combining multiple islands into a single contiguous sequence that is interspersed with regions of high sequence identity. Similarly, genomic islands may be misclassified as core regions simply because of heterogeneous composition with interspersed core and accessory sequences in regions of high heterogeneity. Panseq uses parameters optimized for high-resolution comparison of genomic sequences and we found that it more stringently defined islands unique to a genome sequence.

With default program settings (Nucmer values of b:200, c:50, d:0.12, g:100, l:20) and a minimum novel region size of 50 bp, Panseq identified 214 K- and 304 O-Islands, compared with the 234 K- and 177 O-Islands of Perna et al. [34]; a detailed comparison of the findings can be found in Additional File 2. As a result of differences in the stringency between the two methods, present but heterogeneous islands such as K-Island #4 were not identified as islands by Panseq; K-island #4 matches at 90% sequence identity for 84 bp at positions 111429 - 11511 in K-12 and 116035 - 116118 in EDL933. This represents an example of modest differences in core sequence rather than of "novel" regions that comprise the accessory genome (see also below the discussion of identity thresholds). Additionally, islands identified by Perna et al. with interspersed regions of high sequence identity were split and refined into separate islands by Panseq; for example O-island #7 was split into four unique islands, eliminating the regions of high sequence identity. These results demonstrate that Panseq can correctly identify and rapidly extract novel genomic sequences.

Two-way comparisons similar to those used for E. coli strains K-12 and EDL933 have been used to identify genetic attributes which confer distinct phenotypes on many other taxonomically related strain pairs. This type of analysis has allowed researchers to identify elements unique to each strain and suggest putative functions for these elements in the life cycles of the respective organisms, their virulence and their ability to survive in different niches. Examples of such studies include the identification of genetic differences between the human-restricted Salmonella enterica subspecies enterica serovar Typhi, and Salmonella enterica subspecies enterica serovar Typhimurium, which is a murine pathogen but is not host-restricted [35]; and the differences between the host-restricted causative agent of plague, Yersinia pestis and the enteric pathogen Yersinia pseudotuberculosis[36]. While these types of two strain comparison studies have been extremely insightful, studies using programs such as Panseq will greatly facilitate these comparisons and also allow comparisons between multi-strain groups.

The NRF Module in Multiple Sequence Comparisons

To demonstrate a comparison of recently sequenced genomes to previously completed "reference" genomic sequences, we used the two draft genome sequences of Listeria monocytogenes F6900 and 10403S and compared them to the four complete L. monocytogenes genomes in Genbank: Clip81459, EGD-e, HCC23 and 4bF2365 (Table 1). The NRF module found 45 novel regions ≥500 bp, constituting 126858 bp of genomic DNA not present in the four reference L. monocytogenes genomes. The size distribution of the novel regions from the output file is presented in Figure 1. In addition to the summary file, the novel regions are output in FASTA format with sequence location and size information found in the header. These files are suitable for additional bioinformatic or phylogenetic analysis. In addition to the FASTA file, an SVG graphics file showing the location of the novel regions in each strain is optionally provided. As can be seen in Figure 1, the greatest number of novel regions was found in the 1000 -2000 bp range, with 38 of the 45 regions less than 4000 bp in size.

RefSeq Accession No.	Strain	Genome Status
NC_002973	4b F2365	Complete
NC_003210	EGD-e	Complete
NC_011660	HCC23	Complete
NC_012488	Clip81459	Complete
NZ_AARU02	F6900	Draft
NZ_AARZ02	10403S	Draft

Experimental Confirmation of Outputs from the NRF Module

The CAGF Module in SNP Analysis of the Core Genome

Accessing the core genome is important for phylogenetic studies, which have used "core" gene concatenates ranging in size from a few genes in multi-locus sequence typing (MLST) schemes for Campylobacter jejuni[37] and L. monocytogenes[38] to all known "core" genes in E. coli O157:H7 [39]. In addition to offering the best available data for assessing the phylogenetic reconstruction of the evolutionary history of an organism, a small number of single nucleotide polymorphisms (SNPs) can be useful in defining clusters of epidemiologically related strains [40–42].

We analyzed the six L. monocytogenes strains in Table 1 with the CAGF module of Panseq, using a fragmentation size of 500 bp and a sequence identity cutoff of 85%. The resulting concatenated core data, which include all conserved nucleotides and SNPs, were used to construct a maximum parsimony (MP) tree using Phylip v3.69 [43]. A MP tree was also generated in silico for the same six strains using the L. monocytogenes MLST protocol outlined by Nightingale et al. [38], for comparative purposes (Figure 3). The symmetric distance between the two trees shows that the overall topology of the trees differ between that created from MLST data and that based on the entire concatenated core. This is likely due to the fact that the MLST protocol only considers seven genes, where disproportionate variation among these few loci and the relative paucity of loci compared to that of the entire core genome can reduce the ability of this method to capture the overall relationships among strains. With the continued increase in sequencing throughput, settling for rough approximations of true tree topologies may no longer be necessary.

The CAGF Module in the Analysis of the Accessory Genome

While SNP analysis has proven to be an extremely useful tool, the ideal reconstruction of an evolutionary history would take into account not only the heterogeneity among all core genomic regions, but also the presence or absence of regions in the accessory gene pool, especially since the accessory genes can directly affect phenotype (e.g. niche specificity, antimicrobial resistance, virulence, etc.).

Traditional methods of estimating phylogenies rely on SNPs within the core genome, but it has been shown for C. jejuni[44], E. coli O157:H7 [45] and Streptococcus pneumoniae[46] that the distribution of accessory genes provides a very similar overall tree topology to methods based on variability in the core genome. Although variation in the core genome among a small number of loci may be sufficient for identifying clusters of related strains, discrimination among strains is often more difficult because there is less variation and fewer phylogenetically informative loci, leading to fewer genotypes. Accessory genome content, which can be highly variably among strains, appears in many cases to be consistent with phylogenetic analyses of core genes and as a result of greater variability provides a higher degree of discrimination among strains. Analysis of accessory gene content thus presents an opportunity to integrate valuable links to phenotype into a genetic classification scheme.

To examine the performance of accessory genome information for phylogenetic reconstruction, we used the binary presence/absence data of the accessory genome computed for the set of six L. monocytogenes strains above, and constructed a MP tree (Figure 3). When comparing the tree topology of this accessory-based tree to that of the SNP-based core tree and the MLST tree, we found that the core- and accessory-based MP trees had an identical topology, but differed from the MLST-based tree in the placement of the two strains 10403S and F6900. Further studies will be required to determine the extent of the phylogenetic concordance between the core genome and the accessory genome among other bacterial groups.

The CAGF Module in the Examination of Pan-Genomic Differences

As well as producing a concatenated core, the CAGF module also produces a table listing each SNP and its location and allele within each original sequence. This can be useful in examining individual differences or for hierarchical clustering of data. To demonstrate how the results of Panseq can be visualized, both tabular output files from the Listeria core/accessory analysis were used to create hierarchical clustering dendrograms: one based on the SNP character data from the core regions (Figure 4) and the other based on the binary presence/absence data of the accessory regions (Figure 5). Both dendrograms have the same tree topology, and the underlying data for both core and accessory regions shows only the differences among the genomic sequences, making comparisons between strains clearer than they might be from whole-genome comparisons where conserved as well as variable loci are considered in the analysis.

The LS Module in the Comparison of Randomly Generated Data

The approach of the LS module in Panseq is to iteratively build the final loci set, including only the loci that produce the most unique fingerprints, and offer the most variability among input sequences at each step. This allows it to efficiently examine datasets that would be computationally prohibitive if all possible combinations were considered.

To test the Panseq LS module, we created a random binary dataset of 100 loci among 10 sequences in Microsoft Excel (Additional File 3). We subsequently ran Panseq with the 'best' option and found that 4 loci generated a unique profile for all 10 sequences, and that the four loci: locus84, locus40, locus79 and locus29 provided the maximum possible discrimination of 100 POD for the dataset. While these loci cannot be guaranteed to be the only four loci to generate the same results, they will always be one of the most discriminatory sets.

We then ran the same dataset through a Perl script that generated all possible 4-loci combinations, outputting those that produced 10 unique sequence profiles (Additional File 4). We found that 99132 unique combinations of 4 loci from the dataset yielded 10 unique sequence profiles. Evaluating every possible combination required 449 s on a computer running Ubuntu 9.10 with 3.6 GB of available RAM and AMD Phenom 8450 triple-core processors. Panseq was able to sort through the data to generate a single group of loci that contained not only the most unique profiles, but that provided the most discrimination among the input sequences, completing the task in one s.

The LS Module in the Analysis of SNP Genotyping Data

Any variable character data can be used as inputs in the LS module. To illustrate the functionality of the LS module, we used a set of 96 SNPs identified from E. coli O157:H7 by Manning et al. [47], for which the nucleotide value of each SNP was determined among a set of 17 E. coli O157:H7 genomic sequences [48]. This dataset was first modified to represent any unknown character with the '?' symbol, and to replace any locus described by two characters with a single character (Additional File 5).

Advantages of Panseq over Other Related Programs

While many sequence analysis programs exists, some with overlapping capabilities, no two are identical with respect to the tasks they perform. With enough time and knowledge one can parse the output of a sequence alignment program such as BLAST or MAUVE manually, but there is a considerable time savings and ease of use with a program such as Panseq that automates the process. Panseq is unique in its single-step novel region finding options among groups and specific to individual sequences, which other in silico subtractive programs such as mGenomeSubtractor [22] and nWayComp [23] lack. Panseq also provides a comprehensive analysis of the pan-genome, automatically generating analyses that can only be partially accomplished by any other single program; eg. With MAUVE [49] a list of SNPs can be generated and a display of the similarities/differences between the genomes is produced, but the underlying nucleotide sequence is not automatically extracted. With Panseq, the underlying sequence data is automatically extracted, the segments compared among all sequences and presented in tabular form and input files for phylogenetic programs are automatically created. Further, the SNP table or binary presence/absence table of the accessory genome can be used directly in the LS module, for the selection of the most discriminatory loci.