BACKGROUND
systematic studies have shown that horizontal gene transfer  is prevalent in prokaryotes , where it serves as an important driving force of microbial evolution  <cit> .
hgt challenges the detection of vertical inheritance patterns in prokaryotes, and the application of conventional phylogenetic approaches to infer evolutionary history of microbial clades has seen increased limitations .
in essence, the ubiquitous nature of this process calls for the need to separate the vertical and horizontal patterns in evolutionary history of bacterial genomes.
however, this is not straightforward in practice and is especially difficult for deep historical events because the horizontally acquired genes evolved along with the recipient genomes, gradually losing the signatures of their original hosts   <cit> .
furthermore, hgts between closely related organisms, although common, are difficult to detect because in these cases donor and recipient share common compositional and phylogenetic features.
so far, multiple computational methods have been developed to facilitate hgt detection, which may be loosely categorized into three main strategies based on sequence composition, phylogenetic analysis, or best blast matches .
however, there appears to be poor agreement between outcomes of diverse methods, and comparative studies have repeatedly demonstrated that depending on method different sets of putative hgt-derived genes are identified from the same dataset, reflecting limitations in the current methodology for hgt prediction .
given the rapid increase in available annotated genome data, and the associated computational challenge of analyzing such data, the blast best match method has remained a popular surrogate for first pass discovery analyses of gene histories that differ from the strict vertical pattern  <cit> .
specifically, this strategy is practiced by sorting blast hits by measures such as bit scores, an indicator of sequence similarity, and the best match organism represented by the top hit is identified for each gene  <cit> .
if the best match is a distantly related organism, instead of one expected by vertical inheritance, then the gene is categorized as likely horizontally acquired  <cit>  .
in practice, researchers often identify the best match using the criterion of bidirectional best hits   <cit>  to rule out potential paralogs.
this approach has been applied in numerous studies .
examples of programs featuring this approach include pyphy  <cit> , phylogenie  <cit> , ngibws  <cit> , and darkhorse  <cit> , although the latter also employs a user-definable filter threshold in combination with taxonomic scaling.
as expected, blast-based hgt detection has limitations.
the bit score is based on sequence similarity, and provides only a rough estimate of the phylogenetic history between organisms.
also, best hits do not necessarily represent the nearest neighbors  <cit> .
most importantly, unexpected best hit can also be caused by reasons other than hgt, such as lack of sequence information in related organisms  <cit> , gene loss events  <cit> , stochastic similarity  <cit>  as well as database error  <cit> .
in practice, the predicted hgt candidates are often rejected by downstream phylogenetic analyses  <cit> .
similar to phylogenetic hgt detection methods, blast best match is effective in detecting recent hgt events, but shows reduced sensitivity for ancient events, when donor and recipient sequences have already diverged over the long history of evolution  <cit> .
additionally, merely a best match does not necessarily provide insights into the direction of gene flow.
despite the listed issues, the best blast match methods make use of all sequence data available in genbank.
this sidesteps the need for manual subsampling and curation of comparative sequence data that is a challenging and time-consuming step in phylogenetic methods.
this strategy is therefore likely to remain a feasible solution to explore microbial evolution in a first pass step, utilizing all of the ever increasing genomic data  <cit> .
considering this trend, we introduce a blast-based method to facilitate the detection of horizontal gene histories that aims to remedy some of the above outlined issues.
this approach starts with standard all-against-all blastp, and is followed by an investigation of the weight distribution of all hits grouped by phylogenetically informed user defined categories .
a general pattern of blast hit distributions  of the genomes of interest is computed, and blast hit weights of each single gene are compared to the general fingerprint.
this decreases sensitivity to stochastic disturbances.
because phylogenetic information is incorporated into this process, each resulting distribution is divided by uniquely defined cutoffs into typical and atypical gene populations.
using a combination of rules, a pool of genes that is putatively horizontally derived is reported.
it is recommended that the atypical gene pool be subjected to downstream phylogenetic validation of hgt, which is implemented in this pipeline.figure  <dig>
flowchart illustrating the procedures of the hgtector method.
parallelograms indicate input data or information, rectangles indicate processing steps, diamonds indicate decisions, and rounded rectangles indicate start and end of the work flow.
graphic illustrations of hypothetical phylogenetic tree, blast hit table and fingerprint are drawn on the right side of their corresponding steps.
tree topology and fingerprint  of tests on simulated genomes.
one representative test using either the idealized tree topology  or the randomized tree topology  is depicted .
the kernel density function of close weight distribution for both topologies  shows the distribution of all genes in the input genomes in black, and that of actual positive genes  in red.
genes involved in other simulated evolutionary events are shown in different colors in the lower panels.
locations of these genes in the general distribution are indicated as a rug below each plot.
the scales of x-axes between the upper and the lower panels are identical.
cutoffs computed by the program distinguishing the atypical region from the typical region are represented in dashed  and dotted  lines.
with this approach, the method retains all advantages of blast-based methods, such as efficient and exhaustive searches, which also facilitates re-analysis of data following the addition of new data.
rather than using general filter thresholds, and subsequent refining by taxonomic metrics , this method combines the two steps into one, and defines unique thresholds for each hierarchical level under consideration.
in order to assess the performance and robustness of this pipeline regarding the identification of putatively hgt-derived genes, it was applied towards simulated genomic data with known hgt events under consideration of various evolutionary forces.
the method was also tested on real genomic datasets from multiple organismal groups, as exemplified by rickettsia, whose hgt patterns have been previously studied .
the consistency of these results was compared to those obtained by other methods, including sequence composition and phylogenetic approaches.
overlaps between results were investigated and discussed in the framework of technical and biological challenges behind each method.
methods
method overview
at the core of this approach is an all-against-all blastp search of the protein product of each protein-coding gene  of the genome of interest against the genome database.
for each protein, the blast hits are recorded and sorted by their bit scores from high to low.
the bit scores are then normalized by dividing them by the bit score of the query sequence in order to account for protein length variance, so that every hit has a normalized bit score within the 0– <dig> range  <cit> .
the organism corresponding to each hit and the taxonomic ranks of each organism are identified and recorded from the ncbi taxonomy database.
to proceed with the analyses, hits of each gene have to be divided into self, close and distal groups.
in other words, the pipeline doesn’t use phylogenetic trees or taxonomic lineage information directly, but rather allows the user to define three relational hierarchical categories, with the biological questions of the research in mind .
this approach allows flexibility, because it can be scaled to the level of taxonomic/phylogenetic interest , and it can be adjusted to frequent updates in taxonomy.
the self group is considered the recipient, and always has to include the query genome, and, depending on analytical scale may also include its immediate sister organisms .
the close group will include representatives of the putative vertical inheritance history of the group .
the distal group includes all other organisms, which are considered phylogenetically distant from the query genome .
the method will then aim to identify genes that are likely derived from directional gene flow from groups of organisms within the distal group to members of the self group.
we introduced a measure to quantify quality of blast hits of each group.
this measure  is calculated per group by summing up the normalized blast bit scores of hits.
because three categories were defined a priori, this step will result in three weights per gene.
the three weights of all genes of the query genome are considered as three independent statistical populations.
if multiple genomes are analyzed together, the weight populations can be merged.
the three distributions together are defined as a fingerprint of the input genome.
a cutoff is set up to divide the weight distribution of each group  into typical  and atypical .
if most genes of the genome have a vertical inheritance history, the typical portion should include a majority of the genes, while the atypical portion should only include a small, but significant subset of genes whose hits of this specific group are underrepresented.
the cutoff defines the stringency of prediction: the higher the weight cutoff, the more genes are considered as being atypical.
the pipeline implements statistical approaches to compute the cutoffs, but the user is free to implement and use their own statistics.
because each genome or set of genomes generates a different fingerprint, cutoffs will vary, and are not transferrable across analyses.
for any gene to be predicted as putatively horizontally acquired, the following rules apply, which take into account all three weight distributions for that gene: rule  <dig> the gene is below cutoff  in the close weight distribution, suggesting that the orthologs of the query gene are absent in all or most of the sister groups of the organism of interest.
this means that the blast hits belonging to this group are significantly underrepresented, in terms of number of hits or bit score, or both .
this phenomenon can be explained by:  the gene lacks a vertical inheritance history, or alternatively,  the gene was vertically inherited, but underwent multiple independent gene loss events in the sister groups, a case that is usually less likely to be true  <cit> .
in the ideal scenario of a high likelihood of hgt, the weight should be zero, meaning that there is no close hit .
in real datasets however, there are sometimes sporadic close hits with low bit scores .
these hits may be due to stochastic sequence similarity, secondary hgts, paralogous genes or other mechanisms.
this situation typically causes false negatives in conventional best blast match methods  <cit> .
however, in the present method, these sporadic hits will not significantly alter the overall weight of a gene, thus hardly affecting the prediction results.
on the other hand, for a vertically transmitted gene, its orthologs may not always be present in each and every sister lineage.
occasional gene loss events may take away some of the expected number of close hits .
this situation is a major source of false positives in best blast match methods  <cit> .
however, this problem is effectively overcome in the present method, because sporadic absences of hits do not make the overall weight atypical.
moreover, one or a few high-score distal hits caused by natural  or artificial reasons   can easily deceive conventional methods  <cit> .
however, the present method is immune to this problem, as the close weight remains unchanged in this situation.
rule  <dig> the gene is equal to or above cutoff  in the distal weight distribution, meaning that the hits from distant organisms are not underrepresented .
this criterion sets up a filter at the donor side of potential hgt events: given the gene was transferred from a representative or ancestor of organisms that belong to the distal group, blast hits of the distal group itself should not be underrepresented.
otherwise, it is not convincing to conclude that the gene is hgt-derived.
the goal of this rule is stringency, in order to better distinguish putative hgt events from other scenarios that can also make a gene’s close weight atypical:  de novo gene origination within the self group,  inaccurate genome annotation that considers a non-coding region a gene, or  hgt from an unsequenced organism that is not phylogenetically close to any sequenced groups of organism.
these scenarios usually result in few to no distal hits .
meanwhile,  sequence similarity due to randomness instead of homology may also bring in some distal hits with low bit scores .
rule  <dig>  the gene is below cutoff  in the self weight distribution.
it means that a gene is only sporadically detected, rather than being prevalent across the whole self group of organisms .
this option will restrict the prediction to a subset of putatively hgt-derived genes that were acquired by specific lineages of the self group, instead of the whole self group.
often, these could be recent transfer events.
to assess the source of a putative hgt event, the best match organism from the distal group is reported.
it is important to understand that this best match is likely not the actual physical donor, but may be an extant representative of an ancestral, extinct donor.
we recommend using “donor link” to describe the directionality of transfer, and relationship between these organisms, instead of “donor”  <cit> .
a flowchart of the procedures of this method is illustrated in figure  <dig>
computational pipeline
general procedure
the hgtector pipeline  is written in perl, and is cross-platform supported, running in windows, mac os and linux systems.
it requires only the perl interpreter with its core modules, a default component of most mac os and linux distributions and is very easy to install in windows systems.
in the default mode, the program depends on no additional software or local databases to run.
this characteristic maximizes the ease of installation for users without professional computer background and resources.
optionally, it calls r  <cit>  to perform advanced statistical computing and graphing.
parameters of the program are managed by a central configuration file, which can be created and edited manually or via a graphical user interface .
the program is composed of the following procedures:
first, the program performs batch blastp of protein products of multiple genes supplied by the user.
multiple formats ranging from simple lists of ncbi accession numbers to annotated genomes in genbank format are supported.
it runs blastp either via web connection to the ncbi server, or with a standalone blast program and a local database.
faster alternatives to the original blastp may also be used, as long as their output files are in a compatible format.
it also harvests taxonomic information of each hit for each gene from the ncbi taxonomy database.
hierarchical taxonomic reports  and sequences of hits  can be retrieved optionally.
blast results are saved in nexus format  <cit> , which can be directly viewed by text editors, or opened as multiple sequence alignments by external programs such as seaview  <cit> , for additional analyses.
this characteristic facilitates downstream analyses, and compatibility with other programs.
blast hits are filtered by multiple optional functions to overcome putative taxon sampling bias that may affect blast hit distribution:  when multiple hits are present for one organism , only the best hit is maintained, representing the putative ortholog of the query protein  <cit> .
when multiple genomes are present for one species , only the genome that contains the best hit is maintained.
taxonomic name keywords or ids that represent unwanted blast targets, such as inaccurately defined taxonomic ranks  and biological categories , can be specified so that these hits will be omitted.
by default, the program will exclude genes without any non-self hits from subsequent analyses, because they may represent orfans  <cit> , resulting from de novo gene origination events , or transfer events from unknown sources that are very dissimilar from any sequenced genomes.
alternatively, they may represent genome annotation errors, which have been long recognized as a common and perturbing issue .
while these genes are not considered in the subsequent analysis, the genes are reported as “putative orfans or annotation errors”, or poe, in this pipeline .
this allows the user to check which poes were omitted, and if necessary make further adjustments to the analytical set up.
based on the retrieved taxonomic information, the program can automatically formulate a grouping scenario, in which the lowest rank that includes all input genomes is defined as the self group, and the next higher rank as the close group.
because taxonomic classifications in genbank may not always reflect the most current natural groupings of organisms, users may manually define hierarchies based on current knowledge of phylogeny and the purpose of their research.
with a properly defined grouping scenario, the program then calculates the three weights of each gene, computes a fingerprint of the whole genome, defines proper cutoffs and determines the population of atypical events, and possible hgts based on the selected rules.
basic statistical parameters of the three distributions of weights, as well as the weight populations themselves are reported.
the fingerprint may be visualized by box plots, histograms, density plots and scatter plots .
statistical analyses and graphing of blast hit distributions are automated in the program.
they are performed using perl codes, or by sending commands to r  <cit> .
the communication between perl and r is utilized by the perl module “statistics::r”.
while multiple statistical approaches are available for the user’s choice, the typical procedures, which are used in the tests described in this article, are as follows:figure  <dig>
fingerprint of seven
rickettsia
genomes.
blast hit weights of all protein-coding genes in the seven rickettsia genomes are plotted.
kernel density functions of the self, close and distal weights.
the x-axis represents the weight of each gene.
the y-axis represents the probability density of genes with the corresponding weight in the genomes.
in this example, rule  <dig>  and rule  <dig>  were applied.
the close and distal cutoffs computed under the conservative criterion are indicated by dashed lines.
the values of the cutoffs are denoted in each panel.
a scatter plot of the distal weight against the close weight, showing the clustering pattern of the genes.
each dot represents one gene.
genes predicted to be hgt-derived are framed by a red rectangle.
a zoom-in view of the left part of the previous plot.
genes that fall within the atypical region in the close weight distribution are colored by a blue-red color scheme based on the density-based silhouette , a measure of confidence that this gene belongs to the atypical cluster of genes .
the close cutoff used in the subsequent analyses is indicated by a dashed line.
cutoff determination
the program performs kernel density estimation  <cit>  to obtain a function of probability density distribution of the close weight for all genes.
by default it uses gaussian kernel smoothing with silverman's rule-of-thumb bandwidth selector  <cit> .
the user is allowed to choose a proper bandwidth selection factor that controls the smoothness of the curve.
the function is plotted and made visible to the user in real time.
statistically significant local minima  and maxima  are computed using the “pastecs” package  <cit>  with default parameters following kendall’s information theory  <cit> , and their x-coordinates are recorded and displayed to the user.
the program then automatically identifies a local minimum separating the typical from the atypical proportion of the gene population under consideration.
specifically, the biggest peak or group of continuous peaks in terms of number of genes it covers is identified as the typical region, and the rest is defined as the atypical region .
in addition to the apparent typical peak, there is usually a clearly identifiable peak located close to zero .
this peak usually includes “ideal” putatively hgt-derived genes .
between the zero peak and the typical peak is a transitional zone that likely consists of genes with ambiguous evolutionary history .
the program automatically reports two cutoffs: the x-coordinate of the identified local minimum is naturally chosen as a cutoff .
however, based on results of repeated tests, we recommend using the second cutoff, which is defined as the arithmetic mean of the x-coordinates of the local maximum of the zero peak and the local minimum of the selected pit  .
the choice between the two cutoffs depends on the goal of research, but for the identification of putative hgt events, the conservative cutoff is preferred as it meets a balance between precision and recall .
the program also implements several functions to assess the statistical significance of separating atypical genes from typical ones.
for the whole weight population, the program performs hartigans’ dip test  <cit>  to assess the non-unimodality of the weight distribution, which essentially is the statistical significance that a distribution can be divided into two or more distinct parts.
the test is performed by calling the “diptest” package in r  <cit> .
the dip statistic and the p-value for the test for unimodality are reported.
for each individual gene, the program computes a density-based silhouette  , which is a statistical measure of confidence that a certain data point  is allocated to a cluster   <cit> .
the same procedures apply to the distribution of the distal weight .
in addition to the above described statistical approaches, power users may also perform extra statistical analyses based on the weight data output by the program, and type user-defined cutoffs.
based on the cutoffs, the program reports a population of genes with an atypical, non-vertical history, which in the context of the a priori provided phylogenetic information represents a putative horizontal history.
the results are summarized in a choice of plain text, web page  or excel spreadsheet formats.
the latter allows for convenient downstream statistics of outputs, and includes hyperlinks that allow users to track each of the genes back to their original blast report.
it not only reports the number and percentage of putatively hgt-derived genes, but also optionally categorizes each gene in three contexts:  by putative donor group, which is described by user-designated higher taxonomic rank of the best match organism .
by functional annotation of protein products, which is provided by external sources, such as the output of blast2go  <cit> .
by gene orthology , which is identified by a built-in function of blast hits clustering or from external sources, such as the output of orthomcl  <cit> .
these reports  allow users an intuitive view of the prevalence of hgt-derived genes and the evolutionary, ecological and functional implications of hgts at levels of individual genes, gene families, whole genomes and multiple phylogenetically related genomes.
the whole analysis can be performed on a regular personal computer, as it does not require extensive computing power.
batch-blast is the most time-consuming step, which typically lasts several hours to several days, depending on the number of protein-coding genes in the genome of interest.
the computation time can be reduced if faster blast alternatives or smaller databases are used.
the statistical analysis typically takes only minutes.
data generated by the program can be parsed and re-used by other programs for multiple purposes.
as an additional, and important function, the program provides a complete phylogenetic pipeline, which automates the process of multiple sequence alignment, alignment trimming and phylogenetic tree reconstruction, by calling external local programs such as clustalw  <cit> , gblocks  <cit>  and raxml  <cit> , and parsing their outputs.
reconstructed phylogenetic trees are annotated with organismal names and are attached to blast reports, which in turn can be directly viewed by external programs such as figtree  <cit> .
this function allows users to monitor and validate prediction results by manually checking the evolutionary scenarios of individual genes.
analysis of simulated genomic data
to assess the performance of this method under the impact of various evolutionary scenarios, as well as to compare hgtector to conventional blast methods, we tested the above-described pipeline on simulated genomic data.
simulation of genome evolution
simulated genomic data were generated by alf  version  <dig> , a program that simulates genome evolution  <cit> .
in each simulation,  <dig> species evolved from one randomly generated root genome containing  <dig> protein-coding genes that are no shorter than  <dig> aa.
during the process random inter-genomic hgt events occurred under a pre-defined global rate, which varied between simulations .
the process also simulated the following evolutionary forces in addition to hgt, at random rates: speciation, character substitution, insertion and deletion, gc-content amelioration, rate variation among sites and among genes, gene duplication and loss.
many of these forces are known to affect hgt prediction  <cit> .
two types of simulations were performed.
first, an idealized, pre-defined tree topology was used, in which all representatives of close- and self-designation are grouped in a polytomy, signifying  equal genetic distance, thus eliminating the impact from evolutionary artifacts.
the rest of the taxa are placed relatively distant  in the tree to lower the possibility of stochastic blast matches.
the topology is detailed as follows:  <dig> clades branch from the tree base simultaneously at time point  <dig> .
in each clade,  <dig> species branch off simultaneously at time point  <dig>  the tips are at time point  <dig> for all species.
for each clade, one species was randomly chosen as the input genome , while nine species were defined as the close group.
all  <dig> species in other clades are considered as the distal group .
all ten clades were analyzed in this manner.
this simulation was replicated  <dig> times , with the hgt rate ranging from  <dig>  to  <dig>  with an interval of  <dig> .
second, a randomized birth-death tree was generated in alf per simulation , mimicking a more realistic topological situation.
a random clade was manually chosen from the tree, as long as it met the following criteria:  3– <dig> self species that formed a clear monophyletic group;  10– <dig> close species that were closely related to the self clade;  the self and close species together formed a clear monophyletic group that was independent from all other species  .
this simulation was replicated for  <dig> times, with the hgt rate of each replicate randomly sampled from a range of  <dig> to  <dig> .
evolutionary events  simulated in each analysis were extracted from the alf log file.
the time, species and genes involved in these events were recorded.
hgt events involving the self group were further categorized by their donor, recipient and time into the following groups: target hgt , ancient hgt , outgoing hgt , and secondary hgt .
gene loss and gene duplication events taking place in the close group were singled out as they directly influence the distribution of the close weight.
for each simulation , a blastp database including the protein sequences from all  <dig> genomes was created using the standalone blast program  <cit> .
all-against-all blastp was performed with an e-value cutoff at 1 × 10- <dig> for all genes in the selected self genomes, considering at least  <dig> hits.
in order to demonstrate the effect of gene duplication on the prediction result, the program option of excluding paralogs was turned off.
a modified version of the hgtector pipeline was created to parse the alf outputs and to perform analysis.
both the conservative and the relaxed cutoff criteria were tested.
the fingerprint was plotted along with the distribution of actual positives  as well as other evolutionary scenarios  for manual observation in addition to statistical analysis.
comparison to conventional blast approach
the performance of hgtector in the context of identifying atypical and putatively hgt-derived genes, was compared against that of commonly used conventional blast-based methods by modifying the pipeline to mimic the conventional approach, which does not consider overall hit distribution.
specifically, we considered hgt-events for the scenario where no close hits are present ; or the scenario where there is at least one distal hit with a bit score higher than those of any close hits .
assessing performance
the predicting power in the context of target hgts under each criterion, as well as for conventional blast was assessed by precision and recall.
precision describes the number of true positives over the number of all predicted positives , and recall describes the number of true positives over number of all actual positives , which are commonly used performance markers for binary classifications.
statistical analysis and plotting were conducted using r  <cit> .
application of hgtector to real datasets
hgtector was further tested on a variety of real-world genomic data, covering bacteria , unicellular eukaryotes  and human  , the last of which serves as a realistic negative control since hgts into genomes of higher animals are known to be very rare  <cit> .
the most in depth analysis was conducted on the rickettsia dataset, because it affords comparisons to previous results obtained with a variety of methods  <cit> .table  <dig>
real genomic datasets tested in this study
sfg
rickettsia
r. conorii str.
malish  <dig> 1
streptococcus
helicobacter
erwinia
mycobacterium africanum
galdieria sulphuraria
g. sulphuraria 3
homo sapiens
1the genomes used in this study are identical to those used in  <cit> .
2two independent analyses were conducted on different dates, and similar outcomes were obtained.
the more recent result was reported.
3the genome used in this study is identical to that used in  <cit> .
4for genes with multiple isoforms, the longest cds was extracted using an in-house perl script and used for the analysis.
analysis on the rickettsia dataset
out of all available rickettsia genomes, we selected seven representative rickettsia species with fully-annotated genomes for this analysis .
all of these species belong to the spotted fever group , a traditional classification group of rickettsia.
a grouping scenario was chosen based on the taxonomy and phylogeny of major rickettsia species, which has been well resolved by recent studies  <cit> .
specifically, we defined the self group as sfg , the close group as order rickettsiales , excluding sfg, and the distal group as all non-rickettsiales organisms.
all-against-all blastp was performed against the ncbi non-redundant protein sequences  database with an e-value cutoff at 1 × 10- <dig>  a soft filtering for low sequence complexity regions, which was suggested as the optimal parameter for ortholog identification  <cit> , was used.
hits with organism names containing these words were excluded for the purpose of this analysis: unknown, uncultured, unidentified, unclassified, environmental, plasmid, vector, synthetic, and phage.
up to one hit from each organism was preserved.
a maximum number of  <dig> hits were preserved for each protein.
a global fingerprint was computed and graphed to describe the pattern of blast hit distribution of all seven genomes.
cutoffs for the three groups of weights were computed using the built-in kernel density estimation function under the conservative criterion.
the default rules  <dig> and  <dig> were applied.
assessing stochastic events using real datasets
to test the impact of stochastic events on prediction results, the following simulations were carried out on the rickettsia dataset:  database error ,  gene loss in the close group,  rate variation in the close group,  rate variation in the input genomes,  taxon sampling bias , and  smaller sample size of genes.
multiple degrees of modification intensity  for each type of events were tested, each having  <dig> replicates .
specifically:database error.
for each hit, there is x proportion of chance that its corresponding organism was assigned to an organism randomly sampled from the pool of blast results of all genes.
gene loss in the close group.
for each close hit, there is x proportion of chance of being removed from the blast hit table.
rate variation in the close group.
for each close hit, there is x proportion of chance that its bit score is divided by a factor subject to a gamma distribution with shape parameter k =  <dig> and scale parameter θ = 1:
where s <dig> and s <dig> refer to the bit score before and after manipulation, respectively.
rate variation in the input genomes.
for each query gene, there is x proportion of chance that the bit scores of all its hits are divided by a factor subject to a gamma distribution same as above.
taxon sampling bias.
for selected representative groups of organisms from close and distal groups , all hits belonging to this group were replicated into x copies .
smaller sample size of genes.
x out of all  <dig> genes were randomly selected for fingerprint calculation in each replicate.
the resulting cutoffs were in turn used for prediction on all genes.
replicates that failed to pass the hartigans’ dip test   were counted and excluded from the analysis.
hgtector analysis was conducted on these replicates using the same procedures as in the standard analysis on the unmodified dataset .
the results were compared to the results derived from a conventional blast-based approach under the d > c criterion .
precision and recall of the results were computed using the result of the standard analysis as the reference.
cross-method comparison
in order to evaluate the performance of hgtector on the rickettsia dataset in the context of other available methods, results were compared to two examples from each of the three currently employed strategies: bbh   <cit>  and darkhorse  <cit>  based on best blast match; gist  <cit>  and islandviewer  <cit>  based on sequence composition, and two studies conducted by merhej et al.
<cit>  and le et al.
<cit> , using phylogenetic approaches.
we exemplified this comparison on the r. felis genome, which previously has been demonstrated to have high hgt frequency  <cit> .
specifically, bbh analysis was performed as a built-in function of the present pipeline.
this method is similar to the d > c criterion as described above, except for an additional reverse blast step with the same parameters to confirm that the two sequences are each other’s best match within their host genomes.
the result by darkhorse was downloaded from the darkhorse server .
default parameters were used, in which the blastp e-value cutoff is also 1 × 10- <dig>  all three available phylogenetic granularities, strain, species and genus, were used and the results were merged, in order to maximize the discovery rate.
both gist and islandviewer are targeting large pieces of heterogeneous genomic regions, or genomic islands   <cit> .
gist is a synchronization of five subprograms: alienhunter  <cit> , islandpath  <cit> , sigi-hmm  <cit> , indegenius  <cit>  and pai-ida  <cit> .
the subprograms were executed in a local system and the results were processed using the egid algorithm  <cit>  to get a consensus result.
islandviewer is an integration of three subprograms: islandpick, sigi-hmm and islandpath.
the integrated result was downloaded from the islandviewer server .
results of gist and islandviewer were further processed by an in-house perl script to extract the genes included in the genomic islands.
the putative hgt-derived genes predicted by merhej et al.
<cit>  and by le et al.
<cit>  were extracted from the original publications.
specifically, merhej et al.
identified  <dig> hgt-derived genes in the r. felis genome that are linked with organisms other than sfg rickettsia
<cit> .
le et al.
identified  <dig> instances of hgt from outside rickettsiales into the r. felis genome  <cit> .
the predicted hgt-derived genes or genomic islands by different methods were spatially mapped to the r. felis genome and visualized in geneious  <dig>   <cit> .
an “overlap factor”  was employed as a criterion to compare the outcomes of different methods by assessing the overlap.
this was expressed as the negative logarithm of the likelihood that the overlap was obtained by chance.
to compute an of, the number of the same genes predicted by each method pair was counted, and the of was calculated following the probability mass function of the hypergeometric distribution:
where n is the total number of genes; a and b are the numbers of genes predicted by two methods, respectively; k is the number of genes overlapping by two sets of results.
the larger an of is, the more overlapping, and thus more consistent the two sets of results are, and the more likely it is that the two methods are identifying the same group of genes.
RESULTS
performance on simulated genomic data
testing precision and recall
in all experimental groups under the idealized tree topology, a clear bimodal distribution was observed , which is expected when hgt is present in the data.
meanwhile, none of the negative control groups have an identifiable zero peak, which is equivalent to a vertical history for all genes .
both cutoff criteria achieve high precision and recall simultaneously.
in particular, under the conservative criterion,  <dig> % of the prediction results are true positives.
meanwhile, they cover over  <dig> % of all actual hgt-derived genes.
the more relaxed criterion still achieves a precision of  <dig> % and a recall of  <dig> %.
in tests under randomized tree topologies, a larger transitional zone is present between peaks of the expected bimodal distribution, which is also frequently observed in real datasets.
both, precision and recall are affected compared to the idealized scenario, but in different measures.
the conservative cutoff still maintains reasonably high precision  and recall  simultaneously, relative to the known number of hgt events.
the relaxed criterion keeps equally high recall  but its precision drops significantly  .figure  <dig>
comparison of performance of hgtector and conventional blast-based method on simulated genomes.
the methods were tested on the simulated genomic data under an idealized  or randomized  tree topology.
“con” and “rel” represent conservative and relaxed cutoffs in the hgtector analysis.
“c = 0” and “d > c” are two criteria in the conventional blast-based method.
each experimental group is composed of  <dig> tests.
the distribution of results in terms of precision  and recall  is depicted by box plots.
the mean value of each group is label above the corresponding column.
because the test under randomized topologies is a more realistic representation of datasets, its result serve as a better reference for the practical consideration of the hgtector application.
given the simultaneously high precision and recall, we recommend using the conservative cutoff for both initial hgt candidate screening  and standalone hgt discovery .
the relaxed cutoff may be considered only when the user wants to maximize discovery rate in an initial screening, in spite of its higher false positive rate.
the comparison between idealized and randomized topologies further indicates the positive correlation between prediction success and a properly defined grouping scenario, in which:  the self-close clade is relatively distant from any other organisms; and  there are multiple subclades in the close group, each having similar number of taxa represented in the database.
varying global hgt rates seemingly show little effect on the stability of precision and recall of this method , suggesting that the predicting power is independent of hgt rate.
in comparison, the performance of both conventional blast-based approaches  is notably unbalanced.
atypical genes falling under the c =  <dig> criterion have the highest precision , but very low recall .
the d > c criterion has high recall  but extremely low precision , showing an intolerably high false positive rate .
from a practical perspective, c =  <dig> is too stringent, thus omitting a big portion of true hgt-derived genes affected by stochastic events, while d > c is too relaxed and not capable of differentiating genes that have high-score distal hits merely due to stochastic reasons instead of hgt .
it should be mentioned that c =  <dig> is not applicable to real datasets due to frequent genome annotation errors and orfans, both of which may have a zero close weight.
evaluating other evolutionary scenarios
the impact of other evolutionary events on the fingerprint and the division between typical and atypical gene populations was explored.
outgoing hgt  seemingly does not significantly alter the close weight of a gene.
gene loss decreases close weight and gene duplication increases it, both within an insignificant range .
most importantly, genes within these three categories of evolutionary history still fall within the typical region and were not mistakenly detected as atypical by hgtector.
it is particularly notable that the majority of genes derived from ancient hgt events are located in the transitional zone .
expectedly, they constitute a considerable portion of false positives in our analyses.
in other words, depending on cutoff, more of them are likely to be placed in the typical population, instead of the atypical, and presumably non-vertical population.
a similar pattern was observed for secondary hgt events.
although not frequent in the simulated genomes here, these events are actually very frequent in real datasets, as hgt frequency is higher between closely related organisms  <cit> .
however, in the conventional blast method  most false positives are composed of mainly outgoing hgts but only a few ancient hgts.
application to real datasets
the close weight distributions for real datasets exhibit a bimodal distribution containing a broad typical region and a sharp atypical peak that is located close to zero  .
therefore, the cutoffs that divide genes with atypical blast hit distribution from typical ones can be set accordingly and hgt events can be predicted based on the cutoffs.
as an example, fingerprint plots for rickettsia and galdieria exhibit an apparent overlapping pattern between the atypical peak and the hgt-derived genes identified by previous phylogenetic studies  <cit>  .
in contrast, the human genome does not have an apparent atypical peak, which is also expected .
analysis of the rickettsia genomes
the global fingerprint describing the pattern of blast hit distribution of the rickettsia analysis is illustrated in figure 3a-c. results from hartigans’ dip test strongly support the non-unimodality of all three weight distributions .
the seven rickettsia genomes contain a total number of  <dig> annotated chromosomal protein-coding genes, of which  <dig> genes have an atypical close weight and a typical distal weight, and thus are potentially hgt-derived .
the percentage of putative hgt-derived genes in a genome ranges from  <dig> %  in r. akari to  <dig> %  in r. felis.
the number of putative hgt-derived genes per genome is positively correlated to the size of the genome , implying contribution of hgt to the relative genome expansions in the overall reductive trend of rickettsia genome evolution, confirming previous studies  <cit> .
it is notable that the r. felis genome contains significantly more putative hgt-derived genes  than the rest of the genomes , suggesting a particular prevalence of potential hgt events in r. felis evolution.
the comparison between the fingerprint on r. felis genome alone  and the global fingerprint  clearly reveals that r. felis has a much larger atypical peak in the close weight distribution.
this reinforces outcomes from merhej et al.
<cit> , which found that the frequency of cross-species bacterial recombination in r. felis had been underestimated.table  <dig>
summary of genes predicted to be horizontally acquired in seven
rickettsia
genomes
r. akari
r. felis
r. massiliae
r. slovaca
r. rickettsii
r. africae
r. conorii
mean
to further explore the biological information behind the predicted patterns, the results have been summarized in three separate contexts: by putative donor group, by functional annotation and by gene orthology.
similar to previous in depth analyses by merhej et al.
<cit> , our analyses revealed frequent donor links, such as legionellales, enterobacteriales and burkholderiales, and frequently transferred gene categories, such as genes encoding transposases and genes involved in phage and plasmid activity .
stochastic manipulation of rickettsia data: stability of prediction
the results of these simulations  suggest that hgtector’s performance is notably insensitive to database error.
even under an extremely high proportion  of mislabeled genes, its precision remains close to 100%.
as a comparison, the precision of the conventional blast-based approach d > c drops below 15% .
gene loss and rate variation in the close group also have remarkably weaker deleterious effects on precision of hgtector than that of conventional blast based methods under the d > c scenario .
rate variations in input genomes destabilize prediction results at a certain range , but overall, precision remains high in hgtector when rate increases, as compared to conventional blast .
taxon sampling bias has little effect when it happens to distal organisms.
however, in some close organisms , this effect may severely compromise the predicting power if the bias is strong  .
smaller sample size of input genes has a minor effect on the prediction result , suggesting that the method is still effective on incomplete genomes.
it should be noted that when the sample size is extremely small , the statistical power of gene clustering can be significantly impaired .
the above results suggest that our method is generally unaffected by stochastic events, some of which are major challenges to current methods .
a concern is taxon sampling bias in the close group, an issue that can be alleviated by properly defining grouping scenario and masking redundant taxonomic groups using hgtector’s flexible functionality .
comparison of hgtector with other methods for analysis of the r. felis genome
together, hgtector and six other methods identified  <dig> putative hgt-derived genes, out of all  <dig> chromosomal protein-coding genes in the r. felis genome  .
of these  <dig> genes,  <dig> have been uniquely identified by hgtector.
there’s a considerable portion of overlapping hits between each pair of results, but none of the programs produce completely identical hgt predictions.
as expected, there is relatively more overlap between two methods of the same category.
meanwhile, the overlap factors  between two methods from different categories are significantly lower.
however, hgtector is a notable exception because its overlap with any other method of the three categories is significantly higher than between other methods from different categories .figure  <dig>
predicted hgts by multiple methods mapped onto the
rickettsia felis
genome.
a circular view of the whole chromosome of r. felis, with genomic islands  predicted by islandviewer and gist indicated by boxes, and putative hgt-derived genes predicted by other methods indicated by arrowheads.
by comparing the results of hgtector and the results of the other two blast-based methods, it is noticeable that hgtector can effectively exclude false positives caused by database errors.
for instance, bbh and darkhorse respectively predicted  <dig> and  <dig> genes to be acquired from ixodes scapularis , which is a common host of rickettsia.
therefore, it is likely that the sequenced ixodes samples contained rickettsia endosymbionts, and its dna was also sequenced and indiscriminately labeled as ixodes dna in the unassembled shotgun sequences .
in other words, these are potential database errors.
in contrast, in hgtector’s result, none of these genes were predicted as hgt-derived, because they have considerable amount of close hits , despite the presence of a single high-score ixodes hit.
detection of hgt is fraught with challenges, such as ambiguity in compositional features and phyletic patterns  <cit> , difficulty in phylogenetic reconstruction  <cit> , interference from database error and incompleteness  <cit> .
the observed low consistency between methods is not surprising, and has precedents in multiple previous studies .
given the significant overlap of hgtector with all other previously employed approaches, we suggest that it is effective in producing meaningful prediction results on real datasets.
CONCLUSIONS
in this paper a novel method for genome-wide detection of vertical versus non-vertical gene history  is presented.
it features a statistical analysis of blast hit distribution patterns in the context of a priori defined phylogenetic hierarchies.
the innovation of this method is the systematic consideration of all blast hits of all genes within selected genomes.
this is in contrast to conventional blast-based approaches, which typically rely on a single best hit for each gene .
the three-category grouping scenario is a simplified but effective implementation of prior phylogenetic knowledge into a blast hit distribution analysis.
the set-up allows high flexibility in group assignments that best match the taxonomic level of the user’s interest, as well as immediate response to frequent changes in microbial taxonomy.
the advantage of this systematic approach is that it captures the overall image of gene evolution while being significantly insensitive to stochastic events.
as demonstrated in this study, stochastic events such as gene loss, rate variation and database error may impose serious problem to conventional methods , but have comparatively negligible effects on hgtector.
the core assumption of the method is that the typical and atypical parts of the blast hit distribution are distinguishable.
this assumption is repeatedly supported by tests on both simulated and real-world genomic data .
with the proposed procedures of computing cutoffs and the rules of targeting genes that are likely to be hgt-derived, remarkable prediction success was achieved .
given these results, we suggest that hgtector is a useful addition to conventional blast-based approaches.
the hgtector pipeline has advantages of speed, automation, compatibility and low requirement for computational resources, making this program a generally applicable tool for discovery of vertical, and non-vertical history of genes, as well as initial hgt prediction.
it is especially suitable for gaining a rapid and comprehensive overview of newly sequenced genomes to identify their evolutionary and ecological linkages with other organisms, facilitating further exploration of the functional drivers of the dynamics of genome evolution.
it has to be made clear that an atypical blast hit distribution is an empirical observation, rather than a strict certification of hgt.
since an hgt predominantly reflects a past evolutionary event, it is theoretically impossible to identify exact gene donors and mechanisms, and any analysis is only an approximation to possible scenarios.
therefore, we recommend hgtector as a discovery tool for a detection of potential hgt-derived genes that can be further analyzed with phylogenetic approaches.
this is much more effective approach than a very time-consuming and technically challenging process of a priori phylogenetic analysis of all genes within all target genomes, which becomes decreasingly feasible as more genomic data are present.
electronic supplementary material
additional file 1: figure s1:
figure s <dig>  illustration of patterns of blast hit distribution and possible explanations.
figure s <dig>  precision-recall plot of results on simulated genomes.
figure s <dig>  relationship between global hgt rate and performance of methods.
figure s <dig>  fingerprints of genomes of various organisms.
figure s <dig>  distribution of blast hit weights of the r. felis genome.
figure s <dig>  heat map indicating percentages of predicted hgt-derived genes by putative bacterial donor groups in rickettsia genomes.
figure s <dig>  heat map indicating percentages of predicted hgt-derived genes by functional annotations in rickettsia genomes.
figure s <dig>  stability of results on the rickettsia dataset with various simulated stochastic events.
figure s <dig>  stability of results on the rickettsia dataset with simulated taxon sampling bias.
figure s <dig>  stability of results on the rickettsia dataset with smaller sample size of genes for fingerprint calculation.
figure s <dig>  comparison of prediction results in the r. felis genome by multiple methods.
additional file 2: table s1: protein accession numbers of predicted hgt-derived genes by hgtector.
table s <dig>  predicted hgt-derived genes categorized by gene orthology.
table s <dig>  prediction results of protein-coding genes in the r. felis genome by multiple methods.
table s <dig>  pairwise comparison of hgt-prediction results in the r. felis genome by multiple methods.
abbreviations
bbhbidirectional best match
blastbasic local alignment search tool
blastpprotein blast
dbsdensity-based silhouette
gigenomic island
gogene ontology
guigraphical user interface
hgthorizontal gene transfer
htmlhypertext markup language
lpilineage probability index
mlmaximum likelihood
ncbinational center for biotechnology information
nrnon-redundant protein sequences database
ofoverlap factor
poeputative orfans or annotation errors
sfgspotted fever group.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
qz and kd conceived and designed the study.
qz developed the algorithm, wrote the program, performed the data analyses and drafted the manuscript.
kd contributed to the design of algorithm and program, supervised the data analyses and edited the manuscript.
mk contributed to the discussion and helped to revise the manuscript.
all authors read and approved the final manuscript.

 systematic studies have shown that horizontal gene transfer is prevalent in prokaryotes, where it serves as an important driving force of microbial evolution. hgt challenges the detection of vertical inheritance patterns in prokaryotes, and the application of conventional phylogenetic approaches to infer evolutionary history of microbial clades has seen increased limitations. in essence, the ubiquitous nature of this process calls for the need to separate the vertical and horizontal patterns in evolutionary history of bacterial genomes. however, this is not straightforward in practice and is especially difficult for deep historical events because the horizontally acquired genes evolved along with the recipient genomes, gradually losing the signatures of their original hosts. furthermore, hgts between closely related organisms, although common, are difficult to detect because in these cases donor and recipient share common compositional and phylogenetic features. so far, multiple computational methods have been developed to facilitate hgt detection, which may be loosely categorized into three main strategies based on sequence composition, phylogenetic analysis, or best blast matches. however, there appears to be poor agreement between outcomes of diverse methods, and comparative studies have repeatedly demonstrated that depending on method different sets of putative hgt-derived genes are identified from the same dataset, reflecting limitations in the current methodology for hgt prediction. given the rapid increase in available annotated genome data, and the associated computational challenge of analyzing such data, the blast best match method has remained a popular surrogate for first pass discovery analyses of gene histories that differ from the strict vertical pattern. specifically, this strategy is practiced by sorting blast hits by measures such as bit scores, an indicator of sequence similarity, and the best match organism represented by the top hit is identified for each gene. if the best match is a distantly related organism, instead of one expected by vertical inheritance, then the gene is categorized as likely horizontally acquired. in practice, researchers often identify the best match using the criterion of bidirectional best hits to rule out potential paralogs. this approach has been applied in numerous studies. examples of programs featuring this approach include pyphy, phylogenie, ngibws, and darkhorse, although the latter also employs a user-definable filter threshold in combination with taxonomic scaling. as expected, blast-based hgt detection has limitations. the bit score is based on sequence similarity, and provides only a rough estimate of the phylogenetic history between organisms. also, best hits do not necessarily represent the nearest neighbors. most importantly, unexpected best hit can also be caused by reasons other than hgt, such as lack of sequence information in related organisms, gene loss events, stochastic similarity as well as database error. in practice, the predicted hgt candidates are often rejected by downstream phylogenetic analyses. similar to phylogenetic hgt detection methods, blast best match is effective in detecting recent hgt events, but shows reduced sensitivity for ancient events, when donor and recipient sequences have already diverged over the long history of evolution. additionally, merely a best match does not necessarily provide insights into the direction of gene flow. despite the listed issues, the best blast match methods make use of all sequence data available in genbank. this sidesteps the need for manual subsampling and curation of comparative sequence data that is a challenging and time-consuming step in phylogenetic methods. this strategy is therefore likely to remain a feasible solution to explore microbial evolution in a first pass step, utilizing all of the ever increasing genomic data. considering this trend, we introduce a blast-based method to facilitate the detection of horizontal gene histories that aims to remedy some of the above outlined issues. this approach starts with standard all-against-all blastp, and is followed by an investigation of the weight distribution of all hits grouped by phylogenetically informed user defined categories. a general pattern of blast hit distributions of the genomes of interest is computed, and blast hit weights of each single gene are compared to the general fingerprint. this decreases sensitivity to stochastic disturbances. because phylogenetic information is incorporated into this process, each resulting distribution is divided by uniquely defined cutoffs into typical and atypical gene populations. using a combination of rules, a pool of genes that is putatively horizontally derived is reported. it is recommended that the atypical gene pool be subjected to downstream phylogenetic validation of hgt, which is implemented in this pipeline.figure flowchart illustrating the procedures of the hgtector method. parallelograms indicate input data or information, rectangles indicate processing steps, diamonds indicate decisions, and rounded rectangles indicate start and end of the work flow. graphic illustrations of hypothetical phylogenetic tree, blast hit table and fingerprint are drawn on the right side of their corresponding steps. tree topology and fingerprint of tests on simulated genomes. one representative test using either the idealized tree topology or the randomized tree topology is depicted. the kernel density function of close weight distribution for both topologies shows the distribution of all genes in the input genomes in black, and that of actual positive genes in red. genes involved in other simulated evolutionary events are shown in different colors in the lower panels. locations of these genes in the general distribution are indicated as a rug below each plot. the scales of x-axes between the upper and the lower panels are identical. cutoffs computed by the program distinguishing the atypical region from the typical region are represented in dashed and dotted lines. with this approach, the method retains all advantages of blast-based methods, such as efficient and exhaustive searches, which also facilitates re-analysis of data following the addition of new data. rather than using general filter thresholds, and subsequent refining by taxonomic metrics, this method combines the two steps into one, and defines unique thresholds for each hierarchical level under consideration. in order to assess the performance and robustness of this pipeline regarding the identification of putatively hgt-derived genes, it was applied towards simulated genomic data with known hgt events under consideration of various evolutionary forces. the method was also tested on real genomic datasets from multiple organismal groups, as exemplified by rickettsia, whose hgt patterns have been previously studied. the consistency of these results was compared to those obtained by other methods, including sequence composition and phylogenetic approaches. overlaps between results were investigated and discussed in the framework of technical and biological challenges behind each method. methods method overview at the core of this approach is an all-against-all blastp search of the protein product of each protein-coding gene of the genome of interest against the genome database. for each protein, the blast hits are recorded and sorted by their bit scores from high to low. the bit scores are then normalized by dividing them by the bit score of the query sequence in order to account for protein length variance, so that every hit has a normalized bit score within the 0– range. the organism corresponding to each hit and the taxonomic ranks of each organism are identified and recorded from the ncbi taxonomy database. to proceed with the analyses, hits of each gene have to be divided into self, close and distal groups. in other words, the pipeline doesn’t use phylogenetic trees or taxonomic lineage information directly, but rather allows the user to define three relational hierarchical categories, with the biological questions of the research in mind. this approach allows flexibility, because it can be scaled to the level of taxonomic/phylogenetic interest, and it can be adjusted to frequent updates in taxonomy. the self group is considered the recipient, and always has to include the query genome, and, depending on analytical scale may also include its immediate sister organisms. the close group will include representatives of the putative vertical inheritance history of the group. the distal group includes all other organisms, which are considered phylogenetically distant from the query genome. the method will then aim to identify genes that are likely derived from directional gene flow from groups of organisms within the distal group to members of the self group. we introduced a measure to quantify quality of blast hits of each group. this measure is calculated per group by summing up the normalized blast bit scores of hits. because three categories were defined a priori, this step will result in three weights per gene. the three weights of all genes of the query genome are considered as three independent statistical populations. if multiple genomes are analyzed together, the weight populations can be merged. the three distributions together are defined as a fingerprint of the input genome. a cutoff is set up to divide the weight distribution of each group into typical and atypical. if most genes of the genome have a vertical inheritance history, the typical portion should include a majority of the genes, while the atypical portion should only include a small, but significant subset of genes whose hits of this specific group are underrepresented. the cutoff defines the stringency of prediction: the higher the weight cutoff, the more genes are considered as being atypical. the pipeline implements statistical approaches to compute the cutoffs, but the user is free to implement and use their own statistics. because each genome or set of genomes generates a different fingerprint, cutoffs will vary, and are not transferrable across analyses. for any gene to be predicted as putatively horizontally acquired, the following rules apply, which take into account all three weight distributions for that gene: rule the gene is below cutoff in the close weight distribution, suggesting that the orthologs of the query gene are absent in all or most of the sister groups of the organism of interest. this means that the blast hits belonging to this group are significantly underrepresented, in terms of number of hits or bit score, or both. this phenomenon can be explained by: the gene lacks a vertical inheritance history, or alternatively, the gene was vertically inherited, but underwent multiple independent gene loss events in the sister groups, a case that is usually less likely to be true. in the ideal scenario of a high likelihood of hgt, the weight should be zero, meaning that there is no close hit. in real datasets however, there are sometimes sporadic close hits with low bit scores. these hits may be due to stochastic sequence similarity, secondary hgts, paralogous genes or other mechanisms. this situation typically causes false negatives in conventional best blast match methods. however, in the present method, these sporadic hits will not significantly alter the overall weight of a gene, thus hardly affecting the prediction results. on the other hand, for a vertically transmitted gene, its orthologs may not always be present in each and every sister lineage. occasional gene loss events may take away some of the expected number of close hits. this situation is a major source of false positives in best blast match methods. however, this problem is effectively overcome in the present method, because sporadic absences of hits do not make the overall weight atypical. moreover, one or a few high-score distal hits caused by natural or artificial reasons can easily deceive conventional methods. however, the present method is immune to this problem, as the close weight remains unchanged in this situation. rule the gene is equal to or above cutoff in the distal weight distribution, meaning that the hits from distant organisms are not underrepresented. this criterion sets up a filter at the donor side of potential hgt events: given the gene was transferred from a representative or ancestor of organisms that belong to the distal group, blast hits of the distal group itself should not be underrepresented. otherwise, it is not convincing to conclude that the gene is hgt-derived. the goal of this rule is stringency, in order to better distinguish putative hgt events from other scenarios that can also make a gene’s close weight atypical: de novo gene origination within the self group, inaccurate genome annotation that considers a non-coding region a gene, or hgt from an unsequenced organism that is not phylogenetically close to any sequenced groups of organism. these scenarios usually result in few to no distal hits. meanwhile, sequence similarity due to randomness instead of homology may also bring in some distal hits with low bit scores. rule the gene is below cutoff in the self weight distribution. it means that a gene is only sporadically detected, rather than being prevalent across the whole self group of organisms. this option will restrict the prediction to a subset of putatively hgt-derived genes that were acquired by specific lineages of the self group, instead of the whole self group. often, these could be recent transfer events. to assess the source of a putative hgt event, the best match organism from the distal group is reported. it is important to understand that this best match is likely not the actual physical donor, but may be an extant representative of an ancestral, extinct donor. we recommend using “donor link” to describe the directionality of transfer, and relationship between these organisms, instead of “donor”. a flowchart of the procedures of this method is illustrated in figure  computational pipeline general procedure the hgtector pipeline is written in perl, and is cross-platform supported, running in windows, mac os and linux systems. it requires only the perl interpreter with its core modules, a default component of most mac os and linux distributions and is very easy to install in windows systems. in the default mode, the program depends on no additional software or local databases to run. this characteristic maximizes the ease of installation for users without professional computer background and resources. optionally, it calls r to perform advanced statistical computing and graphing. parameters of the program are managed by a central configuration file, which can be created and edited manually or via a graphical user interface. the program is composed of the following procedures: first, the program performs batch blastp of protein products of multiple genes supplied by the user. multiple formats ranging from simple lists of ncbi accession numbers to annotated genomes in genbank format are supported. it runs blastp either via web connection to the ncbi server, or with a standalone blast program and a local database. faster alternatives to the original blastp may also be used, as long as their output files are in a compatible format. it also harvests taxonomic information of each hit for each gene from the ncbi taxonomy database. hierarchical taxonomic reports and sequences of hits can be retrieved optionally. blast results are saved in nexus format, which can be directly viewed by text editors, or opened as multiple sequence alignments by external programs such as seaview, for additional analyses. this characteristic facilitates downstream analyses, and compatibility with other programs. blast hits are filtered by multiple optional functions to overcome putative taxon sampling bias that may affect blast hit distribution: when multiple hits are present for one organism, only the best hit is maintained, representing the putative ortholog of the query protein. when multiple genomes are present for one species, only the genome that contains the best hit is maintained. taxonomic name keywords or ids that represent unwanted blast targets, such as inaccurately defined taxonomic ranks and biological categories, can be specified so that these hits will be omitted. by default, the program will exclude genes without any non-self hits from subsequent analyses, because they may represent orfans, resulting from de novo gene origination events, or transfer events from unknown sources that are very dissimilar from any sequenced genomes. alternatively, they may represent genome annotation errors, which have been long recognized as a common and perturbing issue. while these genes are not considered in the subsequent analysis, the genes are reported as “putative orfans or annotation errors”, or poe, in this pipeline. this allows the user to check which poes were omitted, and if necessary make further adjustments to the analytical set up. based on the retrieved taxonomic information, the program can automatically formulate a grouping scenario, in which the lowest rank that includes all input genomes is defined as the self group, and the next higher rank as the close group. because taxonomic classifications in genbank may not always reflect the most current natural groupings of organisms, users may manually define hierarchies based on current knowledge of phylogeny and the purpose of their research. with a properly defined grouping scenario, the program then calculates the three weights of each gene, computes a fingerprint of the whole genome, defines proper cutoffs and determines the population of atypical events, and possible hgts based on the selected rules. basic statistical parameters of the three distributions of weights, as well as the weight populations themselves are reported. the fingerprint may be visualized by box plots, histograms, density plots and scatter plots. statistical analyses and graphing of blast hit distributions are automated in the program. they are performed using perl codes, or by sending commands to r. the communication between perl and r is utilized by the perl module “statistics::r”. while multiple statistical approaches are available for the user’s choice, the typical procedures, which are used in the tests described in this article, are as follows:figure fingerprint of seven rickettsia genomes. blast hit weights of all protein-coding genes in the seven rickettsia genomes are plotted. kernel density functions of the self, close and distal weights. the x-axis represents the weight of each gene. the y-axis represents the probability density of genes with the corresponding weight in the genomes. in this example, rule and rule were applied. the close and distal cutoffs computed under the conservative criterion are indicated by dashed lines. the values of the cutoffs are denoted in each panel. a scatter plot of the distal weight against the close weight, showing the clustering pattern of the genes. each dot represents one gene. genes predicted to be hgt-derived are framed by a red rectangle. a zoom-in view of the left part of the previous plot. genes that fall within the atypical region in the close weight distribution are colored by a blue-red color scheme based on the density-based silhouette, a measure of confidence that this gene belongs to the atypical cluster of genes. the close cutoff used in the subsequent analyses is indicated by a dashed line. cutoff determination the program performs kernel density estimation to obtain a function of probability density distribution of the close weight for all genes. by default it uses gaussian kernel smoothing with silverman's rule-of-thumb bandwidth selector. the user is allowed to choose a proper bandwidth selection factor that controls the smoothness of the curve. the function is plotted and made visible to the user in real time. statistically significant local minima and maxima are computed using the “pastecs” package with default parameters following kendall’s information theory, and their x-coordinates are recorded and displayed to the user. the program then automatically identifies a local minimum separating the typical from the atypical proportion of the gene population under consideration. specifically, the biggest peak or group of continuous peaks in terms of number of genes it covers is identified as the typical region, and the rest is defined as the atypical region. in addition to the apparent typical peak, there is usually a clearly identifiable peak located close to zero. this peak usually includes “ideal” putatively hgt-derived genes. between the zero peak and the typical peak is a transitional zone that likely consists of genes with ambiguous evolutionary history. the program automatically reports two cutoffs: the x-coordinate of the identified local minimum is naturally chosen as a cutoff. however, based on results of repeated tests, we recommend using the second cutoff, which is defined as the arithmetic mean of the x-coordinates of the local maximum of the zero peak and the local minimum of the selected pit. the choice between the two cutoffs depends on the goal of research, but for the identification of putative hgt events, the conservative cutoff is preferred as it meets a balance between precision and recall. the program also implements several functions to assess the statistical significance of separating atypical genes from typical ones. for the whole weight population, the program performs hartigans’ dip test to assess the non-unimodality of the weight distribution, which essentially is the statistical significance that a distribution can be divided into two or more distinct parts. the test is performed by calling the “diptest” package in r. the dip statistic and the p-value for the test for unimodality are reported. for each individual gene, the program computes a density-based silhouette, which is a statistical measure of confidence that a certain data point is allocated to a cluster. the same procedures apply to the distribution of the distal weight. in addition to the above described statistical approaches, power users may also perform extra statistical analyses based on the weight data output by the program, and type user-defined cutoffs. based on the cutoffs, the program reports a population of genes with an atypical, non-vertical history, which in the context of the a priori provided phylogenetic information represents a putative horizontal history. the results are summarized in a choice of plain text, web page or excel spreadsheet formats. the latter allows for convenient downstream statistics of outputs, and includes hyperlinks that allow users to track each of the genes back to their original blast report. it not only reports the number and percentage of putatively hgt-derived genes, but also optionally categorizes each gene in three contexts: by putative donor group, which is described by user-designated higher taxonomic rank of the best match organism. by functional annotation of protein products, which is provided by external sources, such as the output of blast2go. by gene orthology, which is identified by a built-in function of blast hits clustering or from external sources, such as the output of orthomcl. these reports allow users an intuitive view of the prevalence of hgt-derived genes and the evolutionary, ecological and functional implications of hgts at levels of individual genes, gene families, whole genomes and multiple phylogenetically related genomes. the whole analysis can be performed on a regular personal computer, as it does not require extensive computing power. batch-blast is the most time-consuming step, which typically lasts several hours to several days, depending on the number of protein-coding genes in the genome of interest. the computation time can be reduced if faster blast alternatives or smaller databases are used. the statistical analysis typically takes only minutes. data generated by the program can be parsed and re-used by other programs for multiple purposes. as an additional, and important function, the program provides a complete phylogenetic pipeline, which automates the process of multiple sequence alignment, alignment trimming and phylogenetic tree reconstruction, by calling external local programs such as clustalw, gblocks and raxml, and parsing their outputs. reconstructed phylogenetic trees are annotated with organismal names and are attached to blast reports, which in turn can be directly viewed by external programs such as figtree. this function allows users to monitor and validate prediction results by manually checking the evolutionary scenarios of individual genes. analysis of simulated genomic data to assess the performance of this method under the impact of various evolutionary scenarios, as well as to compare hgtector to conventional blast methods, we tested the above-described pipeline on simulated genomic data. simulation of genome evolution simulated genomic data were generated by alf version, a program that simulates genome evolution. in each simulation, species evolved from one randomly generated root genome containing protein-coding genes that are no shorter than aa. during the process random inter-genomic hgt events occurred under a pre-defined global rate, which varied between simulations. the process also simulated the following evolutionary forces in addition to hgt, at random rates: speciation, character substitution, insertion and deletion, gc-content amelioration, rate variation among sites and among genes, gene duplication and loss. many of these forces are known to affect hgt prediction. two types of simulations were performed. first, an idealized, pre-defined tree topology was used, in which all representatives of close- and self-designation are grouped in a polytomy, signifying equal genetic distance, thus eliminating the impact from evolutionary artifacts. the rest of the taxa are placed relatively distant in the tree to lower the possibility of stochastic blast matches. the topology is detailed as follows: clades branch from the tree base simultaneously at time point. in each clade, species branch off simultaneously at time point the tips are at time point for all species. for each clade, one species was randomly chosen as the input genome, while nine species were defined as the close group. all species in other clades are considered as the distal group. all ten clades were analyzed in this manner. this simulation was replicated times, with the hgt rate ranging from to with an interval of. second, a randomized birth-death tree was generated in alf per simulation, mimicking a more realistic topological situation. a random clade was manually chosen from the tree, as long as it met the following criteria: 3– self species that formed a clear monophyletic group; 10– close species that were closely related to the self clade; the self and close species together formed a clear monophyletic group that was independent from all other species. this simulation was replicated for times, with the hgt rate of each replicate randomly sampled from a range of to. evolutionary events simulated in each analysis were extracted from the alf log file. the time, species and genes involved in these events were recorded. hgt events involving the self group were further categorized by their donor, recipient and time into the following groups: target hgt, ancient hgt, outgoing hgt, and secondary hgt. gene loss and gene duplication events taking place in the close group were singled out as they directly influence the distribution of the close weight. for each simulation, a blastp database including the protein sequences from all genomes was created using the standalone blast program. all-against-all blastp was performed with an e-value cutoff at 1 × 10- for all genes in the selected self genomes, considering at least hits. in order to demonstrate the effect of gene duplication on the prediction result, the program option of excluding paralogs was turned off. a modified version of the hgtector pipeline was created to parse the alf outputs and to perform analysis. both the conservative and the relaxed cutoff criteria were tested. the fingerprint was plotted along with the distribution of actual positives as well as other evolutionary scenarios for manual observation in addition to statistical analysis. comparison to conventional blast approach the performance of hgtector in the context of identifying atypical and putatively hgt-derived genes, was compared against that of commonly used conventional blast-based methods by modifying the pipeline to mimic the conventional approach, which does not consider overall hit distribution. specifically, we considered hgt-events for the scenario where no close hits are present; or the scenario where there is at least one distal hit with a bit score higher than those of any close hits. assessing performance the predicting power in the context of target hgts under each criterion, as well as for conventional blast was assessed by precision and recall. precision describes the number of true positives over the number of all predicted positives, and recall describes the number of true positives over number of all actual positives, which are commonly used performance markers for binary classifications. statistical analysis and plotting were conducted using r. application of hgtector to real datasets hgtector was further tested on a variety of real-world genomic data, covering bacteria, unicellular eukaryotes and human, the last of which serves as a realistic negative control since hgts into genomes of higher animals are known to be very rare. the most in depth analysis was conducted on the rickettsia dataset, because it affords comparisons to previous results obtained with a variety of methods.table real genomic datasets tested in this study sfg rickettsia r. conorii str. malish 1 streptococcus helicobacter erwinia mycobacterium africanum galdieria sulphuraria g. sulphuraria 3 homo sapiens 1the genomes used in this study are identical to those used in. 2two independent analyses were conducted on different dates, and similar outcomes were obtained. the more recent result was reported. 3the genome used in this study is identical to that used in. 4for genes with multiple isoforms, the longest cds was extracted using an in-house perl script and used for the analysis. analysis on the rickettsia dataset out of all available rickettsia genomes, we selected seven representative rickettsia species with fully-annotated genomes for this analysis. all of these species belong to the spotted fever group, a traditional classification group of rickettsia. a grouping scenario was chosen based on the taxonomy and phylogeny of major rickettsia species, which has been well resolved by recent studies. specifically, we defined the self group as sfg, the close group as order rickettsiales, excluding sfg, and the distal group as all non-rickettsiales organisms. all-against-all blastp was performed against the ncbi non-redundant protein sequences database with an e-value cutoff at 1 × 10- a soft filtering for low sequence complexity regions, which was suggested as the optimal parameter for ortholog identification, was used. hits with organism names containing these words were excluded for the purpose of this analysis: unknown, uncultured, unidentified, unclassified, environmental, plasmid, vector, synthetic, and phage. up to one hit from each organism was preserved. a maximum number of hits were preserved for each protein. a global fingerprint was computed and graphed to describe the pattern of blast hit distribution of all seven genomes. cutoffs for the three groups of weights were computed using the built-in kernel density estimation function under the conservative criterion. the default rules and were applied. assessing stochastic events using real datasets to test the impact of stochastic events on prediction results, the following simulations were carried out on the rickettsia dataset: database error, gene loss in the close group, rate variation in the close group, rate variation in the input genomes, taxon sampling bias, and smaller sample size of genes. multiple degrees of modification intensity for each type of events were tested, each having replicates. specifically:database error. for each hit, there is x proportion of chance that its corresponding organism was assigned to an organism randomly sampled from the pool of blast results of all genes. gene loss in the close group. for each close hit, there is x proportion of chance of being removed from the blast hit table. rate variation in the close group. for each close hit, there is x proportion of chance that its bit score is divided by a factor subject to a gamma distribution with shape parameter k =  and scale parameter θ = 1: where s and s refer to the bit score before and after manipulation, respectively. rate variation in the input genomes. for each query gene, there is x proportion of chance that the bit scores of all its hits are divided by a factor subject to a gamma distribution same as above. taxon sampling bias. for selected representative groups of organisms from close and distal groups, all hits belonging to this group were replicated into x copies. smaller sample size of genes. x out of all genes were randomly selected for fingerprint calculation in each replicate. the resulting cutoffs were in turn used for prediction on all genes. replicates that failed to pass the hartigans’ dip test were counted and excluded from the analysis. hgtector analysis was conducted on these replicates using the same procedures as in the standard analysis on the unmodified dataset. the results were compared to the results derived from a conventional blast-based approach under the d > c criterion. precision and recall of the results were computed using the result of the standard analysis as the reference. cross-method comparison in order to evaluate the performance of hgtector on the rickettsia dataset in the context of other available methods, results were compared to two examples from each of the three currently employed strategies: bbh and darkhorse based on best blast match; gist and islandviewer based on sequence composition, and two studies conducted by merhej et al. and le et al., using phylogenetic approaches. we exemplified this comparison on the r. felis genome, which previously has been demonstrated to have high hgt frequency. specifically, bbh analysis was performed as a built-in function of the present pipeline. this method is similar to the d > c criterion as described above, except for an additional reverse blast step with the same parameters to confirm that the two sequences are each other’s best match within their host genomes. the result by darkhorse was downloaded from the darkhorse server. default parameters were used, in which the blastp e-value cutoff is also 1 × 10- all three available phylogenetic granularities, strain, species and genus, were used and the results were merged, in order to maximize the discovery rate. both gist and islandviewer are targeting large pieces of heterogeneous genomic regions, or genomic islands. gist is a synchronization of five subprograms: alienhunter, islandpath, sigi-hmm, indegenius and pai-ida. the subprograms were executed in a local system and the results were processed using the egid algorithm to get a consensus result. islandviewer is an integration of three subprograms: islandpick, sigi-hmm and islandpath. the integrated result was downloaded from the islandviewer server. results of gist and islandviewer were further processed by an in-house perl script to extract the genes included in the genomic islands. the putative hgt-derived genes predicted by merhej et al. and by le et al. were extracted from the original publications. specifically, merhej et al. identified hgt-derived genes in the r. felis genome that are linked with organisms other than sfg rickettsia. le et al. identified instances of hgt from outside rickettsiales into the r. felis genome. the predicted hgt-derived genes or genomic islands by different methods were spatially mapped to the r. felis genome and visualized in geneious. an “overlap factor” was employed as a criterion to compare the outcomes of different methods by assessing the overlap. this was expressed as the negative logarithm of the likelihood that the overlap was obtained by chance. to compute an of, the number of the same genes predicted by each method pair was counted, and the of was calculated following the probability mass function of the hypergeometric distribution: where n is the total number of genes; a and b are the numbers of genes predicted by two methods, respectively; k is the number of genes overlapping by two sets of results. the larger an of is, the more overlapping, and thus more consistent the two sets of results are, and the more likely it is that the two methods are identifying the same group of genes.  performance on simulated genomic data testing precision and recall in all experimental groups under the idealized tree topology, a clear bimodal distribution was observed, which is expected when hgt is present in the data. meanwhile, none of the negative control groups have an identifiable zero peak, which is equivalent to a vertical history for all genes. both cutoff criteria achieve high precision and recall simultaneously. in particular, under the conservative criterion, % of the prediction results are true positives. meanwhile, they cover over % of all actual hgt-derived genes. the more relaxed criterion still achieves a precision of % and a recall of %. in tests under randomized tree topologies, a larger transitional zone is present between peaks of the expected bimodal distribution, which is also frequently observed in real datasets. both, precision and recall are affected compared to the idealized scenario, but in different measures. the conservative cutoff still maintains reasonably high precision and recall simultaneously, relative to the known number of hgt events. the relaxed criterion keeps equally high recall but its precision drops significantly.figure comparison of performance of hgtector and conventional blast-based method on simulated genomes. the methods were tested on the simulated genomic data under an idealized or randomized tree topology. “con” and “rel” represent conservative and relaxed cutoffs in the hgtector analysis. “c = 0” and “d > c” are two criteria in the conventional blast-based method. each experimental group is composed of tests. the distribution of results in terms of precision and recall is depicted by box plots. the mean value of each group is label above the corresponding column. because the test under randomized topologies is a more realistic representation of datasets, its result serve as a better reference for the practical consideration of the hgtector application. given the simultaneously high precision and recall, we recommend using the conservative cutoff for both initial hgt candidate screening and standalone hgt discovery. the relaxed cutoff may be considered only when the user wants to maximize discovery rate in an initial screening, in spite of its higher false positive rate. the comparison between idealized and randomized topologies further indicates the positive correlation between prediction success and a properly defined grouping scenario, in which: the self-close clade is relatively distant from any other organisms; and there are multiple subclades in the close group, each having similar number of taxa represented in the database. varying global hgt rates seemingly show little effect on the stability of precision and recall of this method, suggesting that the predicting power is independent of hgt rate. in comparison, the performance of both conventional blast-based approaches is notably unbalanced. atypical genes falling under the c =  criterion have the highest precision, but very low recall. the d > c criterion has high recall but extremely low precision, showing an intolerably high false positive rate. from a practical perspective, c =  is too stringent, thus omitting a big portion of true hgt-derived genes affected by stochastic events, while d > c is too relaxed and not capable of differentiating genes that have high-score distal hits merely due to stochastic reasons instead of hgt. it should be mentioned that c =  is not applicable to real datasets due to frequent genome annotation errors and orfans, both of which may have a zero close weight. evaluating other evolutionary scenarios the impact of other evolutionary events on the fingerprint and the division between typical and atypical gene populations was explored. outgoing hgt seemingly does not significantly alter the close weight of a gene. gene loss decreases close weight and gene duplication increases it, both within an insignificant range. most importantly, genes within these three categories of evolutionary history still fall within the typical region and were not mistakenly detected as atypical by hgtector. it is particularly notable that the majority of genes derived from ancient hgt events are located in the transitional zone. expectedly, they constitute a considerable portion of false positives in our analyses. in other words, depending on cutoff, more of them are likely to be placed in the typical population, instead of the atypical, and presumably non-vertical population. a similar pattern was observed for secondary hgt events. although not frequent in the simulated genomes here, these events are actually very frequent in real datasets, as hgt frequency is higher between closely related organisms. however, in the conventional blast method most false positives are composed of mainly outgoing hgts but only a few ancient hgts. application to real datasets the close weight distributions for real datasets exhibit a bimodal distribution containing a broad typical region and a sharp atypical peak that is located close to zero. therefore, the cutoffs that divide genes with atypical blast hit distribution from typical ones can be set accordingly and hgt events can be predicted based on the cutoffs. as an example, fingerprint plots for rickettsia and galdieria exhibit an apparent overlapping pattern between the atypical peak and the hgt-derived genes identified by previous phylogenetic studies. in contrast, the human genome does not have an apparent atypical peak, which is also expected. analysis of the rickettsia genomes the global fingerprint describing the pattern of blast hit distribution of the rickettsia analysis is illustrated in figure 3a-c. results from hartigans’ dip test strongly support the non-unimodality of all three weight distributions. the seven rickettsia genomes contain a total number of annotated chromosomal protein-coding genes, of which genes have an atypical close weight and a typical distal weight, and thus are potentially hgt-derived. the percentage of putative hgt-derived genes in a genome ranges from % in r. akari to % in r. felis. the number of putative hgt-derived genes per genome is positively correlated to the size of the genome, implying contribution of hgt to the relative genome expansions in the overall reductive trend of rickettsia genome evolution, confirming previous studies. it is notable that the r. felis genome contains significantly more putative hgt-derived genes than the rest of the genomes, suggesting a particular prevalence of potential hgt events in r. felis evolution. the comparison between the fingerprint on r. felis genome alone and the global fingerprint clearly reveals that r. felis has a much larger atypical peak in the close weight distribution. this reinforces outcomes from merhej et al., which found that the frequency of cross-species bacterial recombination in r. felis had been underestimated.table summary of genes predicted to be horizontally acquired in seven rickettsia genomes r. akari r. felis r. massiliae r. slovaca r. rickettsii r. africae r. conorii mean to further explore the biological information behind the predicted patterns, the results have been summarized in three separate contexts: by putative donor group, by functional annotation and by gene orthology. similar to previous in depth analyses by merhej et al., our analyses revealed frequent donor links, such as legionellales, enterobacteriales and burkholderiales, and frequently transferred gene categories, such as genes encoding transposases and genes involved in phage and plasmid activity. stochastic manipulation of rickettsia data: stability of prediction the results of these simulations suggest that hgtector’s performance is notably insensitive to database error. even under an extremely high proportion of mislabeled genes, its precision remains close to 100%. as a comparison, the precision of the conventional blast-based approach d > c drops below 15%. gene loss and rate variation in the close group also have remarkably weaker deleterious effects on precision of hgtector than that of conventional blast based methods under the d > c scenario. rate variations in input genomes destabilize prediction results at a certain range, but overall, precision remains high in hgtector when rate increases, as compared to conventional blast. taxon sampling bias has little effect when it happens to distal organisms. however, in some close organisms, this effect may severely compromise the predicting power if the bias is strong. smaller sample size of input genes has a minor effect on the prediction result, suggesting that the method is still effective on incomplete genomes. it should be noted that when the sample size is extremely small, the statistical power of gene clustering can be significantly impaired. the above results suggest that our method is generally unaffected by stochastic events, some of which are major challenges to current methods. a concern is taxon sampling bias in the close group, an issue that can be alleviated by properly defining grouping scenario and masking redundant taxonomic groups using hgtector’s flexible functionality. comparison of hgtector with other methods for analysis of the r. felis genome together, hgtector and six other methods identified putative hgt-derived genes, out of all chromosomal protein-coding genes in the r. felis genome. of these genes, have been uniquely identified by hgtector. there’s a considerable portion of overlapping hits between each pair of results, but none of the programs produce completely identical hgt predictions. as expected, there is relatively more overlap between two methods of the same category. meanwhile, the overlap factors between two methods from different categories are significantly lower. however, hgtector is a notable exception because its overlap with any other method of the three categories is significantly higher than between other methods from different categories.figure predicted hgts by multiple methods mapped onto the rickettsia felis genome. a circular view of the whole chromosome of r. felis, with genomic islands predicted by islandviewer and gist indicated by boxes, and putative hgt-derived genes predicted by other methods indicated by arrowheads. by comparing the results of hgtector and the results of the other two blast-based methods, it is noticeable that hgtector can effectively exclude false positives caused by database errors. for instance, bbh and darkhorse respectively predicted and genes to be acquired from ixodes scapularis, which is a common host of rickettsia. therefore, it is likely that the sequenced ixodes samples contained rickettsia endosymbionts, and its dna was also sequenced and indiscriminately labeled as ixodes dna in the unassembled shotgun sequences. in other words, these are potential database errors. in contrast, in hgtector’s result, none of these genes were predicted as hgt-derived, because they have considerable amount of close hits, despite the presence of a single high-score ixodes hit. detection of hgt is fraught with challenges, such as ambiguity in compositional features and phyletic patterns, difficulty in phylogenetic reconstruction, interference from database error and incompleteness. the observed low consistency between methods is not surprising, and has precedents in multiple previous studies. given the significant overlap of hgtector with all other previously employed approaches, we suggest that it is effective in producing meaningful prediction results on real datasets.  in this paper a novel method for genome-wide detection of vertical versus non-vertical gene history is presented. it features a statistical analysis of blast hit distribution patterns in the context of a priori defined phylogenetic hierarchies. the innovation of this method is the systematic consideration of all blast hits of all genes within selected genomes. this is in contrast to conventional blast-based approaches, which typically rely on a single best hit for each gene. the three-category grouping scenario is a simplified but effective implementation of prior phylogenetic knowledge into a blast hit distribution analysis. the set-up allows high flexibility in group assignments that best match the taxonomic level of the user’s interest, as well as immediate response to frequent changes in microbial taxonomy. the advantage of this systematic approach is that it captures the overall image of gene evolution while being significantly insensitive to stochastic events. as demonstrated in this study, stochastic events such as gene loss, rate variation and database error may impose serious problem to conventional methods, but have comparatively negligible effects on hgtector. the core assumption of the method is that the typical and atypical parts of the blast hit distribution are distinguishable. this assumption is repeatedly supported by tests on both simulated and real-world genomic data. with the proposed procedures of computing cutoffs and the rules of targeting genes that are likely to be hgt-derived, remarkable prediction success was achieved. given these results, we suggest that hgtector is a useful addition to conventional blast-based approaches. the hgtector pipeline has advantages of speed, automation, compatibility and low requirement for computational resources, making this program a generally applicable tool for discovery of vertical, and non-vertical history of genes, as well as initial hgt prediction. it is especially suitable for gaining a rapid and comprehensive overview of newly sequenced genomes to identify their evolutionary and ecological linkages with other organisms, facilitating further exploration of the functional drivers of the dynamics of genome evolution. it has to be made clear that an atypical blast hit distribution is an empirical observation, rather than a strict certification of hgt. since an hgt predominantly reflects a past evolutionary event, it is theoretically impossible to identify exact gene donors and mechanisms, and any analysis is only an approximation to possible scenarios. therefore, we recommend hgtector as a discovery tool for a detection of potential hgt-derived genes that can be further analyzed with phylogenetic approaches. this is much more effective approach than a very time-consuming and technically challenging process of a priori phylogenetic analysis of all genes within all target genomes, which becomes decreasingly feasible as more genomic data are present. electronic supplementary material additional file 1: figure s1: figure s illustration of patterns of blast hit distribution and possible explanations. figure s precision-recall plot of results on simulated genomes. figure s relationship between global hgt rate and performance of methods. figure s fingerprints of genomes of various organisms. figure s distribution of blast hit weights of the r. felis genome. figure s heat map indicating percentages of predicted hgt-derived genes by putative bacterial donor groups in rickettsia genomes. figure s heat map indicating percentages of predicted hgt-derived genes by functional annotations in rickettsia genomes. figure s stability of results on the rickettsia dataset with various simulated stochastic events. figure s stability of results on the rickettsia dataset with simulated taxon sampling bias. figure s stability of results on the rickettsia dataset with smaller sample size of genes for fingerprint calculation. figure s comparison of prediction results in the r. felis genome by multiple methods. additional file 2: table s1: protein accession numbers of predicted hgt-derived genes by hgtector. table s predicted hgt-derived genes categorized by gene orthology. table s prediction results of protein-coding genes in the r. felis genome by multiple methods. table s pairwise comparison of hgt-prediction results in the r. felis genome by multiple methods. abbreviations bbhbidirectional best match blastbasic local alignment search tool blastpprotein blast dbsdensity-based silhouette gigenomic island gogene ontology guigraphical user interface hgthorizontal gene transfer htmlhypertext markup language lpilineage probability index mlmaximum likelihood ncbinational center for biotechnology information nrnon-redundant protein sequences database ofoverlap factor poeputative orfans or annotation errors sfgspotted fever group. competing interests the authors declare that they have no competing interests. authors’ contributions qz and kd conceived and designed the study. qz developed the algorithm, wrote the program, performed the data analyses and drafted the manuscript. kd contributed to the design of algorithm and program, supervised the data analyses and edited the manuscript. mk contributed to the discussion and helped to revise the manuscript. all authors read and approved the final manuscript. 
