Applied - Sequence Composition Analysis

In one line: GC content, k-mers, codon usage and sliding windows — where Biopython's objects meet NumPy's speed.

Uses: Seq · Bio.SeqUtils · ndarray · Vectorization · Memory Layout and Strides · Broadcasting

Encoding sequence as an array

The move that unlocks everything else:

import numpy as np
from Bio import SeqIO

BASES = np.array(list("ACGT"))

def encode(seq):
    """DNA string → int8 codes; N and other letters → 4."""
    arr = np.frombuffer(str(seq).upper().encode(), dtype=np.uint8)
    lut = np.full(256, 4, dtype=np.int8)
    for i, b in enumerate(b"ACGT"):
        lut[b] = i
    return lut[arr]

codes = encode(record.seq)              # (L,) of 0..4

np.frombuffer on the encoded bytes plus a 256-entry lookup table converts a 250 Mb chromosome in well under a second. A Python loop over the same string takes minutes. See Vectorization.

onehot = (codes[:, None] == np.arange(4)[None, :]).astype(np.int8)   # (L, 4)

That is Broadcasting: (L,1) against (1,4) gives (L,4).

GC content

from Bio.SeqUtils import gc_fraction

gc_fraction(record.seq)                        # 0.0–1.0 — a FRACTION
gc_fraction(record.seq, ambiguous="ignore")    # exclude N from the denominator

Remember SeqUtils.GC() was removed in 1.82 and returned a percentage. A naive swap introduces a 100× error. See Biopython Deprecations.

Vectorised, over an array:

is_gc = (codes == 1) | (codes == 2)            # C or G
is_n = codes == 4
gc = is_gc.sum() / (~is_n).sum()               # exclude N

Sliding windows — free, via strides

from numpy.lib.stride_tricks import sliding_window_view

W = 1000
gc_windows = sliding_window_view(is_gc, W).mean(axis=1)
n_windows = sliding_window_view(is_n, W).mean(axis=1)
gc_windows[n_windows > 0.5] = np.nan           # mask N-rich windows

fig, ax = plt.subplots(figsize=(11, 3))
ax.plot(np.arange(len(gc_windows)) + W // 2, gc_windows, lw=0.5)
ax.axhline(np.nanmean(gc_windows), ls="--", c="crimson", lw=0.8)
ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda v, _: f"{v/1e6:.0f}"))
ax.set(xlabel="Position (Mb)", ylabel="GC fraction")

sliding_window_view builds the (L-W+1, W) window matrix by stride manipulation — it allocates nothing. The memory cost appears only at the reduction. See Memory Layout and Strides.

Masking N-rich windows matters: a centromeric gap of Ns otherwise reads as GC = 0, producing a dramatic and meaningless dip.

GC skew

from Bio.SeqUtils import gc_skew

is_g = codes == 2
is_c = codes == 1
W = 10_000
g = sliding_window_view(is_g, W).sum(axis=1)[::W]      # non-overlapping
c = sliding_window_view(is_c, W).sum(axis=1)[::W]
skew = (g - c) / np.maximum(g + c, 1)
cum = np.cumsum(skew)

origin = int(np.argmin(cum) * W)
terminus = int(np.argmax(cum) * W)

In bacteria, cumulative GC skew has a minimum at the replication origin and a maximum at the terminus — leading and lagging strands have asymmetric mutational biases. A clean V-shape in a bacterial genome is a good sign; a noisy one suggests a misassembly.

k-mer counting

def kmer_counts(codes, k=6):
    valid = codes < 4
    windows = sliding_window_view(codes, k)
    ok = sliding_window_view(valid, k).all(axis=1)      # skip windows with N
    powers = 4 ** np.arange(k - 1, -1, -1)
    idx = (windows[ok] * powers).sum(axis=1)             # base-4 encoding
    return np.bincount(idx, minlength=4**k)

counts = kmer_counts(codes, k=6)                         # 4096 entries

def kmer_name(i, k):
    return "".join("ACGT"[(i >> (2*(k-1-j))) & 3] for j in range(k))

top = np.argsort(counts)[::-1][:20]
for i in top:
    print(kmer_name(i, 6), counts[i])

Encoding each k-mer as a base-4 integer turns counting into a single np.bincount. For k=6 that is 4,096 bins; for k=12 it is 16 million (still fine); beyond k≈14 you need a hash-based counter (jellyfish, KMC).

k-mer spectra are how you estimate genome size and heterozygosity from raw reads before assembling (GenomeScope), and how you detect contamination.

Codon usage

from Bio import SeqIO
from collections import Counter
import pandas as pd

def codon_table(cds_records):
    c = Counter()
    for rec in cds_records:
        s = str(rec.seq).upper()
        s = s[:len(s) - len(s) % 3]
        c.update(s[i:i+3] for i in range(0, len(s), 3))
    return pd.Series(c).sort_values(ascending=False)

from Bio.Data import CodonTable
std = CodonTable.unambiguous_dna_by_id[1]
aa_of = {**std.forward_table, **{s: "*" for s in std.stop_codons}}

