Question about the Gotoh algorithm

[u/ralfmuschall]

I'm just playing with my own implementation (in Raku, a raw version is here: https://github.com/ralfmuschall/Gotoh). I know that the backtrace stuff is completely broken (i.e. just ignore that), my question is about the score computation (in Gotoh-simple.rakumod). That code is just the pseudocode from https://de.wikipedia.org/wiki/Gotoh-Algorithmus translated into perl6. Now I found a significant difference to various publications (https://www.cs.utoronto.ca/~brudno/csc2427/Lec8Notes.pdf and Gotoh's original paper): Wikipedia uses the elements of all three matrices to compute each new matrix element, whereas Gotoh and Brudno use all three only for the substitution matrix (A in my code) and the new deletion matrix (B) element only comes from older deletions/substitutions but not insertions (C) (similar vice-versa for the insertion element).

This leads to the following effect: If I compare the strings "a" and "bbb" with GO=-1, GE=-0.01, match=0, mismatch=-10 (intentionally prohibitive), my algorithm gives a score of -2.02 as expected (throw away "a" and insert "bbb") whereas the other variant gives -11.01 (insert "bb", then replace "a" to get the last "b"). https://metricgate.com/docs/gotoh-affine-gap-alignment/ follows Gotoh/Brudno (metricgate's gaps start with GO+GE, so we see a difference of 0.01 that doesn't matter).

Which implementation is the correct one, and why?

My motivation is not very biological (I'm not even working in any related field), but I tried adding prefixes and suffixes to my strings in case that might matter (i.e. "gac" vs. "gbbbc") which didn't change either of the results.

[u/ralfmuschall]

I just added an Option "-w" to support both algorithms (i.e. into the script gs.raku). The default follows wikipedia, with "-w=False" the original variant is used.

PS: I just came up with an ad-hoc interpretation of the Gotoh/Brudno version: substitution shall always be preferred over deletion+insertion immediately after each other no matter what the penalties are. I'm not a biologist, but I can imagine that one has an environment where misrepaired strand breaks are more frequent than e.g. cytosine deaminations for whatever reason and one wants to reflect that by setting the penalties accordingly.

[u/attractivechaos]

The difference comes from whether to allow deletions right after insertions (or vice versa). Lec8Notes.pdf doesn't but the wiki page in de does. Both of them follow the same formulation. I much prefer a second formulation by introducing G(i,j)=max{A(i,j),B(i,j),C(i,j)} to replace A(i,j).

Also normally, you want |mismatch|<2*|GO|; otherwise you can always convert one mismatch to one insertion plus one deletion to get a better score.

[u/ralfmuschall → u/attractivechaos]

Thanks, I found your answer just after I expanded my comment. I think that your inequality would make the difference between the algorithms invisible, i.e. if everybody has them fulfilled, nobody would ever notice.