# Bio Core Biopython Essentials

> Manipulate Seq/SeqRecord objects, parse FASTA/FASTQ/GenBank with SeqIO, query NCBI via Entrez, and run PairwiseAligner in Biopython. Use for sequence I/O, translation, reverse complement, GC content, or NCBI fetch in Python.

- Skill: `pavel-kravchenko/bio-core-biopython-essentials` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-core-biopython-essentials`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-core-biopython-essentials/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-core-biopython-essentials

---


# BioPython Essentials

## When to Use

- Reading/writing FASTA, FASTQ, GenBank, EMBL, or Stockholm files with `SeqIO`.
- Computing reverse complement, transcription, translation, GC content, or molecular weight of a `Seq`.
- Fetching sequences or GenBank records from NCBI via `Bio.Entrez` (esearch/efetch/elink).
- Extracting a CDS from a `SeqRecord`'s `features` and translating it to protein.
- Running pairwise global/local alignment (Needleman-Wunsch / Smith-Waterman) with `PairwiseAligner`.

## Version Compatibility

Biopython >= 1.79 (current stable ~1.84), Python >= 3.9. `Bio.pairwise2` is deprecated since 1.80 — use `Bio.Align.PairwiseAligner` instead. In Biopython >= 1.79, `Seq` and `str` interoperate freely in comparisons/concatenation; older versions required explicit `str()` conversion.

## Prerequisites

```bash
pip install biopython
```
Concepts: FASTA/FASTQ/GenBank file formats, genetic code tables, Phred quality scores.

**Goal:** Manipulate a DNA sequence and translate it to protein.
**Approach:** Build a `Seq`, use its built-in methods for transcription/translation, and always call `to_stop=True` for CDS translation so the stop codon isn't included as `*`.

```python
from Bio.Seq import Seq, MutableSeq
from Bio.SeqUtils import gc_fraction, molecular_weight

def translate_cds(dna_str: str, table: int = 1) -> str:
    """Translate a coding DNA sequence to protein, stopping at the first stop codon.

    table=1 standard, table=2 vertebrate mitochondrial (TGA=Trp), table=11 bacterial.
    """
    dna = Seq(dna_str)
    return str(dna.translate(table=table, to_stop=True))

dna = Seq("ATGCGATCGATCGTAA")
dna.complement()            # 3'->5' complement
dna.reverse_complement()    # 5'->3' reverse complement
dna.transcribe()            # DNA -> RNA (T->U)
dna.transcribe().back_transcribe()  # RNA -> DNA round trip

gc_fraction(dna)             # 0.0-1.0
molecular_weight(dna)        # Da (DNA)
molecular_weight(Seq("MKPG"), seq_type="protein")  # Da (protein)

mutable = MutableSeq("ATGCGATCG")
mutable[3] = "T"              # in-place point mutation, G->T
```

**Goal:** Annotate a sequence with features and extract/translate a CDS.
**Approach:** `SeqRecord` holds the sequence plus metadata; `SeqFeature.location.extract()` slices out the feature's subsequence directly from the parent record.

```python
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation

record = SeqRecord(
    Seq("ATGCGATCGATCGATCGATCGATCGTAA"),
    id="BRCA1_001", name="BRCA1",
    description="BRCA1 partial CDS",
)
record.annotations["organism"] = "Homo sapiens"

cds = SeqFeature(FeatureLocation(0, 27), type="CDS",
                 qualifiers={"gene": ["BRCA1"]})
record.features.append(cds)

cds_seq = cds.location.extract(record.seq)
protein = cds_seq.translate(to_stop=True)

# Per-letter annotations (e.g. FASTQ quality)
record.letter_annotations["phred_quality"] = [30, 30, 28, 35]
```

**Goal:** Read, filter, and write sequence files in bulk.
**Approach:** `SeqIO.parse()` always returns an iterator (safe for multi-record files); `SeqIO.read()` requires exactly one record. Filter FASTQ reads by mean Phred quality before writing back out.

```python
from Bio import SeqIO

def filter_fastq_by_quality(in_path: str, out_path: str, min_mean_q: float = 25.0) -> int:
    """Keep only reads whose mean Phred quality >= min_mean_q; return count kept."""
    good = [r for r in SeqIO.parse(in_path, "fastq")
            if sum(r.letter_annotations["phred_quality"]) / len(r) >= min_mean_q]
    return SeqIO.write(good, out_path, "fastq")

# Multi-record read
for rec in SeqIO.parse("sequences.fasta", "fasta"):
    print(rec.id, len(rec))

