Write-up - Sequence Aligner

←Back to the aligner Group A1: Sophia, Thomas, Brady, and Alex

The question we set out to answer

If the alignment you get depends on the scoring numbers you pick, how do you know which alignment is actually right?

We built our own tunable aligner in Python, starting from the global and local alignment code we wrote this module and generalizing it to proteins. Then we used it on real DNA and protein, hemoglobin across species and the SARS spike gene, to figure out what makes one alignment better than another. The short version is that the raw score cannot be the judge, because we chose the numbers that produce it. This page walks through how we know that, and what we trust instead.

Good settings versus careless settings

Here is the whole problem in miniature. These are the same two short DNA sequences, aligned twice, changing only the gap penalty.

Toy DNA aligned with good and bad gap settings — Our aligner on a toy pair. Good settings (gap = 3) give one clean insertion: score 4, identity 92.3 percent, one gap. Careless settings (gap = 0.1) score 23.2 with a perfect 100 percent identity, but the alignment shatters into five scattered gaps to dodge every mismatch. The higher score and higher identity are the worse answer.

You might think that was an artifact of a tiny made-up example, so we tried it on a real gene: hemoglobin beta, human versus zebrafish, scored with the BLOSUM62 protein matrix.

Hemoglobin beta aligned with good and bad gap settings — Human versus zebrafish hemoglobin beta. Good settings (gap = 10) give score 419 with one clean gap. Free gaps (gap = 0) raise the score to 497 and even raise the identity, but the gap count explodes from 1 to 29 separate insertions and deletions. Twenty-nine independent events is a far less believable history than one.

So even on real biology, the settings that score higher and read as more identical are the ones that are biologically wrong. That is the trap this whole project is about. Open the aligner, load one of these examples, and drag the gap slider down to zero to watch it happen live.

Why one gap penalty is not enough

The reason careless gap settings break things is that a single flat gap penalty is asked to do two opposite jobs at once. It has to allow one long real insertion, but block many scattered fake ones. With a linear penalty a gap of length nine costs exactly nine times the penalty, so a real block insertion is charged the same per letter as nine unrelated single-letter gaps. The model literally cannot tell a chunky indel from confetti.

Linear versus affine gap penalties — The fix is affine gaps: charge a lot to open a gap and a little to extend it. Then one long gap is cheap to grow but scattered gaps each pay the opening cost. The affine setting in the aligner turns this on for global mode.

Why we align protein, not just DNA

We took the SARS-CoV and SARS-CoV-2 spike genes, aligned them as DNA and again as the protein they code for, then sorted every change.

Synonymous versus nonsynonymous changes in the SARS spike — DNA identity is 77.3 percent and protein identity is 78.0 percent, almost the same headline. But of the 713 DNA changes, 438 are synonymous: the DNA letter changed but the amino acid did not. Only 275 actually change the protein. The DNA sees hundreds of differences the protein never notices, which is the redundancy of the genetic code.

That is the difference between a mutation that exists and a mutation that matters, and it is why the amino acid alignment tells you more about what the virus is actually doing.

Seeing the changes in 3D

Those 275 protein changes are not abstract. We drew them onto the real SARS-CoV-2 spike structure from the Protein Data Bank.

Changed residues mapped onto the spike structure — Gray is the conserved backbone. Red marks the 297 amino acids that changed, about 23 percent of the 1273-residue spike. Our own aligner decides which residues are red. They cluster on the outer surface, the part that grabs human cells and that antibodies target.

So how do we judge a good alignment?

Since the score cannot be trusted on its own, here are the six things we actually look at. What they have in common is that they do not move when we change the reward knobs.

Percent identity. It does not depend on the reward you set, so it ranks pairs honestly. Human and gorilla come out at 99.3 percent, human and zebrafish at 51.7 percent, which is the right order.

Gap structure. A few long gaps beat a spray of tiny ones, even when the tiny ones score higher. One insertion is one event; twenty-nine is not.

Conserved regions line up. The parts we already know matter, like the two histidines that grip the heme in hemoglobin, should land on top of each other.

Better than chance. We shuffle one sequence thousands of times and align those. If the real score is not far out past the shuffled ones, it is not meaningful.

Independent evidence agrees. A known structure, or a professional tool, should tell the same story. Ours does (see below).

The history is believable. The mutations the alignment implies should be a plausible way for one sequence to have become the other.

The single best tell is stability. When the same alignment survives as we wiggle every parameter, that steadiness is the real answer, not the biggest number.

Further directions we explored

Once the core aligner worked, we kept pushing it. The interactive page runs all of the alignment methods; the rest we built and tested alongside them.

Semi-global alignment Affine gaps Banded whole-genome alignment GPU acceleration On-device AI Dot plots and conservation maps

Whole genomes. A full alignment table for two 30,000-letter SARS genomes would be billions of cells and run out of memory, so we added a banded method that only fills the stripe near the diagonal. It aligns two whole genomes in about a second and drops the memory from quadratic to linear in the band width.

Speed. For jobs that repeat, like the shuffle test that aligns against thousands of random sequences, we moved the work onto the graphics card with a diagonal wavefront method. On an Intel Arc GPU it runs about 28 times faster than the processor and gives the exact same scores, and it falls back to the processor automatically when there is no GPU.

On-device AI. In the desktop version we trained two small neural networks from scratch that run locally: one predicts secondary structure (helix, sheet, or coil) at about 85 percent accuracy, and a transformer language model turns a protein into a fingerprint that sorts proteins into the right family about 90 percent of the time.

A sanity check. To be sure our aligner is correct and not just plausible, we aligned human versus gorilla hemoglobin beta with our engine and with Biopython, a tool professionals use, under identical settings. Both returned the exact same optimal score, 777.0, and the exact same 99.32 percent identity. Not close, identical.

Try it yourself

Everything above comes out of our own aligner. Head back to the tool, load an example, and move the sliders. The one line we would leave you with: the score rewards exactly what you tell it to reward, and biology decides what is actually true.

←Back to the aligner

Group A1: Sophia, Thomas, Brady, and Alex. CMU Pre-College Computational Biology, Module 2.