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.
- 原文作者:春江暮客
- 原文链接:https://www.bobobk.com/en/protein-fasta-audit-python.html
- 版权声明:本作品采用 知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议 进行许可,非商业转载请注明出处(作者,原文链接),商业转载请联系作者获得授权。
相关文章
- Protein ML in Python: Build an Amino Acid Composition Baseline
- nanoBERT and VHHBERT: A Practical Guide to Nanobody Language Models
- ESM Model Family Explained: ESM-2, ESM C, ESMFold2, and ESM3
- From Protein Language Models to CoFoldArena: How to Evaluate Predictions
- Protein Language Models for Antibodies: Choose Models and Test Them Fairly