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

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 *.

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.

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.

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.

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.

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