<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>

CC BY-SA 4.0
2022-02-15T16:21:53.487Z
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.
user3754317
CC BY-SA 4.0
2022-02-15T18:05:23.890Z
Okay, the implication is there is a bug in the code, I had a brief look.
M__
CC BY-SA 4.0
2022-02-15T21:12:51.433Z
My mistake, I removed the duplicated code.
user3754317
CC BY-SA 4.0
2022-02-15T22:03:07.767Z
Okay, I looked through the code again and `in` and `dele` would a key source of bugs
M__
CC BY-SA 4.0
2022-02-15T15:01:49.853Z
2022-02-15T22:02:16.810Z
2022-02-15T22:02:16.810Z
M__
<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>

CC BY-SA 4.0
2022-02-16T14:26:29.817Z
2022-02-16T14:26:29.817Z
Daniel Standage
<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>

CC BY-SA 4.0
2022-02-15T14:49:38.177Z
The good thing is, this sure ain't a contaminant :-) Is this the untranslated locus of mtDNA?
M__
CC BY-SA 4.0
2022-02-15T14:54:28.780Z
I believe so. The fasta file states that it is a complete human mitochondrial dna sequence.
user3754317
CC BY-SA 4.0
2022-02-15T15:11:29.057Z
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.
M__
<phylogenetics><phylogeny>
Hirschberg's Alignment Algorithm Implementation. Works on the wiki example but not on large sequences of dissimilar size
# 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.


---

