PairwiseAligner

In one line: Biopython's modern pairwise alignment engine — configurable, fast, and the replacement for the deprecated Bio.pairwise2.

Basic use

from Bio import Align

aligner = Align.PairwiseAligner()
aligner.mode = "global"                 # or "local"
aligner.match_score = 2
aligner.mismatch_score = -1
aligner.open_gap_score = -10
aligner.extend_gap_score = -0.5

score = aligner.score("ACGTACGT", "ACGTTCGT")     # score only — FAST
alignments = aligner.align("ACGTACGT", "ACGTTCGT")
print(len(alignments))                             # how many optimal alignments
best = alignments[0]
print(best)
print(best.score)

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

mode	Algorithm	Use for
"global"	Needleman–Wunsch	full-length sequences of similar length
"local"	Smith–Waterman	finding a subsequence, or sequences of very different length

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:

from Bio.Align import substitution_matrices

substitution_matrices.load()                     # list what's available
aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
aligner.open_gap_score = -11
aligner.extend_gap_score = -1

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

aligner.target_end_gap_score = 0.0      # free gaps at the ends of the target
aligner.query_end_gap_score = 0.0       # free at the ends of the query
aligner.target_internal_open_gap_score = -10
aligner.query_internal_extend_gap_score = -0.5

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.

print(aligner)      # prints every parameter — useful for checking your setup

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:

Need	Tool
Search a large database	BLAST, DIAMOND, MMseqs2
Align reads to a genome	minimap2, bwa, bowtie2
Fast pairwise in Python	parasail, edlib, mappy (minimap2 bindings)
Multiple sequence alignment	MAFFT, MUSCLE via subprocess → AlignIO
Whole-genome alignment	minimap2, LAST, MUMmer

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