in-depth analysis of both simulated and real genomic data suggests that the method is notably insensitive to stochastic events such as gene loss, rate variation and database error, which are common challenges to the current methodology. a new computational method of rapid, exhaustive and genome-wide detection of hgt was developed, featuring the systematic analysis of blast hit distribution patterns in the context of a priori defined hierarchical evolutionary categories. however, these methods often produce misleading results due to their inability to resolve indirect phylogenetic links and their vulnerability to stochastic events. this method was further tested on real genomic datasets, including rickettsia genomes, and was compared to previous studies. results show consistency with currently employed categories of hgt prediction methods. the program is versatile, easily deployed, has a low requirement for computational resources. this will continue to be necessary given the rapid growth of genomic data and the technical difficulties in conducting large-scale explicit phylogenetic analyses. tests on simulated genomic data suggest that this approach effectively targets atypically distributed genes that are highly likely to be hgt-derived, and exhibits robust performance compared to conventional blast-based approaches.
 first pass methods based on blast match are commonly used as an initial step to separate the different phylogenetic histories of genes in microbial genomes, and target putative horizontal gene transfer events. this will continue to be necessary given the rapid growth of genomic data and the technical difficulties in conducting large-scale explicit phylogenetic analyses. however, these methods often produce misleading results due to their inability to resolve indirect phylogenetic links and their vulnerability to stochastic events.  a new computational method of rapid, exhaustive and genome-wide detection of hgt was developed, featuring the systematic analysis of blast hit distribution patterns in the context of a priori defined hierarchical evolutionary categories. genes that fall beyond a series of statistically determined thresholds are identified as not adhering to the typical vertical history of the organisms in question, but instead having a putative horizontal origin. tests on simulated genomic data suggest that this approach effectively targets atypically distributed genes that are highly likely to be hgt-derived, and exhibits robust performance compared to conventional blast-based approaches. this method was further tested on real genomic datasets, including rickettsia genomes, and was compared to previous studies. results show consistency with currently employed categories of hgt prediction methods. in-depth analysis of both simulated and real genomic data suggests that the method is notably insensitive to stochastic events such as gene loss, rate variation and database error, which are common challenges to the current methodology. an automated pipeline was created to implement this approach and was made publicly available at: https://github.com/dittmarlab/hgtector. the program is versatile, easily deployed, has a low requirement for computational resources.  hgtector is an effective tool for initial or standalone large-scale discovery of candidate hgt-derived genes. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords horizontal gene transferbacterial genomesblastissue-copyright-statement© the author 2014 