Skillversitas

Genome Analysis

Genome-scale analysis with Biopython means parsing annotated genome files (GenBank/EMBL), iterating features, and computing genome-wide statistics: gene counts, GC profiles, intergenic distances, and motif scans. The SeqRecord model makes this a matter of loops over records and features.

Parsing an Annotated Genome

from Bio import SeqIO

rec = SeqIO.read("genome.gb", "genbank")
print(rec.id, len(rec), "bp")
print("source:", rec.annotations.get("source"))

from collections import Counter
types = Counter(f.type for f in rec.features)
print(types.most_common(5))
# e.g. [('CDS', 1740), ('gene', 1740), ('mobile_element', 40), ...]

A Gene Census

cds_features = [f for f in rec.features if f.type == "CDS"]
print("CDS count:", len(cds_features))

total_len = sum(f.location.end - f.location.start for f in cds_features)
print(f"coding fraction: {total_len / len(rec.seq):.2%}")

# genes by product
products = Counter(
    f.qualifiers.get("product", ["unknown"])[0] for f in cds_features
)
for p, n in products.most_common(5):
    print(n, p[:50])

GC Content in Windows

seq = str(rec.seq)
window = 10000
gc_values = [
    (seq[i:i + window].count("G") + seq[i:i + window].count("C")) / window
    for i in range(0, len(seq) - window, window)
]
print("min GC", round(min(gc_values), 3), "max GC", round(max(gc_values), 3))

# plot with matplotlib
import matplotlib.pyplot as plt
plt.plot(range(len(gc_values)), gc_values)
plt.xlabel("window index"); plt.ylabel("GC")
plt.show()

Intergenic Gaps

cds = sorted(
    (int(f.location.start), int(f.location.end)) for f in cds_features
)
gaps = []
for (s1, e1), (s2, e2) in zip(cds, cds[1:]):
    if s2 - e1 > 1000:
        gaps.append(s2 - e1)
print("large intergenic gaps:", len(gaps), "max", max(gaps) if gaps else 0)

The Reusable Pattern

  1. Read the record once (SeqIO.read for a single genome, parse for many).
  2. Filter features by type (CDS, gene, repeat_region) into working lists.
  3. Compute per-feature or per-window statistics with plain Python + NumPy.
  4. Aggregate into tables (pandas) and plots (matplotlib).

Summary

Genome analysis in Biopython is: parse the annotated record, filter features by type, and compute statistics - gene counts and coding fraction from CDS features, GC content in sliding windows, intergenic gap sizes - then aggregate with pandas and plot with matplotlib. The SeqRecord feature model makes whole-genome summaries a few dozen lines of Python.