DABDGE

TP2.1 Sequence handling with Biopython

Thursday 09:15-10:30, roughly 1 h 15 min. No submission. The habits you build here are used in every block that follows.

Yesterday we said that FASTA has nowhere to put metadata, so people put it in the definition line, and that every downstream script therefore has to parse that line. This morning you write those scripts.

Set up your working directory

Everything you do today happens in one place: inside the downloaded repository, with the environment activated. Consult the repository’s root README.md for how to do that.

When all is set up, check that Biopython is available:

python -c "import Bio; print(Bio.__version__)"

If that fails, create and activate the conda environment described in the repository README.

1. Core: Get some sequences

We will use a FASTA file as BOLD used to produce. The reason why we use this older version is that it demonstrates how definition lines in FASTA files are sometimes overloaded with metadata in idiosyncratic ways.

cp data/Danaus.v3.fas data/Danaus.fas

BOLD v3 was retired in July 2026 and no longer serves this FASTA endpoint. The file in this repository is an archived teaching artefact: keep the retrieval note with it, because an API that moves is a reproducibility problem and you have just met one.

Look at what you have before you write any code:

grep '>' data/Danaus.fas

You should see a few hundred records, with BOLD v3-style definition lines.

2. Core: Read the definition line with shell tools

The definition line is structured as >ID|Scientific binomial|marker|accession. That structure is a convention, not a standard, but it is enough for a first pass:

grep '>' data/Danaus.fas | cut -f 3 -d '|' | sort | uniq -c | sort -rn

You should see about twenty of marker groups, one group being the dominant marker in BOLD

How many markers are in the file? Did you expect more than one? Could anything have told you in advance?

3. Core: Hit the wall

Now try to pull out only the COI-5P sequences with the same tools.

Write the command. Then explain why it does not work.

Intentional failure. This command cannot work with a line-oriented tool. The failure is the answer. It looks like headers without their full sequences, or sequences detached from the marker you filtered on. A FASTA record spans an unpredictable number of lines, so you must use a parser.

4. Core: Filter with Biopython

Copy the code below, paste it into a text editor, and save it as scripts/filter_marker.py. (You will find that the script already exists as a fallback in the repository - we do this simply to get comfortable with saving plain text code and running it.)

import sys
from Bio import SeqIO

infile, marker = sys.argv[1], sys.argv[2]

for record in SeqIO.parse(infile, "fasta"):
    fields = record.description.split('|')
    if len(fields) > 2 and fields[2] == marker:
        print('>' + record.description)
        print(record.seq)

Run it:

python3 scripts/filter_marker.py data/Danaus.fas COI-5P > out/Danaus.COI-5P.fas
grep -c '>' out/Danaus.COI-5P.fas

You should see the number 379, i.e. the count of filtered records.

5. Core: Describe what you have

Create scripts/summarise.py:

import sys
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

lengths = []
for record in SeqIO.parse(sys.argv[1], "fasta"):
    lengths.append(len(record.seq))
    print(f"{record.id}\t{len(record.seq)}\t{gc_fraction(record.seq):.3f}")

print(f"\n# {len(lengths)} records, "
      f"min {min(lengths)}, max {max(lengths)}, "
      f"mean {sum(lengths) / len(lengths):.1f}", file=sys.stderr)

And run it:

python3 scripts/summarise.py out/Danaus.COI-5P.fas > out/Danaus.COI-5P.tsv
more out/Danaus.COI-5P.tsv # or open in something that views TSV tables, e.g. RStudio or excel

You should see the 379 records, with their ID, sequence lengths and GC content.

Why is the ID formatted the way it is? Could this be improved? How?

6. Stretch: Translate, and find out what the reading frame is

Create scripts/tp1_translate_frames.py

import sys
from Bio import SeqIO
from Bio.Data.CodonTable import TranslationError

print("id\tframe\tstops")
for record in SeqIO.parse(sys.argv[1], "fasta"):
    counts = []
    for frame in range(3):
        seq = record.seq[frame:]
        stops = 0
        for i in range(0, len(seq) - len(seq) % 3, 3):
            try:
                if seq[i:i + 3].translate(table=5) == "*":
                    stops += 1
            except TranslationError:
                pass
        counts.append((stops, frame))
    stops, frame = min(counts)
    print(f"{record.id}\t{frame}\t{stops}")

And run it:

python3 scripts/tp1_translate_frames.py out/Danaus.COI-5P.fas

Notice how the optimal frame varies a bit from one record to the next but that the optimum always has 0 stop codons. Now change table=5 to the default table=1. What happens, and why?

Intentional failure. Changing to table=1 is expected to fail biologically for this case. The failure is the answer. It looks like many extra stop codons in otherwise plausible coding sequence.

What to keep

Leave scripts/ and out/ where they are. You will retrieve records into the same tree this morning, annotate them this afternoon, and package the lot at the end of the day.

Next: TP2.2 Retrieval

This site is open source. Improve this page.