df = codon_table(cds).rename("count").to_frame()
df["aa"] = df.index.map(aa_of)
df["rscu"] = df.groupby("aa")["count"].transform(
    lambda s: s / s.mean()                      # relative synonymous codon usage
)

fig, ax = plt.subplots(figsize=(11, 3.5))
sns.barplot(data=df.reset_index().sort_values(["aa", "index"]),
            x="index", y="rscu", hue="aa", dodge=False, ax=ax, legend=False)
ax.tick_params(axis="x", rotation=90, labelsize=5)
ax.axhline(1, ls="--", c="grey", lw=0.8)

RSCU (relative synonymous codon usage) normalises within each amino acid, so RSCU = 1 means "used exactly as often as its synonyms". Values far from 1 indicate codon bias, which correlates with expression level and tRNA abundance.

The groupby().transform() is the per-group normalisation pattern from groupby — normalising within amino acid without a merge.

GC3 and codon bias

from Bio.SeqUtils import GC123

stats = pd.DataFrame([
    {"gene": r.id, "gc": GC123(r.seq)[0]/100, "gc1": GC123(r.seq)[1]/100,
     "gc2": GC123(r.seq)[2]/100, "gc3": GC123(r.seq)[3]/100, "length": len(r)}
    for r in SeqIO.parse("cds.fasta", "fasta")
])

sns.scatterplot(data=stats, x="gc3", y="gc", alpha=0.3)

GC3 — GC at the third codon position — is the standard codon-bias measure, because third positions are largely synonymous and therefore reflect mutational and translational pressure rather than protein constraint. A "GC3 plot" (GC3 against overall GC) is a classic comparative-genomics figure.

Sequence complexity

def shannon_entropy(codes, k=1):
    valid = codes[codes < 4]
    if k == 1:
        p = np.bincount(valid, minlength=4) / len(valid)
    else:
        p = kmer_counts(codes, k); p = p / p.sum()
    p = p[p > 0]
    return -(p * np.log2(p)).sum()

def dust_score(codes, window=64):
    """Low-complexity detection, roughly like the DUST algorithm."""
    tri = sliding_window_view(codes, 3)
    ok = (tri < 4).all(axis=1)
    idx = (tri[ok] * np.array([16, 4, 1])).sum(axis=1)
    counts = np.bincount(idx, minlength=64)
    return (counts * (counts - 1)).sum() / (2 * max(len(idx) - 1, 1))

Low-complexity regions (homopolymers, simple repeats) produce spurious alignments and inflated motif scores. Most search tools mask them by default; if you write your own scan, you need to.

CpG islands

def cpg_stats(codes, window=200, step=1):
    is_c, is_g = codes == 1, codes == 2
    cg = np.zeros(len(codes), dtype=bool)
    cg[:-1] = is_c[:-1] & is_g[1:]                     # CpG dinucleotide

    W = window
    c_n = sliding_window_view(is_c, W).sum(axis=1)
    g_n = sliding_window_view(is_g, W).sum(axis=1)
    cg_n = sliding_window_view(cg, W).sum(axis=1)

    gc_frac = (c_n + g_n) / W
    expected = np.maximum(c_n * g_n / W, 1e-9)
    obs_exp = cg_n / expected
    return gc_frac, obs_exp

gc_frac, obs_exp = cpg_stats(codes)
island = (gc_frac > 0.55) & (obs_exp > 0.65)           # Takai & Jones criteria

The observed/expected CpG ratio is the informative quantity: vertebrate genomes are globally CpG-depleted (methylated CpG deaminates to TpG), so an undepleted region marks an unmethylated CpG island, which usually marks a promoter.

Bringing it into pandas

summary = pd.DataFrame([
    {"id": r.id, "length": len(r),
     "gc": gc_fraction(r.seq),
     "n_frac": str(r.seq).upper().count("N") / len(r),
     "entropy": shannon_entropy(encode(r.seq)),
     "gc3": GC123(r.seq)[3] / 100 if len(r) % 3 == 0 else np.nan}
    for r in SeqIO.parse("sequences.fasta", "fasta")
])

sns.pairplot(summary[["length", "gc", "entropy", "n_frac"]],
             diag_kind="kde", corner=True, plot_kws=dict(s=8, alpha=0.4))

Common mistakes

SeqUtils.GC() — removed; 100× different from gc_fraction.
Not upper-casing soft-masked sequence.
Counting N in the GC denominator without deciding to.
Windows spanning assembly gaps, producing meaningless dips.
Python loops over long sequences. Encode to an array first.
as_strided instead of sliding_window_view. The former reads past the buffer.
Not filtering low-complexity regions before motif or alignment work.
Comparing GC between species without controlling for gene set.
Codon usage on sequences whose length is not a multiple of 3.

Applied - Sequence Composition Analysis

Applied - Sequence Composition Analysis

Encoding sequence as an array

GC content

Sliding windows — free, via strides

GC skew

k-mer counting

Codon usage

GC3 and codon bias

Sequence complexity

CpG islands

Bringing it into pandas

Common mistakes

See also