春江暮客

春江暮客的个人学习分享网站

Audit Protein FASTA Files with Python Before Generating Embeddings

2026-09-23 Technology
Audit Protein FASTA Files with Python Before Generating Embeddings

Before running a protein encoder, check what is actually in the FASTA file. Two different record IDs can contain the same sequence. One ID can also appear twice with different sequences. Both cases make it easy to join an embedding to the wrong measurement or place repeated inputs on opposite sides of an evaluation split.

This tutorial builds a small audit tool using only Python’s standard library. It produces a CSV report without changing the input. Use it before the amino acid composition baseline or the ESM-2 embeddings workflow.

Define the input policy first

A FASTA record begins with a description line starting with >, followed by sequence lines. The NCBI format guide describes that structure and the symbols accepted by BLAST. A model pipeline may impose a narrower policy.

This example accepts the 20 standard amino acid letters, converts lowercase to uppercase for comparison, and joins wrapped sequence lines. It skips blank lines and accepts a UTF-8 byte-order mark. These are explicit parser choices, not a claim that every FASTA consumer behaves the same way.

Check Reported behavior
Identifier The first whitespace-separated field after >; it must be unique within the file.
Empty sequence Flag it for review.
Symbols outside the 20-letter alphabet Flag them; do not delete or substitute residues.
Spaces or tabs inside sequence lines Flag them, including leading and trailing spaces.
Length above the supplied limit Flag it; do not truncate.
Identical canonical sequences Point later records to the first matching record number.

X, B, Z, U, O, gaps, and stop symbols fail this deliberately narrow policy. That does not make every such sequence biologically invalid. Check the dataset and the selected model’s tokenizer before deciding how to handle them. Even an all-ACGT string passes this alphabet check; the tool cannot establish that the input is a protein rather than a nucleotide sequence.

The length limit is required because it belongs to your chosen workflow. A tokenizer’s token budget may include special tokens, and some encoders use different tokenization schemes. Do not treat a model’s advertised token limit as an automatic residue limit.

Run a small, inspectable example

Download audit_fasta.py and demo.fasta into the same directory. No third-party packages are needed. The script was tested with Python 3.14.7.

The demonstration contains artificial strings for software testing:

>sample_a synthetic example
ACDEFGHIKL
MNPQRSTVWY
>sample_b same sequence in lowercase
acdefghiklmnpqrstvwy
>sample_c requires a residue policy
ACDXEFG
>sample_a reused identifier
MKT
>sample_empty

Run:

python3 audit_fasta.py demo.fasta --max-residues 20 --out audit.csv

The limit of 20 keeps this example small; it is not a recommended model setting. Output:

records=5 flagged=4 exact_repeats=1
report=audit.csv

Both occurrences of sample_a are flagged because their IDs collide. sample_c contains X, and sample_empty has no residues. sample_b repeats record 1 after uppercase conversion; the duplicate reference appears in a separate column.

The exit status is 0 when there are no flagged records or exact repeats, 1 when a completed report needs review, and 2 for a parsing, argument, or file error. This demonstration therefore exits with status 1 even though it successfully creates the report. An existing report is never overwritten; choose a fresh --out path for another run.

The complete audit script

"""Audit protein FASTA under an explicit 20-residue policy; never rewrite input."""
import argparse
import csv
import hashlib
import sys
from collections import Counter
from pathlib import Path

AA = set("ACDEFGHIKLMNPQRSTVWY")


def records(path):
    header, chunks = None, []
    with path.open(encoding="utf-8-sig") as handle:
        for line_no, raw in enumerate(handle, 1):
            line = raw.rstrip("\r\n")
            if not line.strip():
                continue
            if line.startswith(">"):
                if header is not None:
                    yield header, "".join(chunks)
                header, chunks = line[1:].strip(), []
                if not header:
                    raise ValueError(f"line {line_no}: empty FASTA header")
            else:
                if header is None:
                    raise ValueError(f"line {line_no}: sequence before first header")
                chunks.append(line)
    if header is not None:
        yield header, "".join(chunks)


