Biopython Tutorial: Read, Write, and Analyze Sequences in Python

Quick answer: Biopython is the standard Python library for biological data. It parses FASTA, FASTQ, GenBank and dozens of other formats through one API (SeqIO), represents sequences as Seq objects that complement, transcribe and translate correctly, scores alignments with PairwiseAligner, and reaches NCBI’s databases through Entrez. Install it with pip install biopython and the day-one tasks — read a file, filter records, convert a format, fetch a sequence by accession — collapse into three or four lines of readable code instead of hand-written string parsing that breaks on the first wrapped FASTA header.

The code companion to Python for bioinformatics and the FASTA file explainer on this blog: those teach the language and the format, this one teaches the library that reads them.

What Biopython actually gives you

Biopython is a library, not a framework or a pipeline tool. You import the pieces you need, and nothing runs in the background. The modules you will touch in the first month:

• Bio.SeqIO — the workhorse. One parser interface for FASTA, FASTQ, GenBank, FASTQ-solexa, and roughly eighty other formats, plus the matching writers for converting between them.
• Bio.Seq — the Seq object: a sequence with methods for complement(), reverse_complement(), transcribe(), translate() and slicing, so you stop treating DNA as a bare string and re-checking your codon table by hand.
• Bio.Align — PairwiseAligner, the maintained replacement for the old pairwise2 module, with local and global modes and configurable gap scores.
• Bio.Entrez — programmatic access to NCBI’s E-utilities: search nucleotide/protein/pubmed databases, fetch records, page through results.
• Bio.Phylo — reading and writing Newick trees, so a tree file becomes a traversable object instead of nested parentheses.
• Bio.PDB — parsing PDB structure files: residues, atoms, and coordinates for structural work.

Everything else in the package (restriction enzymes, motifs, CodonTables, SeqUtils) is there when you need it. The point is that these modules agree on the same record objects, so a file read with SeqIO can flow straight into an alignment, a translation, or a tree without conversion glue.

Install it in your project environment:

python -m pip install biopython
python -c "import Bio; print(Bio.__version__)"

Biopython releases are frequent and the API has moved in the last few years — pairwise2 is gone from recent versions, and this tutorial uses the current Align API throughout.

Reading files with SeqIO

Every format shares one entry point: SeqIO.parse(handle, format) yields records, and each record carries an id, a description, a Seq and — for FASTQ — its per-base quality scores.

from Bio import SeqIO

for rec in SeqIO.parse("sample.fastq", "fastq"):
    quals = rec.letter_annotations["phred_quality"]
    if len(rec) >= 150 and sum(quals) / len(quals) >= 35:
        print(rec.id, len(rec), "good")

Three habits worth forming immediately:

• parse is a generator. It streams the file, so a 20 GB FASTQ costs the memory of one record at a time. Only call list(...) when you genuinely need everything at once.
• SeqIO.read is for single-record files (one fasta with one sequence). Multi-record files must use parse or you get a ValueError.
• Qualities live in letter_annotations. They are integers, Phred-scaled, aligned base-for-base with the sequence — the same numbers the sequencer wrote into the file, which is why FASTQ exists as a format at all (see the FASTA/FASTQ comparison for why one carries them and the other does not).

Converting formats is one line, because reading and writing share the same record type:

records = SeqIO.parse("sample.fastq", "fastq")
n = SeqIO.write(records, "sample.fasta", "fasta")
print("wrote", n, "records")

The quality scores do not survive the conversion — FASTA has nowhere to put them — which is exactly the kind of silent information loss you want to decide on deliberately rather than discover later.

Seq objects instead of raw strings

A Seq object behaves like a string but carries operations that are correct by construction:

from Bio.Seq import Seq

cds = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
print(cds.translate())              # MAIVMGR*KGAR*
print(cds.reverse_complement())     # correct strand flip
print(cds.translate(table=2, to_stop=True))   # alternative genetic code, stop at *

This matters more than it sounds. Reverse complementing a string by hand is the classic off-by-one-and-wrong-strand bug; translate() silently gives you a stop codon in the middle of a protein when your reading frame is wrong, which is exactly the diagnostic you wanted. Slicing works as it does on strings (cds[10:40]), so extracting a region from a GenBank feature is ordinary Python:

from Bio import SeqIO

gb = SeqIO.read("entry.gb", "genbank")
for feature in gb.features:
    if feature.type == "CDS":
        gene = feature.qualifiers.get("gene", ["?"])[0]
        print(gene, feature.location, len(feature.location.extract(gb.seq)))

feature.location.extract() is the detail beginners miss: features can sit on the minus strand or span a join, and slicing the record directly does not reverse-complement for you.

Pairwise alignment the modern way

