The Biopython basics I use in almost every script — Part 01

Most of my early bioinformatics code did things that Biopython already does in one line. I just didn’t know it yet. I was splitting strings by hand, writing my own reverse-complement function (and getting it subtly wrong), counting nucleotides with loops. It worked, mostly. It was also a waste of an afternoon.

So this is the stuff I wish someone had put in front of me on day one. Nothing fancy. Just the handful of tools I end up reaching for in almost every script, with the small traps that cost me time so they don’t cost you any.

I’m on Biopython 1.87 as I write this. If you’re on something older, most of this still holds, but a few things have changed over the years, so check your version if a line misbehaves.

Install it first

pip install biopython

Then confirm it worked:

import Bio
print(Bio.__version__)

One thing that threw me at the start: the import is Bio, not biopython. If import Bio runs without complaining, you're set.

1. Seq — the thing everything else is built on

A Seq holds a single sequence. That's it. But almost every other tool hands you one of these, so it's worth getting comfortable with.

>>> from Bio.Seq import Seq
>>> my_dna = Seq("ATGCGTACGTTAG")
>>> print(my_dna)
ATGCGTACGTTAG
>>> len(my_dna)
13

It behaves a lot like a normal string, which is the point. You can measure it, slice it, search it. The difference is that it knows it’s biology, so it can do things a plain string can’t.

2. Reverse complement, and other things you’ll do constantly

The reverse complement is probably the function I call most. If you’ve ever written this yourself, you know there are two easy ways to mess it up: forgetting to reverse, or fumbling the base pairing. Biopython just does it.

>>> my_dna.reverse_complement()
Seq('CTAACGTACGCAT')

There’s a plain complement() too, but in real work you almost always want the reverse complement, because genes on the opposite strand are read the other direction.

Counting a base:

>>> my_dna.count("A")
3

Small trap here: count does not count overlaps. In AAAA, counting AA gives you 2, not 3. I've been bitten by that.

Finding where something starts:

>>> my_dna.find("TAC")
3
>>> my_dna.find("GGG")
-1

That -1 is the "not found" signal. Get used to checking for it, because it's easy to treat -1 as a real position and end up slicing from the wrong end of your sequence.

And slicing works how you’d expect, with the usual Python rule that the start is included and the end is not:

>>> my_dna[0:3]
Seq('ATG')

3. From DNA to protein in three lines

This is the part that made me actually like the library. The whole central-dogma path, DNA to RNA to protein, is three method calls.

dna = Seq("ATGGCCATTGTAATGGGCCGCTGA")

print(dna.transcribe())         # DNA -> RNA, T becomes U
print(dna.translate())          # -> protein
print(dna.translate(to_stop=True))

Output:

AUGGCCAUUGUAAUGGGCCGCUGA
MAIVMGR*
MAIVMGR

A quick read of that protein: each letter is one amino acid, and the * at the end is a stop codon, not an amino acid. Pass to_stop=True and it cuts the protein at the first stop and drops the star, which is usually what you want.

One thing worth knowing before it burns you: the genetic code isn’t universal. Mitochondria and some organisms read certain codons differently. By default translate uses the standard code, but you can pass a different table:

>>> dna.translate(table="Vertebrate Mitochondrial")
Seq('MAIVMGRW')

Same DNA, different last letter, because TGA means "stop" in the standard code but codes for tryptophan in the vertebrate mitochondrial one. If you're working on mitochondrial genes and forget this, your protein will be quietly wrong.

GC content comes from a separate module:

>>> from Bio.SeqUtils import gc_fraction
>>> round(gc_fraction(Seq("ATGCGTACGGGGCCCC")) * 100, 1)
75.0

It returns a fraction from 0 to 1, so multiply by 100 if you want a percentage.

4. SeqRecord — because a sequence without a name is useless

Here’s the thing nobody explained to me clearly, so I’ll try to.

A Seq is just letters. If you load a thousand sequences and they're all bare Seq objects, you have a thousand strings and no idea which is which. Which one is the human gene? Which is the mouse one? You can't tell. The letters don't carry a name.

A SeqRecord fixes that. It's the sequence plus its identity: an id, a name, a description, and room for extra notes. When you read a real file, this is what you get back. Not a bare Seq, a SeqRecord.

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord

record = SeqRecord(
    Seq("ATGCGTACGTTAG"),
    id="GENE001",
    name="MyGene",
    description="An example gene",
)

print(record.id, "|", record.description)
print(len(record))

GENE001 | An example gene
13

The id is the unique identifier (in real data, an accession number). The description is the human-readable line. And you can stash extra facts in annotations, which is just a dictionary:

record.annotations["organism"] = "Homo sapiens"

