Answers.Body: <p>So it turns out the <a href="https://academic.oup.com/nar/article/25/17/3389/1061651" rel="nofollow noreferrer">Gapped Blast paper</a> doesn't actually give the extension algorithm. It gives a high level overview and refers to an article in preparation at the time:</p>
  <blockquote>
  <p>The central idea is to consider only cells for which the optimal local alignment score falls no more than <span class="math-container">$X_g$</span> 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 <span class="math-container">$X_g$</span> 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.</p>
  </blockquote>
  <p>The motivation for this approach is to be a &quot;smarter&quot; 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 &quot;<a href="https://www.bx.psu.edu/miller_lab/dist/lowscore.pdf" rel="nofollow noreferrer">Alignments Without Low-Scoring Regions</a>&quot; comes into play. This paper is the one cited as &quot;in-prep&quot; 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 <a href="http://pipmaker.bx.psu.edu/dist/greedy.pdf" rel="nofollow noreferrer">here</a> 2 years later by the same lead author (the former paper in '98 having 67 citations and the latter in 2000 having 4700+).</p>
  <p>For completeness, I've pasted some pseudocode from the more popular x-drop paper from 2000. Note that <span class="math-container">$X$</span> is a user specified parameter, with larger <span class="math-container">$X$</span> allowing for a larger search space but slower speed.</p>
  <p><a href="https://i.sstatic.net/craWo.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/craWo.png" alt="enter image description here" /></a></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-05T12:11:42.287Z
Answers.Text: 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?
Answers.UserDisplayName: Daniel
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-05T13:31:27.103Z
Answers.Text: 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.
Answers.UserDisplayName: Throckmorton
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-05T13:32:23.967Z
Answers.Text: 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.
Answers.UserDisplayName: Throckmorton
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-05T16:43:00.910Z
Answers.Text: 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!
Answers.UserDisplayName: Daniel
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-07-03T21:10:32.483Z
Answers.LastActivityDate: 2021-07-05T13:32:02.417Z
Answers.LastEditDate: 2021-07-05T13:32:02.417Z
Answers.UserDisplayName: Throckmorton
Body: <p>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.</p>
  <p>Can anyone help me with this?</p>
  <p>Thank you all. Best,</p>
  <p>Daniel</p>
  
Tags: <blast>
Title: Pseudocode for gapped BLAST seed extension
ThreadText: # Pseudocode for gapped BLAST seed extension
  
  (Asked by: username_0 on 2021-07-02)
  
  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,
  
  username_0
  
  **Tags:** blast
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2021-07-03. Score: 3 | Accepted answer
  
  So it turns out the [Gapped Blast paper](https://academic.oup.com/nar/article/25/17/3389/1061651) 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](https://www.bx.psu.edu/miller_lab/dist/lowscore.pdf)" 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](http://pipmaker.bx.psu.edu/dist/greedy.pdf) 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](https://i.sstatic.net/craWo.png)](https://i.sstatic.net/craWo.png)
  
  > **Comments:**
  > **username_0**: 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?  
  > **username_1**: 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.  
  > **username_1**: 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.  
  > **username_0**: 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!  
  
  ---
  
  