From FASTA to a gene and its protein
Trace real HBB and HBD transcripts through their annotated coding regions, then verify translation against RefSeq.
Goal and how to run
This lesson assumes basic Python knowledge. Read each method and caption before changing parameters. AI exercises are optional; no AI account is needed for the core workflow. Data lives in notebook metadata, so there is no companion data download. Keep the notebook saved in the active working folder.
Your goal is to understand the biological question, execute transparent calculations, check the results, and state the limits of the conclusion.
Setup and parameters
from pathlib import Path
import base64, gzip, hashlib, json, platform
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display, HTML
RECIPE_ID='01-fasta-to-protein'
OUTPUT_DIR=Path.cwd()/"results"/RECIPE_ID
plt.rcParams.update({"figure.dpi":120,"font.size":10,"axes.spines.top":False,"axes.spines.right":False,
"axes.prop_cycle":plt.cycler(color=["#18785f","#416aa6","#bc6541","#9974a5"])})
print("Python:",platform.python_version(),"; NumPy:",np.__version__,"; Pandas:",pd.__version__)
print("Working folder:",Path.cwd())
print("Core data mode: embedded snapshot; no network requests")
Python: 3.13.13 ; NumPy: 2.5.3 ; Pandas: 3.0.6 Working folder: /Users/zhiluo/Developer/vibeit-cookbook/cookbook/downloads Core data mode: embedded snapshot; no network requests
Load the embedded snapshot
def load_snapshot(recipe_id, expected_sha256, snapshot_key):
"""Read embedded data from this notebook, including a renamed copy."""
for path in sorted(Path.cwd().glob("*.ipynb")):
try:
notebook = json.loads(path.read_text(encoding="utf-8"))
metadata = notebook.get("metadata", {}).get("vibeit_cookbook", {})
blob = metadata.get("snapshots", {}).get(snapshot_key)
if metadata.get("recipe_id") != recipe_id or not blob:
continue
if blob.get("sha256") != expected_sha256:
continue
raw = gzip.decompress(base64.b64decode(blob["gzip_base64"], validate=True))
if hashlib.sha256(raw).hexdigest() != expected_sha256:
raise ValueError("Embedded snapshot checksum mismatch: " + snapshot_key)
return json.loads(raw)
except (OSError, json.JSONDecodeError):
continue
raise FileNotFoundError(
"Save/open the complete .ipynb in the active working folder. "
"The embedded snapshot could not be located for " + recipe_id)
SNAPSHOT_HASHES={'sequences': 'a8ae02a43c841e2932efafaf82a0bf78199baa8921a821fcbd774c8d7649e8b7'}
def snapshot(key):
return load_snapshot(RECIPE_ID,SNAPSHOT_HASHES[key],key)
print("Embedded snapshots:",list(SNAPSHOT_HASHES))
Embedded snapshots: ['sequences']
Read the method functions
def parse_fasta(text):
records, label, chunks = [], None, []
for line in text.splitlines():
line = line.strip()
if not line:
continue
if line.startswith(">"):
if label is not None:
records.append((label, "".join(chunks).upper()))
label, chunks = line[1:].strip(), []
if not label:
raise ValueError("FASTA header is empty")
else:
if label is None:
raise ValueError("FASTA requires a header before sequence lines")
chunks.append(line)
if label is not None:
records.append((label, "".join(chunks).upper()))
if not records:
raise ValueError("No FASTA records")
if len({label.split()[0] for label, seq in records}) != len(records):
raise ValueError("Duplicate FASTA identifiers")
for label, seq in records:
if not seq or set(seq) - set("ACGTRYSWKMBDHVN"):
raise ValueError("Empty or invalid DNA sequence: " + label)
return records
print("Method ready: parse_fasta")
Method ready: parse_fasta
def translate_cds(sequence):
sequence = sequence.upper()
if len(sequence) % 3:
raise ValueError("Complete CDS length must be divisible by three")
bases = "TCAG"
amino_acids = "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"
codons = [a+b+c for a in bases for b in bases for c in bases]
table = dict(zip(codons, amino_acids))
protein = "".join(table.get(sequence[i:i+3], "X") for i in range(0,len(sequence),3))
if "*" in protein[:-1]:
raise ValueError("Internal stop codon: check the reading frame or genetic code")
return protein.removesuffix("*")
print("Method ready: translate_cds")
Method ready: translate_cds
def gc_fraction(sequence):
sequence = sequence.upper()
denominator = sum(sequence.count(base) for base in "ACGT")
if denominator == 0:
return float("nan")
return (sequence.count("G") + sequence.count("C")) / denominator
print("Method ready: gc_fraction")
Method ready: gc_fraction
Steps
1. Read the record before reading the bases
HBB and HBD are human paralogs, not two patients. RefSeq NM_000518.5 and NM_000519.4 describe mRNA. Their sequences include untranslated regions (UTRs). The annotated CDS is 51..494 in both records. GenBank uses 1-based inclusive coordinates; Python slicing is 0-based, end-exclusive. Therefore the correct slice is [start-1:end]. Keep the accession version: a sequence without its version is harder to reproduce.
records = snapshot("sequences")["records"]
fasta = "\n".join(">"+r["accession"]+" "+r["gene"]+"\n"+r["sequence"] for r in records)
parsed = parse_fasta(fasta)
display(pd.DataFrame([dict(gene=r["gene"], accession=r["accession"],
mrna_bases=len(r["sequence"]), cds=f'{r["cds_start"]}..{r["cds_end"]}',
protein=r["protein_id"]) for r in records]))
print(fasta[:180] + "…")
| gene | accession | mrna_bases | cds | protein | |
|---|---|---|---|---|---|
| 0 | HBB | NM_000518.5 | 628 | 51..494 | NP_000509.1 |
| 1 | HBD | NM_000519.4 | 620 | 51..494 | NP_000510.1 |
>NM_000518.5 HBB ACATTTGCTTCTGACACAACTGTGTTCACTAGCAACCTCAAACAGACACCATGGTGCATCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATGAAGTTGGTGGTGAGGCCCTGGGCAGGCTGCTGGTGGTCTACCCTTG…
2. Draw the transcript map
A transcript diagram checks our coordinates before translation. Each rectangle spans annotated bases; the green region is the CDS and the gray regions are UTRs. The x-axis is the mRNA coordinate, not the chromosome coordinate. This drawing comes from the downloaded record, not a hand-drawn gene model.
fig, ax = plt.subplots(figsize=(8,3))
for i, r in enumerate(records):
start, end, length = r["cds_start"], r["cds_end"], len(r["sequence"])
ax.broken_barh([(1,start-1),(end+1,length-end)],(i-.22,.44),facecolors="#d9dedb")
ax.broken_barh([(start,end-start+1)],(i-.22,.44),facecolors="#18785f")
ax.text(start,i+.3,f"CDS {start}–{end}",fontsize=10)
ax.set(yticks=range(len(records)),yticklabels=[r["gene"] for r in records],
xlabel="mRNA position (1-based)",title="Annotated coding regions and untranslated regions")
ax.spines[["top","right","left"]].set_visible(False)
plt.tight_layout(); plt.show()
3. Inspect composition and ambiguous bases
GC content is a descriptive sequence statistic. Count only A/C/G/T in its denominator and report ambiguous bases separately. A GC difference between these two transcripts does not establish expression, stability or evolutionary selection. Change the window later to see how smoothing changes the display.
from collections import Counter
composition = pd.DataFrame({r["gene"]:Counter(r["sequence"]) for r in records}).reindex(list("ACGT")).fillna(0)
composition.plot.bar(figsize=(7,3),color=["#18785f","#416aa6"],rot=0)
plt.ylabel("Base count"); plt.title("Composition of the complete mRNA sequences")
plt.tight_layout(); plt.show()
display(pd.DataFrame([dict(gene=r["gene"],gc_fraction=gc_fraction(r["sequence"]),
ambiguous_bases=sum(b not in "ACGT" for b in r["sequence"])) for r in records]))
| gene | gc_fraction | ambiguous_bases | |
|---|---|---|---|
| 0 | HBB | 0.512739 | 0 |
| 1 | HBD | 0.506452 | 0 |
4. Compare local GC content
Each point uses an 80-base window and a 20-base step. The position is the window midpoint; terminal partial windows are omitted. A short transcript supports only a few windows, so avoid treating this curve as a genome-wide pattern.
WINDOW, STEP = 80, 20
fig, ax = plt.subplots(figsize=(8,3))
for r in records:
starts = range(0,len(r["sequence"])-WINDOW+1,STEP)
ax.plot([s+(WINDOW+1)/2 for s in starts],
[100*gc_fraction(r["sequence"][s:s+WINDOW]) for s in starts],marker="o",label=r["gene"])
ax.set(xlabel="Window midpoint (mRNA coordinate)",ylabel="GC (%)",title=f"Local GC: {WINDOW}-base windows")
ax.legend(); plt.tight_layout(); plt.show()
5. Translate the annotated coding sequence
Use the standard genetic code, not a search for the first ATG in the whole transcript. A complete CDS includes its terminal stop codon; the RefSeq protein does not. Here 444 coding bases correspond to 147 amino acids plus the terminal stop. The first methionine is present in the translated precursor. Mature hemoglobin chain numbering may remove that methionine.
proteins, cds_records = {}, {}
for r in records:
cds = r["sequence"][r["cds_start"]-1:r["cds_end"]]
protein = translate_cds(cds)
assert protein == r["translation"], "Translation differs from RefSeq"
proteins[r["gene"]], cds_records[r["gene"]] = protein, cds
print(r["gene"], "CDS:",len(cds), "bases; protein:",len(protein),"aa; RefSeq match: True")
print(protein)
HBB CDS: 444 bases; protein: 147 aa; RefSeq match: True MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVANALAHKYH HBD CDS: 444 bases; protein: 147 aa; RefSeq match: True MVHLTPEEKTAVNALWGKVNVDAVGGEALGRLLVVYPWTQRFFESFGDLSSPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFSQLSELHCDKLHVDPENFRLLGNVLVCVLARNFGKEFTPQMQAAYQKVVAGVANALAHKYH
6. Check and export the result
Assertions make coordinate or translation mistakes visible. Export both nucleotide and protein sequences so the analysis can move to another tool. The output directory is created next to the notebook. Re-running replaces this lesson's generated files, not your input notebook.
assert len(parsed) == 2
assert all(len(cds) == 444 for cds in cds_records.values())
assert all(len(p) == 147 for p in proteins.values())
OUTPUT_DIR.mkdir(parents=True,exist_ok=True)
(OUTPUT_DIR/"transcripts.fasta").write_text(fasta+"\n")
(OUTPUT_DIR/"proteins.fasta").write_text("".join(f">{g}\n{p}\n" for g,p in proteins.items()))
pd.DataFrame([dict(gene=g,cds_bases=len(cds_records[g]),protein_aa=len(p),gc=gc_fraction(cds_records[g]))
for g,p in proteins.items()]).to_csv(OUTPUT_DIR/"sequence-summary.csv",index=False)
print("Checks passed. Exported:",[p.name for p in OUTPUT_DIR.iterdir()])
Checks passed. Exported: ['proteins.fasta', 'transcripts.fasta', 'sequence-summary.csv']
Exercises and reference explanation
- Explain why
[50:494]is correct for CDS 51..494, and what changes if the start index is 51. - Change WINDOW from 80 to 40. Keep the same transcript data and label the new window size.
Reference explanation. The correct CDS has 444 bases. Starting at index 51 removes one base: the 443-base length fails the complete-CDS check. A smaller window shows more local variation; it cannot establish a regulatory mechanism.
Keep a copy before changing parameters. Compare old and new figures and record which inputs, calculations and conclusions changed.
Optional AI coding exercises
Prompt 1
Add an ambiguous-base count beside each GC fraction using only the standard library and the current packages. Keep the A/C/G/T denominator explicit.
Prompt 2
Add a reverse-complement function for IUPAC DNA symbols and a double-reverse-complement assertion. Explain why reverse-complementing an annotated forward CDS is not a replacement for strand annotation.
Human acceptance. Check that the original two RefSeq translations still match; test N-only GC as NaN, invalid symbols as an error, and reverse complement twice as identity.
Review the edited source, then manually run affected cells and subsequent checks. Ask the assistant to explain the purpose, packages and data provenance of its changes, and reconcile that explanation with actual output.
Optional online extension
RUN_ONLINE_UPDATE=False
if RUN_ONLINE_UPDATE:
import requests
try:
response=requests.get('https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NM_000518.5&rettype=fasta&retmode=text',timeout=20)
response.raise_for_status()
print("Source:",response.url)
print(response.text[:700])
except requests.RequestException as error:
print("Online extension unavailable:",type(error).__name__,str(error)[:160])
else:
print("Online extension skipped; embedded core data unchanged")
Online extension skipped; embedded core data unchanged
Sources, snapshots and next steps
- RefSeq HBB/HBD — NCBI public sequence data; retrieved 2026-10-09. SHA-256
ea64e38e8e5ffad06f6d8770eaee066ef556e1488c4b9a6b340cbd445abe04ad. GenBank XML CDS uses 1-based inclusive coordinates
Next steps. Transfer these methods to your own public or authorized data while retaining input checks, recorded parameters and result validation. Document the new experimental design and method limits; this example does not validate a new experiment. Full data licenses and generation instructions accompany the cookbook.