用 Biopython 计算蛋白质序列一致性与覆盖度

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

>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

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

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

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

"""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())