2
$\begingroup$

I am trying to understand in detail how does the gapped BLAST seed extension DP algorithm works. Despite searching for days now, I could not find anywhere a pseudocode that would allow implementing this in python. Not being a C programmer, it would be impossible for me to grasp this from the BLAST source code.

Can anyone help me with this?

Thank you all. Best,

Daniel

$\endgroup$
0

1 Answer 1

3
$\begingroup$

So it turns out the Gapped Blast paper doesn't actually give the extension algorithm. It gives a high level overview and refers to an article in preparation at the time:

The central idea is to consider only cells for which the optimal local alignment score falls no more than $X_g$ below the best alignment score yet found. Starting from a single aligned pair of residues, called the seed, the dynamic programming proceeds both forward and backward through the path graph (Zheng Zhang et al., manuscript in preparation) (Figs 3a and 4). The advantage of this approach is that the region of the path graph explored adapts to the alignment being constructed. The alignment can wander arbitrarily many diagonals away from the seed, but the number of cells expanded on each row tends to remain limited, and may even shrink to zero before a boundary of the path graph is encountered (Fig. 4). The $X_g$ parameter serves a similar function to the band-width parameter of the earlier heuristic, but the region of the path graph it implicitly specifies be explored is in general more productively chosen.

The motivation for this approach is to be a "smarter" version of the banded alignment approach. The issue with banded alignment approaches is that they implicitly won't allow for large insertions/deletions within the alignment, even if the areas around the alignment are high scoring. Increasing the band with to allow for this reduces the speedup gained by the band constraint. This is where the idea of "Alignments Without Low-Scoring Regions" comes into play. This paper is the one cited as "in-prep" by the gapped blast paper, and it describes an early version of the gapped extension approach. However, that approach seems to have been superseded by the one described here 2 years later by the same lead author (the former paper in '98 having 67 citations and the latter in 2000 having 4700+).

For completeness, I've pasted some pseudocode from the more popular x-drop paper from 2000. Note that $X$ is a user specified parameter, with larger $X$ allowing for a larger search space but slower speed.

enter image description here

$\endgroup$
4
  • $\begingroup$ Thanks a lot. From what I understand, though, the seed extension is nothing but a "semi-global" alignment, in which the position 0;0 of the matrix is initialized to the score of the kmer match, and the first row and column are initialized to the that score, minus the gap open + gap extend penalties. Right? $\endgroup$
    – Daniel
    Jul 5, 2021 at 12:11
  • $\begingroup$ Yes, the algorithm is semi-global in the sense that it seeks two prefixes of the input sequences that maximizes the alignment score. However, the initial value of the origin is arbitrary. The algorithm only cares that no cell in the search path falls $X_g$ below the best value seen thus far. The novelty of the X-drop algorithm in Blast is (1) the search termination bounds (i.e. the x-drop) and (2) the greedy algorithm which is faster than typical DP. $\endgroup$ Jul 5, 2021 at 13:31
  • $\begingroup$ I've updated the links to be to pdfs that don't have a paywall. Note that the DP algorithm pictured in the answer is not the main one presented by the paper. The greedy algorithm (which is more efficient) is presented in section 3. $\endgroup$ Jul 5, 2021 at 13:32
  • $\begingroup$ Thanks a lot. It's fun that nowhere it is specified that the seed-and-extend algorithm is a just a case of "semi-global" alignment. You have been superhelpful! $\endgroup$
    – Daniel
    Jul 5, 2021 at 16:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.