Answers.Body: <p><strong>Comparisons</strong> First of all I would simply say pipe it through <code>Muscle</code> and then compare the results.</p>
  <p>I forget the algorithm advances in the 'new' (now quite old) generation of aligners, but it is worth reading up on the latest Clustal omega, possibly MAFFT and how they view <code>Muscle</code>. I do remember is 2018 Clustal Omega upgrade eulogising the historic innovation of Muscle around accuracy and describing how this was implemented. If that works for you, you simply replicate it in your code.</p>
  <p>The modern development in aligning algorithms is how they compare in speed versus number of sequences, but accuracy is always measured.</p>
  <p><strong>Amino acids -&gt; nucleotides</strong> The reason I don't follow the 'accuracy' issues is because its easier to align via amino acids and then switch back to codons. If you are sticking with this gap, penalty system using amino acids is a good approach and then simply back-translating, but if its a UTR then that will not work.</p>
  <p><strong>Brute-forcing</strong> the gap penalties is certainly an option (which is the advantages of coding it). Proving why it has worked might be a bit tricky, but its doable. You could simply fire parameters onto the data and use regex (<code>import re</code>) or straight character matching <code>==</code> to assess the improvements. Regex is better to assess improvements, whilst <code>==</code> approach is much faster. It depends on the parameter space (I think regex). Moving to an coded 'optical analysis' is an easy modification and could be useful.</p>
  <hr />
  <p>The implication from the comment is there's suboptimal code ...</p>
  <ul>
  <li>The raw default output to Biopython's pairwise should be used to debug this. Then parameters changed stepwise to understand the anomalous output.</li>
  </ul>
  <p>Skimming the code I don't spot anything obvious. What I don't understand is why the code has been duplicated</p>
  <pre><code>from Bio import SeqIO
  from Bio import Align
  from Bio import pairwise2
  insertionPenalty = -3
  ambigSub = 0
  match = 5
  deletePenalty = -3
  mismatch = -1
  
  #Ambiguous nucleotide meaning
  ambNucleotideDict = {&quot;N&quot;:[&quot;A&quot;,&quot;C&quot;,&quot;G&quot;,&quot;T&quot;,&quot;U&quot;],
                       &quot;R&quot;:[&quot;A&quot;,&quot;G&quot;],
                       &quot;Y&quot;:[&quot;T&quot;,&quot;C&quot;],
                       &quot;K&quot;:[&quot;G&quot;,&quot;T&quot;],
                       &quot;M&quot;:[&quot;A&quot;,&quot;C&quot;],
                       &quot;S&quot;:[&quot;G&quot;,&quot;C&quot;],
                       &quot;W&quot;:[&quot;A&quot;,&quot;T&quot;],
                       &quot;B&quot;:[&quot;C&quot;,&quot;G&quot;,&quot;T&quot;],
                       &quot;D&quot;:[&quot;A&quot;,&quot;G&quot;,&quot;T&quot;],
                       &quot;H&quot;:[&quot;A&quot;,&quot;C&quot;,&quot;T&quot;],
                       &quot;V&quot;:[&quot;A&quot;,&quot;C&quot;,&quot;G&quot;]}
  toPrint = False
  #Scoring functions
  def ins(y):
      return insertionPenalty
  
  
  def dele(x):
      return deletePenalty
  
  
  def sub(x, y):
      if x == y:
          return match
      elif x in ambNucleotideDict:
          if y in ambNucleotideDict[x]:
              return ambigSub
      elif y in ambNucleotideDict:
          if x in ambNucleotideDict[y]:
              return ambigSub
      else:
          return mismatch
  </code></pre>
  <p>This is present twice. I think Python would complain about the duplication of functions.
  The duplication of variables is fine and probably imports.</p>
  <hr />
  <p>A looked through the code a final time and yes there are bugs</p>
  <pre><code>def ins(y):
      return insertionPenalty
  
  
  def dele(x):
      return deletePenalty
  </code></pre>
  <p>That can't be right because whatever you call the functions (which you do), each time the output is <code>-3</code></p>
  <p>No matter what you input the output is always <code>-3</code>. If this was OOP and these were class methods in some sort of method chain that  might be different. These are functions which are operating not different to calling the variable direct. Thus if you substitute</p>
  <pre><code>... + dele(&quot;a&quot;)
  </code></pre>
  <p>with</p>
  <pre><code>... + deletePenalty
  </code></pre>
  <p>or even</p>
  <pre><code>.... + in(&quot;a&quot;) # I know this isn't supposed to be here but its to make a point
  </code></pre>
  <p>The answer in all cases is <code>-3</code>.</p>
  <p>If you asking for bugs, thats got to be involved.</p>
  <p>Out of interest are you using an IDE with breakpoints etc...?</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-02-15T16:21:53.487Z