Bio.Align.PairwiseAligner replaces the retired pairwise2 module and is the alignment engine inside Biopython — for comparing two sequences, not for mapping reads to a genome.

from Bio.Align import PairwiseAligner

aligner = PairwiseAligner()
aligner.mode = "local"
aligner.match_score = 2
aligner.mismatch_score = -1
aligner.open_gap_score = -5
aligner.extend_gap_score = -2

seq1 = "GATCACAGGTCTATCACCCT"
seq2 = "GATCACAGGTCTATCACCCTATC"
print(aligner.score(seq1, seq2))
aln = aligner.align(seq1, seq2)
print(aln[0])

aligner.score() gives you the number without constructing alignments; aligner.align() returns a lazy container whose first element renders as a formatted alignment. Setting mode = "global" forces end-to-end alignment (Needleman-Wunsch semantics), "local" finds the best matching region (Smith-Waterman) — the same distinction covered from first principles in the alignment algorithms guide. What Biopython gives you is the scoring machinery with sane defaults; what it does not give you is a genome-scale index, so mapping millions of reads to a reference stays the job of BWA or Bowtie2, never this.

Querying NCBI with Entrez

Bio.Entrez wraps NCBI’s E-utilities so a record by accession is two lines:

from Bio import Entrez, SeqIO

Entrez.email = "you@example.com"     # required by NCBI; identify yourself
with Entrez.efetch(db="nuccore", id="MN905015", rettype="fasta", retmode="text") as h:
    record = SeqIO.read(h, "fasta")

print(record.id, len(record))

The rules that keep you off NCBI’s block list: always set Entrez.email, stay under three requests per second for anonymous use, use esearch to resolve names into accession lists before efetch, and never pull a thousand records in a loop when retmax paging exists. For anything bigger than a few hundred records — entire projects, raw runs — drop down to proper download tooling, which is a different workflow entirely and covered step by step in the public sequencing data guide.

Writing results and chaining steps

Scripts earn their keep when one step’s output is the next step’s input, and SeqIO is the interchange layer:

from Bio import SeqIO

clean = (r for r in SeqIO.parse("raw.fastq", "fastq")
         if min(r.letter_annotations["phred_quality"]) >= 20)
SeqIO.write(clean, "clean.fastq", "fastq")

sizes = [len(r) for r in SeqIO.parse("clean.fastq", "fastq")]
print("reads:", len(sizes), "mean len:", sum(sizes) / len(sizes))

A generator expression filters a FASTQ of any size with constant memory; SeqIO.write returns the count of records written, which belongs in your log. For pipelines that must survive re-runs and failures — everything a real analysis needs — the file-chaining pattern graduates into a workflow engine such as Nextflow or Snakemake, as described in the workflows article.

Where Biopython stops

Being honest about the boundary saves months: Biopython is the format-and-sequences layer, not the analysis layer.

• It reads and writes BAM/VCF-adjacent formats’ simpler cousins and GenBank/FASTQ well, but it is not a variant caller, an aligner for genomes, or a differential-expression engine.
• It has no statistical models — no multiple-testing correction, no negative binomial fitting. That toolkit lives in R/Bioconductor, not here.
• For raw sequencing files at scale, samtools, bcftools and fastp will always be faster than Python equivalents; Biopython is for glue, inspection, and small transformations.
• For BLAST, Bio.Blast.NCBIWWW.qblast exists but is slow and rate-limited against the public server — local BLAST+ with the Bio.Blast parser is the usable combination.

The division of labour that works: Biopython reads, converts and inspects; the specialist tools analyze; your workflow engine sequences them.

Key takeaways

• SeqIO.parse is one API for nearly every format, and it streams — read is only for single-record files.
• Seq objects make complement, translation and slicing correct by construction; hand-rolled string surgery is where strand and frame bugs come from.
• PairwiseAligner (local or global) is the current alignment API; pairwise2 is retired.
• Entrez gives scripted NCBI access — set your email, respect the rate limit, and hand bulk downloads to the SRA/ENA tooling.
• Conversion between formats is one line, but information (quality scores, annotations) only survives if the target format has somewhere to put it.
• Biopython is the glue layer: pair it with samtools/BLAST+ for heavy lifting and a workflow engine to make the chain reproducible.

Challenge

Take the GenBank record for any gene you have worked with, load it with SeqIO, and print every CDS with its strand, its extracted length, and the protein its translate() produces — then compare that translation to the /translation qualifier the submitters stored. When the two disagree, work out why (frame, pseudo-gene, alternative table, or a genuinely different annotation). That discrepancy hunt is exactly the debugging loop of annotation work, and it takes ten lines of Biopython.

Want guided practice turning library scripts into defensible bioinformatics analysis? Ampersand Academy offers one-to-one mentoring in bioinformatics — part of Mahadhi, a development and digital marketing company.

FAQ