Waterman-Smith-Beyer implementation in Python

Question

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

I am using Freiburg University Tools 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 @gringer, 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 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:

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

Answer 1

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.