Protein Sequence Identity and Coverage with Biopython

A protein alignment can report 100% identity while matching less than half of the longer sequence. That can happen when a short query matches a region inside a larger protein. Without the alignment length, denominator, and coverage of both sequences, the percentage is easy to misread.

This tutorial uses Biopython to align two FASTA records and report exact residue matches, two identity definitions, and coverage for both proteins. The examples are synthetic and need no GPU or sequence database. It follows the FASTA audit workflow and helps make the similarity checks behind grouped protein dataset splits explicit.

1 Choose the alignment question

Use local alignment when you want the highest-scoring matching region. Use global alignment when you want an alignment spanning both complete inputs, with gaps allowed. These modes optimize a scoring function; neither is designed to maximize the identity percentage you report afterward. The Biopython alignment tutorial explains both modes and protein substitution matrices.

Here, both modes use BLOSUM62, a gap opening score of -10, and a gap extension score of -0.5. The settings apply to internal and terminal gaps. Under this affine convention, a gap of length k contributes -10 - 0.5 * (k - 1). These are explicit example settings, not a claim that one scoring policy suits every protein comparison.

Before comparing files from different tools, agree on the denominators:

The two coverage fields deliberately count paired residues. They are not alignment-span coverage. A sequence insertion inside the alignment span is excluded from paired coverage if it sits opposite gaps. The script also reports query_span and target_span using zero-based, half-open coordinates so you can inspect that difference.

2 Install Biopython and download the files

Download compare_proteins.py, motif.fasta, and insertion.fasta into a new directory. Use Python 3.10 or newer. The example was tested with Python 3.14.7 and Biopython 1.88.

python3 -m venv .venv
source .venv/bin/activate
python -m pip install "biopython==1.88"
python compare_proteins.py --help

The input must contain exactly two unaligned FASTA records: query first, target second. IDs must be nonempty and distinct. This example accepts only the 20 standard amino acid letters, converts lowercase to uppercase, and rejects ambiguous symbols, stop symbols, and input gaps. It accepts a UTF-8 byte order mark. Normal FASTA sequence whitespace is handled by the parser; for a broader input audit, use the earlier FASTA tutorial.

3 Run the local motif example

The first fixture contains an eight-residue query inside an eighteen-residue target:

>query synthetic short sequence
ACDEFGHI
>target synthetic sequence with unmatched flanks
LLLLLACDEFGHILLLLL

python compare_proteins.py motif.fasta --mode local --show > motif-local.json
python -m json.tool motif-local.json

The script writes JSON to standard output and, with --show, the readable alignment to standard error. The command therefore saves only JSON in motif-local.json; use a new output filename to preserve an earlier report. The displayed alignment is:

target            5 ACDEFGHI 13
                  0 ||||||||  8
query             0 ACDEFGHI  8

The target has ten residues outside the selected local alignment. The 100% identity describes the eight paired positions. It does not describe the whole target, establish a shared function, or provide a significance estimate. This script reports a raw alignment score, not a BLAST bit score or E-value.

4 Compare global alignment and an insertion

python compare_proteins.py motif.fasta --mode global --show > motif-global.json
python compare_proteins.py insertion.fasta --mode global --show > insertion-global.json

The second fixture inserts four alanines into a twenty-residue synthetic sequence:

>query synthetic reference
ACDEFGHIKLMNPQRSTVWY
>target synthetic sequence with a four-residue insertion
ACDEFGHIAAAAKLMNPQRSTVWY

Global alignment includes the motif target’s flanks as gap columns. It therefore changes the column-based identity while leaving paired-residue identity at 100%. In the insertion example, the target span is [0, 24], but only twenty target residues face query residues. Its span covers the entire target while its paired coverage is 83.33%.

A conservative substitution can receive a positive BLOSUM62 score and still be a mismatch for exact identity. The score helps choose the alignment; the script counts identity by comparing the aligned letters exactly.

5 Use the complete comparison script

The following code matches the downloadable file:

"""Compare exactly two protein FASTA records with explicit identity denominators."""
import argparse
from itertools import islice
import json
from pathlib import Path
import sys

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

ALPHABET = set("ACDEFGHIKLMNPQRSTVWY")


