Working with sequences

Why it matters

Sequencing produces text: long strings of A, C, G and T. Before and after every experiment you will handle sequences: checking a primer site, finding a restriction site before cloning, measuring GC content of a region that is hard to amplify, or confirming that a read contains an expected motif. Doing this by eye is slow and error-prone. A few lines of Python do it exactly, every time, and the same ideas power the big bioinformatics tools.

Sequences are strings

In Python, a DNA sequence is just a string. Real sequences arrive messy, with lowercase letters, spaces and line breaks (for example when copied from a FASTA file), so the first step is always to clean them:

def clean(seq):
    return "".join(seq.split()).upper()

clean("gaat tc\nGGA")   # 'GAATTCGGA'

seq.split() splits on any whitespace and "".join(...) glues the pieces back together.

The FASTA format

Most sequence files and databases use FASTA: a header line starting with >, then sequence lines.

>sample_01 partial gene region
GATATATGCATATACTT

Reading a FASTA file means skipping the > line and joining the rest.

GC content

GC content is the percentage of bases that are G or C. It matters because G-C pairs have three hydrogen bonds and A-T pairs two, so GC-rich DNA melts at higher temperature, is harder to amplify and can form stable secondary structures.

def gc_content(seq):
    s = clean(seq)
    return 100 * (s.count("G") + s.count("C")) / len(s)

Worked example. GATATATGCATATACTT has 17 bases: G 2, C 2. GC content = 4 / 17 x 100 = 23.5 percent: AT-rich.

Reverse complement

DNA is double-stranded and the two strands run in opposite directions. The other strand, read 5' to 3', is the reverse complement: swap A with T and G with C, then reverse.

pairs = {"A": "T", "T": "A", "G": "C", "C": "G"}
"".join(pairs[b] for b in reversed("ATGCC"))   # 'GGCAT'

Many restriction sites are palindromes: GAATTC (EcoRI) is its own reverse complement, so the enzyme recognises it on both strands.

Finding a motif

A motif is a short sequence of interest: a restriction site, a primer binding site, a promoter element. You often need every position where it occurs, including overlapping matches. Python's str.find and str.count skip overlaps, so a sliding window is safer:

for i in range(len(s) - len(m) + 1):
    if s[i:i + len(m)] == m:
        positions.append(i + 1)   # biologists count from 1

Worked example. In GATATATGCATATACTT, the motif ATAT starts at positions 2, 4 and 10 (1-based). The matches at 2 and 4 overlap. "GATATATGCATATACTT".count("ATAT") returns only 2, which would be wrong for this task.

Notice the off-by-one detail: Python counts from 0, but sequence coordinates in biology (and in tools such as BLAST and genome browsers) count from 1. Always state which you use.

Libraries for real work

For bigger jobs, bioinformaticians use libraries such as Biopython, which reads FASTA and FASTQ files, translates sequences, runs alignments and talks to NCBI. The ideas are the ones above; the library handles edge cases and file formats. Learning the plain-Python version first means you understand what the library is doing. (Biopython is not available in this site's activity runner, so the activity uses plain Python.)

Remember also to treat sequence data with the same care as any lab record: keep the original file, note where it came from, and never edit a sequence by hand without recording why.

Common mistakes

• Forgetting to remove whitespace and line breaks, which distorts length and GC content.
• Mixing uppercase and lowercase, so "atg" == "ATG" is False.
• Counting non-overlapping matches when you need all of them.
• Reporting 0-based positions as if they were 1-based.
• Writing the complement but forgetting to reverse it.

Quick recap

• Clean sequences first: remove whitespace, uppercase.
• GC content = (G + C) / length x 100.
• Reverse complement = complement each base, then reverse.
• Find motifs with a sliding window to catch overlaps, and report 1-based positions.

Try it

Open the Do tab. The GC and reverse-complement functions are written; you fill in the motif check.