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

Question

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)

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

Answer 2

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.