> For the complete documentation index, see [llms.txt](https://pop-lab.gitbook.io/bioinformatics-tutorials/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://pop-lab.gitbook.io/bioinformatics-tutorials/readme/how-to-interpret-database-searches.md). # How to interpret database searches Database search is one of the most common operations in bioinformatics. Despite (or perhaps because) the widespread use of database search tools, there are common misconceptions about the tools used to search databases and the interpretation of their output. Here, I try to describe some of this nuance. ### Basic Local Alignment Search Tool (BLAST) BLAST ([Altschul et al., 1990](https://pubmed.ncbi.nlm.nih.gov/2231712/)) is one of the oldest and most widely used database search tools and it underlies many bioinformatics workflows. A key innovation introduced by BLAST was a statistical model that helps distinguish "real" alignments from noise. Each database "hit" found by BLAST is accompanied by an E-value—a statistical estimation of the number of random alignments that have as high a score as the alignment between the query and the corresponding database entry. Before delving a bit deeper into this concept, I want to provide a high-level overview of how BLAST operates. More details are available on the BLAST [Wikipedia page](https://en.wikipedia.org/wiki/BLAST_\(biotechnology\)) and at [NCBI](https://blast.ncbi.nlm.nih.gov/doc/blast-help/). The goal of this tool is to find potentially imperfect matches between a query sequence and a collection of sequences organized in a database. In order to perform this task efficiently even for large databases, BLAST must take shortcuts since finding an imperfect alignment between two sequences is computationally expensive (aligning two sequences of length $$L$$ requires $$L^2$$ time in the worst case). While eventually BLAST does compute an alignment between the query sequence and the database entries it matches, the search is accelerated by initially focusing on matches between short "words" shared by the query and database entries. The matches sought are imperfect (to allow for differences between the query and database entry, particularly important for amino-acid sequences) but ungapped since most of the cost of inexact sequence alignment derives from the handling of gaps. Words that match sufficiently well are then extended to longer ungapped segments shared by the query and the database entry called "high scoring pairs" (HSPs). It is, at this point, that E-values are used to focus on just the most promising HSPs, those that suggest that the corresponding alignment does not occur by chance. A brief aside is useful at this point. Search/alignment algorithms will typically find matches whether or not they are the ones you are looking for. You have certainly experienced that when using a web search engine only to find that many of the pages returned were not relevant to your query. The E-values are intended as a way of deciding which alignments are meaningful. "Meaning", in this case, is defined in terms of the statistical estimate of the expected number of **random** sequences that would align to the query sequence with an equal or greater alignment "score". The score is an estimate of the likelihood that the query and the database hit could be transformed into each other by the process of evolution, i.e., the likelihood that the two sequences are evolutionarily related. Thus, the biological "meaning" encoded in BLAST is evolutionary relatedness between sequences. The statistical framework for defining the alignment score, however, only works for ungapped alignments and it is for this reason that BLAST performs the E-value-based filtering at the HSP stage, before gapped alignments are computed. There are several underappreciated implications of this approach: * Good alignments may be excluded because none of their HSPs have a good enough E-value/score. This could limit the sensitivity of BLAST-based searches when sequences are distantly related to each other, since we expect more gaps in alignments across long evolutionary distances. * The order in which BLAST reports results favors ungapped alignments. Thus, in certain cases, the top hit returned may not actually be the best hit. * When there are multiple equivalent alignments (at least according to BLAST's scoring) the order of the hits is related to the order of the corresponding entries in the database. Thus, the output may change as the database changes. These points are important to remember since, frequently, code that interprets the BLAST output tends to only pick the first hit for a query sequence. BLAST even provides a shortcut, the command-line parameter -max-target-seqs which can limit the number of database hits returned. One can mistakenly assume that -max-target-seqs 1 returns just the best match, and therefore miss better alignments that are reported later in the output. A long discussion of this topic is described here, and here , as well as in several papers:\ Shah, Nidhi, et al. "[Misunderstood parameter of NCBI BLAST impacts the correctness of bioinformatics workflows](https://academic.oup.com/bioinformatics/article-abstract/35/9/1613/5106166)." *Bioinformatics* 35.9 (2019): 1613-1614. Madden, Thomas L., Ben Busby, and Jian Ye. "[Reply to the paper: Misunderstood parameters of NCBI BLAST impacts the correctness of bioinformatics workflows](https://academic.oup.com/bioinformatics/article-abstract/35/15/2699/5259186)." *Bioinformatics* 35.15 (2019): 2699-2700. González-Pech, Raúl A., Timothy G. Stephens, and Cheong Xin Chan. "[Commonly misunderstood parameters of NCBI BLAST and important considerations for users.](https://academic.oup.com/bioinformatics/article-abstract/35/15/2697/5239655)" *Bioinformatics* 35.15 (2019): 2697-2698 ### The random null model is not always appropriate Above we have referred to both alignment scores and their statistical interpretation through E-values by the BLAST tool. These concepts are commonly mis-interpreted to have biological meaning, in part because the alignment scores are based on information inferred from sequences that are presumed to be related evolutionarily. While the whole point of the machinery embedded in BLAST is to make its output more biologically meaningful, it is critical to recognize that there is an imperfect and biased link between BLAST's output and actual biology. It is important to understand the assumptions that BLAST makes and how they differ from the setting in which you are working. First of all, the alignment scores are, at a high level, intended to capture evolutionary relationships between sequences. The reality is a bit more nuanced - the BLOSUM matrices from which these scores are computed, are inferred from alignments between sequences that have a certain level of sequence similarity with each other. The pairs of sequences within this "training" data set are not necessarily related to each other; they simply have a certain level of amino-acid similarity. If the goal of your study is to capture evolutionarily-related sequences, then you need to look at other tools, beyond BLAST. One should also be mindful of the fact that software and data sets contain errors: it was fairly recently discovered that the BLOSUM matrices had been mis-calculated: Styczynski, Mark P., et al. "[BLOSUM62 miscalculations improve search performance." *Nature biotechnology*](https://citeseerx.ist.psu.edu/document?repid=rep1\&type=pdf\&doi=bd6f2abba2055c18e70063ac5ee8d27b707303d0) 26.3 (2008): 274-275. Hess, Martin, et al. "[Addressing inaccuracies in BLOSUM computation improves homology search performance](https://link.springer.com/article/10.1186/s12859-016-1060-3)." *Bmc Bioinformatics* 17 (2016): 1-10. Second, as described above, the E-values represent the expected number of **random** sequences that would match the query with the same score, or higher, than that of the hit being evaluated. In many cases, we're not interested whether a particular hit is different from that to a random other sequence. For example, assume we are looking at a sequence derived from the 16S rRNA gene of a bacterium and trying to find matches to sequences related to it (e.g., from the same species). Any match between our query sequence and another 16S rRNA sequence, irrespective of the organism from which it is derived, will look to be significant, since it's highly unlikely that a random sequence would match equally well. Thus, all that E-values will tell us is that our query sequence "looks like" a 16S rRNA sequence, which is something we knew in the first place. To capture more subtle relationships, such as the relatedness between 16S rRNA sequences, one needs to use more sophisticated statistical tests. An approach, that we described in Shah, N., Altschul, S.F. & Pop, M. Outlier detection in BLAST hits. *Algorithms Mol Biol* 13, 7 (2018). , relies on an idea similar to a likelihood ratio test, comparing the likelihood that all the "top" hits to the query are similar to each other to the likelihood that a subset of the hits are more related to each other and to the query than the rest. If the latter hypothesis has a higher likelihood, then we can infer that the corresponding set of hits is related to the query in contrast to the remaining hits. ### Be careful what you search for with profile Hidden Markov Models The end of the previous section brought up the question of model comparison. In the BLAST example, the E-value is a typical statistical test, whereas the outlier detection approach used for 16S rRNA sequences compared two different "models" of the data. This latter formulation is important to keep in mind when performing searches with profile hidden Markov models (pHMMs). As a quick refresher, pHMMs are statistical models that are typically built from the multiple alignment of proteins from the same family. They are a form of a generative model, and, thus, can be used to generate sequences that "look like" their training data (the sequences in the multiple alignment). More commonly, though, they are used to interrogate whether a query sequence is consistent with the model, indicating that the query sequence may be from the same protein family as the sequences in the training set. Just as in the case of BLAST, hmmer, a search package for pHMMs, also provides estimates of the statistical significance of a match, again comparing against a random background model. Usually, users of pHMMs, also (frequently incorrectly) consider this statistical test as sufficient. A more accurate way of thinking about this is to consider each pHMM (e.g., each different protein family) as a different model from which the query sequence may have been drawn. The query sequence may match multiple profiles in a statistically significant way, and finding the correct fit for the sequence requires one to compare the likelihood of the "fit" between the sequence and the different models. A nice example of the kind of pitfalls one encounters when assuming that a statistically-significant hit to a pHMM is sufficient, is shown in the figure below:

Length of genes matching a pHMM constructed for the recA gene. The different segments of the length distribution indicate that the model captures several distinct "phenomena".

In this figure, you can see the length distribution of gene sequences that had statistically-significant matches to a pHMM constructed from the bacterial gene recA. As you can see, there are three sections of the distribution. The middle, where most of the density of the distribution is located, are the genes that are likely bacterial variants of recA. A second, lower peak to the right of this distribution turns out to represent variants of radA, an archeal gene that is related to the bacterial recA gene. Software that simply defines a match based on an E-value cutoff would miss this nuance. A likelhood ratio test comparing matches to the recA and radA pHMMs, is necessary to distinguish between these two genes. Though, in this particular case, the different length of the two genes offers another (potentially easier to implement) strategy. A third segment of the distribution is the section to the left of the main peak. These are short or fragmentary genes that happen to have significant hits to the recA gene. These could be actual fragments of the gene itself (which is quite possible in a metagenomic setting) but could also represent artifactual hits to domains of the recA gene. Again, the simple statistical estimate provided by the hmmer package is not sufficient to distinguish these from actual matches to recA, and additional checks must be implemented to avoid analysis errors.