# Hirschberg's Alignment Algorithm Implementation. Works on the wiki example but not on large sequences of dissimilar size

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[j] = score[j-1] + ins("a")
for i in range(1,len(X)+1):

# Delete Penalty Row
score[i] = score[i-1] + 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.seqA
W = a.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])+Hirschberg(X[xmid:xlen],Y[ymid:ylen])
W = W+Hirschberg(X[:xmid],Y[:ymid]) +  Hirschberg(X[xmid:xlen],Y[ymid:ylen])

return Z,W

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)

• The good thing is, this sure ain't a contaminant :-) Is this the untranslated locus of mtDNA?
– M__
Feb 15, 2022 at 14:49
• I believe so. The fasta file states that it is a complete human mitochondrial dna sequence. Feb 15, 2022 at 14:54
• 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__
Feb 15, 2022 at 15:11

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

• 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. Feb 15, 2022 at 16:21
• Okay, the implication is there is a bug in the code, I had a brief look.
– M__
Feb 15, 2022 at 18:05
• My mistake, I removed the duplicated code. Feb 15, 2022 at 21:12
• Okay, I looked through the code again and in and dele would a key source of bugs
– M__
Feb 15, 2022 at 22:03

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.