Python for Biologists

Lesson 13 of 14 · 11 min

Mini Project: Sequence Analysis Report

This mini project ties together everything from the module into one script you can run. You will read DNA sequences from a file, and for each one measure its GC content, build its reverse complement, and translate it into protein. The result is a short, tidy report printed to your screen.

Create the FASTA file

Biological sequences are usually stored in FASTA format, a plain text file where each entry starts with a header line beginning with a greater-than sign, followed by the sequence letters on the lines below. Save the two sequences below in a file named sequences.fasta in the same folder where your script will live.

text
>seq1 sample gene
ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG
>seq2 short peptide
ATGTTTAAAGGGCCCTAA

Write the script

python
from Bio import SeqIO                       # tools for reading sequence files
from Bio.SeqUtils import gc_fraction        # helper that measures GC content

# SeqIO.parse reads the file one record at a time; "fasta" is the format
for record in SeqIO.parse("sequences.fasta", "fasta"):
    seq = record.seq                        # the DNA letters for this record
    gc = gc_fraction(seq) * 100             # fraction of G and C, turned into a percent
    rev_comp = seq.reverse_complement()     # the opposite strand, read 5' to 3'
    protein = seq.translate(to_stop=True)   # DNA codons into amino acids, stop at first stop

    print(f"ID: {record.id}")               # the sequence name from the FASTA header
    print(f"Length: {len(seq)} bp")         # how many bases long
    print(f"GC content: {gc:.1f}%")         # one decimal place
    print(f"Reverse complement: {rev_comp}")
    print(f"Protein: {protein}")
    print("-" * 40)                         # a divider line of 40 dashes

Save this file as analyze.py. It imports two Biopython tools, then SeqIO.parse walks through the FASTA file one record at a time. For each record it stores the DNA in seq, then computes the GC content, the reverse complement (the opposite strand read 5' to 3'), and the protein by reading the DNA three letters at a time. Each group of three bases is called a codon, and to_stop tells translate to stop at the first stop codon, which is why the protein ends with a real amino acid and no asterisk.

Tip: gc_fraction returns a value between 0 and 1, not a percentage, so we multiply by 100. GC content matters because sequences rich in G and C bind more tightly and melt at higher temperatures, which affects experiments like PCR.

Run it

bash
pip install biopython
python analyze.py
ID: seq1
Length: 39 bp
GC content: 56.4%
Reverse complement: CTATCGGGCACCCTTTCAGCGGCCCATTACAATGGCCAT
Protein: MAIVMGR
----------------------------------------
ID: seq2
Length: 18 bp
GC content: 38.9%
Reverse complement: TTAGGGCCCTTTAAACAT
Protein: MFKGP
----------------------------------------

Try it yourself: Add a third DNA sequence to sequences.fasta and run the script again. Then add one more line inside the loop that prints the protein length, print(f"Protein length: {len(protein)} aa"), and see which sequence codes for the longest protein.