BACKGROUND
bacteriophages represent a numerically vast, highly dynamic, evolutionarily ancient, and genetically highly diverse population  <cit> .
phage genomes are typically small compared to those of their bacterial hosts  and no longer present significant technical challenges to sequence determination  <cit> .
as genomic sequencing approaches get simpler and cheaper, the availability of individual phage isolates for characterization becomes limiting, a need that can be effectively addressed through integrated research-education programs involving undergraduate and high school student investigators  <cit> .
in spite of their relatively small size, phages present significant challenges to accurate genome annotation including gene identification.
two principal issues arise.
first, phage genes tend on average to be small , approximately two-thirds the average size of bacterial genes  <cit> .
many of the genes required for virion structure and assembly are relatively large , but those in the non-structural genomic segments are small, often shorter than  <dig> codons.
secondly, phage genomes are replete with genes of unknown function for which no homologues have been described  <cit> .
mosaic architectures are hallmarks of phage genomes, and individual phages can be considered as particular combinations of interchangeable modules, each of which can be present in two or more different genomic contexts  <cit> .
in some cases, where the recombination events giving rise to these mosaic structures occurred relatively recently in evolutionary time, mosaicism is apparent through nucleotide sequence comparisons  <cit> .
when the events occurred in more remote evolutionary times the evidence of common ancestry is usually no longer apparent at the nucleotide level, but often can be revealed from comparison of the predicted amino acid sequences  <cit> .
such comparisons reveal that individual phage genomes are typically constructed from multiple modules - often corresponding to single genes - each of which has a distinctly different phylogeny  <cit> .
as such, accurate compilations of whole genome phylogenies that reflect the evolutionary history of the entire genome are not possible, and reticulate-based representations are needed to capture this evolutionary complexity  <cit> .
the mechanisms giving rise to genome mosaicism are unclear but must accommodate the striking observation that module boundaries correspond closely with gene boundaries, and in some cases, domain boundaries  <cit> .
one model invokes homologous recombination events targeted to short conserved boundary sequences between genes, and there is evidence for this in some phage genomes  <cit> .
however, there are numerous examples where no conserved boundary sequences are evident, raising the possibility that mosaicism results largely from illegitimate recombination events between randomly chosen partners sharing little or no sequence identity  <cit> .
in this second model, correspondence between module and gene boundaries results from the selection for gene function, not from targeting of the recombination events  <cit> .
comparison of genomes from phages that infect taxonomically diverse hosts typically provides little information into their evolution because of only very limited similarity at either the nucleotide or amino acid sequence level  <cit> .
phages of a common host, however, have the advantage that they are more likely to have been in recent genetic communication with each other and to have exchanged modules in recent evolutionary times  <cit> .
large sets of phage genomes are now available for several hosts including burkholderia  <cit> , bacillus  <cit> , enterobacteriaceae  <cit> , mycobacteria  <cit> , prochlorococcus and synechococcus  <cit> , pseudomonas  <cit> , and staphylococcus  <cit> , although even these can span enormous genetic diversity  <cit> .
a large number of phages that infect mycobacterium smegmatis mc <dig> have been isolated and a comparative analysis of  <dig> has been described  <cit> .
although these are genetically diverse, the diversity is heterogeneous, and phages can be grouped into 'clusters' according to their overall nucleotide sequence relationships  <cit> .
of the  <dig> published completely sequenced mycobacteriophage genomes,  <dig> can be grouped into ten major clusters, seven of which can be further subdivided into subclusters according to the extent of the nucleotide similarities  <cit> .
five of the genomes have no close relatives and are referred to as 'singletons'  <cit> .
because the currently sequenced mycobacteriophage genomes under-represent the mycobacteriophage population-at-large, these cluster designations will undergo modifications as new genomes are sequenced  <cit> .
there are, however, numerous examples of genes that are shared between phages of different clusters and whose common ancestry is only apparent from amino acid sequence similarity  <cit> .
we have proposed previously  <cit>  that mycobacteriophage genes related to each other can be grouped into phamilies  and that mosaic relationships can be analyzed and represented using pham-annotated genomes maps and phamily circles that show the patterns of which phages contain members of particular phams.
although manual or semi-automated approaches are applicable when only small numbers of genomes are analyzed  <cit> , this becomes an impossible task as the number of genomes expands.
we describe here a software program 'phamerator' that provides bioinformatic tools for both analyzing and representing phage genome mosaicism.
the core functionality of phamerator performs pair wise similarity searches between predicted protein products of a set of phage genomes, and assorts them into phamilies  of related sequences.
genome maps can be displayed that illustrate the relationships between phages at both the nucleotide and amino acid sequence level.
moreover, the evolutionary histories of specific genes can be displayed by phamily circles in which all gene members of particular phams are represented, and for which multiple phams can be compared.
we illustrate the utilities of phamerator using a set of  <dig> completely sequenced mycobacteriophage genomes, but the program is applicable to any set of phage genomes for which comparative analysis is desired.
methods
phamerator database architecture
phamerator is written entirely in the python computer programming language and makes use of a number of modules, including the biopython framework for computational biology  <cit> .
biopython provides a programmatic interface for sequence manipulation, the construction and parsing of files in relevant formats, and access to external command line applications such as those used for sequence alignment.
for the latter, phamerator uses biopython to interact with local instances of blastp and clustalw.
due to the significant computational time required for performing large numbers of sequence alignments, phamerator employs a distributed processing model that is implemented using python remote objects .
phamerator was developed on ubuntu linux but should be able to run on any modern unix-based operating system.
phamerator uses the mysql database software with a simple, custom database schema that incorporates and extends the relevant information found in genbank records.
the phage and gene tables are populated with data from genbank files, while the remainder store data relevant to phamerator analysis or imported data from external databases such as the ncbi conserved domain database .
additional tables are used to store current and historical pham assignments, and records are maintained of the splitting or joining of phams that can occur as new sequences are added to the database.
in the event that a new mycobacteriophage protein is added to the database that has similarity to members of more than one existing pham, the phams are merged, their pham names  retired, and a new pham created with a new number.
the new pham contains each of the members of the joined phams in addition to the new protein.
conversely, the addition of new proteins to the database can also invalidate an existing pham because the blastp e values used for determination of pham membership are in part dependent on the size of the data set.
thus, after adding a new genome to the database, if a protein in a pham is no longer related to any members of the pham it is removed from the pham, and it is either placed into an existing pham if that pham includes a related protein, or it becomes an orpham .
use of distributed computing resources
adding genomes to a large phamerator database is time-consuming, with the length of time required being proportional to the size of the existing database and the number of genes in the genome to be added.
however, the length of time required to perform this operation can be reduced by distributing the required calculations to several computers.
the reduction in time required scales approximately linearly with the number of available computers.
communication between computers is achieved using the python remote objects  library, and can be distributed to computers on the local network or across the internet.
availability and distribution
phamerator operates using a client/server model.
a relational database resides on the server computer and is distributed to client computers when they run the phamerator client program.
this enables all clients to have the performance benefits associated with interacting with data stored on the local computer along with a benefit normally associated with processing data on a server--assurance that the data being used is up-to-date.
because the data files being disseminated are relatively small , a modest server can handle the load of hundreds of users without performance deterioration.
phamerator can be downloaded from: http://phamerator.csm.jmu.edu/files/phamerator-current.tar.gz.
further information and installation instructions are provided in additional files  <dig> and  <dig>
RESULTS
rationale for phamerator construction and operations
the pervasive mosaicism of bacteriophage genomes requires bioinformatic tools that can organize and display their complex relationships.
two key questions arise in phage comparative analysis: what are the relationships between genes that are evolutionarily mobile within a given set of phage genomes, and how are they related to genes found in other genomes.
both approaches are complicated by the presence of intragenic mosaicism reflecting distinct evolutionary histories of gene segments  <cit> .
phamerator is a computational tool designed to sort phage genes into phamilies of related sequences using pairwise amino acid sequence comparisons of predicted gene products.
rather than using ortholog identification programs such as hmmer or pfam  <cit>  we have employed blastp and clustalw to perform pairwise comparisons that are then assembled into phamilies of related proteins.
both of these programs are relatively fast computationally, a critical factor when large numbers of computations are required.
for example, in the dataset of  <dig> phage genomes used here, a total of  <dig>  ×  <dig> comparisons must be computed.
we note that the use of efficient pairwise blastp comparisons to generate families of protein sequences has been described previously  <cit> .
to identify homologues of previously identified proteins phamerator performs automated searches of genbank non-redundant protein sequences, as well as searches for conserved domains in the ncbi conserved domain database using the rps-blast tool.
this information can be exported in tabular form, or represented in a whole-genomic context.
details of the phamerator program are described in the methods section and an overview of database structure is shown in figure  <dig>
pham-building parameters
the building of phams is strongly influenced by the specific parameters used for amino acid comparisons.
in early studies initiated prior to phamerator development we used a blastp cutoff value of  <dig>  and a clustalw cutoff of 25% amino acid sequence identity for manual pham building  <cit> .
however, as the number of mycobacteriophage genomes increased and computational processing became essential, it became clear that these largely arbitrarily chosen parameters promoted assembly of many large phams that require time-consuming manual deconvolution  <cit> .
we therefore explored the impact of varying the threshold values for blastp and clustalw on pham assembly.
we first evaluated the effect of changing the threshold for clustalw comparisons.
we varied the threshold between 50% amino acid sequence identity and  <dig> %, and for each level determined the number of phams generated, the size of the largest pham, the number of orphams , the percentage of orphams, and the mean pham size .
these data are informative and provide guidance as to the optimal parameters to use for routine database construction.
in particular, we note that as the threshold for amino acid sequence similarity is made less stringent  there is a reduction in the total number of phams  reflecting the process of pham assembly.
interestingly, this relationship is linear between 50% and  <dig> %, with a reduction in the total pham number of about  <dig> for every percentage of identity that is reduced .
as the percent identity falls below  <dig> % the relationship becomes non-linear, with a progressively greater reduction in the number of phams as the percent identity threshold falls from  <dig> % to  <dig> % .
the number of orphams also reduces as stringency is relaxed, while mean pham size increases as stringency is relaxed, and there are notable changes between values above and below  <dig> % identity .
there is also a dramatic change in the size of the largest pham as the threshold level varies from  <dig> % to  <dig> % .
at  <dig> % the largest pham contains  <dig> members, but increases to  <dig>   <dig> at  <dig> % .
the size of the largest pham is more stable between  <dig> % and 50% identity thresholds and varies from  <dig> to  <dig>  we interpret these data as indicating that between 50% and  <dig> % identity, pham assembly proceeds in a manner that simply reflects the variation in the overall relationships between genes.
however, at levels below  <dig> %, there is an increasing proportion of phams that are more complex, such that not all pairwise matches within the pham are above the threshold level.
one example might be where two genes  have been fused into a single open reading frame , such that although genes a and c, and gene b and c, both surpass the threshold, the unrelated genes a and b do not .
for routine database building purposes where we wish to avoid the assembly of large phams that then warrant subsequent deconvolution, we have chosen to use a  <dig> % clustalw threshold, but note that comparison of phams generated with these parameters and those with lower stringencies should be useful in analyzing intragenic mosaic relationships.
one advantage of clustalw as an assembly program is that the threshold values are independent of gene length.
nonetheless, we predict there are instances where large genes may not exceed the clustalw threshold but are evidently homologues because of statistically informative blastp scores.
we therefore examined the impact of including a blastp search along with the clustalw comparison  and varying the blastp cut off value  .
as threshold values are made less stringent we observe a reduction in the total number of phams and the number of orphams, and a corresponding increase in mean pham size and the size of the largest pham .
between blastp thresholds of 10- <dig> and 10- <dig> these is only a modest change in the total pham number  and number of orphams , but the size of the largest pham changes from  <dig> to  <dig> .
upon manual inspection of all phams with  <dig> or more members at the blastp thresholds between 10- <dig> and 10- <dig> we see 'false' pham assembly occur as illustrated by joining of a phamily of tape measure proteins with a phamily of lysin proteins, through sharing of small but closely related domains.
nonetheless, inclusion of the blastp comparison with a 10- <dig> threshold joins several lysin phams that are separate when using just the clustalw comparison.
we conclude that inclusion of blastp contributes rather little to the pham assembly process, but that a combination of a clustalw threshold of  <dig> % and a blastp cut off of 10- <dig> offers optimal parameters for this dataset, building phams of homologous proteins while minimizing construction of complex phams in which only segments of the proteins are related to each other.
when these parameters are applied to this dataset, the  <dig> mycobacteriophage genomes contain a total of  <dig>   <dig> genes that assemble into  <dig>   <dig> phamilies with a mean size of  <dig>  genes/pham;  <dig>   <dig> phams are orphams  and the largest pham contains  <dig> members.
identification of known homologues and conserved domains
once a novel genome has been sequenced and annotated, questions about the functions of individual genes encoded within the genome can be addressed.
this process is facilitated by analyzing the predicted gene products for the presence of conserved domains.
numerous tools already exist for this purpose, but the ncbi conserved domain database  aggregates many of them into a single, searchable dataset.
these domain databases often use different, complementary techniques such as hidden markov models or position-specific scoring matrices to define domains and for the matching of novel sequences to existing domain models.
to leverage the power of each of these approaches, we have implemented a system in phamerator whereby phage proteins are used to query cdd, and the results are presented in a searchable database browser in addition to being displayed on genome maps.
the latter affords both a global view of the genomes and provides a way to visualize the conservation of specific sequences in the phamerator dataset with those in other organisms.
an analysis of the distribution of cdd hits among the mycobacteriophage protein sequences reveals that a current search produces a total of  <dig>   <dig> matches among the  <dig>   <dig> predicted proteins, for an average of  <dig>  matches per protein.
however, as expected, the distribution of matches is non-random, with only  <dig>   <dig> proteins  having at least one cdd entry match, with the average number of matches for these proteins being  <dig> .
while  <dig> proteins each match a single domain, a single protein, myrna gp <dig> , matches  <dig> domain models, most of which are aminotransferases.
generally, when numerous matches of a single protein to domains in the cdd are found, it reflects the redundant nature of the cdd dataset, although in some cases it reflects the presence of multiple domains within a single protein.
computation of nucleotide sequence similarities
phamerator uses the blast "align two sequences" program  to perform pairwise local alignments of whole genome sequences.
an e value of 1e- <dig> was chosen as the blast threshold.
the alignments are performed between adjacent genomes on the linear genome maps and the results are overlaid on the maps using shading between the genomes to depict the aligned regions.
this shading is color-coded according to the e value with violet representing an e value of zero and red an e value equal to the threshold used.
representation of genome maps
perhaps the key functionality of phamerator is the construction of phage genome maps that incorporate the nucleotide similarity and pham assignment information.
an example of this is shown in figure  <dig>  when the six genomes that currently constitute cluster d are displayed, each is represented by a horizontal bar with coordinate markers with putative genes shown as colored boxes either above or below, corresponding to rightwards or leftwards transcription respectively .
the name of each gene is shown within the gene box, and the pham number is shown above with the total number of pham members shown in parentheses.
each pham has a designated color, with the exception of orphams that are shown as white boxes.
because the genomes shown in figure  <dig> are all members of the same cluster they share substantial nucleotide sequence similarity, which is reflected by the extensive violet shading between adjacent genomes in the stack of maps.
genomes can be easily re-positioned both vertically and horizontally within the display such that different pairwise relationships can be captured.
for genomes within a designated cluster - such as those in figure  <dig> - interruptions in the nucleotide sequence similarity are readily apparent, seen as either a reduced level of similarity  or by no shading  and pham <dig> ].
dotplot analysis clearly shows that gumball gene  <dig> and troll <dig> gene  <dig> segments are unrelated at the nucleotide level  and sequence alignment reveals that the discontinuities occur at the start codons of these genes, and those of the downstream genes .
the map function of phamerator provides a tool for readily identifying and analyzing these module boundaries.
phamerator-generated maps optionally can also display conserved domains identified with the automated cdd function .
domain hits are shown as yellow boxes or lines  within each gene box.
hovering the mouse over any domain pops up a description of that domain hit .
phamily circle representations of gene phylogenies
phamily circles provide a graphic way to illustrate the relationship between proteins in a phamily, and to display which genomes within the dataset contain members of that pham .
the strength of the pairwise relationships for a particular phamily can be simply accessed from the 'phams' function in the left hand panel .
the 'phams' window displays two separate panels with the upper one showing a numerical list of phams, the number of phamily members, and the clusters and subclusters that are represented.
selecting a pham directs a display in the lower window of each of the phamily members by gene name and phage.
when a gene name is selected, phamerator reports the clustalw and blastp score of each of the other members relative to the selected gene .
phamily circle diagrams can be generated for individual phams and include the name for each phage in the database positioned around the circumference of a circle, ordered and colored according to cluster and subcluster designation .
if a given phage has a gene that is a member of the phamily represented in the diagram, the protein name is included with the phage name.
arcs are drawn between pairs of genomes that contain a gene member of that phamily; relationships derived from clustalw analyses are represented in blue, and blastp in red.
in the pham <dig> example shown in figure  <dig>  the phamily of small proteins is present in  <dig> of the genomes and distributed among several cluster and subclusters.
some of the relationships are shown as blue-only arcs, indicating that the relationships exceed the threshold of  <dig> % amino acid sequence identity of the clustalw comparisons but does not meet the e value of 10- <dig> for the blastp comparisons .
the gene sequences can be readily exported for each pham and used to construct neighbor-joining trees for comparison with the phamily circles.
an abundance of orphams
the great genetic diversity of the mycobacteriophage population is reflected in the large number of orphams , the relatively low mean and median pham sizes , and the observation that 91% of the phams contain ten or fewer members .
the question arises as to whether the orpham designation is useful, because one gene member might be deemed insufficient to form a phamily.
we think this is a useful designation because this is an especially informative group when attempting to identify those genes that have been most recently acquired in evolutionary time  <cit> , but also note that orphams are especially abundant in singleton genomes for which no close relatives have yet been isolated .
the abundance therefore reflects the current poor sampling of the mycobacteriophage population, and as relatives of the current singleton genomes are discovered, many of the orphams are anticipated to be assembled into larger phams  <cit> .
we anticipate a substantial reduction in the proportion of orphams as we approach saturation of the phage population, but it is not yet simple to predict when that will occur.
discussion
we have explored the use of phamerator with several other sets of phage genomes including spo1-like phages of bacillus subtilis  <cit>  and a group of streptomyces phages, but we have recently successfully generated a database from  <dig> genomes, substantially larger than the 111-genomes described here.
we recognize that as the number of complete phage genome sequences increases that the computational time required increases as the square of the number of genes, and this could impose considerable limitations.
for example, increasing the number of genomes to  <dig>   <dig> - not an unreasonable expectation given the advances in dna sequencing technology - increases the number of pairwise computations to ~ <dig>  a 1000-fold increase in time over the current dataset.
however, with recent advances in cloud computing and the availability of massively parallel and multi-core computing systems we anticipate that these demands can be readily met.
for example, cloud-computing systems can provide more than a 1000-fold increase in the number of processors at minimal cost.
phamerator will remain a useful tool for comparative phage genome analysis for the next few years.
we also note that recent developments in alternative profile-based similarity searches such as hmmerhead and hhmer <dig> http://hmmer.org/ that greatly increase their efficiencies should provide additional phamerator components  <cit> .
CONCLUSIONS
phamerator provides a simple but useful computational tool for dissecting the genetic relationships among bacteriophage genomes, and displaying them in informative representations.
phamerator is especially useful for analysis of particular sets of phages such as the mycobacteriophages described here, but can be readily expanded to include a broader phage set, in which it is desired to map the horizontal exchange of genes between phage populations .
because of its computational intensity, it is less well suited to mapping global genome-scale relationships among large phage genome sets, but other programs have been described for this purpose  <cit> .
the use of a commonplace biopython framework and mysql database software should facilitate interaction of the phamerator database components with other web-based utilities to make this a broadly accessible utility.
competing interests
the authors declare that they have no competing interests.
authors' contributions
sgc.
wrote the phamerator code with assistance from mb and nd djs, rwh, and gfh contributed to the phamerator design and its applications.
sgc, rwh, djs and gfh wrote the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
phamerator program.
phamerator program.
click here for file
additional file 2
phamerator installation instructions.
this file contains installation instructions for phamerator.
click here for file
acknowledgements
this work was supported in part by a grant to the university of pittsburgh by the howard hughes medical institute  in support of g.f.h.
under hhmi's professors program.
support was also provided by grants from nih to r.w.h.
, to g.f.h.
and to s.g.c.
.

 bacteriophages represent a numerically vast, highly dynamic, evolutionarily ancient, and genetically highly diverse population. phage genomes are typically small compared to those of their bacterial hosts and no longer present significant technical challenges to sequence determination. as genomic sequencing approaches get simpler and cheaper, the availability of individual phage isolates for characterization becomes limiting, a need that can be effectively addressed through integrated research-education programs involving undergraduate and high school student investigators. in spite of their relatively small size, phages present significant challenges to accurate genome annotation including gene identification. two principal issues arise. first, phage genes tend on average to be small, approximately two-thirds the average size of bacterial genes. many of the genes required for virion structure and assembly are relatively large, but those in the non-structural genomic segments are small, often shorter than codons. secondly, phage genomes are replete with genes of unknown function for which no homologues have been described. mosaic architectures are hallmarks of phage genomes, and individual phages can be considered as particular combinations of interchangeable modules, each of which can be present in two or more different genomic contexts. in some cases, where the recombination events giving rise to these mosaic structures occurred relatively recently in evolutionary time, mosaicism is apparent through nucleotide sequence comparisons. when the events occurred in more remote evolutionary times the evidence of common ancestry is usually no longer apparent at the nucleotide level, but often can be revealed from comparison of the predicted amino acid sequences. such comparisons reveal that individual phage genomes are typically constructed from multiple modules - often corresponding to single genes - each of which has a distinctly different phylogeny. as such, accurate compilations of whole genome phylogenies that reflect the evolutionary history of the entire genome are not possible, and reticulate-based representations are needed to capture this evolutionary complexity. the mechanisms giving rise to genome mosaicism are unclear but must accommodate the striking observation that module boundaries correspond closely with gene boundaries, and in some cases, domain boundaries. one model invokes homologous recombination events targeted to short conserved boundary sequences between genes, and there is evidence for this in some phage genomes. however, there are numerous examples where no conserved boundary sequences are evident, raising the possibility that mosaicism results largely from illegitimate recombination events between randomly chosen partners sharing little or no sequence identity. in this second model, correspondence between module and gene boundaries results from the selection for gene function, not from targeting of the recombination events. comparison of genomes from phages that infect taxonomically diverse hosts typically provides little information into their evolution because of only very limited similarity at either the nucleotide or amino acid sequence level. phages of a common host, however, have the advantage that they are more likely to have been in recent genetic communication with each other and to have exchanged modules in recent evolutionary times. large sets of phage genomes are now available for several hosts including burkholderia, bacillus, enterobacteriaceae, mycobacteria, prochlorococcus and synechococcus, pseudomonas, and staphylococcus, although even these can span enormous genetic diversity. a large number of phages that infect mycobacterium smegmatis mc have been isolated and a comparative analysis of has been described. although these are genetically diverse, the diversity is heterogeneous, and phages can be grouped into'clusters' according to their overall nucleotide sequence relationships. of the published completely sequenced mycobacteriophage genomes, can be grouped into ten major clusters, seven of which can be further subdivided into subclusters according to the extent of the nucleotide similarities. five of the genomes have no close relatives and are referred to as'singletons'. because the currently sequenced mycobacteriophage genomes under-represent the mycobacteriophage population-at-large, these cluster designations will undergo modifications as new genomes are sequenced. there are, however, numerous examples of genes that are shared between phages of different clusters and whose common ancestry is only apparent from amino acid sequence similarity. we have proposed previously that mycobacteriophage genes related to each other can be grouped into phamilies and that mosaic relationships can be analyzed and represented using pham-annotated genomes maps and phamily circles that show the patterns of which phages contain members of particular phams. although manual or semi-automated approaches are applicable when only small numbers of genomes are analyzed, this becomes an impossible task as the number of genomes expands. we describe here a software program'phamerator' that provides bioinformatic tools for both analyzing and representing phage genome mosaicism. the core functionality of phamerator performs pair wise similarity searches between predicted protein products of a set of phage genomes, and assorts them into phamilies of related sequences. genome maps can be displayed that illustrate the relationships between phages at both the nucleotide and amino acid sequence level. moreover, the evolutionary histories of specific genes can be displayed by phamily circles in which all gene members of particular phams are represented, and for which multiple phams can be compared. we illustrate the utilities of phamerator using a set of completely sequenced mycobacteriophage genomes, but the program is applicable to any set of phage genomes for which comparative analysis is desired. methods phamerator database architecture phamerator is written entirely in the python computer programming language and makes use of a number of modules, including the biopython framework for computational biology. biopython provides a programmatic interface for sequence manipulation, the construction and parsing of files in relevant formats, and access to external command line applications such as those used for sequence alignment. for the latter, phamerator uses biopython to interact with local instances of blastp and clustalw. due to the significant computational time required for performing large numbers of sequence alignments, phamerator employs a distributed processing model that is implemented using python remote objects. phamerator was developed on ubuntu linux but should be able to run on any modern unix-based operating system. phamerator uses the mysql database software with a simple, custom database schema that incorporates and extends the relevant information found in genbank records. the phage and gene tables are populated with data from genbank files, while the remainder store data relevant to phamerator analysis or imported data from external databases such as the ncbi conserved domain database. additional tables are used to store current and historical pham assignments, and records are maintained of the splitting or joining of phams that can occur as new sequences are added to the database. in the event that a new mycobacteriophage protein is added to the database that has similarity to members of more than one existing pham, the phams are merged, their pham names retired, and a new pham created with a new number. the new pham contains each of the members of the joined phams in addition to the new protein. conversely, the addition of new proteins to the database can also invalidate an existing pham because the blastp e values used for determination of pham membership are in part dependent on the size of the data set. thus, after adding a new genome to the database, if a protein in a pham is no longer related to any members of the pham it is removed from the pham, and it is either placed into an existing pham if that pham includes a related protein, or it becomes an orpham. use of distributed computing resources adding genomes to a large phamerator database is time-consuming, with the length of time required being proportional to the size of the existing database and the number of genes in the genome to be added. however, the length of time required to perform this operation can be reduced by distributing the required calculations to several computers. the reduction in time required scales approximately linearly with the number of available computers. communication between computers is achieved using the python remote objects library, and can be distributed to computers on the local network or across the internet. availability and distribution phamerator operates using a client/server model. a relational database resides on the server computer and is distributed to client computers when they run the phamerator client program. this enables all clients to have the performance benefits associated with interacting with data stored on the local computer along with a benefit normally associated with processing data on a server--assurance that the data being used is up-to-date. because the data files being disseminated are relatively small, a modest server can handle the load of hundreds of users without performance deterioration. phamerator can be downloaded from: http://phamerator.csm.jmu.edu/files/phamerator-current.tar.gz. further information and installation instructions are provided in additional files and  rationale for phamerator construction and operations the pervasive mosaicism of bacteriophage genomes requires bioinformatic tools that can organize and display their complex relationships. two key questions arise in phage comparative analysis: what are the relationships between genes that are evolutionarily mobile within a given set of phage genomes, and how are they related to genes found in other genomes. both approaches are complicated by the presence of intragenic mosaicism reflecting distinct evolutionary histories of gene segments. phamerator is a computational tool designed to sort phage genes into phamilies of related sequences using pairwise amino acid sequence comparisons of predicted gene products. rather than using ortholog identification programs such as hmmer or pfam we have employed blastp and clustalw to perform pairwise comparisons that are then assembled into phamilies of related proteins. both of these programs are relatively fast computationally, a critical factor when large numbers of computations are required. for example, in the dataset of phage genomes used here, a total of × comparisons must be computed. we note that the use of efficient pairwise blastp comparisons to generate families of protein sequences has been described previously. to identify homologues of previously identified proteins phamerator performs automated searches of genbank non-redundant protein sequences, as well as searches for conserved domains in the ncbi conserved domain database using the rps-blast tool. this information can be exported in tabular form, or represented in a whole-genomic context. details of the phamerator program are described in the methods section and an overview of database structure is shown in figure pham-building parameters the building of phams is strongly influenced by the specific parameters used for amino acid comparisons. in early studies initiated prior to phamerator development we used a blastp cutoff value of and a clustalw cutoff of 25% amino acid sequence identity for manual pham building. however, as the number of mycobacteriophage genomes increased and computational processing became essential, it became clear that these largely arbitrarily chosen parameters promoted assembly of many large phams that require time-consuming manual deconvolution. we therefore explored the impact of varying the threshold values for blastp and clustalw on pham assembly. we first evaluated the effect of changing the threshold for clustalw comparisons. we varied the threshold between 50% amino acid sequence identity and %, and for each level determined the number of phams generated, the size of the largest pham, the number of orphams, the percentage of orphams, and the mean pham size. these data are informative and provide guidance as to the optimal parameters to use for routine database construction. in particular, we note that as the threshold for amino acid sequence similarity is made less stringent there is a reduction in the total number of phams reflecting the process of pham assembly. interestingly, this relationship is linear between 50% and %, with a reduction in the total pham number of about for every percentage of identity that is reduced. as the percent identity falls below % the relationship becomes non-linear, with a progressively greater reduction in the number of phams as the percent identity threshold falls from % to %. the number of orphams also reduces as stringency is relaxed, while mean pham size increases as stringency is relaxed, and there are notable changes between values above and below % identity. there is also a dramatic change in the size of the largest pham as the threshold level varies from % to %. at % the largest pham contains members, but increases to at %. the size of the largest pham is more stable between % and 50% identity thresholds and varies from to we interpret these data as indicating that between 50% and % identity, pham assembly proceeds in a manner that simply reflects the variation in the overall relationships between genes. however, at levels below %, there is an increasing proportion of phams that are more complex, such that not all pairwise matches within the pham are above the threshold level. one example might be where two genes have been fused into a single open reading frame, such that although genes a and c, and gene b and c, both surpass the threshold, the unrelated genes a and b do not. for routine database building purposes where we wish to avoid the assembly of large phams that then warrant subsequent deconvolution, we have chosen to use a % clustalw threshold, but note that comparison of phams generated with these parameters and those with lower stringencies should be useful in analyzing intragenic mosaic relationships. one advantage of clustalw as an assembly program is that the threshold values are independent of gene length. nonetheless, we predict there are instances where large genes may not exceed the clustalw threshold but are evidently homologues because of statistically informative blastp scores. we therefore examined the impact of including a blastp search along with the clustalw comparison and varying the blastp cut off value. as threshold values are made less stringent we observe a reduction in the total number of phams and the number of orphams, and a corresponding increase in mean pham size and the size of the largest pham. between blastp thresholds of 10- and 10- these is only a modest change in the total pham number and number of orphams, but the size of the largest pham changes from to. upon manual inspection of all phams with or more members at the blastp thresholds between 10- and 10- we see'false' pham assembly occur as illustrated by joining of a phamily of tape measure proteins with a phamily of lysin proteins, through sharing of small but closely related domains. nonetheless, inclusion of the blastp comparison with a 10- threshold joins several lysin phams that are separate when using just the clustalw comparison. we conclude that inclusion of blastp contributes rather little to the pham assembly process, but that a combination of a clustalw threshold of % and a blastp cut off of 10- offers optimal parameters for this dataset, building phams of homologous proteins while minimizing construction of complex phams in which only segments of the proteins are related to each other. when these parameters are applied to this dataset, the mycobacteriophage genomes contain a total of genes that assemble into phamilies with a mean size of genes/pham; phams are orphams and the largest pham contains members. identification of known homologues and conserved domains once a novel genome has been sequenced and annotated, questions about the functions of individual genes encoded within the genome can be addressed. this process is facilitated by analyzing the predicted gene products for the presence of conserved domains. numerous tools already exist for this purpose, but the ncbi conserved domain database aggregates many of them into a single, searchable dataset. these domain databases often use different, complementary techniques such as hidden markov models or position-specific scoring matrices to define domains and for the matching of novel sequences to existing domain models. to leverage the power of each of these approaches, we have implemented a system in phamerator whereby phage proteins are used to query cdd, and the results are presented in a searchable database browser in addition to being displayed on genome maps. the latter affords both a global view of the genomes and provides a way to visualize the conservation of specific sequences in the phamerator dataset with those in other organisms. an analysis of the distribution of cdd hits among the mycobacteriophage protein sequences reveals that a current search produces a total of matches among the predicted proteins, for an average of matches per protein. however, as expected, the distribution of matches is non-random, with only proteins having at least one cdd entry match, with the average number of matches for these proteins being. while proteins each match a single domain, a single protein, myrna gp, matches domain models, most of which are aminotransferases. generally, when numerous matches of a single protein to domains in the cdd are found, it reflects the redundant nature of the cdd dataset, although in some cases it reflects the presence of multiple domains within a single protein. computation of nucleotide sequence similarities phamerator uses the blast "align two sequences" program to perform pairwise local alignments of whole genome sequences. an e value of 1e- was chosen as the blast threshold. the alignments are performed between adjacent genomes on the linear genome maps and the results are overlaid on the maps using shading between the genomes to depict the aligned regions. this shading is color-coded according to the e value with violet representing an e value of zero and red an e value equal to the threshold used. representation of genome maps perhaps the key functionality of phamerator is the construction of phage genome maps that incorporate the nucleotide similarity and pham assignment information. an example of this is shown in figure when the six genomes that currently constitute cluster d are displayed, each is represented by a horizontal bar with coordinate markers with putative genes shown as colored boxes either above or below, corresponding to rightwards or leftwards transcription respectively. the name of each gene is shown within the gene box, and the pham number is shown above with the total number of pham members shown in parentheses. each pham has a designated color, with the exception of orphams that are shown as white boxes. because the genomes shown in figure are all members of the same cluster they share substantial nucleotide sequence similarity, which is reflected by the extensive violet shading between adjacent genomes in the stack of maps. genomes can be easily re-positioned both vertically and horizontally within the display such that different pairwise relationships can be captured. for genomes within a designated cluster - such as those in figure - interruptions in the nucleotide sequence similarity are readily apparent, seen as either a reduced level of similarity or by no shading and pham ]. dotplot analysis clearly shows that gumball gene and troll gene segments are unrelated at the nucleotide level and sequence alignment reveals that the discontinuities occur at the start codons of these genes, and those of the downstream genes. the map function of phamerator provides a tool for readily identifying and analyzing these module boundaries. phamerator-generated maps optionally can also display conserved domains identified with the automated cdd function. domain hits are shown as yellow boxes or lines within each gene box. hovering the mouse over any domain pops up a description of that domain hit. phamily circle representations of gene phylogenies phamily circles provide a graphic way to illustrate the relationship between proteins in a phamily, and to display which genomes within the dataset contain members of that pham. the strength of the pairwise relationships for a particular phamily can be simply accessed from the'phams' function in the left hand panel. the'phams' window displays two separate panels with the upper one showing a numerical list of phams, the number of phamily members, and the clusters and subclusters that are represented. selecting a pham directs a display in the lower window of each of the phamily members by gene name and phage. when a gene name is selected, phamerator reports the clustalw and blastp score of each of the other members relative to the selected gene. phamily circle diagrams can be generated for individual phams and include the name for each phage in the database positioned around the circumference of a circle, ordered and colored according to cluster and subcluster designation. if a given phage has a gene that is a member of the phamily represented in the diagram, the protein name is included with the phage name. arcs are drawn between pairs of genomes that contain a gene member of that phamily; relationships derived from clustalw analyses are represented in blue, and blastp in red. in the pham example shown in figure the phamily of small proteins is present in of the genomes and distributed among several cluster and subclusters. some of the relationships are shown as blue-only arcs, indicating that the relationships exceed the threshold of % amino acid sequence identity of the clustalw comparisons but does not meet the e value of 10- for the blastp comparisons. the gene sequences can be readily exported for each pham and used to construct neighbor-joining trees for comparison with the phamily circles. an abundance of orphams the great genetic diversity of the mycobacteriophage population is reflected in the large number of orphams, the relatively low mean and median pham sizes, and the observation that 91% of the phams contain ten or fewer members. the question arises as to whether the orpham designation is useful, because one gene member might be deemed insufficient to form a phamily. we think this is a useful designation because this is an especially informative group when attempting to identify those genes that have been most recently acquired in evolutionary time, but also note that orphams are especially abundant in singleton genomes for which no close relatives have yet been isolated. the abundance therefore reflects the current poor sampling of the mycobacteriophage population, and as relatives of the current singleton genomes are discovered, many of the orphams are anticipated to be assembled into larger phams. we anticipate a substantial reduction in the proportion of orphams as we approach saturation of the phage population, but it is not yet simple to predict when that will occur. discussion we have explored the use of phamerator with several other sets of phage genomes including spo1-like phages of bacillus subtilis and a group of streptomyces phages, but we have recently successfully generated a database from genomes, substantially larger than the 111-genomes described here. we recognize that as the number of complete phage genome sequences increases that the computational time required increases as the square of the number of genes, and this could impose considerable limitations. for example, increasing the number of genomes to - not an unreasonable expectation given the advances in dna sequencing technology - increases the number of pairwise computations to ~ a 1000-fold increase in time over the current dataset. however, with recent advances in cloud computing and the availability of massively parallel and multi-core computing systems we anticipate that these demands can be readily met. for example, cloud-computing systems can provide more than a 1000-fold increase in the number of processors at minimal cost. phamerator will remain a useful tool for comparative phage genome analysis for the next few years. we also note that recent developments in alternative profile-based similarity searches such as hmmerhead and hhmer http://hmmer.org/ that greatly increase their efficiencies should provide additional phamerator components.  phamerator provides a simple but useful computational tool for dissecting the genetic relationships among bacteriophage genomes, and displaying them in informative representations. phamerator is especially useful for analysis of particular sets of phages such as the mycobacteriophages described here, but can be readily expanded to include a broader phage set, in which it is desired to map the horizontal exchange of genes between phage populations. because of its computational intensity, it is less well suited to mapping global genome-scale relationships among large phage genome sets, but other programs have been described for this purpose. the use of a commonplace biopython framework and mysql database software should facilitate interaction of the phamerator database components with other web-based utilities to make this a broadly accessible utility. competing interests the authors declare that they have no competing interests. authors' contributions sgc. wrote the phamerator code with assistance from mb and nd djs, rwh, and gfh contributed to the phamerator design and its applications. sgc, rwh, djs and gfh wrote the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 phamerator program. phamerator program. click here for file additional file 2 phamerator installation instructions. this file contains installation instructions for phamerator. click here for file acknowledgements this work was supported in part by a grant to the university of pittsburgh by the howard hughes medical institute in support of g.f.h. under hhmi's professors program. support was also provided by grants from nih to r.w.h., to g.f.h. and to s.g.c..
bacteriophage genomes have mosaic architectures and are replete with small open reading frames of unknown function, presenting challenges in their annotation, comparative analysis, and representation. phamerator also generates phamily circle representations of gene phamilies, facilitating analysis of the different evolutionary histories of individual genes that migrate through phage populations by horizontal genetic exchange. we describe here a bioinformatic tool, phamerator, that assorts protein-coding genes into phamilies of related sequences using pairwise comparisons to generate a database of gene relationships. phamerator represents a useful tool for comparative genomic analysis and comparative representations of bacteriophage genomes.
 bacteriophage genomes have mosaic architectures and are replete with small open reading frames of unknown function, presenting challenges in their annotation, comparative analysis, and representation.  we describe here a bioinformatic tool, phamerator, that assorts protein-coding genes into phamilies of related sequences using pairwise comparisons to generate a database of gene relationships. this database is used to generate genome maps of multiple phages that incorporate nucleotide and amino acid sequence relationships, as well as genes containing conserved domains. phamerator also generates phamily circle representations of gene phamilies, facilitating analysis of the different evolutionary histories of individual genes that migrate through phage populations by horizontal genetic exchange.  phamerator represents a useful tool for comparative genomic analysis and comparative representations of bacteriophage genomes. 