ML for Biology

Lesson 11 of 13 · 12 min

Deep Learning Teaser: CNNs on Sequence

Transcription-factor binding sites are short sequence motifs, and for decades we modeled them with position weight matrices (PWMs) that score how well each position matches a preferred base. A convolutional filter over one-hot-encoded DNA is a strict generalization of that idea: it is a learnable weight matrix that scans the sequence for a pattern, except the network fits the pattern from labeled data instead of you specifying it. In this lesson you will one-hot a sequence, build a small 1D CNN for binding-site classification, and weigh it against a k-mer-plus-tree baseline so you know when the extra machinery is worth it.

One-hot sequence tensor

python
import numpy as np
import torch

BASES = "ACGT"
base_to_idx = {b: i for i, b in enumerate(BASES)}

def one_hot(seq: str) -> np.ndarray:
    # returns (4, L): channels-first, the layout nn.Conv1d expects
    x = np.zeros((4, len(seq)), dtype=np.float32)
    for j, base in enumerate(seq.upper()):
        i = base_to_idx.get(base)
        if i is not None:      # leave N / ambiguous bases as an all-zero column
            x[i, j] = 1.0
    return x

seqs = ["ACGTACGTNN" * 20, "TTGGCCAAGG" * 20]   # two 200 bp windows
batch = torch.from_numpy(np.stack([one_hot(s) for s in seqs]))
print(batch.shape, batch.dtype)
torch.Size([2, 4, 200]) torch.float32

Filters as motif detectors

A 1D convolution slides a bank of filters along the length axis. Each filter is a weight matrix of shape 4 by kernel_size, so a width-12 filter is a PWM-like scanner that produces a match score at every position by cross-correlating its weights with the one-hot window. ReLU zeroes out weak matches, keeping only positions where the motif is clearly present, and a global max pool collapses each filter to its single strongest hit anywhere in the sequence. That pooling gives position invariance: the classifier learns whether a motif occurs, not exactly where, which is what a binding-site call needs. After training you can visualize each filter as a motif: a quick approximation softmaxes its weights across the four base channels into a position probability matrix, while the more faithful DeepBind/Basset approach stacks the subsequences that most strongly activate the filter and builds a PWM or sequence logo from those aligned hits.

python
import torch.nn as nn
from sklearn.metrics import roc_auc_score

class SeqCNN(nn.Module):
    def __init__(self, n_filters=64, kernel=12):
        super().__init__()
        # 64 filters of width 12 = 64 learnable, PWM-like motif scanners over 4 channels
        self.conv = nn.Conv1d(4, n_filters, kernel_size=kernel)
        self.act  = nn.ReLU()
        self.pool = nn.AdaptiveMaxPool1d(1)   # global max: best match anywhere
        self.head = nn.Linear(n_filters, 1)

    def forward(self, x):                       # x: (N, 4, L)
        z = self.pool(self.act(self.conv(x)))   # (N, n_filters, 1)
        return self.head(z.squeeze(-1)).squeeze(-1)   # (N,) logits

model = SeqCNN()
opt = torch.optim.Adam(model.parameters(), lr=1e-3)
loss_fn = nn.BCEWithLogitsLoss()
print("trainable params:", sum(p.numel() for p in model.parameters()))

for epoch in range(1, 6):
    model.train()
    for xb, yb in train_loader:            # xb: (N,4,200), yb: (N,) float in {0.,1.}
        opt.zero_grad()
        loss_fn(model(xb), yb).backward()
        opt.step()
    model.eval()
    with torch.no_grad():
        p = torch.sigmoid(model(X_val)).numpy()
    print(f"epoch {epoch}  val AUROC {roc_auc_score(y_val, p):.3f}")
trainable params: 3201
epoch 1  val AUROC 0.742
epoch 2  val AUROC 0.836
epoch 3  val AUROC 0.881
epoch 4  val AUROC 0.902
epoch 5  val AUROC 0.911

Versus k-mer plus trees

Before committing to a GPU, benchmark the classic baseline. Count k-mers (for example every 6-mer with scikit-learn's CountVectorizer using analyzer='char' and ngram_range=(6, 6), which yields overlapping 6-mer counts) and feed the counts to gradient-boosted trees such as HistGradientBoostingClassifier or XGBoost. On a few thousand labeled windows this pipeline routinely reaches AUROC around 0.86, trains in seconds on a CPU, and its feature importances are directly interpretable as informative k-mers. The CNN above edged ahead to roughly 0.91, but that gain comes only because convolutions can compose motifs and tolerate variable spacing that a fixed k-mer vocabulary cannot represent.

Deep learning is not free. Below about 10^4 labeled sequences the k-mer-plus-tree model usually matches or beats a CNN, trains on CPU in seconds, and ships more easily. CNNs typically start to pull clearly ahead in the 10^4 to 10^5-plus range, where you want a GPU: production-scale sequence models (deep, wide, long inputs, large batches) train roughly 10 to 50 times faster on a GPU than a CPU, though a toy like the one above is small enough that kernel-launch overhead can wipe out any gain. Budget for hyperparameter search, early stopping on a held-out set, and a fixed seed before you trust a 0.05 AUROC gap over the baseline.

Try it yourself: swap AdaptiveMaxPool1d(1) for AdaptiveAvgPool1d(1) and retrain. Max pooling reports each filter's single strongest match (motif presence), while average pooling reports its mean activation across the window (motif density); compare val AUROC on your binding-site set. Then extract the conv weights, softmax across the 4 base channels for a quick position-probability view, and plot the top filters as sequence logos with logomaker. Do the recovered motifs match known JASPAR PWMs for your transcription factor?