PairwiseAligner

Core Principle: Bio.Align.PairwiseAligner provides a high-performance pairwise alignment engine supporting configurable global and local alignment modes, affine gap penalties, and standard substitution scoring matrices.

Basic use

Use .score() when you only need the number. It skips the traceback entirely and is substantially faster: important when scoring thousands of pairs.

Global vs local

There are also "fogsaa" and end-gap-free variants configurable through the gap-score parameters.

Substitution matrices

For protein alignment, use a real matrix rather than match/mismatch scores:

BLOSUM62 with gap open −11 / extend −1 is the BLASTP default and a sensible starting point. See Substitution Matrices.

Setting a substitution matrix overrides match_score/mismatch_score.

Fine-grained gap scoring

Setting end-gap scores to 0 gives semi-global (glocal) alignment: the whole query is aligned but the target ends are free. This is what you want for aligning a read to a reference, or a primer to a template.

Reading an alignment

aln = aligner.align(seq1, seq2)[0]

aln.score
aln.aligned            # ((t_start, t_end), ...), ((q_start, q_end), ...): blocks
aln[0]; aln[1]         # the two gapped sequence strings
str(aln)               # the pretty-printed alignment
aln.counts()           # identities, mismatches, gaps
aln.format("fasta")    # or "clustal", "phylip", "bed", "psl", "sam"
aln.substitutions      # substitution count matrix

aln.aligned is the useful programmatic output: it gives you the coordinate blocks of the ungapped segments in both sequences, which is what you need to map a position from one sequence to the other.

Percent identity is a derived quantity, and the denominator is a choice:

c = aln.counts()
pid_over_aligned = c.identities / (c.identities + c.mismatches)
pid_over_length  = c.identities / aln.length          # includes gaps

Always say which denominator you used. The two numbers can differ by 30 points on a gappy alignment, and "percent identity" alone in a methods section is ambiguous.

Many optimal alignments

alignments = aligner.align(a, b)
len(alignments)                        # can be enormous
for aln in alignments[:5]:             # SLICE: do not iterate them all
    print(aln)

alignments is a lazy object. The number of co-optimal alignments can be astronomically large for repetitive sequences, so always slice or take [0].

Behaviour change in 1.86

Biopython 1.86 changed the default gap scores to reduce the number of redundant co-optimal alignments returned. If you upgraded and your alignment counts or chosen alignment changed, that is why. Set your gap scores explicitly and the behaviour is stable across versions.

Bio.pairwise2 is deprecated

from Bio import pairwise2                    # deprecated since 1.80
pairwise2.align.globalms(a, b, 2, -1, -10, -0.5)

Slower, a less flexible API, and slated for removal. Migrate:

# old
pairwise2.align.globalms(a, b, 2, -1, -10, -0.5)

# new
aligner = Align.PairwiseAligner(match_score=2, mismatch_score=-1,
                                open_gap_score=-10, extend_gap_score=-0.5)
aligner.align(a, b)

See Biopython Deprecations.

Bioinformatics examples

from Bio import Align, SeqIO
from Bio.Align import substitution_matrices

# protein similarity search against a small database
aligner = Align.PairwiseAligner(mode="local")
aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
aligner.open_gap_score, aligner.extend_gap_score = -11, -1

scores = [(rec.id, aligner.score(query, rec.seq))
          for rec in SeqIO.parse("db.faa", "fasta")]
top = sorted(scores, key=lambda x: -x[1])[:10]

# locate a primer in a template (semi-global)
pa = Align.PairwiseAligner(mode="local", match_score=2, mismatch_score=-3,
                           open_gap_score=-5, extend_gap_score=-2)
hit = pa.align(template, primer)[0]
(t_start, t_end), = hit.aligned[0]
print(f"primer binds at {t_start}-{t_end}")

# check a variant call by aligning ref and alt context
aligner = Align.PairwiseAligner(mode="global", match_score=1, mismatch_score=-1,
                                open_gap_score=-2, extend_gap_score=-0.5)
print(aligner.align(ref_context, alt_context)[0])

# all-vs-all identity matrix for a small set of sequences
recs = list(SeqIO.parse("orthologs.faa", "fasta"))
n = len(recs)
pid = np.zeros((n, n))
for i in range(n):
    for j in range(i, n):
        c = aligner.align(recs[i].seq, recs[j].seq)[0].counts()
        pid[i, j] = pid[j, i] = c.identities / (c.identities + c.mismatches)

sns.clustermap(pd.DataFrame(pid, index=[r.id for r in recs],
                            columns=[r.id for r in recs]),
               cmap="rocket", vmin=0, vmax=1)

When to use something else

Pairwise alignment is O(mn) in time and memory. Biopython is fine for hundreds of pairs of protein-length sequences. Beyond that:

Biopython no longer wraps these: Bio.Align.Applications was removed in 1.86. Call them with subprocess. See Biopython Deprecations.

Common mistakes

• Bio.pairwise2 in new code.
• Iterating all co-optimal alignments. Slice.
• Default gap scores for protein alignment: far too permissive without a substitution matrix.
• Ambiguous percent identity. State the denominator.
• Global alignment of very different lengths. Use local or free end gaps.
• Not using .score() when the traceback is not needed.
• Pairwise alignment against a big database. Use BLAST/DIAMOND.
• Bio.Align.Applications: removed in 1.86.

See also

Substitution Matrices · AlignIO · Seq · Bio.Blast · Biopython Deprecations · Applied - Multiple Sequence Alignment Analysis

Exercise

PairwiseAligner in global mode is Needleman-Wunsch, and Biopython is not installed here, so implement it. Fill a matrix H with H[0][0] = 0, H[i][0] = i * gap, H[0][j] = j * gap and H[i][j] = max(H[i-1][j-1] + s, H[i-1][j] + gap, H[i][j-1] + gap), where s is match when the characters agree and mismatch when they do not. All three scores are integers, so the comparisons are exact.

Write align(a, b, match, mismatch, gap) returning (score, top, bottom, identity). Trace back from H[n][m] to H[0][0], preferring the diagonal step whenever it reproduces the cell's value, then the step from above (which consumes a character of a against a gap), then the step from the left. Co-optimal alignments are the rule rather than the exception, and this tie-break is what makes the answer reproducible. Identity is the number of identical columns divided by the number of columns in which neither sequence has a gap, and is 0.0 when there are none. Print the score, then the two gapped strings on their own lines, then identity=<x> rounded to three decimals. The scores arrive as match mismatch gap on the first line, the target on the second and the query on the third; an empty query line means an empty sequence.