How to align integer (non-DNA/protein) sequences?

Question

I am looking for an algorithm to find the "best" alignment between two sequences of integers similar to how one aligns nucleic acids or amino acids for homology comparisons. For example, the best alignment for the two sequences below is:

                10 06 13 91 22 16 06
                |  |  |  |  |  |  |
       01 01 02 10 04 11 92 22 17 05 

or an example with gaps:

                10 _  _  91 22 16 06 22 12
                |  |  |  |  |  |  |
          11 02 10 04 11 92 22 17 05 

The alignment would be based on the distance between each number similar to aligning amino acid sequences where one would prefer to align a hydrophobic amino acid with another hydrophobic amino acid because they are more "similar". I hope this explanation makes sense.

Another way to think about this problem is if you took the 20 amino acids and replaced them with a number from 1-20 instead of their single letter abbreviation and then tried to align 2 sequenced based on how close the numbers where, how would someone align them? In this example the alignment is based on how close two integers are from one another whereas when aligning amino acids the alignment is based on the chemical similarity of the amino acids, ie leucine is very similar to isoleucine whereas valine is very dissimilar to tryptophan.

I've been experimenting with using the biopython pairwise library and trying to utilize it for my problem. There is a similar approach to what I am looking for here:

https://stackoverflow.com/a/45512865/1255817

This solution solves my problem, but only for the single character digits [0-9] and my numbers are more than one character.

Does anyone know of an approach (algorithm, package, library) to align sequences of integers? Ideally, I am looking for something in Python with a minimal working example.

Answer 1

The tricky part about your request is the size of the alphabet. If you are happy with alphabets of up to about 250 "characters", then the programs in the FASTA package (github.com/wrpearson/fasta36) can be modified to do what you need. The are designed to work with unsigned char strings, and could work with matrices of up to 250 x 250 or so.

There is a non-vectorized version of Smith-Waterman that can calculate optimal local alignments with alphabets extended to 250 characters (with the appropriate scoring matrix).

However, I suspect that pairwise could do the same thing (with modifications), again using an expanded scoring matrix. The affine-gap Smith-Waterman algorithm is pretty simple, the trick is to get your multi-digit "characters" into a single byte (or unsigned char).

Answer 2

To parse data into pairwise2 so you can follow the SO post here fairly easy by simply converting your output values into a higher base.

from numpy import base_repr
base36 = []
for count in range(1,20):
    base36.append( base_repr(count, 36))
myseq36i = str(''.join(base36))

Output

123456789ABCDEFGHIJ

This represents 1 to 20 as single digits (of course there's 36 single digit integers) and then implement the pairwise2 solution. The advantage of this solution is you're not stuck on 20 amino acids. str may not be needed.

The custom matrix described in that post you would supply can be based around the exact base 10 differences under the scoring matrix - to avoid subtracting in the higher base.

---

Please do upvote the answer on SO as it is relevant to your problem. The first solution they presented is cool and am sure there's broader application.