Analyzing DNA Sequences with Python

After Completing This Topic

You'll be able to calculate GC content of DNA sequences in Python, parse FASTA files, and convert codons to amino acids.

Variables: Containers for Data

Just as you write "sample concentration: 2.5 mg/mL" in your lab notebook, Python stores values in variables.

python

gene_name = "BRCA1"
sequence = "ATGCGATCGATCGATCG"
gc_ratio = 0.529

print(f"Gene: {gene_name}")
print(f"Sequence length: {len(sequence)}bp")
print(f"GC ratio: {gc_ratio:.1%}")

assert gene_name == "BRCA1"
assert len(sequence) == 17

Strings: Working with DNA Sequences

DNA sequences are strings made of four characters: A, T, G, C. Python's string methods let you analyze sequences.

python

sequence = "ATGCGATCGATCGATCG"

# Count specific nucleotides
g_count = sequence.count("G")
c_count = sequence.count("C")
print(f"G: {g_count}, C: {c_count}")

# Create the complement (A↔T, G↔C)
complement_table = str.maketrans("ATGC", "TACG")
complement = sequence.translate(complement_table)
print(f"Original:   {sequence}")
print(f"Complement: {complement}")
print(f"RevComp:    {complement[::-1]}")

assert complement == "TACGCTAGCTAGCTAGC"
assert complement[::-1] == "CGATCGATCGATCGCAT"

Functions: A GC Content Calculator

A function is like packaging one step of an experiment protocol to be reusable. Build it once, apply it to any sequence.

python

def calculate_gc_content(sequence: str) -> float:
    sequence = sequence.upper()
    gc_count = sequence.count("G") + sequence.count("C")
    return (gc_count / len(sequence)) * 100

# Test
seq1 = "ATGCGATCGATCGATCG"
seq2 = "AAAAAAAAAA"
seq3 = "GGGGGGGGGG"

print(f"seq1 GC: {calculate_gc_content(seq1):.1f}%")
print(f"seq2 GC: {calculate_gc_content(seq2):.1f}%")
print(f"seq3 GC: {calculate_gc_content(seq3):.1f}%")

assert abs(calculate_gc_content(seq1) - 52.9) < 0.1
assert calculate_gc_content(seq2) == 0.0
assert calculate_gc_content(seq3) == 100.0

Lists and Loops: Processing Multiple Sequences at Once

Just as you apply the same treatment to every well in a 96-well plate, a for loop repeats the same operation on multiple data points.

python

def calculate_gc_content(sequence: str) -> float:
    sequence = sequence.upper()
    gc_count = sequence.count("G") + sequence.count("C")
    return (gc_count / len(sequence)) * 100

genes = {
    "BRCA1": "ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC",
    "TP53":  "ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA",
    "EGFR":  "ATGCGACCCTCCGGGACGGCCGGGGCAGCGCTCCTGGCGCTG",
}

results = []
for name, seq in genes.items():
    gc = calculate_gc_content(seq)
    results.append(gc)
    print(f"{name}: GC={gc:.1f}%, length={len(seq)}bp")

assert len(results) == 3
assert all(0 <= gc <= 100 for gc in results)

Dictionaries: Parsing FASTA Files

A dictionary is like labeling experiment samples. You can look up a sequence (value) directly by gene name (key).

python

def parse_fasta(fasta_text: str) -> dict[str, str]:
    sequences: dict[str, str] = {}
    current_header = ""
    for line in fasta_text.strip().split("\n"):
        if line.startswith(">"):
            current_header = line[1:].strip()
            sequences[current_header] = ""
        else:
            sequences[current_header] += line.strip()
    return sequences

sample_fasta = """>BRCA1_human
ATGGATTTATCTGCTCTTCG
CGTTGAAGAAGTACAAAATGTC
>TP53_human
ATGGAGGAGCCGCAGTCAG
ATCCTAGCGTGAGTTTGCTGTGA"""

result = parse_fasta(sample_fasta)
print(f"Sequences parsed: {len(result)}")
for name, seq in result.items():
    print(f"  {name}: {len(seq)}bp")

assert len(result) == 2
assert result["BRCA1_human"] == "ATGGATTTATCTGCTCTTCGCGTTGAAGAAGTACAAAATGTC"
assert result["TP53_human"] == "ATGGAGGAGCCGCAGTCAGATCCTAGCGTGAGTTTGCTGTGA"

Codon → Amino Acid Translation

Every 3 letters of a DNA sequence (a codon) specify one amino acid. You can build this translation table as a dictionary.

python

CODON_TABLE = {
    "ATG": "M",  # Methionine (start codon)
    "TTT": "F", "TTC": "F",  # Phenylalanine
    "TTA": "L", "TTG": "L", "CTT": "L", "CTC": "L",  # Leucine
    "GAT": "D", "GAC": "D",  # Aspartic acid
    "GAA": "E", "GAG": "E",  # Glutamic acid
    "GCT": "A", "GCC": "A",  # Alanine
    "TAA": "*", "TAG": "*", "TGA": "*",  # Stop codons
}

def translate_sequence(dna: str) -> str:
    protein = []
    for i in range(0, len(dna) - 2, 3):
        codon = dna[i:i+3]
        amino_acid = CODON_TABLE.get(codon, "?")
        if amino_acid == "*":
            break
        protein.append(amino_acid)
    return "".join(protein)

test_seq = "ATGGATTTTGAA"
protein = translate_sequence(test_seq)
print(f"DNA:     {test_seq}")
print(f"Protein: {protein}")

assert protein == "MDFE"

Try It Yourself (Faded Example)

Fill in the blanks to complete the GC content calculation function.

Fill in the Blankspython

def gc_content(seq):
  seq = seq.upper()
  gc = seq.count('G') + seq.count('')
  return  / len(seq) * 100

Common Errors & Solutions

Q: IndentationError: expected an indented block

Python uses indentation to define code blocks. Lines after def, for, if must be indented with 4 spaces.

Q: KeyError: 'BRCA1'

This occurs when a key doesn't exist in the dictionary. Use dict.get("BRCA1", "not found") to return a default value without an error.

Q: Counts are off because the sequence has mixed lowercase

Convert to uppercase first with .upper(). In actual FASTA files, lowercase letters often indicate repeat-masked regions.

In the next article, we'll learn to version control this analysis code with Git.

BioPlayground

Analyzing DNA Sequences with Python

Analyzing DNA Sequences with Python

After Completing This Topic

Variables: Containers for Data

Strings: Working with DNA Sequences

Functions: A GC Content Calculator

Lists and Loops: Processing Multiple Sequences at Once

Dictionaries: Parsing FASTA Files

Codon → Amino Acid Translation

Try It Yourself (Faded Example)

Common Errors & Solutions

💬 Questions & Comments