Applied - Sequence Composition Analysis

Core Principle: Sequence composition analysis quantifies GC content, cumulative GC skew, oligonucleotide frequencies, and codon usage bias using vectorized sliding window views.

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.

See also

Seq · Bio.SeqUtils · Vectorization · Memory Layout and Strides · Bio.motifs · Applied - FASTA and FASTQ Workflows · Broadcasting

Exercise

Biopython is not installed here, so parse the FASTA yourself; the composition statistics themselves are the exercise, and NumPy is available for the windowing.

Write composition(text), returning a list of output lines. The first input line holds the window size and the step, whitespace separated. What follows is a FASTA file whose sequences may be wrapped over several lines and may contain lower-case letters and ambiguity codes such as N, which are not A, C, G or T and must be excluded from every count. Process the records in input order.

Print records=<n> first. Then print one line per record, tab separated, holding the record name -- the first whitespace-separated token of the header, without the > -- the full sequence length including ambiguous bases, the GC fraction over the unambiguous bases rounded to four decimals, the 1-based start of the minimum-GC-skew window, and the most frequent 3-mer. GC skew for a window is (G - C) / (G + C), and a window whose G and C counts are both zero has a skew of zero. Take windows of the given size at the given step, starting at offset zero, and keep only whole windows that fit inside the sequence; break a tie for the minimum in favour of the earliest window, and print - for the window field when the sequence is shorter than one window. The most frequent 3-mer counts overlapping 3-mers made only of A, C, G and T, breaking ties in favour of the lexicographically first, and is - when there is none. A record with no unambiguous base has a GC fraction of 0.0.