You pulled a genome off NCBI, parsed it into a dict of contigs, and now you need two things per contig: the GC content, and the reverse complement so you can scan the minus strand. Both look like two-line problems. The first time you calculate GC content in Python you almost certainly write (seq.count("G") + seq.count("C")) / len(seq), get a plausible number, and move on. That line is wrong on a lot of real sequence data, and the reverse complement most people write is slow in a way that shows up the moment your sequences get big.
Neither bug is hard to fix. They are just quiet.
The two-line version, and when it is fine
Start with the honest baseline:
seq = "ATGGCGCTAGCTAGCTAGCGCGCTA"
gc = (seq.count("G") + seq.count("C")) / len(seq)
print(f"{gc:.3f}")Code language: Python (python)0.600Code language: plaintext (plaintext)str.count is implemented in C and does a single pass over the string, so this is genuinely fast. If your sequence is clean uppercase ACGT with no gaps, this is the right answer and you should not overthink it.
The problem is that clean uppercase ACGT is a thing that exists in tutorials.
Lowercase is not decoration
Reference genomes from UCSC and Ensembl ship soft-masked. RepeatMasker and Dustmasker lowercase the repetitive regions instead of deleting them, so the sequence is intact and the case tells you which bases fell inside a repeat. That is useful information, and it silently destroys the count above.
seq = "ATGGCGCTAGCTAGCTAGCGCGCTA"
soft = "ATGGCGctagctagctagCGCGCTA"
print(soft.upper() == seq)
print((soft.count("G") + soft.count("C")) / len(soft))Code language: Python (python)True
0.36Code language: plaintext (plaintext)Same sequence. 0.600 becomes 0.360, because count never matched the lowercase g and c. Nothing raises. You get a float, the float is plausible, and if you are plotting GC across a chromosome you will see a beautiful signal that is actually just a repeat annotation.
The fix is seq.upper() before counting. Worth knowing what that costs: .upper() allocates a whole new string, so on human chromosome 1 you are asking for another 250 MB. If that matters, count both cases instead and skip the copy. If it does not matter, and usually it does not, just call .upper() and stop thinking about it.
The fast way to calculate GC content in Python
The second problem is the denominator. Sequences contain N. Assemblies use runs of N for gap padding between scaffolds, and base callers emit N when they are not confident. An N is not an A, T, G, or C. Dividing by len(seq) treats every N as a base that is definitely not G or C, which drags GC content toward zero in proportion to how bad your assembly is.
The convention almost everyone uses is to compute GC over the bases you actually called:
def gc_content(seq):
"""GC fraction over unambiguous bases. Returns nan if there are none."""
seq = seq.upper()
gc = seq.count("G") + seq.count("C") + seq.count("S")
at = seq.count("A") + seq.count("T") + seq.count("W")
total = gc + at
if total == 0:
return float("nan")
return gc / totalCode language: Python (python)S and W are the two IUPAC ambiguity codes that are unambiguous about GC even though they are ambiguous about the base. S means “strong”, G or C. W means “weak”, A or T. Every other code (R, Y, K, M, B, D, H, V, N) straddles the GC boundary, so it goes in neither bucket and drops out of the denominator.
print(gc_content("ACGTNNNNACGTACGT"))
print((("ACGTNNNNACGTACGT").count("G") + ("ACGTNNNNACGTACGT").count("C")) / 16)Code language: Python (python)0.5
0.375Code language: plaintext (plaintext)Six of the twelve real bases are G or C, so 0.5 is the number you want. The len-based version reports 0.375 and the difference is entirely gap padding.
The nan return is deliberate. An all-N window has no GC content, and nan propagates through numpy and plots as a hole instead of pretending the answer is zero. If you would rather it raise, raise. What you should not do is return 0.0.
Reverse complements with str.translate
Here is the version almost everyone writes first:
COMP = {"A": "T", "C": "G", "G": "C", "T": "A", "N": "N"}
def reverse_complement_slow(seq):
return "".join(COMP[base] for base in reversed(seq))Code language: Python (python)It works, it reads clearly, and it has the same lowercase bug as before except this one is loud: reverse_complement_slow("atgc") raises KeyError: 'a'. That is arguably an improvement. Any ambiguity code beyond N blows up too.
The real issue is that it runs one Python-level dict lookup per base. On a 5 Mb bacterial genome that is five million trips through the interpreter. Python has a built-in for exactly this shape of problem:
COMPLEMENT = str.maketrans("ACGTRYSWKMBDHVNacgtryswkmbdhvn",
"TGCAYRSWMKVHDBNtgcayrswmkvhdbn")
def reverse_complement(seq):
return seq.translate(COMPLEMENT)[::-1]Code language: Python (python)str.maketrans builds a mapping from character codepoints to codepoints once, at import time. str.translate then walks the string in C and does the whole substitution in one pass, and the slice reverses it in another. No interpreter loop.
The mapping is the full IUPAC set, which is less exotic than it looks. R (A or G) complements to Y (C or T), K (G or T) to M (A or C), B (not A) to V (not T), D (not C) to H (not G). S, W, and N are their own complements. Case is handled by just listing the lowercase letters too, so soft-masking survives the round trip:
print(reverse_complement("ATGGCGctagctagctagCGCGCTA"))Code language: Python (python)TAGCGCGctagctagctagCGCCATCode language: plaintext (plaintext)Worth actually measuring rather than taking my word for it:
import timeit, random
seq = "".join(random.choices("ACGT", k=1_000_000))
print(timeit.timeit(lambda: reverse_complement(seq), number=10))
print(timeit.timeit(lambda: reverse_complement_slow(seq), number=10))Code language: Python (python)The exact numbers depend on your machine, but this is not a few percent. It is a large multiple, and it grows with sequence length.
If you are reading sequence straight off disk and never need it as text, bytes has the same API and skips the unicode overhead entirely:
COMPLEMENT_BYTES = bytes.maketrans(b"ACGTRYSWKMBDHVNacgtryswkmbdhvn",
b"TGCAYRSWMKVHDBNtgcayrswmkvhdbn")
def reverse_complement_bytes(seq: bytes) -> bytes:
return seq.translate(COMPLEMENT_BYTES)[::-1]Code language: Python (python)translate never raises
One catch you need to know about. str.translate leaves any character that is not in the table exactly where it found it. That means garbage passes through untouched:
print(reverse_complement("ATGX"))Code language: Python (python)XCATCode language: plaintext (plaintext)No error. The X is still there, now sitting at the front of a string you are about to treat as sequence. If your input comes from anywhere you do not control, validate once and let the fast path stay fast:
VALID = set("ACGTRYSWKMBDHVNacgtryswkmbdhvn")
def check(seq):
bad = set(seq) - VALID
if bad:
raise ValueError(f"unexpected characters: {sorted(bad)}")
return seqCode language: Python (python)set(seq) on a long string is one pass and the resulting set is at most a couple dozen elements, so this is cheap even on a chromosome. The usual culprits are - from a gapped alignment you forgot to degap, * from a stop codon in a file that turned out to be protein, and \r from a FASTA written on Windows.
Where this actually gets used
Whole-sequence GC content is a summary statistic and mostly tells you which organism you are looking at. The interesting version is GC along the sequence, which is how you find CpG islands, horizontally transferred regions, and coverage bias in your library prep:
def sliding_gc(seq, window, step):
for start in range(0, len(seq) - window + 1, step):
yield start, gc_content(seq[start:start + window])
demo = "AAAA" * 5 + "GCGC" * 5 + "AAAA" * 5
for start, gc in sliding_gc(demo, window=20, step=20):
print(start, f"{gc:.2f}")Code language: Python (python)0 0.00
20 1.00
40 0.00Code language: plaintext (plaintext)That is where the nan decision pays off. Slide a window across a real assembly and some windows land entirely in a gap. Getting nan back means those positions disappear from your plot, which is true, instead of drawing a cliff to zero, which is a lie you will spend an afternoon investigating.
If you would rather not calculate GC content in Python yourself, Biopython’s Bio.SeqUtils.gc_fraction handles ambiguity codes and has an argument controlling how they are counted. Read its docstring before you trust the default, because “what do we do with N” is a choice and different tools choose differently. Knowing which choice you made is the part that matters.