The sequence itself is still right there in record.seq, so anything from earlier still works: record.seq.translate(), record.seq.reverse_complement(), all of it.

5. SeqIO — reading and writing files

Real sequences live in files, and SeqIO is how you read and write them. There are two functions to keep straight:

• parse for a file with many records. It gives you an iterator, so you loop.
• read for a file with exactly one record. It hands you that record directly.

Reading a FASTA:

from Bio import SeqIO

for record in SeqIO.parse("example.fasta", "fasta"):
    print(record.id, len(record))

Note that parse needs two things: the filename and the format string. read will throw an error if the file has more than one record, so don't reach for it unless you're sure there's only one.

Now the mistake I made more than once, and watched other people make too. You can’t read a file that doesn’t exist yet. Obvious when you say it out loud, but in a notebook it’s easy to write the reading code, run it, and get a FileNotFoundError because you never actually created or downloaded the file. So, create it first:

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO

records = [
    SeqRecord(Seq("ATGAAATAG"), id="g1", description="gene one"),
    SeqRecord(Seq("ATGGGGTAA"), id="g2", description="gene two"),
]

SeqIO.write(records, "out.fasta", "fasta")

That writes a proper FASTA:

>g1 gene one
ATGAAATAG
>g2 gene two
ATGGGGTAA

FASTQ files are similar, except each read also carries a quality score for every single base (how sure the sequencer was). Biopython decodes those funny quality symbols into real numbers for you:

for record in SeqIO.parse("reads.fastq", "fastq"):
    print(record.id, record.letter_annotations["phred_quality"])

You’ll get a list like [40, 40, 40, ...], where higher means more reliable.

And converting between formats is a one-liner, which is genuinely one of my favorite small conveniences:

count = SeqIO.convert("demo.gb", "genbank", "demo.fasta", "fasta")
print("Converted", count, "record(s)")

Four arguments: input file, input format, output file, output format. Done.

6. Entrez — getting real data off NCBI

Everything above works on files you already have. But a lot of the time you want data from NCBI, the huge public database, and you don’t want to click through a website to get it. That’s Entrez.

Two rules before anything else. First, this needs an internet connection. Second, you have to give NCBI your email, every time. It’s not optional, it’s how they contact you if your script accidentally hammers their servers. So set it first:

from Bio import Entrez, SeqIO

Entrez.email = "your_real_email@example.com"

If you already know the accession, efetch pulls it straight down. Here's the original SARS-CoV-2 genome:

handle = Entrez.efetch(db="nucleotide", id="MN908947", rettype="fasta", retmode="text")
record = SeqIO.read(handle, "fasta")
handle.close()

print(record.id, len(record))

You’ll get back MN908947.3 and a length of 29903. A whole viral genome, in four lines. efetch gives you a handle, which you can think of as an open pipe to NCBI, and you read from it just like a file. Close it when you're done.

Save what you download, by the way. Don’t re-fetch the same thing over and over. Pull it once, write it to disk, read the local copy next time.

If you don’t know the accession and just have keywords, esearch finds matching ids:

handle = Entrez.esearch(db="nucleotide", term="human insulin", retmax=5)
results = Entrez.read(handle)
handle.close()

print(results["Count"], results["IdList"])

results["Count"] is how many total matches exist, and IdList is the actual ids (up to retmax). Notice you read search results with Entrez.read, not SeqIO.read. They're different things.

Now the lesson that actually cost me time. A plain keyword search like "human insulin" matches a lot of garbage, because NCBI matches words, not meaning. I once searched for human insulin and the top hit was a giant fish genome, twenty-four million bases long, that happened to mention "insulin" and "human" somewhere in its notes. So be specific, and use field tags:

term = "INS[Gene] AND Homo sapiens[Organism] AND refseq[filter] AND biomol_mrna[PROP]"

That says: the INS gene, in humans, from the clean RefSeq set, mRNA only. Way better. And whatever you get, always sanity-check the description and the length before you trust it. If you asked for a gene and got back something twenty million bases long, that’s a chromosome, not a gene.

That’s part one

If you take one thing from this, let it be the boring one: create or download your file before you try to read it. That single mistake accounts for a surprising number of the errors I see, my own included.

Next time I’ll get into the more interesting stuff: aligning sequences to see how similar they are, running BLAST to identify a mystery sequence.

Link to part 2 is as follows:

[embed]Biopython: Comparing sequences and finding out what they are — Part 02 Part 1 was all about handling one sequence at a time: read it, translate it, save it, download it. Useful, but a bit…minhazulhasansohan.medium.com

If you spot something I got wrong, or you have a Biopython trick you swear by, tell me. I’m still learning this too.