Answers.Body: <p>This is simply the need for a gap penalty with the algorithm being used.</p>
  <p>Biopython has I believe a Waterman based aligner via Bio.pairwise2 and this should include a gap penalty. Lowering the gap penalty will sort out the artefacts you are observing, all that happening is the gap penalty is too high. It's really that simple. In bioinformatics we call them &quot;indels&quot;. It's like water off a ducks back.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-07-11T20:54:26.833Z
Answers.Text: @dawnandrew100 given this, can you try setting the default gap open penalty to 1, and see if it resolves your issue? (i.e. in the `__init__` method of the `waterman_smith_beyer` class)
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-07-10T02:12:33.060Z
Answers.LastActivityDate: 2024-07-10T02:12:33.060Z
Answers.UserDisplayName: M__
Body: <p>I am working on a text aligner to help me get a better understanding of specific steps necessary to perform sequence alignment. So far, things have been going great but I noticed yesterday that my Waterman-Smith-Beyer algorithm doesn't align as expected.</p>
  <p>The specific issue that I'm having is that when inputs are longer, all of the text is aligned to the far right, instead of in the correct place.</p>
  <p>For example - For the two input sequences &quot;HOLYWATERISABLESSING&quot; and &quot;WATERBLESSING&quot; the correct alignment should be:</p>
  <pre><code>HOLYWATERISABLESSING
  ----WATER---BLESSING
  </code></pre>
  <p>However, the alignment that I get is:</p>
  <pre><code>HOLYWATERISABLESSING
  -------WATERBLESSING
  </code></pre>
  <p>This doesn't seem to be a problem when the the sequence is short. The input sequences &quot;ABA&quot; and &quot;BA&quot; or &quot;ABA&quot; and &quot;AA&quot; all produce the correct results of</p>
  <pre><code>ABA
  -BA
  
  ABA
  A-A
  </code></pre>
  <p>Here is a link to the whole project if you want to see the full directory layout and tests I've run: <a href="https://github.com/dawnandrew100/limestone/blob/master/limestone/editdistance.py" rel="nofollow noreferrer">Limestone</a></p>
  <p>I am using <a href="http://rna.informatik.uni-freiburg.de/Teaching/index.jsp?toolName=Waterman-Smith-Beyer" rel="nofollow noreferrer">Freiburg University Tools</a> as a validation for my alignments.</p>
  <p>I should note that when the &quot;new gap penalty&quot; is removed from the SWB algorithm (basically just turning it into a Needleman-Wunsch algorithm), everything works exactly as expected. I should also note that the way the initialisation of the matrix is currently configured produces an initial output of the first column and first row which are exactly as expected.</p>
  <p>And here is the specific relevant code from my project:</p>
  <pre><code>from __future__ import annotations
  try:
      # external dependency
      import numpy
      from numpy import float64
      from numpy._typing import NDArray
  except ImportError:
      numpy = None  
  
  def main():
      qqs = &quot;HOLYWATERISABLESSING&quot;
      sss = &quot;WATERBLESSING&quot;
  
      print(watermanSmithBeyer.matrix(qqs,sss))
      print(watermanSmithBeyer.align(qqs,sss))
  
  class _GLOBALBASE():
    def matrix(self, querySequence: str, subjectSequence: str)-&gt;list[list[float]]:
      matrix, _ = self(querySequence, subjectSequence)
      return matrix
  
    def distance(self, querySequence: str, subjectSequence: str)-&gt;float:
      matrix, _ = self(querySequence, subjectSequence)
      return matrix[matrix.shape[0]-1,matrix.shape[1]-1]
  
    def similarity(self, querySequence: str, subjectSequence: str)-&gt;float:
      return max(len(querySequence),len(subjectSequence)) - self.distance(querySequence, subjectSequence)
  
    def normalized_distance(self, querySequence: str, subjectSequence: str)-&gt;float:
      dist = self.distance(querySequence, subjectSequence)
      return dist/max(map(len, [querySequence,subjectSequence]))
  
    def normalized_similarity(self, querySequence: str, subjectSequence: str)-&gt;float:
      sim = self.similarity(querySequence, subjectSequence)
      return sim/max(map(len, [querySequence,subjectSequence]))
  
    def align(self, querySequence: str, subjectSequence: str)-&gt;str: 
        qs,ss= map(lambda x: x.upper(), [querySequence,subjectSequence])
        _, pointerMatrix = self(qs, ss)
        
        qs = [x for x in qs]
        ss = [x for x in ss]
        i = len(qs)
        j = len(ss)
        queryAlign= &quot;&quot;
        subjectAlign = &quot;&quot;
  
        while i &gt; 0 or j &gt; 0: #looks for match/mismatch/gap starting from bottom right of matrix
          if pointerMatrix[i,j] in [2, 5, 6, 9]:
              #appends match/mismatch then moves to the cell diagonally up and to the left
              queryAlign = qs[i-1] + queryAlign
              subjectAlign = ss[j-1] + subjectAlign
              i -= 1
              j -= 1
          elif pointerMatrix[i,j] in [3, 5, 7, 9]:
            #appends gap and accompanying nucleotide, then moves to the cell above
              subjectAlign = '-' + subjectAlign
              queryAlign = qs[i-1] + queryAlign
              i -= 1
          elif pointerMatrix[i,j] in [4, 6, 7, 9]:
              #appends gap and accompanying nucleotide, then moves to the cell to the left
              subjectAlign = ss[j-1] + subjectAlign
              queryAlign = '-' +queryAlign
              j -= 1
  
        return f&quot;{queryAlign}\n{subjectAlign}&quot;
  
  class waterman_smith_beyer(_GLOBALBASE):
    def __init__(self, match_score:int = 0, mismatch_penalty:int = 1, new_gap_penalty:int = 3, continue_gap_penalty:int = 1)-&gt;None:
        self.match_score = match_score
        self.mismatch_penalty = mismatch_penalty
        self.new_gap_penalty = new_gap_penalty
        self.continue_gap_penalty = continue_gap_penalty
  
    def __call__(self, querySequence: str, subjectSequence: str)-&gt;tuple[NDArray[float64], NDArray[float64]]:
        if not numpy:
            raise ImportError('Please pip install numpy!')
  
        qs,ss= map(lambda x: x.upper(), [querySequence,subjectSequence])
        qs = [x for x in qs]
        ss = [x for x in ss]
        qs, ss = frontWhiteSpace(qs, ss) 
        
        #matrix initialisation
        self.alignment_score = numpy.zeros((len(qs),len(ss)))
        #pointer matrix to trace optimal alignment
        self.pointer = numpy.zeros((len(qs), len(ss))) 
        self.pointer[:,0] = 3
        self.pointer[0,:] = 4
        self.alignment_score[:,0] = [self.new_gap_penalty + n * self.continue_gap_penalty for n in range(len(qs))]
        self.alignment_score[0,:] = [self.new_gap_penalty + n * self.continue_gap_penalty for n in range(len(ss))] 
        self.alignment_score[0][0] = 0
        
        for i, subject in enumerate(qs):
          for j, query in enumerate(ss):
            if i == 0 or j == 0:
                #keeps first row and column consistent throughout all calculations
                continue
            if subject == query: 
              matchScore = self.alignment_score[i-1][j-1] - self.match_score
            else:
              matchScore = self.alignment_score[i-1][j-1] + self.mismatch_penalty
            #both gaps defaulted to continue gap penalty
            ugapScore = self.alignment_score[i-1][j] + self.continue_gap_penalty
            lgapScore = self.alignment_score[i][j-1] + self.continue_gap_penalty
            #if cell before i-1 or j-1 is gap, then this is a gap continuation
            if self.alignment_score[i-1][j] != (self.alignment_score[i-2][j]) + self.new_gap_penalty + self.continue_gap_penalty:
              ugapScore += self.new_gap_penalty
            if self.alignment_score[i][j-1] != (self.alignment_score[i][j-2]) + self.new_gap_penalty + self.continue_gap_penalty:
              lgapScore += self.new_gap_penalty
            tmin = min(matchScore, lgapScore, ugapScore)
  
            self.alignment_score[i][j] = tmin #lowest value is best choice
  
            #matrix for traceback based on results from scoring matrix
            if matchScore == tmin: 
              self.pointer[i,j] += 2
            elif ugapScore == tmin:
              self.pointer[i,j] += 3
            elif lgapScore == tmin:
              self.pointer[i,j] += 4
  
        return self.alignment_score, self.pointer
  
  watermanSmithBeyer = waterman_smith_beyer()
  if __name__ == &quot;__main__&quot;:
      main()
  </code></pre>
  <p>So far I've tried editing the order of initialisation for the matrix, changing the match score, inverting how the score is calculated (+ -&gt; -, - -&gt; +, min -&gt; max). As previously mentioned, when the match score is set arbitrarily high, the expected output displays BUT this is not a permanent fix because this then produces bugs with other inputs. Initialising the matrix differently doesn't change anything about the alignment. I've also tried adjusting the pointer matrix. I've mainly changed around the 3 and 4 between ugap and lgap but that just produces a completely wrong result.</p>
  <p>I THINK that there is some problem with the calculation of the score matrix when the new gap penalty is implemented.</p>
  <p><strong>EDITS:</strong></p>
  <p>In response to @terdon, I'm not 100% sure if there's an issue. I've been using the Freiburg tool to validate my alignments, but it's true that there may be an issue in their code, not mine. Currently, I'm treating all mismatches as the same but in the future I want users to be able to toggle matrices like the BLOSOM62 matrix. And the <code>HOLYWATERBLESSING</code> to <code>HOLYWATERISABLESSING</code> alignment does align exactly how I expect it to:</p>
  <pre><code>HOLYWATERISABLESSING
  HOLYWATER---BLESSING
  </code></pre>
  <p>And in response to @gringer, <code>HOLYWATERISABLESSING</code> and <code>HOLYERISSING</code> still have the same end alignment:</p>
  <pre><code>HOLYWATERISABLESSING
  --------HOLYERISSING
  </code></pre>
  <p>But when I add more matches to the front of the sequence, it aligns with multiple spaces:</p>
  <pre><code>HOLYWATERISABLESSING
  HOLYWAT--IS----SSING
  </code></pre>
  <p>and when I remove from the very end (to a point) it'll align with multiple spaces</p>
  <pre><code>HOLYWATERISABLESSING
  HOLYWAT--IS----SSI--
  </code></pre>
  <p>But if I remove an S or remove too many letters from the back then it all aligns to the right:</p>
  <pre><code>HOLYWATERISABLESSING
  -------HOLYWATISSING
  
  and
  
  HOLYWATERISABLESSING
  ---------HOLYWATISSS
  </code></pre>
  <p>In reply to M_'s comment to lower the gap penalty and gringer's suggestion to lower the gap penalty by the new gap penalty to 1 (instead of 3):
  This does help in most scenarios and is indeed better than a gap penalty of 3.
  However, there are some scenarios that I'm not entirely sure if it provides an accurate result such as:</p>
  <pre><code>HOLYWATERISABLESSING
  ---HWAT--IS-BLESSING
  </code></pre>
  <p>which I believe should be</p>
  <pre><code>HOLYWATERISABLESSING
  H---WAT--IS-BLESSING
  </code></pre>
  <p>and</p>
  <pre><code>HOLYWATERISABLESSING
  HOL-W----IS-BLESSING
  </code></pre>
  <p>which I believe should be</p>
  <pre><code>HOLYWATERISABLESSING
  HOLW-----IS-BLESSING
  
  OR
  
  HOLYWATERISABLESSING
  HOL-----WIS-BLESSING
  </code></pre>
  <p>I have tried playing around with the match, mismatch, opening gap, and continuing gap penalty values but it seems like most solve the first problem of the H not being aligned but don't solve the second problem of an unnecessary gap between L and W.</p>
  <p>Currently, all scores are set to 1; so matches give the cell a better score by 1, mismatches and continuing a gap worsen the score by 1 and a new gap worsens the score by 2</p>
  
