Going Further: GenBank, Big Files & Joining Tables

Before you start

You've finished the five core Track 1 lessons - FASTA/FASTQ, formats, Biopython, databases, and pandas. In Colab, run !pip install biopython pandas first.

Learning objectives

By the end of this lesson you will be able to: parse a GenBank record to read its annotations, handle compressed and large sequence files without running out of memory, and join two tables together to combine their data.

The core lessons got you reading sequences and tables. This optional lesson adds the three things that trip people up the moment they touch real, full-size data.

1. GenBank - sequences with annotations

A FASTA file gives you the sequence and a name. A GenBank record gives you the whole story: the organism, the gene names, where the coding regions are, and more - all the annotation a FASTA leaves out. You fetch it the same way, but ask for rettype="gb" and parse it as "genbank".

▶ Code along

Run !pip install biopython pandas first.

Open a blank Colab →

from Bio import Entrez, SeqIO
Entrez.email = "your.email@example.com"

# note rettype="gb" - that's the difference from fetching FASTA
h = Entrez.efetch(db="nucleotide", id="NM_000207", rettype="gb", retmode="text")
rec = SeqIO.read(h, "genbank")
h.close()

print("Name:", rec.description)
print("Length:", len(rec.seq))

The real value is in rec.features - a list of annotated regions. Each feature has a .type (like "gene" or "CDS") and a dictionary of .qualifiers.

for feat in rec.features:
    if feat.type in ("source", "gene", "CDS"):
        print(feat.type, "→", dict(feat.qualifiers))

Decode the jargon: CDS

A CDS (coding sequence) feature marks the exact part of the record that gets translated into protein. Pull its /translation qualifier and you have the official protein - no guessing where the gene starts. Qualifier values come back as lists, so you read them like feat.qualifiers["gene"][0].

2. Compressed & large files

Real sequencing files are huge and almost always gzipped - you'll see names ending in .fastq.gz. You don't need to unzip them; Python's gzip module reads them on the fly.

import gzip
from Bio import SeqIO

# make a gzipped FASTA, then read it WITHOUT unzipping
with gzip.open("demo.fasta.gz", "wt") as f:
    f.write(">geneA\nATGCGTACGT\n>geneB\nATGGGGCCCC\n")

# "rt" = read text; SeqIO reads straight from the gzip handle
for rec in SeqIO.parse(gzip.open("demo.fasta.gz", "rt"), "fasta"):
    print(rec.id, len(rec.seq))

The golden rule for big files: don't load it all

A FASTQ file can be tens of gigabytes - far more than your computer's memory. The trick is to stream: loop over records one at a time with SeqIO.parse (which yields them lazily) instead of building a giant list. The loop above already does this - it never holds more than one record in memory, so it works on a 10-line file or a 100-million-read file unchanged.

# count millions of reads using almost no memory
total = 0
for rec in SeqIO.parse(gzip.open("demo.fasta.gz", "rt"), "fasta"):
    total += 1
print("reads:", total)

3. Joining two tables

This is the pandas skill that unlocks real analysis. Your measurements live in one table; the information about each sample (which is control, which is treated) lives in another. To analyze them together, you merge on a shared column.

import pandas as pd

counts = pd.DataFrame({
    "sample": ["s1", "s2", "s3"],
    "gene_count": [1200, 980, 1500],
})
samples = pd.DataFrame({
    "sample": ["s1", "s2", "s3"],
    "condition": ["control", "control", "treated"],
})

# glue them together on the shared "sample" column
merged = counts.merge(samples, on="sample")
print(merged)

# now you can ask cross-table questions
merged.groupby("condition")["gene_count"].mean()

Why merging is everywhere

Nearly every real analysis is "table A + table B." An RNA-seq count matrix means nothing until you join it to the sample sheet that says which columns are treated vs. control - that join is literally step one of the flagship project. Learn merge and you've learned the connective tissue of data analysis.

🚀 Put it together

Fetch a gene as GenBank and print its organism and every gene name in the record.
Extract the /translation from the CDS feature and compare it to rec.seq.translate().
Download a real .fasta.gz (e.g., a small genome) and count its sequences by streaming - never building a list.
Make a fake counts table and sample sheet, merge them, and compute the mean per condition.

Check your understanding

What does GenBank give you that FASTA doesn't?

Annotation - the organism, gene names, coding regions (CDS), and more. FASTA is just a header plus the sequence.

How do you read a 20 GB FASTQ without running out of memory?

Stream it: loop over records one at a time with SeqIO.parse (and gzip.open if it's compressed), never building a full list. Each iteration holds just one record.

You have a counts table and a separate sample sheet. How do you analyze them together?

Merge them on their shared column (e.g., the sample ID) with pandas.merge, then group and summarize across the combined table.

Why can you loop over a 100-million-read FASTQ file without running out of memory?

SeqIO.parse yields records one at a time rather than building a giant list, so the loop never holds more than a single record in memory.

In a GenBank record, what does a CDS feature mark?

A CDS (coding sequence) feature marks the region that becomes protein. Its /translation qualifier gives the official protein with no guessing where the gene starts.

That's Track 1, in full