Answers.Text: I know there are better ways to do it but I successfully used blastn to align the sequences which operates using scoring and global alignment. So it should be possible. I believe I need to create a system to allow for gaps at the beginning and the end of the smaller sequence.
Answers.UserDisplayName: user3754317
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-02-15T18:05:23.890Z
Answers.Text: Okay, the implication is there is a bug in the code, I had a brief look.
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-02-15T21:12:51.433Z
Answers.Text: My mistake, I removed the duplicated code.
Answers.UserDisplayName: user3754317
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-02-15T22:03:07.767Z
Answers.Text: Okay, I looked through the code again and `in` and `dele` would a key source of bugs
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-02-15T15:01:49.853Z
Answers.LastActivityDate: 2022-02-15T22:02:16.810Z
Answers.LastEditDate: 2022-02-15T22:02:16.810Z
Answers.UserDisplayName: M__
Answers.Body: <p>Large sequences of dissimilar size? Sounds like you want a local alignment algorithm. Hirschberg is a global alignment algorithm, which will often perform unexpectedly on large sequences of dissimilar size.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2022-02-16T14:26:29.817Z
Answers.LastActivityDate: 2022-02-16T14:26:29.817Z
Answers.UserDisplayName: Daniel Standage
Body: <p>I implemented Hirschberg's algorithm in python and used the wiki example to verify correct implementation given the scoring parameters and sequences:</p>
  <pre><code>Let X = AGTACGCA
      Y = TATGC
  Del(X) = -2
  Ins(Y) = -2
  Sub(x,y) = -1 for mismatch +2 for match
  
  This aligns to: 
  
  AGTACGCA
  --TATGC-
  </code></pre>
  <p>Which matches the wiki output. I also printed the recursion and it matched the wiki recursive partitioning too. Then I tried it on a harder problem and got a terrible result. It was a mitochondrial Neanderthal DNA sample of about 500 nucleotides and a complete human DNA mitochondrial sample of 16504 nucleotides. The output of trying it on BlastN is here:</p>
  <p><a href="https://blast.ncbi.nlm.nih.gov/Blast.cgi#alnHdr_Query_40035" rel="nofollow noreferrer">https://blast.ncbi.nlm.nih.gov/Blast.cgi#alnHdr_Query_40035</a></p>
  <p>My results look something like this where the top is neanderthal and the bottom is the full human genome. It appears though it just matches each character as soon as it finds on.</p>
  <pre><code>---C-CA-----A--------------------G-----T----AT---T--G--A----C-TC--------A--C-C-C-AT--CA-----A--C-----A-----A-C-C-----G-C-CA-T----GT-----ATT--T-C------------------G----TAC----A-T-T-AC---------T----------G------------------C--------CA--------------G-----C--CAC--C-A-T---------GA-AT-AT---------T-----------------------G-T-------ACAG---T-A--C-CAT-------AA--------------------T----------TA-C----TT----------------G-----A--C--T---A-C---C----------T--------------------------------------------------------G----T---A---A---------------T---A-C---ATA-----AA--AACC------------------------T--A---A--T--CC-----A---CA-T-CA---A-----------C---C-C-C--C-CCCC-----C-------------CC-A---T--------GCT-T-A---------CA-A-GCAAGCA-C--------AG------CA------A--TCA--AC---C--------------T--------------T-CA----A---C--T--G--T--C-A---------------------T-A----C----A----T---C-AA-------CTA--C-A-ACT--CC--A---------AA----G----A-C-ACC-C--T-------T-A-C--A---C--------CC-----A---------------C------T----A-----G-------------GA-T---A------TC-A----ACA-AA------C--C-----T------A-C----C----CAC-C-----CT-----------TG--A-------C---A----G-T-A-C---A----T-A---G------CA-CA----T-------A-AA------G---------T-----C--A-------T-----T--T--A---CC-----G-------T-------A--C-ATA------G--CA-C---A-----T-T-----A---TA------G---TC---A--AA----T--------C--C--------C-----T-------T-----C-T-----C--G------CCC---------C----------C-A--T-----------G--------GAT-G-----A-----C--------C--------------C------C-------------C---------C-------T----C----A------G------A---T----A-G-----------GG---------------------G--T--------------C----C----C-T--T--G---A---------------------------------------------------------------------------------------------------------------
  GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATCCTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACTTACTAAAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATAACAATTGAATGTCTGCACAGCCACTTTCCACACAGACATCATAACAAAAAATTTCCACCAAACCCCCCCTCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGCCAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAATTTTATCTTTTGGCGGTATGCACTTTTAACAGTCACCCCCCAACTAACACATTATTTTCCCCTCCCACTCCCATACTACTAATCTCATCAATACAACCCCCGCCCATCCTACCCAGCACACACACACCGCTGCTAACCCCATACCCCGAACCAACCAAACCCCAAAGACACCCCCCACAGTTTATGTAGCTTACCTCCTCAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCAAGCACGCAGCAATGCAGCTCAAAACGCTTAGCCTAGCCACACCCCCACGGGAAACAGCAGTGATTAACCTTTAGCAATAAACGAAAGTTTAACTAAGCTATACTAACCCCAGGGTTGGTCAATTTCGTGCCAGCCACCGCGGTCACACGATTAACCCAAGTCAATAGAAGCCGGCGTAAAGAGTGTTTTAGATCACCCCCTCCCCAATAAAGCTAAAACTCACCTGAGTTGTAAAAAACTCCAGTTGACACAAAATAGACTACGAAAGTGGCTTTAACATATCTGAACACACAATAGCTAAGACCCAAACTGGGATTAGATACCCCACTATGCTTAGCCCTAAACCTCAACAGTTAAATCAACAAAACTGCTCGCCAGAACACTACGAGCCACAGCTTAAAACTCAAAGGACCTGGCGGTGCTTCATATCCCTCTAGAGGAGCCTGTTCTGTAATCGATAAACCCCGATCAACCTCACCACCTCTTGCTCAGCCTATATACCGCCATCTTCAGCAAACCCTGATGAAGGCTACAAAGTAAGCGCAAGTACCCACGTAAAGACGTTAGGTCAAGGTGTAGCCCATGAGGTGGCAAGAAATGGGCTACATTTTCTACCCCAGAAAACTACGATAGCCCTTATGAAACTTAAGGGTCGAAGGTGGATTTAGCAGTAAACTAAGAGTAGAGTGCTTAGTTGAACAGGGCCCTGAAGCGCGTACACACCGCCCGTCACCCTCCTCAAGTATACTTCAAAGGACATTTAACTAAAACCCCTACGCATTTATATAGAGGAGACAAGTCGTAACATGGTAAGTGTACTGGAAAGTGCACTTGGACGAACCAGAGTGTAGCTTAACACAAAGCACCCAACTTACACTTAGGAGATTTCAACTTAACTTGACCGCTCTGAGCTAAACCTAGCCCCAAAC
  </code></pre>
  <p>The neanderthal sample matches at the end of the human mitochondrial dna. I suspect the gap score penalty is far too costly to insert enough gaps to reach the end of the alignment but I set blast to a linear gap penalty. How do I overcome this? Here is the implementation code for my program:</p>
  <pre><code>from Bio import SeqIO
  from Bio import Align
  from Bio import pairwise2
  insertionPenalty = -3
  ambigSub = 0
  match = 5
  deletePenalty = -3
  mismatch = -1
  
  #Ambiguous nucleotide meaning
  ambNucleotideDict = {&quot;N&quot;:[&quot;A&quot;,&quot;C&quot;,&quot;G&quot;,&quot;T&quot;,&quot;U&quot;],
                       &quot;R&quot;:[&quot;A&quot;,&quot;G&quot;],
                       &quot;Y&quot;:[&quot;T&quot;,&quot;C&quot;],
                       &quot;K&quot;:[&quot;G&quot;,&quot;T&quot;],
                       &quot;M&quot;:[&quot;A&quot;,&quot;C&quot;],
                       &quot;S&quot;:[&quot;G&quot;,&quot;C&quot;],
                       &quot;W&quot;:[&quot;A&quot;,&quot;T&quot;],
                       &quot;B&quot;:[&quot;C&quot;,&quot;G&quot;,&quot;T&quot;],
                       &quot;D&quot;:[&quot;A&quot;,&quot;G&quot;,&quot;T&quot;],
                       &quot;H&quot;:[&quot;A&quot;,&quot;C&quot;,&quot;T&quot;],
                       &quot;V&quot;:[&quot;A&quot;,&quot;C&quot;,&quot;G&quot;]}
  toPrint = False
  #Scoring functions
  def ins(y):
      return insertionPenalty
  
  
  def dele(x):
      return deletePenalty
  
  
  def sub(x, y):
      if x == y:
          return match
      elif x in ambNucleotideDict:
          if y in ambNucleotideDict[x]:
              return ambigSub
      elif y in ambNucleotideDict:
          if x in ambNucleotideDict[y]:
              return ambigSub
      else:
          return mismatch
  
  
  
  
  # NWScore
  def NWScore(X, Y):
      # Initialize NM Scoring Matrix. Array needs to be the length of X +1 and Y + 1
      score = [[0 for i in range(len(Y)+1)] for j in range(len(X)+1)]
      for j in range(1,len(Y)+1):
          # Insertion scoring penalty column
          score[0][j] = score[0][j-1] + ins(&quot;a&quot;)
      for i in range(1,len(X)+1):
  
          # Delete Penalty Row
          score[i][0] = score[i-1][0] + dele(&quot;a&quot;)
      for i in range(1,len(X)+1):
          for j in range(1,len(Y)+1):
              # Scoring Matrix
              scoreSub = score[i-1][j - 1] + sub(X[i-1], Y[j-1])
              scoreDel = score[i-1][j] + dele(X[i-1])
              scoreIns = score[i][j - 1] + ins(Y[j-1])
              score[i][j] = max(scoreSub, scoreDel, scoreIns)
      if toPrint:
          print(X,Y)
          print(&quot;Score&quot;)
          for row in score:
              print(row)
      lastLine = score[-1]
      return lastLine
  
  
  # Hirschberg
  def Hirschberg(X, Y):
      Z = &quot;&quot;
      W = &quot;&quot;
      # If sequence X has no more char in the seq
      if len(X) == 0:
          for i in range(len(Y)):
              # Concatenate a gap(?) to Z and Y[i] to W
              Z = Z + '-'
              W = W + Y[i]
      # If sequence Y has no more char in the seq
      elif len(Y) == 0:
          for i in range(len(X)):
              Z = Z + X[i]
              W = W + '-'
      elif len(X) == 1 or len(Y) == 1:
          #If one seq has a single use full global alignment. For now there is a biopython implementation of global alignment
  
          a = pairwise2.align.globalms(X,Y,2,-1,-2,-2)
          Z = a[0].seqA
          W = a[0].seqB
  
  
  
      else:
          xlen = int(len(X))
          xmid = round(len(X)/2)
          ylen = int(len(Y))
          # Ok so this is the meat of the algorithm. The last few conditionals dealt with the base cases.
          # Scoreleft takes NWScore for the first half of sequence X, and Sequence Y
  
  
          ScoreL = NWScore(X[:xmid], Y)
          # I imagine this is to start from the right side.
          ScoreR = NWScore(X[xmid:][::-1], Y[::-1])[::-1]
  
          ymid = numpy.argmax([x+y for (x,y) in zip(ScoreL,ScoreR)])
  
          # Recursion. This is divide and conquer. Creating partitions that somehow dynamically reassemble even
          # though portions of it are getting reversed?!
          Z = Z+ Hirschberg(X[:xmid],Y[:ymid])[0]+Hirschberg(X[xmid:xlen],Y[ymid:ylen])[0]
          W = W+Hirschberg(X[:xmid],Y[:ymid])[1] +  Hirschberg(X[xmid:xlen],Y[ymid:ylen])[1]
  
      return Z,W
  seq1 = SeqIO.read(&quot;seq1.fasta&quot;,&quot;fasta&quot;)
  seq2 = SeqIO.read(&quot;seq2.fasta&quot;,&quot;fasta&quot;)
  
  Z,W = Hirschberg(seq1.seq,seq2.seq)
  print(Z)
  print(W)
  score = 0
  for i in range(len(Z)):
      if Z[i] == '-' or W[i] == '-':
          score+= insertionPenalty
      elif Z[i]  == W[i]:
          score+=match
      else:
          score+=mismatch
  
  print(&quot;Score: &quot;,score)
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-02-15T14:49:38.177Z
Comments.Text: The good thing is, this sure ain't a contaminant :-) Is this the untranslated locus of mtDNA?
Comments.UserDisplayName: M__
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-02-15T14:54:28.780Z
Comments.Text: I believe so. The fasta file states that it is a complete human mitochondrial dna sequence.
Comments.UserDisplayName: user3754317
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2022-02-15T15:11:29.057Z
Comments.Text: Okay thanks, aligning via amino acids will not work then (which is what I woudl have done). The coding looks fine and its a nicely focused question.
Comments.UserDisplayName: M__
Tags: <phylogenetics><phylogeny>
Title: Hirschberg's Alignment Algorithm Implementation. Works on the wiki example but not on large sequences of dissimilar size
ThreadText: # Hirschberg's Alignment Algorithm Implementation. Works on the wiki example but not on large sequences of dissimilar size
  
  (Asked by: username_0 on 2022-02-15)
  
  I implemented Hirschberg's algorithm in python and used the wiki example to verify correct implementation given the scoring parameters and sequences:
  
  ```
  Let X = AGTACGCA
      Y = TATGC
  Del(X) = -2
  Ins(Y) = -2
  Sub(x,y) = -1 for mismatch +2 for match
  
  This aligns to: 
  
  AGTACGCA
  --TATGC-
  
  ```
  
  Which matches the wiki output. I also printed the recursion and it matched the wiki recursive partitioning too. Then I tried it on a harder problem and got a terrible result. It was a mitochondrial Neanderthal DNA sample of about 500 nucleotides and a complete human DNA mitochondrial sample of 16504 nucleotides. The output of trying it on BlastN is here:
  
  <https://blast.ncbi.nlm.nih.gov/Blast.cgi#alnHdr_Query_40035>
  
  My results look something like this where the top is neanderthal and the bottom is the full human genome. It appears though it just matches each character as soon as it finds on.
  
  ```
  ---C-CA-----A--------------------G-----T----AT---T--G--A----C-TC--------A--C-C-C-AT--CA-----A--C-----A-----A-C-C-----G-C-CA-T----GT-----ATT--T-C------------------G----TAC----A-T-T-AC---------T----------G------------------C--------CA--------------G-----C--CAC--C-A-T---------GA-AT-AT---------T-----------------------G-T-------ACAG---T-A--C-CAT-------AA--------------------T----------TA-C----TT----------------G-----A--C--T---A-C---C----------T--------------------------------------------------------G----T---A---A---------------T---A-C---ATA-----AA--AACC------------------------T--A---A--T--CC-----A---CA-T-CA---A-----------C---C-C-C--C-CCCC-----C-------------CC-A---T--------GCT-T-A---------CA-A-GCAAGCA-C--------AG------CA------A--TCA--AC---C--------------T--------------T-CA----A---C--T--G--T--C-A---------------------T-A----C----A----T---C-AA-------CTA--C-A-ACT--CC--A---------AA----G----A-C-ACC-C--T-------T-A-C--A---C--------CC-----A---------------C------T----A-----G-------------GA-T---A------TC-A----ACA-AA------C--C-----T------A-C----C----CAC-C-----CT-----------TG--A-------C---A----G-T-A-C---A----T-A---G------CA-CA----T-------A-AA------G---------T-----C--A-------T-----T--T--A---CC-----G-------T-------A--C-ATA------G--CA-C---A-----T-T-----A---TA------G---TC---A--AA----T--------C--C--------C-----T-------T-----C-T-----C--G------CCC---------C----------C-A--T-----------G--------GAT-G-----A-----C--------C--------------C------C-------------C---------C-------T----C----A------G------A---T----A-G-----------GG---------------------G--T--------------C----C----C-T--T--G---A---------------------------------------------------------------------------------------------------------------
  GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATCCTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACTTACTAAAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATAACAATTGAATGTCTGCACAGCCACTTTCCACACAGACATCATAACAAAAAATTTCCACCAAACCCCCCCTCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGCCAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAATTTTATCTTTTGGCGGTATGCACTTTTAACAGTCACCCCCCAACTAACACATTATTTTCCCCTCCCACTCCCATACTACTAATCTCATCAATACAACCCCCGCCCATCCTACCCAGCACACACACACCGCTGCTAACCCCATACCCCGAACCAACCAAACCCCAAAGACACCCCCCACAGTTTATGTAGCTTACCTCCTCAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCAAGCACGCAGCAATGCAGCTCAAAACGCTTAGCCTAGCCACACCCCCACGGGAAACAGCAGTGATTAACCTTTAGCAATAAACGAAAGTTTAACTAAGCTATACTAACCCCAGGGTTGGTCAATTTCGTGCCAGCCACCGCGGTCACACGATTAACCCAAGTCAATAGAAGCCGGCGTAAAGAGTGTTTTAGATCACCCCCTCCCCAATAAAGCTAAAACTCACCTGAGTTGTAAAAAACTCCAGTTGACACAAAATAGACTACGAAAGTGGCTTTAACATATCTGAACACACAATAGCTAAGACCCAAACTGGGATTAGATACCCCACTATGCTTAGCCCTAAACCTCAACAGTTAAATCAACAAAACTGCTCGCCAGAACACTACGAGCCACAGCTTAAAACTCAAAGGACCTGGCGGTGCTTCATATCCCTCTAGAGGAGCCTGTTCTGTAATCGATAAACCCCGATCAACCTCACCACCTCTTGCTCAGCCTATATACCGCCATCTTCAGCAAACCCTGATGAAGGCTACAAAGTAAGCGCAAGTACCCACGTAAAGACGTTAGGTCAAGGTGTAGCCCATGAGGTGGCAAGAAATGGGCTACATTTTCTACCCCAGAAAACTACGATAGCCCTTATGAAACTTAAGGGTCGAAGGTGGATTTAGCAGTAAACTAAGAGTAGAGTGCTTAGTTGAACAGGGCCCTGAAGCGCGTACACACCGCCCGTCACCCTCCTCAAGTATACTTCAAAGGACATTTAACTAAAACCCCTACGCATTTATATAGAGGAGACAAGTCGTAACATGGTAAGTGTACTGGAAAGTGCACTTGGACGAACCAGAGTGTAGCTTAACACAAAGCACCCAACTTACACTTAGGAGATTTCAACTTAACTTGACCGCTCTGAGCTAAACCTAGCCCCAAAC
  
  ```
  
  The neanderthal sample matches at the end of the human mitochondrial dna. I suspect the gap score penalty is far too costly to insert enough gaps to reach the end of the alignment but I set blast to a linear gap penalty. How do I overcome this? Here is the implementation code for my program:
  
  ```
  from Bio import SeqIO
  from Bio import Align
  from Bio import pairwise2
  insertionPenalty = -3
  ambigSub = 0
  match = 5
  deletePenalty = -3
  mismatch = -1
  
  #Ambiguous nucleotide meaning
  ambNucleotideDict = {"N":["A","C","G","T","U"],
                       "R":["A","G"],
                       "Y":["T","C"],
                       "K":["G","T"],
                       "M":["A","C"],
                       "S":["G","C"],
                       "W":["A","T"],
                       "B":["C","G","T"],
                       "D":["A","G","T"],
                       "H":["A","C","T"],
                       "V":["A","C","G"]}
  toPrint = False
  #Scoring functions
  def ins(y):
      return insertionPenalty
  
  def dele(x):
      return deletePenalty
  
  def sub(x, y):
      if x == y:
          return match
      elif x in ambNucleotideDict:
          if y in ambNucleotideDict[x]:
              return ambigSub
      elif y in ambNucleotideDict:
          if x in ambNucleotideDict[y]:
              return ambigSub
      else:
          return mismatch
  
  # NWScore
  def NWScore(X, Y):
      # Initialize NM Scoring Matrix. Array needs to be the length of X +1 and Y + 1
      score = [[0 for i in range(len(Y)+1)] for j in range(len(X)+1)]
      for j in range(1,len(Y)+1):
          # Insertion scoring penalty column
          score[0][j] = score[0][j-1] + ins("a")
      for i in range(1,len(X)+1):
  
          # Delete Penalty Row
          score[i][0] = score[i-1][0] + dele("a")
      for i in range(1,len(X)+1):
          for j in range(1,len(Y)+1):
              # Scoring Matrix
              scoreSub = score[i-1][j - 1] + sub(X[i-1], Y[j-1])
              scoreDel = score[i-1][j] + dele(X[i-1])
              scoreIns = score[i][j - 1] + ins(Y[j-1])
              score[i][j] = max(scoreSub, scoreDel, scoreIns)
      if toPrint:
          print(X,Y)
          print("Score")
          for row in score:
              print(row)
      lastLine = score[-1]
      return lastLine
  
  # Hirschberg
  def Hirschberg(X, Y):
      Z = ""
      W = ""
      # If sequence X has no more char in the seq
      if len(X) == 0:
          for i in range(len(Y)):
              # Concatenate a gap(?) to Z and Y[i] to W
              Z = Z + '-'
              W = W + Y[i]
      # If sequence Y has no more char in the seq
      elif len(Y) == 0:
          for i in range(len(X)):
              Z = Z + X[i]
              W = W + '-'
      elif len(X) == 1 or len(Y) == 1:
          #If one seq has a single use full global alignment. For now there is a biopython implementation of global alignment
  
          a = pairwise2.align.globalms(X,Y,2,-1,-2,-2)
          Z = a[0].seqA
          W = a[0].seqB
  
      else:
          xlen = int(len(X))
          xmid = round(len(X)/2)
          ylen = int(len(Y))
          # Ok so this is the meat of the algorithm. The last few conditionals dealt with the base cases.
          # Scoreleft takes NWScore for the first half of sequence X, and Sequence Y
  
          ScoreL = NWScore(X[:xmid], Y)
          # I imagine this is to start from the right side.
          ScoreR = NWScore(X[xmid:][::-1], Y[::-1])[::-1]
  
          ymid = numpy.argmax([x+y for (x,y) in zip(ScoreL,ScoreR)])
  
          # Recursion. This is divide and conquer. Creating partitions that somehow dynamically reassemble even
          # though portions of it are getting reversed?!
          Z = Z+ Hirschberg(X[:xmid],Y[:ymid])[0]+Hirschberg(X[xmid:xlen],Y[ymid:ylen])[0]
          W = W+Hirschberg(X[:xmid],Y[:ymid])[1] +  Hirschberg(X[xmid:xlen],Y[ymid:ylen])[1]
  
      return Z,W
  seq1 = SeqIO.read("seq1.fasta","fasta")
  seq2 = SeqIO.read("seq2.fasta","fasta")
  
  Z,W = Hirschberg(seq1.seq,seq2.seq)
  print(Z)
  print(W)
  score = 0
  for i in range(len(Z)):
      if Z[i] == '-' or W[i] == '-':
          score+= insertionPenalty
      elif Z[i]  == W[i]:
          score+=match
      else:
          score+=mismatch
  
  print("Score: ",score)
  
  ```
  
  **Tags:** phylogenetics, phylogeny
  
  > **Comments:**
  > **username_1**: The good thing is, this sure ain't a contaminant :-) Is this the untranslated locus of mtDNA?  
  > **username_0**: I believe so. The fasta file states that it is a complete human mitochondrial dna sequence.  
  > **username_1**: Okay thanks, aligning via amino acids will not work then (which is what I woudl have done). The coding looks fine and its a nicely focused question.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2022-02-15. Score: 1
  
  **Comparisons** First of all I would simply say pipe it through `Muscle` and then compare the results.
  
  I forget the algorithm advances in the 'new' (now quite old) generation of aligners, but it is worth reading up on the latest Clustal omega, possibly MAFFT and how they view `Muscle`. I do remember is 2018 Clustal Omega upgrade eulogising the historic innovation of Muscle around accuracy and describing how this was implemented. If that works for you, you simply replicate it in your code.
  
  The modern development in aligning algorithms is how they compare in speed versus number of sequences, but accuracy is always measured.
  
  **Amino acids -> nucleotides** The reason I don't follow the 'accuracy' issues is because its easier to align via amino acids and then switch back to codons. If you are sticking with this gap, penalty system using amino acids is a good approach and then simply back-translating, but if its a UTR then that will not work.
  
  **Brute-forcing** the gap penalties is certainly an option (which is the advantages of coding it). Proving why it has worked might be a bit tricky, but its doable. You could simply fire parameters onto the data and use regex (`import re`) or straight character matching `==` to assess the improvements. Regex is better to assess improvements, whilst `==` approach is much faster. It depends on the parameter space (I think regex). Moving to an coded 'optical analysis' is an easy modification and could be useful.
  
  ---
  
  The implication from the comment is there's suboptimal code ...
  
  * The raw default output to Biopython's pairwise should be used to debug this. Then parameters changed stepwise to understand the anomalous output.
  
  Skimming the code I don't spot anything obvious. What I don't understand is why the code has been duplicated
  
  ```
  from Bio import SeqIO
  from Bio import Align
  from Bio import pairwise2
  insertionPenalty = -3
  ambigSub = 0
  match = 5
  deletePenalty = -3
  mismatch = -1
  
  #Ambiguous nucleotide meaning
  ambNucleotideDict = {"N":["A","C","G","T","U"],
                       "R":["A","G"],
                       "Y":["T","C"],
                       "K":["G","T"],
                       "M":["A","C"],
                       "S":["G","C"],
                       "W":["A","T"],
                       "B":["C","G","T"],
                       "D":["A","G","T"],
                       "H":["A","C","T"],
                       "V":["A","C","G"]}
  toPrint = False
  #Scoring functions
  def ins(y):
      return insertionPenalty
  
  def dele(x):
      return deletePenalty
  
  def sub(x, y):
      if x == y:
          return match
      elif x in ambNucleotideDict:
          if y in ambNucleotideDict[x]:
              return ambigSub
      elif y in ambNucleotideDict:
          if x in ambNucleotideDict[y]:
              return ambigSub
      else:
          return mismatch
  
  ```
  
  This is present twice. I think Python would complain about the duplication of functions.
  The duplication of variables is fine and probably imports.
  
  ---
  
  A looked through the code a final time and yes there are bugs
  
  ```
  def ins(y):
      return insertionPenalty
  
  def dele(x):
      return deletePenalty
  
  ```
  
  That can't be right because whatever you call the functions (which you do), each time the output is `-3`
  
  No matter what you input the output is always `-3`. If this was OOP and these were class methods in some sort of method chain that might be different. These are functions which are operating not different to calling the variable direct. Thus if you substitute
  
  ```
  ... + dele("a")
  
  ```
  
  with
  
  ```
  ... + deletePenalty
  
  ```
  
  or even
  
  ```
  .... + in("a") # I know this isn't supposed to be here but its to make a point
  
  ```
  
  The answer in all cases is `-3`.
  
  If you asking for bugs, thats got to be involved.
  
  Out of interest are you using an IDE with breakpoints etc...?
  
  > **Comments:**
  > **username_0**: I know there are better ways to do it but I successfully used blastn to align the sequences which operates using scoring and global alignment. So it should be possible. I believe I need to create a system to allow for gaps at the beginning and the end of the smaller sequence.  
  > **username_1**: Okay, the implication is there is a bug in the code, I had a brief look.  
  > **username_0**: My mistake, I removed the duplicated code.  
  > **username_1**: Okay, I looked through the code again and `in` and `dele` would a key source of bugs  
  
  ---
  
  ### Answer by username_2 on 2022-02-16. Score: 0
  
  Large sequences of dissimilar size? Sounds like you want a local alignment algorithm. Hirschberg is a global alignment algorithm, which will often perform unexpectedly on large sequences of dissimilar size.
  
  
  ---
  
  