Pairwise sequence alignment and scoring

Pairwise alignment lines up two sequences so that matching letters sit in the same column. It is how we measure how alike two sequences are, even when some letters changed or went missing.

You will score an alignment by hand, then let a program find the best one.

Write one sequence above the other. Every column is one of three things. A match has the same letter in both rows. A mismatch has different letters. A gap is a dash, where one sequence has a letter and the other has none.

A gap stands for a letter that was added or lost. Gaps let sequences of different length line up.

A scoring rule gives each kind of column a number. A common choice is +1 for a match, -1 for a mismatch and -1 for a gap. A different rule gives a different score for the same alignment.

The score of the whole alignment is the total of its columns. A higher score means a better alignment.

# Score an alignment column by column.
top = "A-GT"
bottom = "ACGT"
match = 1
mismatch = -1
gap = -1
score = 0
for i in range(len(top)):
    if top[i] == "-" or bottom[i] == "-":
        score = score + gap
    elif top[i] == bottom[i]:
        score = score + match
    else:
        score = score + mismatch
print("score:", score)
# Change the top row to "AAGT" and predict the score first.

docker exec jupyter python3 -c "print(sum(1 if a == b else -1 for a, b in zip('AC-GT', 'A-AGT')))"

The command walks down the five columns of AC-GT over A-AGT. The columns score +1, -1, -1, +1 and +1.

It prints 1, the total. Here a mismatch and a gap cost the same, so one test is enough.

Many alignments are possible for two sequences, and trying every one is far too slow. Dynamic programming is the way out. It builds a table where each cell holds the best score for a pair of prefixes.

Each cell comes from one of three neighbours. A diagonal step is a match or mismatch. A step down or a step right is a gap. Take the highest.

# Fill the table: the best score for each pair of prefixes.
a = "ACGT"
b = "AGT"
match = 1
mismatch = -1
gap = -1
table = []
for i in range(len(a) + 1):
    table.append([0] * (len(b) + 1))
for i in range(1, len(a) + 1):
    table[i][0] = i * gap
for j in range(1, len(b) + 1):
    table[0][j] = j * gap
for i in range(1, len(a) + 1):
    for j in range(1, len(b) + 1):
        if a[i - 1] == b[j - 1]:
            diagonal = table[i - 1][j - 1] + match
        else:
            diagonal = table[i - 1][j - 1] + mismatch
        up = table[i - 1][j] + gap
        left = table[i][j - 1] + gap
        table[i][j] = max(diagonal, up, left)
for row in table:
    print(row)
# The bottom right cell is the best score. Change b to "ACT".