# Single-record read (raises if 0 or >1 records)
record = SeqIO.read("single.gb", "genbank")

records_dict = SeqIO.to_dict(SeqIO.parse("seqs.fasta", "fasta"))
SeqIO.write(list(records_dict.values()), "output.fasta", "fasta")
SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")  # quality info is lost
```

Supported formats: `fasta`, `fastq`, `genbank` (or `gb`), `embl`, `stockholm`, `clustal`, `phylip`.

**Goal:** Fetch a gene's mRNA from NCBI and translate its annotated CDS.
**Approach:** `esearch` for the ID, `efetch` for the GenBank record (has `features`), extract the CDS feature, translate, and compute basic stats — always close handles and set `Entrez.email`.

```python
from Bio import Entrez

def gene_to_protein(gene_name: str, organism: str = "Homo sapiens", email: str = "you@example.com"):
    """Fetch a gene's RefSeq mRNA from NCBI, extract its CDS, and translate to protein."""
    Entrez.email = email  # required by NCBI
    handle = Entrez.esearch(
        db="nucleotide",
        term=f"{gene_name}[Gene] AND {organism}[Organism] AND RefSeq[Filter] AND mRNA[Filter]",
        retmax=1,
    )
    ids = Entrez.read(handle)["IdList"]
    handle.close()
    if not ids:
        return None

    handle = Entrez.efetch(db="nucleotide", id=ids[0], rettype="gb", retmode="text")
    gb_record = SeqIO.read(handle, "genbank")
    handle.close()

    cds_seq = next(
        (f.location.extract(gb_record.seq) for f in gb_record.features if f.type == "CDS"),
        None,
    )
    if cds_seq is None:
        return None
    protein = cds_seq.translate(to_stop=True)
    return {"accession": gb_record.id, "mrna_bp": len(gb_record),
            "cds_bp": len(cds_seq), "protein_aa": len(protein)}

# Cross-database links (nucleotide -> protein)
handle = Entrez.elink(dbfrom="nucleotide", db="protein", id="NM_000518.5")
link_record = Entrez.read(handle); handle.close()
```

**Goal:** Align two sequences globally or locally.
**Approach:** Configure `PairwiseAligner` with a substitution matrix and gap scores once, then reuse it; switch `mode` between calls.

```python
from Bio.Align import PairwiseAligner, substitution_matrices

aligner = PairwiseAligner()
aligner.substitution_matrix = substitution_matrices.load("BLOSUM62")
aligner.open_gap_score = -11
aligner.extend_gap_score = -1

aligner.mode = "global"  # Needleman-Wunsch
alns = aligner.align(Seq("MVHLTPEEKSAVTALWGKVN"), Seq("MVHLTDAEKAAVNGLWGKVN"))
print(alns[0].score, alns[0])

aligner.mode = "local"   # Smith-Waterman
alns = aligner.align(Seq("XXXXMVHLTPEEKXXXXXX"), Seq("YYYMVHLTDAEKYYYY"))
print(alns[0].score, alns[0])
```

## Pitfalls

- **`translate()` stop codons**: default `dna.translate()` renders stops as `*` and continues past them; use `to_stop=True` for CDS protein extraction so trailing junk/`*` isn't included.
- **`SeqIO.parse()` vs `SeqIO.read()`**: `parse()` is an iterator, safe for any file; `read()` raises `ValueError` unless the file has exactly one record.
- **Entrez etiquette**: always call `handle.close()`; NCBI allows 3 req/s without a key, 10 req/s with `Entrez.api_key` set. Missing `Entrez.email` will get requests blocked/rate-limited harder.
- **Genetic code table**: mitochondrial code (`table=2`) reads TGA as Trp, not stop; bacterial (`table=11`) differs subtly from standard (`table=1`). Always pass `table=` explicitly for non-standard organisms.
- **`Bio.pairwise2` is deprecated** (removed in future releases) — use `Bio.Align.PairwiseAligner`, which is faster and vectorized.
- **`FeatureLocation.extract()` handles strand/joins**: for features on the minus strand or spanning multiple exons (`CompoundLocation`), `.extract()` automatically reverse-complements/concatenates correctly — don't manually slice `record.seq`.

## See Also

- `bio-database-access-entrez-search`, `bio-database-access-entrez-fetch` — deeper Entrez query patterns
- `bio-sequence-io-read-sequences`, `bio-sequence-io-write-sequences` — dedicated SeqIO I/O patterns
- `bio-sequence-manipulation-transcription-translation` — translation edge cases
- `bio-alignment-pairwise-alignment` — advanced PairwiseAligner scoring and multi-alignment workflows