def read_pair(path, max_length):
    if max_length < 1:
        raise ValueError("max length must be positive")
    with Path(path).open(encoding="utf-8-sig") as handle:
        if handle.read(1) != ">":
            raise ValueError("FASTA must begin with a header")
        handle.seek(0)
        records = list(islice(SeqIO.parse(handle, "fasta"), 3))
    if len(records) != 2:
        raise ValueError("supply exactly two FASTA records: query, then target")
    if not all(record.id for record in records) or records[0].id == records[1].id:
        raise ValueError("record IDs must be nonempty and distinct")
    sequences = [str(record.seq).upper() for record in records]
    for record, sequence in zip(records, sequences):
        if not sequence or len(sequence) > max_length:
            raise ValueError(f"{record.id}: sequence length must be 1..{max_length}")
        invalid = set(sequence) - ALPHABET
        if invalid:
            raise ValueError(f"{record.id}: unsupported symbols {''.join(sorted(invalid))!r}")
    return records, sequences


def compare(path, mode="local", max_length=5000):
    records, (query, target) = read_pair(path, max_length)
    aligner = PairwiseAligner(mode=mode)
    aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
    aligner.open_gap_score = -10.0
    aligner.extend_gap_score = -0.5
    alignment = next(iter(aligner.align(target, query)), None)
    if alignment is None:
        raise ValueError("no positive-scoring local alignment")
    paired = matches = 0
    for (t0, t1), (q0, q1) in zip(*alignment.aligned):
        left, right = target[t0:t1], query[q0:q1]
        paired += len(left)
        matches += sum(a == b for a, b in zip(left, right))
    if paired == 0:
        raise ValueError("alignment has no residue-to-residue columns")
    columns = int(alignment.length)
    result = {
        "query_id": records[0].id, "target_id": records[1].id,
        "query_length": len(query), "target_length": len(target),
        "mode": mode, "matrix": "BLOSUM62", "gap_open": -10.0, "gap_extend": -0.5,
        "biopython": Bio.__version__, "score": float(alignment.score),
        "matches": matches, "paired_residues": paired, "alignment_columns": columns,
        "identity_paired_percent": round(100 * matches / paired, 2),
        "identity_columns_percent": round(100 * matches / columns, 2),
        "query_paired_coverage_percent": round(100 * paired / len(query), 2),
        "target_paired_coverage_percent": round(100 * paired / len(target), 2),
        "query_span": [int(x) for x in alignment.coordinates[1, [0, -1]]],
        "target_span": [int(x) for x in alignment.coordinates[0, [0, -1]]],
    }
    return result, str(alignment)


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("fasta", type=Path, help="two records: query, then target")
    parser.add_argument("--mode", choices=["local", "global"], default="local")
    parser.add_argument("--max-length", type=int, default=5000)
    parser.add_argument("--show", action="store_true", help="print alignment to stderr")
    args = parser.parse_args()
    try:
        result, alignment = compare(args.fasta, args.mode, args.max_length)
    except (OSError, ValueError) as error:
        print(f"ERROR: {error}", file=sys.stderr)
        return 1
    if args.show:
        print(alignment, file=sys.stderr)
    print(json.dumps(result, indent=2))
    return 0


if __name__ == "__main__":
    raise SystemExit(main())

Bio.SeqIO.parse reads FASTA records. The script reads at most three records to reject files that contain more than the required pair. Biopython’s Alignment.aligned property identifies paired sequence blocks, and alignment.length counts the printed alignment columns. The reported matches and coverage follow the definitions in the table above.

The script selects the first returned optimal alignment without enumerating alternatives. Ties can differ in gap placement or selected coordinates; keep the Biopython version and scoring settings with the report. Percentages are rounded to two decimal places, and the integer counts remain available for recalculation.

6 Fix invalid inputs and interpret results

The script uses dynamic programming for one sequence pair. Large inputs can consume substantial time and memory, even below the configurable limit. For database searches, use a search workflow designed for that scale. This example does not cluster a dataset, verify every train/test sequence pair, or supply a universal threshold for biological function.

If you use sequence similarity to design a machine-learning split, document the alignment mode, matrix, gap policy, identity denominator, and minimum coverage of both sequences before looking at model scores. Apply the same definitions to the cross-subset comparisons you actually perform. Distinct cluster labels alone do not establish an identity ceiling.

Summary

Report exact matches, the identity denominator, and coverage for both proteins together. Local matches and gapped global alignments can share the same paired-residue identity while describing different amounts of sequence. Preserve the input records and scoring settings so another reader can reproduce the comparison.

The fixtures and cover are illustrative; they do not represent experimental protein structures or a biological benchmark.