def audit(path, max_residues):
    rows, first_sequence = [], {}
    for number, (header, raw) in enumerate(records(path), 1):
        ident = header.split()[0]
        sequence = raw.upper()
        problems = []
        if not sequence:
            problems.append("empty_sequence")
        if any(c.isspace() for c in raw):
            problems.append("whitespace_in_sequence")
        invalid = sorted(set(raw) - AA - set("acdefghiklmnpqrstvwy"))
        if invalid:
            problems.append("unsupported=" + repr("".join(invalid)))
        if len(sequence) > max_residues:
            problems.append("over_limit")
        # Hash and duplicate comparison only for canonical, nonempty sequences.
        canonical = bool(sequence) and not invalid
        digest = hashlib.sha256(sequence.encode("ascii")).hexdigest() if canonical else ""
        duplicate = ""
        if canonical:
            duplicate = first_sequence.get(sequence, "")
            first_sequence.setdefault(sequence, number)
        rows.append(dict(record=number, id=ident, length=len(sequence),
                         lowercase=raw != sequence, sha256=digest,
                         duplicate_of_record=duplicate, issues=";".join(problems)))
    if not rows:
        raise ValueError("no FASTA records")
    counts = Counter(row["id"] for row in rows)
    for row in rows:
        if counts[row["id"]] > 1:
            row["issues"] = ";".join(filter(None, [row["issues"], "duplicate_id"]))
    return rows


def main():
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument("fasta", type=Path)
    parser.add_argument("--max-residues", type=int, required=True)
    parser.add_argument("--out", type=Path, default=Path("audit.csv"))
    args = parser.parse_args()
    if args.max_residues < 1:
        parser.error("--max-residues must be positive")
    try:
        rows = audit(args.fasta, args.max_residues)
        with args.out.open("x", newline="", encoding="utf-8") as handle:
            writer = csv.DictWriter(handle, fieldnames=list(rows[0]))
            writer.writeheader()
            writer.writerows(rows)
    except (OSError, ValueError) as error:
        parser.exit(2, f"error: {error}\n")
    flagged = sum(bool(row["issues"]) for row in rows)
    repeats = sum(row["duplicate_of_record"] != "" for row in rows)
    print(f"records={len(rows)} flagged={flagged} exact_repeats={repeats}")
    print(f"report={args.out}")
    return 1 if flagged or repeats else 0


if __name__ == "__main__":
    sys.exit(main())

The script uses csv.DictWriter with newline="" to write the report. Sequence fingerprints come from hashlib.sha256. They describe the uppercase, joined canonical sequence, not the original file bytes, header, or measurement.

Duplicate detection compares the actual normalized sequence strings. Unsupported and empty sequences receive no fingerprint and are excluded from that comparison. A canonical sequence above the chosen length limit still receives a fingerprint and can be marked as a repeat.

The tool reads the input line by line but retains report rows and distinct canonical sequences in memory. Use it for datasets that fit comfortably in RAM. For a large corpus, move the ID and sequence index into a disk-backed store and preserve the same validation rules.

Turn the report into a reproducible dataset

Review duplicate IDs before joining labels or embeddings. The script keeps the entire first header field, including any pipe characters; it does not infer a database accession. If your labels use a different identifier scheme, write and check that mapping explicitly.

Repeated sequences need context. They may represent redundant downloads, repeated measurements, or the same protein tested under different conditions. Keep the measurement provenance. When evaluating generalization to unseen sequences, assign identical sequences to the same split, and also account for related sequences. Exact matching cannot detect homology, shared antibody lineages, or antigen overlap; the antibody evaluation guide discusses those broader grouping choices.

After resolving the report, save the source FASTA, the chosen normalization policy, the final ID mapping, and the split assignments. If you change a sequence, regenerate its fingerprint and its embedding. Join by a checked identifier and sequence fingerprint rather than relying on array row order.

Passing this audit does not validate labels or prove that a dataset is free of leakage. Learned transformations such as scaling and feature selection still need to be fitted on training data only, as described in scikit-learn’s leakage guidance.

Troubleshooting

Message or observation What to check
sequence before first header The file must begin with a FASTA header before any sequence data. Headers must start in column one.
empty FASTA header Add a meaningful identifier after >.
duplicate_id Inspect both records; do not silently keep the last one.
unsupported= Check the raw symbols and document a dataset-specific policy.
over_limit Review the encoder’s usable residue budget and choose an explicit long-sequence strategy.
Exit status 1 Open the completed CSV; flagged rows or repeated sequences need review.

Once the report is resolved, run the composition baseline and save its split. That gives you a traceable input dataset for the next embedding experiment.

Sources checked September 23, 2026. The cover is AI-generated conceptual artwork.

友情链接

其它