Tags: <phylogenetics><python><multiple-sequence-alignment>
Title: Waterman-Smith-Beyer implementation in Python
ThreadText: # Waterman-Smith-Beyer implementation in Python
  
  (Asked by: username_0 on 2024-07-08)
  
  I am working on a text aligner to help me get a better understanding of specific steps necessary to perform sequence alignment. So far, things have been going great but I noticed yesterday that my Waterman-Smith-Beyer algorithm doesn't align as expected.
  
  The specific issue that I'm having is that when inputs are longer, all of the text is aligned to the far right, instead of in the correct place.
  
  For example - For the two input sequences "HOLYWATERISABLESSING" and "WATERBLESSING" the correct alignment should be:
  
  ```
  HOLYWATERISABLESSING
  ----WATER---BLESSING
  
  ```
  
  However, the alignment that I get is:
  
  ```
  HOLYWATERISABLESSING
  -------WATERBLESSING
  
  ```
  
  This doesn't seem to be a problem when the the sequence is short. The input sequences "ABA" and "BA" or "ABA" and "AA" all produce the correct results of
  
  ```
  ABA
  -BA
  
  ABA
  A-A
  
  ```
  
  Here is a link to the whole project if you want to see the full directory layout and tests I've run: [Limestone](https://github.com/username_0/limestone/blob/master/limestone/editdistance.py)
  
  I am using [Freiburg University Tools](http://rna.informatik.uni-freiburg.de/Teaching/index.jsp?toolName=Waterman-Smith-Beyer) as a validation for my alignments.
  
  I should note that when the "new gap penalty" is removed from the SWB algorithm (basically just turning it into a Needleman-Wunsch algorithm), everything works exactly as expected. I should also note that the way the initialisation of the matrix is currently configured produces an initial output of the first column and first row which are exactly as expected.
  
  And here is the specific relevant code from my project:
  
  ```
  from __future__ import annotations
  try:
      # external dependency
      import numpy
      from numpy import float64
      from numpy._typing import NDArray
  except ImportError:
      numpy = None  
  
  def main():
      qqs = "HOLYWATERISABLESSING"
      sss = "WATERBLESSING"
  
      print(watermanSmithBeyer.matrix(qqs,sss))
      print(watermanSmithBeyer.align(qqs,sss))
  
  class _GLOBALBASE():
    def matrix(self, querySequence: str, subjectSequence: str)->list[list[float]]:
      matrix, _ = self(querySequence, subjectSequence)
      return matrix
  
    def distance(self, querySequence: str, subjectSequence: str)->float:
      matrix, _ = self(querySequence, subjectSequence)
      return matrix[matrix.shape[0]-1,matrix.shape[1]-1]
  
    def similarity(self, querySequence: str, subjectSequence: str)->float:
      return max(len(querySequence),len(subjectSequence)) - self.distance(querySequence, subjectSequence)
  
    def normalized_distance(self, querySequence: str, subjectSequence: str)->float:
      dist = self.distance(querySequence, subjectSequence)
      return dist/max(map(len, [querySequence,subjectSequence]))
  
    def normalized_similarity(self, querySequence: str, subjectSequence: str)->float:
      sim = self.similarity(querySequence, subjectSequence)
      return sim/max(map(len, [querySequence,subjectSequence]))
  
    def align(self, querySequence: str, subjectSequence: str)->str: 
        qs,ss= map(lambda x: x.upper(), [querySequence,subjectSequence])
        _, pointerMatrix = self(qs, ss)
        
        qs = [x for x in qs]
        ss = [x for x in ss]
        i = len(qs)
        j = len(ss)
        queryAlign= ""
        subjectAlign = ""
  
        while i > 0 or j > 0: #looks for match/mismatch/gap starting from bottom right of matrix
          if pointerMatrix[i,j] in [2, 5, 6, 9]:
              #appends match/mismatch then moves to the cell diagonally up and to the left
              queryAlign = qs[i-1] + queryAlign
              subjectAlign = ss[j-1] + subjectAlign
              i -= 1
              j -= 1
          elif pointerMatrix[i,j] in [3, 5, 7, 9]:
            #appends gap and accompanying nucleotide, then moves to the cell above
              subjectAlign = '-' + subjectAlign
              queryAlign = qs[i-1] + queryAlign
              i -= 1
          elif pointerMatrix[i,j] in [4, 6, 7, 9]:
              #appends gap and accompanying nucleotide, then moves to the cell to the left
              subjectAlign = ss[j-1] + subjectAlign
              queryAlign = '-' +queryAlign
              j -= 1
  
        return f"{queryAlign}\n{subjectAlign}"
  
  class waterman_smith_beyer(_GLOBALBASE):
    def __init__(self, match_score:int = 0, mismatch_penalty:int = 1, new_gap_penalty:int = 3, continue_gap_penalty:int = 1)->None:
        self.match_score = match_score
        self.mismatch_penalty = mismatch_penalty
        self.new_gap_penalty = new_gap_penalty
        self.continue_gap_penalty = continue_gap_penalty
  
    def __call__(self, querySequence: str, subjectSequence: str)->tuple[NDArray[float64], NDArray[float64]]:
        if not numpy:
            raise ImportError('Please pip install numpy!')
  
        qs,ss= map(lambda x: x.upper(), [querySequence,subjectSequence])
        qs = [x for x in qs]
        ss = [x for x in ss]
        qs, ss = frontWhiteSpace(qs, ss) 
        
        #matrix initialisation
        self.alignment_score = numpy.zeros((len(qs),len(ss)))
        #pointer matrix to trace optimal alignment
        self.pointer = numpy.zeros((len(qs), len(ss))) 
        self.pointer[:,0] = 3
        self.pointer[0,:] = 4
        self.alignment_score[:,0] = [self.new_gap_penalty + n * self.continue_gap_penalty for n in range(len(qs))]
        self.alignment_score[0,:] = [self.new_gap_penalty + n * self.continue_gap_penalty for n in range(len(ss))] 
        self.alignment_score[0][0] = 0
        
        for i, subject in enumerate(qs):
          for j, query in enumerate(ss):
            if i == 0 or j == 0:
                #keeps first row and column consistent throughout all calculations
                continue
            if subject == query: 
              matchScore = self.alignment_score[i-1][j-1] - self.match_score
            else:
              matchScore = self.alignment_score[i-1][j-1] + self.mismatch_penalty
            #both gaps defaulted to continue gap penalty
            ugapScore = self.alignment_score[i-1][j] + self.continue_gap_penalty
            lgapScore = self.alignment_score[i][j-1] + self.continue_gap_penalty
            #if cell before i-1 or j-1 is gap, then this is a gap continuation
            if self.alignment_score[i-1][j] != (self.alignment_score[i-2][j]) + self.new_gap_penalty + self.continue_gap_penalty:
              ugapScore += self.new_gap_penalty
            if self.alignment_score[i][j-1] != (self.alignment_score[i][j-2]) + self.new_gap_penalty + self.continue_gap_penalty:
              lgapScore += self.new_gap_penalty
            tmin = min(matchScore, lgapScore, ugapScore)
  
            self.alignment_score[i][j] = tmin #lowest value is best choice
  
            #matrix for traceback based on results from scoring matrix
            if matchScore == tmin: 
              self.pointer[i,j] += 2
            elif ugapScore == tmin:
              self.pointer[i,j] += 3
            elif lgapScore == tmin:
              self.pointer[i,j] += 4
  
        return self.alignment_score, self.pointer
  
  watermanSmithBeyer = waterman_smith_beyer()
  if __name__ == "__main__":
      main()
  
  ```
  
  So far I've tried editing the order of initialisation for the matrix, changing the match score, inverting how the score is calculated (+ -> -, - -> +, min -> max). As previously mentioned, when the match score is set arbitrarily high, the expected output displays BUT this is not a permanent fix because this then produces bugs with other inputs. Initialising the matrix differently doesn't change anything about the alignment. I've also tried adjusting the pointer matrix. I've mainly changed around the 3 and 4 between ugap and lgap but that just produces a completely wrong result.
  
  I THINK that there is some problem with the calculation of the score matrix when the new gap penalty is implemented.
  
  **EDITS:**
  
  In response to @terdon, I'm not 100% sure if there's an issue. I've been using the Freiburg tool to validate my alignments, but it's true that there may be an issue in their code, not mine. Currently, I'm treating all mismatches as the same but in the future I want users to be able to toggle matrices like the BLOSOM62 matrix. And the `HOLYWATERBLESSING` to `HOLYWATERISABLESSING` alignment does align exactly how I expect it to:
  
  ```
  HOLYWATERISABLESSING
  HOLYWATER---BLESSING
  
  ```
  
  And in response to @username_2, `HOLYWATERISABLESSING` and `HOLYERISSING` still have the same end alignment:
  
  ```
  HOLYWATERISABLESSING
  --------HOLYERISSING
  
  ```
  
  But when I add more matches to the front of the sequence, it aligns with multiple spaces:
  
  ```
  HOLYWATERISABLESSING
  HOLYWAT--IS----SSING
  
  ```
  
  and when I remove from the very end (to a point) it'll align with multiple spaces
  
  ```
  HOLYWATERISABLESSING
  HOLYWAT--IS----SSI--
  
  ```
  
  But if I remove an S or remove too many letters from the back then it all aligns to the right:
  
  ```
  HOLYWATERISABLESSING
  -------HOLYWATISSING
  
  and
  
  HOLYWATERISABLESSING
  ---------HOLYWATISSS
  
  ```
  
  In reply to M\_'s comment to lower the gap penalty and username_2's suggestion to lower the gap penalty by the new gap penalty to 1 (instead of 3):
  This does help in most scenarios and is indeed better than a gap penalty of 3.
  However, there are some scenarios that I'm not entirely sure if it provides an accurate result such as:
  
  ```
  HOLYWATERISABLESSING
  ---HWAT--IS-BLESSING
  
  ```
  
  which I believe should be
  
  ```
  HOLYWATERISABLESSING
  H---WAT--IS-BLESSING
  
  ```
  
  and
  
  ```
  HOLYWATERISABLESSING
  HOL-W----IS-BLESSING
  
  ```
  
  which I believe should be
  
  ```
  HOLYWATERISABLESSING
  HOLW-----IS-BLESSING
  
  OR
  
  HOLYWATERISABLESSING
  HOL-----WIS-BLESSING
  
  ```
  
  I have tried playing around with the match, mismatch, opening gap, and continuing gap penalty values but it seems like most solve the first problem of the H not being aligned but don't solve the second problem of an unnecessary gap between L and W.
  
  Currently, all scores are set to 1; so matches give the cell a better score by 1, mismatches and continuing a gap worsen the score by 1 and a new gap worsens the score by 2
  
  **Tags:** phylogenetics, python, multiple-sequence-alignment
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2024-07-10. Score: 0
  
  This is simply the need for a gap penalty with the algorithm being used.
  
  Biopython has I believe a Waterman based aligner via Bio.pairwise2 and this should include a gap penalty. Lowering the gap penalty will sort out the artefacts you are observing, all that happening is the gap penalty is too high. It's really that simple. In bioinformatics we call them "indels". It's like water off a ducks back.
  
  > **Comments:**
  > **username_2**: @username_0 given this, can you try setting the default gap open penalty to 1, and see if it resolves your issue? (i.e. in the `__init__` method of the `waterman_smith_beyer` class)  
  
  ---
  
  