PAM, BLOSUM and scoring matrices from 1,000 human mitochondrial genomes
Two sequences can agree at most positions and still leave us with a few awkward decisions. Is a leucine opposite an isoleucine a reasonable match? Should an [alignment](<en.wikipedia.org>) insert a gap here, or accept three mismatches and keep going? A scoring scheme puts numbers on those choices, and changing the numbers can change the alignment you get.
PAM and BLOSUM are families of substitution matrices for [proteins](<genome.gov>). They describe how amino-acid pairs occur in aligned related sequences, relative to the pairs we would expect from their individual abundance. We will work through those counts, then use BioseqDB to calculate [nucleotide](<genome.gov>) and protein comparisons from 1,000 human mitochondrial [genome](<genome.gov>) records. [SQL](<postgresql.org>) will do the counting and matrix calculations. Python will draw the results.
A protein is a chain of [amino acids](<genome.gov>). A coding region of [DNA](<genome.gov>) supplies its instructions in three-letter groups called [codons](<genome.gov>), each specifying an amino acid or a stop signal.
Mitochondria carry a small genome of their own, separate from the nuclear [chromosomes](<genome.gov>). That gives us complete mitochondrial records small enough to inspect and align in a tutorial, along with protein-coding [genes](<genome.gov>) that connect the DNA analysis to BLOSUM. Keep the two alphabets in mind: the A in a DNA sequence is adenine, while the A in BLOSUM is alanine, an amino acid.
What a scoring matrix does
An alignment places two sequences into corresponding columns. A column can contain the same letter twice, two different letters, or a letter opposite a gap. Those gaps let us represent insertions and deletions, usually shortened to indels.
The scoring system has two parts. A substitution matrix assigns a score to each aligned letter pair, and a gap model charges for inserting or extending gaps. If we call the alignment 𝒜, its score is
!S(\mathcal A)= \sum_{(a,b)\in\text{paired columns of }\mathcal A}s(a,b) -\sum_{g\in\text{gap runs of }\mathcal A}G(|g|).
Here s(a,b) is the substitution score and G(|g|) is the cost of a gap run of length |g|. The algorithm searches for a high-scoring arrangement under that rule. A global alignment consumes both full sequences, whereas a local alignment can select matching regions within them. The [alignment foundations chapter](<bioseqdb.com>) gives the formal definitions and shows the different boundary conditions.
Why use biological scores at all? Related proteins can accumulate substitutions while retaining a similar structure and function. A scoring matrix lets an alignment retain evidence from those changes instead of relying entirely on identical letters. We can then assess the resulting score against a chance model, which helps identify similarities consistent with common ancestry. The significance also depends on sequence length, composition and the size of the database being searched.
Start with an identity matrix
The smallest useful example rewards a match with +1 and charges -1 for a mismatch:
!s(a,b)= \begin{cases} +1,&a=b,\ -1,&a\ne b. \end{cases}
Some versions give mismatches zero instead. Either way, all mismatches receive the same treatment. For DNA this is a useful starting point, and simple match/mismatch schemes appear in many nucleotide aligners. It can also score proteins, although we lose the distinction between chemically similar and very different amino acids.
Leucine and isoleucine, written L and I, have similar hydrophobic side chains, which tend to avoid water. Replacing one with the other can preserve the local environment of a protein. Glycine and tryptophan, G and W, are a different pairing: glycine is small and flexible, while tryptophan has a bulky aromatic side chain. Their exchange can demand a much larger structural adjustment.
We call the first kind of replacement conservative and the second radical. The surrounding protein still matters. A substitution that is tolerated at one position can disrupt another, especially in a binding site or tightly packed core. A general substitution matrix averages observations across many positions and protein families.
!An identity matrix gives all mismatches the same score, while BLOSUM62 distinguishes amino-acid pairs.
An identity matrix gives all mismatches the same score, while BLOSUM62 distinguishes amino-acid pairs.
The letters in the protein examples are amino-acid codes. BLOSUM62 rewards L/I with +2, while G/W receives -2.
From observed pairs to log odds
Suppose Q_ab is the probability of an ordered pair (a,b) in our aligned related sequences. Let pₐ and p_b be the background frequencies of those two residues. Under an independent chance model, the expected pair frequency is p_a p_b.
The odds ratio compares the two:
!O_{ab}=\frac{Q_{ab}}{p_ap_b}.
An odds ratio above one means the pair appears more frequently in the related-sequence alignments than the background frequencies predict. Below one means it appears less frequently. With a model that combines evidence from individual columns, taking a logarithm turns multiplication into addition:
!s_{ab}=\log_2\left(\frac{Q_{ab}}{p_ap_b}\right), \qquad \log_2\prod_r O_{a_rb_r}=\sum_r\log_2 O_{a_rb_r}.
This is convenient for alignment algorithms: they can add a score as they advance through each paired column. Positive entries favor an alignment, negative entries reduce its score, and zero means that the observed and expected frequencies agree.
BLOSUM62 stores integer scores in half-bit units. Before rounding, its scoring equation has the form used in the [ForensicMCQ BLOSUM guide](<forensicmcq.com>):
!S_{ij}=2\log_2\left(\frac{\text{observed frequency}}{\text{expected frequency}}\right).
For a hypothetical pair observed with frequency 0.04 and expected with frequency 0.01, the calculation is
!2\log_2(0.04/0.01)=2\log_2(4)=2\times2=+4.
Those numbers illustrate the arithmetic. Actual BLOSUM62 entries come from its measured pair frequencies. For example, the aligned protein fragments ALK and AID score
!S(\mathrm{ALK},\mathrm{AID})= S(\mathrm A,\mathrm A)+S(\mathrm L,\mathrm I)+S(\mathrm K,\mathrm D) =4+2+(-1)=5.
There are no gaps in that example, so we can simply add the three matrix entries. Longer alignments use the same operation repeatedly. Here are both calculations in SQL:
SELECT round(2 * log(2::numeric, 0.04 / 0.01)) AS half_bit_score;
SELECT sum(score) AS alignment_score
FROM (VALUES ('A', 'A', 4), ('L', 'I', 2), ('K', 'D', -1))
AS aligned(first_residue, second_residue, score);
Result 1
half_bit_score / numeric
------------------------
4
Result 2
alignment_score / bigint
------------------------
5
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
PAM: estimating change over evolutionary distance
Margaret Dayhoff and her colleagues developed the PAM approach during the 1960s and 1970s. Their original data came from 71 groups of closely related proteins, with alignments above roughly 85% identity. At that distance, ancestral reconstruction and counting changes along a tree are easier than they would be for highly divergent sequences.
PAM stands for Point Accepted [Mutation](<genome.gov>). An accepted amino-acid replacement is one retained in an evolving lineage. It may preserve protein function well enough to persist, without giving the organism any particular advantage.
One PAM unit represents an evolutionary interval with an expected one accepted replacement per 100 amino-acid positions. It is a measure of accumulated change. The amount of elapsed time needed to accumulate it depends on the proteins and lineages under study.
The construction starts with global alignments of related proteins. From these alignments, a phylogenetic tree and likely ancestral sequences let us infer replacements along branches. Counts of those replacements, together with residue abundance, provide estimates of relative mutability, mₐ: how readily each amino acid changes in the fitted model. The probabilities are then normalized to the small amount of change represented by one PAM unit.
This gives a PAM1 transition matrix, M(1). An entry describes the probability of moving from one amino acid to another over that unit of evolutionary distance. The diagonal contains the probabilities of retaining the same residue.
To obtain a model for a longer distance, multiply the transition matrix by itself:
!M(n)=M(1)^n.
For PAM250, that is M(1)²⁵⁰. Each multiplication sums over possible intermediate residues, so a site can change several times and can return to its starting state. After calculating the transition probabilities at the chosen distance, we convert them into log-odds alignment scores using the background frequencies.
The distinction between these two matrices matters: the matrix being raised to a power contains probabilities. The familiar printed PAM scoring table contains logarithms of odds.
This also explains how sequences can be more than 100 PAM units apart. A site can undergo several replacements. PAM250 means 250 expected replacements per 100 positions, and the original model gives approximately 20% identity at that distance because repeated changes and reversions hide some of that history. Counting different letters at the endpoints only reveals the differences that remain.
Lower-numbered matrices, such as PAM40, describe relatively little change. Higher-numbered matrices, such as PAM250, are useful when searching for more distant protein relationships.
!PAM250 amino-acid scores, divided by three to display log odds in bits.
PAM250 amino-acid scores, divided by three to display log odds in bits.
The PAM250 plot puts the amino acids on both axes. Each cell scores one pair. Dividing its integer entries by three expresses them in bits, which lets us compare their colors with the BLOSUM62 plot below.
What PAM assumes
The original observations cover a small collection of closely related protein families. Extrapolating them over longer distances assumes that the same replacement process continues to apply. An error in a rare replacement estimate can therefore affect many later calculations.
Protein positions also have different constraints. A residue buried in the core, one exposed to solvent and one involved in catalysis can tolerate different changes, even when they start with the same amino acid. A single general matrix averages across those environments. Rate differences between proteins and changes in constraints over time add further variation.
These limitations [explain](<postgresql.org>) the interest in larger datasets and alternative models. They also matter when interpreting an alignment that extends beyond a credible homologous region: the matrix and gap costs may keep adding weakly favorable columns after the biological correspondence has ended.
BLOSUM: count replacements in conserved blocks
Steven and Jorja Henikoff introduced BLOSUM, short for BLOcks SUbstitution Matrix, in 1992. Their starting point was a collection of conserved, ungapped protein alignment blocks, drawn from hundreds of protein families. The [original paper](<pmc.ncbi.nlm.nih.gov>) describes about 2,000 blocks representing more than 500 groups.
A block is a short region for which the correspondence between residues is sufficiently reliable to count pairs column by column. Conserved motifs help locate these regions. Selecting ungapped columns also keeps the substitution estimate separate from the problem of deciding where insertions and deletions belong.
The family can contain both close and distant relatives. That makes it possible to observe replacement patterns over a broad range of divergence directly, instead of deriving every distance from a PAM1-style transition model.
First reduce the weight of near duplicates
Imagine a protein family represented by 90 nearly identical mammalian sequences and a handful of more distant sequences. Counting every pair equally would let those near duplicates dominate the result.
BLOSUM groups sufficiently similar sequence segments into clusters and weights each cluster as one sequence when counting pairs between clusters. For BLOSUM80, the threshold is 80% identity within the block. BLOSUM62 uses 62%. Sequences at or above the threshold are grouped together, and their contributions are averaged.
The number therefore describes how redundancy was handled during construction. A BLOSUM62 alignment can have much less than 62% identity. It also follows that lowering the threshold groups more similar sequences together and shifts the comparison toward more distant relationships.
For the human cohort later in this tutorial, we will retain the sample counts because we want to describe variation within that cohort. The resulting log-odds matrix uses the same observed-versus-expected calculation. Its sampling rule and evolutionary scale will be visible alongside the results.
Count the pairs in one column
The [Columbia course presentation](<cs.columbia.edu>) walks through this construction with small examples. We can express the unweighted counting step directly in terms of residue counts, which is also convenient in SQL.
Take one alignment column containing A A A A I L. There are four alanines, one isoleucine and one leucine. Every distinct pair of sequences contributes one pair of letters.
The number of A/A pairs is 4 × 3 / 2 = 6. There are four A/I pairs, four A/L pairs and one I/L pair. The total is 15, exactly the number of ways to choose two of the six sequences. Grouping the letters first gives the counts directly:
WITH residues(letter) AS (
VALUES ('A'), ('A'), ('A'), ('A'), ('I'), ('L')
), counts AS (
SELECT letter, count(*) AS n FROM residues GROUP BY letter
), pairs AS (
SELECT a.letter AS first_residue, b.letter AS second_residue,
CASE WHEN a.letter = b.letter THEN a.n * (a.n - 1) / 2
ELSE a.n * b.n END AS observations
FROM counts a JOIN counts b ON a.letter <= b.letter
)
SELECT * FROM pairs WHERE observations > 0
ORDER BY first_residue, second_residue;
Result · 4 rows
first_residue / text | second_residue / text | observations / bigint
---------------------+-----------------------+----------------------
A | A | 6
A | I | 4
A | L | 4
I | L | 1
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
!Six letters in one column produce 15 unordered sequence pairs, grouped by letter combination.
Six letters in one column produce 15 unordered sequence pairs, grouped by letter combination.
For residue a occurring n_(a,r) times in column r, the general formulas are
!c^{(r)}_{aa}=\binom{n_{a,r}}{2} =\frac{n_{a,r}(n_{a,r}-1)}{2}, \qquad c^{(r)}_{ab}=n_{a,r}n_{b,r}\quad(a<b).
Add the counts across all usable columns and normalize:
!c_{ab}=\sum_r c^{(r)}_{ab},\qquad T=\sum_r\binom{n_r}{2},\qquad q_{ab}=\frac{c_{ab}}{T},
where nᵣ is the number of usable residues in column r. If every column contains the same n sequences and there are w columns, this becomes T=wn(n-1)/2. The column-specific version lets us exclude unresolved bases without pretending that every position has the same coverage.
Keep ordered and unordered pairs straight
We counted A/B and B/A together. For these unordered frequencies, the residue background and chance expectations are
!p_a=q_{aa}+\frac12\sum_{b\ne a}q_{ab},
!e_{aa}=p_a^2,\qquad e_{ab}=2p_ap_b\quad(a\ne b).
The factor of two accounts for the two possible orders of an unlike pair. With this convention, the score calculation is
!s_{ab}=\log_2\frac{q_{ab}}{e_{ab}},\qquad S_{ab}=\operatorname{round}(2s_{ab}).
The result agrees with the earlier ordered-pair equation: an off-diagonal ordered frequency is half its unordered counterpart. Keeping that factor consistent prevents a systematic error in every mismatch score.
A pair that never occurs despite a positive expected frequency has an unsmoothed score of negative infinity. If a residue itself is absent, its background frequency is zero and the corresponding odds are undefined. Small datasets can have both kinds of cells. We will preserve their counts and mark unobserved cells in the plots, rather than displaying them as a neutral score.
What the BLOSUM family gives us
BLOSUM80 emphasizes closer relationships, while BLOSUM45 is commonly used for more distant ones. Each matrix is estimated from observations under its clustering threshold. BLOSUM62 became a widely used default for protein searches, including BLASTP.
The larger and more varied training collection was an advantage over the original PAM dataset, and the Henikoffs reported improved alignment and search results in their comparisons. Performance still depends on the sequences and the scoring scheme used with the matrix.
The blocks emphasize conserved regions. They provide reliable alignment columns, but they also select a particular part of protein sequence space. A matrix learned from those regions averages over many functions and structural environments. Short domains can present another difficulty: a permissive matrix aimed at distant relationships may accumulate too little information over a short alignment, or extend a good domain match into neighboring unrelated sequence.
!BLOSUM62 amino-acid scores, divided by two to display log odds in bits.
BLOSUM62 amino-acid scores, divided by two to display log odds in bits.
The diagonal rewards identical residues. Off the diagonal, pairs such as L/I retain a positive score because they occur more frequently than their background abundances predict. Both reference matrix plots use the same bit scale.
Comparing PAM and BLOSUM
<table><tr><th>Property</th><th>PAM</th><th>BLOSUM</th></tr><tr><td>Starting observations</td><td>Closely related protein families in the original dataset</td><td>Conserved ungapped blocks from a broader set of families</td></tr><tr><td>Alignment context</td><td>Global alignments and reconstructed evolutionary changes</td><td>Local blocks and weighted residue-pair counts</td></tr><tr><td>Construction</td><td>An explicit transition model</td><td>Direct observed pair frequencies</td></tr><tr><td>Other family members</td><td>Extrapolate transition probabilities from PAM1</td><td>Repeat estimation with another clustering threshold</td></tr><tr><td>Closer relationships</td><td>Lower numbers, such as PAM40</td><td>Higher numbers, such as BLOSUM80</td></tr><tr><td>More distant relationships</td><td>Higher numbers, such as PAM250</td><td>Lower numbers, such as BLOSUM45</td></tr></table>
One way to compare matrices is their average information per aligned pair, also called relative entropy:
!H=\sum_{a\le b}q_{ab}\log_2\frac{q_{ab}}{e_{ab}}.
Approximate comparisons based on information content put PAM120 near BLOSUM80, PAM160 near BLOSUM60-62, and PAM250 near BLOSUM45. These pairings describe a similar overall scale of discrimination. Individual entries and the alignments they favor can differ.
Score units need attention too. Standard BLOSUM62 entries are half bits, while the standard PAM250 table is in thirds of a bit. Divide by two or three respectively when placing their values on a common color scale.
Choose for the sequences and the alignment length
A shallow matrix targets relatively similar sequences. Low PAM numbers and high BLOSUM numbers belong toward this end of their respective families. A deep matrix allows more divergence, as with PAM250 or BLOSUM45. For protein searches, shallow models can target roughly 50-90% identity, while deep models can target similarities around 20-30%, depending on the matrix and gap settings.
If we expect close orthologs, a shallow matrix can focus the search on that range. Orthologs are genes whose relationship follows a speciation event. If we are looking for distant relatives across very different organisms, a deeper matrix can retain useful evidence from conservative replacements.
Length changes the choice. A short domain has fewer positions over which to accumulate a score. [Pearson's discussion of matrix selection](<pmc.ncbi.nlm.nih.gov>) explains why shallow matrices can work better for short comparisons and how gap penalties affect the useful range of a matrix.
An initial alignment and its identity can help establish the scale of the problem, especially when close relatives are expected. Inspect the aligned region as well: a high overall identity does not reveal whether an endpoint has extended into the wrong domain. BLOSUM62 is a practical starting point for many protein searches, with the matrix changed when the expected divergence or region length gives a reason to do so.
Specialized matrices offer another option when the sequence environment differs from the mixture represented in general protein datasets. JTTtm was developed for transmembrane segments. [PHAT](<academic.oup.com>), short for predicted hydrophobic and transmembrane, uses observations from hydrophobic and transmembrane regions. Their relevance comes from the different residue composition and structural constraints of membrane-spanning sequence.
Gap costs change what the matrix can do
Suppose an alignment can avoid several unfavorable substitutions by inserting a gap. Whether that helps its total score depends on both the substitution matrix and the gap cost. Low costs make it cheap to add gaps. High costs can keep unlike residues opposite one another even when an insertion or deletion is a reasonable explanation.
A constant penalty charges the same amount for every gap run:
!G_{\mathrm{constant}}(L)=c.
A linear penalty charges for every missing residue:
!G_{\mathrm{linear}}(L)=eL.
That is easy to compute, but it gives a five-residue deletion the same total cost as five separate one-residue gaps. An affine model distinguishes opening a gap from extending it:
!G_{\mathrm{affine}}(L)=o+eL.
Here o is the opening component and e the extension cost. With o=5 and e=1, a one-residue gap costs 6 and a five-residue gap costs 10. Five separate one-residue gaps cost 30. This favors a contiguous indel over several scattered gap openings, reflecting the possibility that a whole stretch was inserted or deleted in one event.
Some software writes the same family of functions as o'+e(L-1). Its opening parameter includes the first extension cost, so o'=o+e makes the conventions agree. Check which convention an implementation uses before transferring parameter values. BioseqDB's [alignment scoring definition](<bioseqdb.com>) uses o+eL.
For long indels, a two-piece affine model can use a cheaper extension slope once a gap grows beyond a crossover point:
!G_{\mathrm{dual}}(L)=\min(o_1+e_1L,\ o_2+e_2L).
!Constant, linear, affine and two-piece affine curves charge differently for gaps of increasing length.
Constant, linear, affine and two-piece affine curves charge differently for gaps of increasing length.
The second affine branch can charge more to open a gap but less for each additional base. It becomes the cheaper branch for long gaps.
Affine scoring remains a common choice because [dynamic programming](<xlinux.nist.gov>) can evaluate it with a small set of states per matrix cell. Two-piece affine costs and other models extend that idea when the alignment problem calls for different treatment of long gaps.
Where these scores appear in seed-and-extend aligners
Modern nucleotide aligners frequently begin with exact matches. In [my explanation on Bioinformatics Stack Exchange](<bioinformatics.stackexchange.com>), I separated the process into finding short matching pieces, joining them into chains and optionally calculating a detailed alignment. It helps to keep the scores at those stages separate.
First, an aligner extracts seeds, such as [k-mers](<sourmash.readthedocs.io>) or selected combinations of short substrings, and finds matching seeds in the reference. Each matching seed provides an anchor with coordinates in both sequences. The [Python read-mapper tutorial](<bioseqdb.com>) works through this construction with code.
[Chaining](<lh3/minimap2>) then finds anchors that can plausibly belong to the same alignment. It rewards added matching coverage and penalizes awkward jumps between anchor coordinates. If the query advances 100 bases while the reference advances 120, that difference suggests a gap that the chaining score should account for. The [chaining chapter](<bioseqdb.com>) describes these scores and the algorithms used to evaluate them.
!Seeds provide exact anchors, chaining joins their coordinates, and base alignment resolves the intervening letters.
Seeds provide exact anchors, chaining joins their coordinates, and base alignment resolves the intervening letters.
Anchor chaining can locate a corresponding region before a base-by-base alignment has been calculated.
[Minimap2](<lh3/minimap2>) uses a concave gap cost in its chaining calculation. Its long-gap refinements can use the [minigraph chaining approach with range-minimum-query data structures](<bioseqdb.com>), which helps reconnect anchors across long indels. This is a geometric calculation on anchors. The two-piece affine cost described above belongs to its later base-level alignment.
At the end of chaining, an aligner already has approximate mapping coordinates. Minimap2 can output those without calculating a [CIGAR string](<samtools.github.io>). Requesting base-level alignment makes it resolve the letters and gaps between anchors, using its [KSW2](<lh3/ksw2>) dynamic-programming implementation. The [minimap2 manual](<lh3.github.io>) distinguishes the chaining and alignment controls, and the [base-level alignment chapter](<bioseqdb.com>) covers the corresponding BioseqDB implementations.
Several optimizations reduce that work. Banding restricts which part of the dynamic-programming grid is evaluated. X-drop and Z-drop stop an extension when its score falls sufficiently far below a previous best. Z-drop allows a larger drop when the path moves farther from that best-scoring cell's diagonal, which helps it continue across a gap. [Block Aligner](<pmc.ncbi.nlm.nih.gov>) evaluates selected regions of the grid as the alignment advances. Its [source code](<Daniel-Liu-c0deb0t/block-aligner>) shows how those blocks shift, grow and shrink.
These nucleotide alignments use nucleotide scoring rules. Protein searches have a different alphabet and can use BLOSUM during both candidate selection and detailed alignment. For example, protein word-neighborhood seeding can accept a high-scoring word that differs from the query. The [MMseqs2 user guide](<mmseqs.com>) describes its protein-search stages. The [seed-and-extend](<pubmed.ncbi.nlm.nih.gov>) structure tells us how the search is organized, while the alphabet and scoring choices tell us what each stage considers similar.
Count substitutions in 1,000 human mitochondrial genomes
We can carry the same calculation through a real dataset. I chose 1,000 complete human mitochondrial genomes from the Danish cohort published by [Li and colleagues](<pmc.ncbi.nlm.nih.gov>). Mitochondria have their own small genome, about 16,600 bases long, which makes the download manageable while giving us enough samples to see repeated differences at individual positions.
The study deposited 2,000 mitochondrial sequences. We take records in accession order and retain the first 1,000 whose [GenBank](<ncbi.nlm.nih.gov>) title explicitly says "complete genome". This excludes 11 partial records and gives a list beginning at KF161060.1 and ending at KF162070.1. Each has a different isolate identifier. The examples repository includes that list and its metadata, so the sample selection stays the same when somebody reruns the tutorial.
These are mitochondrial genomes from one Danish study. Their frequencies describe that cohort, including its sampling and shared ancestry. A common variant carried by 400 people contributes 400 observations relative to the reference, even if those people inherited it from the same ancestor. That distinction will matter when we interpret the matrix.
Our reference is the [revised Cambridge Reference Sequence, NC_012920.1](<ncbi.nlm.nih.gov>). We will align each sample to it, walk the CIGAR operations to recover corresponding positions, then count the letters in each column. The reference's gene annotations also tell us which groups of three bases form codons in its 13 protein-coding genes.
<table><tr><th>Input</th><th>Role in the analysis</th></tr><tr><td>1,000 complete mitochondrial sample records</td><td>The observations used to estimate the cohort frequencies</td></tr><tr><td><code>NC_012920.1</code></td><td>A coordinate system for alignment and coding-region annotations</td></tr><tr><td>13 annotated coding regions</td><td>Corresponding codons for the amino-acid comparison</td></tr><tr><td>NCBI BLOSUM62 joint frequencies and scores</td><td>The protein reference matrix</td></tr><tr><td>NCBI PAM250 scores</td><td>A second historical matrix to inspect on the same bit scale</td></tr></table>
The DNA calculation uses the four letters A, C, G and T. For the comparison with BLOSUM62 we translate the coding regions and use the 20 amino-acid alphabet. In BLOSUM, A means alanine, C cysteine, G glycine and T threonine. Keeping those two alphabets in separate tables makes the comparisons much easier to read.
Install BioseqDB and get the example files
Follow the [installation guide](<bioseqdb.com>) to create a [PostgreSQL](<postgresql.org>) database with BioseqDB. This tutorial uses its NCBI integration to download sequences through SQL.
Warning
Availability on 16 October 2026. The published 0.1.0-alpha.1 [Docker](<docs.docker.com>) images and native package builds omit the external-database integration used here (due to my oversight). It is planned for the next alpha release, but you can obviously compile source from Github directly. Currently using runtime SET, an API key, or a GRANT won't work as the code isn't present in the binary - when building one must use BIOSEQ_BUILD_EXTERNAL_API_CLIENTS=ON to enable this feature.
The [source-build instructions in the butterfly example](<github.com/BioseqDB/examples/tree/main/butterfies…>) use the same external module. Check that it is present before starting the download:
SELECT to_regnamespace('bioseq_external') IS NOT NULL AS external_schema,
to_regclass('bioseq_external.sequences') IS NOT NULL AS remote_sequences;
Result · 1 row
external_schema / boolean | remote_sequences / boolean
--------------------------+---------------------------
t | t
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
Both columns should be true. Your database role also needs USAGE on [bioseq_external](<bioseqdb.com>) and SELECT on [bioseq_external.sequences](<bioseqdb.com>), which the database owner can grant.
Get the SQL files, accession metadata and plotting code from [BioseqDB/examples/scoring_matrices](<github.com/BioseqDB/examples/tree/main/scoring_ma…>), or download the [complete example bundle](<bioseqdb.com>):
git clone https://github.com/BioseqDB/examples.git
cd examples/scoring_matrices
export DATABASE_URL=bioseq_examples
Set DATABASE_URL to your database name. For a remote server, set PGHOST, PGPORT and PGUSER for your connection. Use a fresh database or one without a scoring schema. The pipeline keeps its intermediate tables in that schema so we can inspect each step.
Download and inspect the sequences
Our table keeps the accession, its role and the sequence itself:
CREATE SCHEMA scoring;
CREATE TABLE scoring.genomes (
accession text NOT NULL,
role text NOT NULL,
seq bioseq.seq('dna5', 64) NOT NULL
) USING bioseq;
CREATE SCHEMA
CREATE TABLE
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
USING [bioseq](<bioseqdb.com>) selects BioseqDB's table storage. [bioseq.seq](<bioseqdb.com>)('dna5', 64) stores DNA with the extra symbol N for an unresolved base. The 64 selects 64-bit packed storage units, and a sequence can occupy many of them. The [sequence representation chapter](<bioseqdb.com>) goes through the details.
Here are the requests for the reference and one sample:
INSERT INTO scoring.genomes (accession, role, seq)
SELECT 'NC_012920.1', 'reference', sequence::bioseq.seq('dna5', 64)
FROM bioseq_external.sequences
WHERE provider = 'ncbi' AND accession = 'NC_012920.1';
INSERT INTO scoring.genomes (accession, role, seq)
SELECT 'KF161060.1', 'sample', sequence::bioseq.seq('dna5', 64)
FROM bioseq_external.sequences
WHERE provider = 'ncbi' AND accession = 'KF161060.1';
SELECT accession, role, bioseq.len(seq) AS bases
FROM scoring.genomes ORDER BY role, accession;
INSERT 0 1 — 1 row affected
INSERT 0 1 — 1 row affected
accession / text | role / text | bases / bigint
-----------------+-------------+---------------
NC_012920.1 | reference | 16569
KF161060.1 | sample | 16570
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
The external table reads the NCBI response, and INSERT ... SELECT saves the sequence locally. The complete import repeats that request for the accession list and the reference. It pauses between requests and skips accessions already stored, so an interrupted download can resume.
Run the full SQL workflow with:
bash run.sh
The script downloads the records, checks the sample count, aligns the genomes and writes the result tables to results/. All the sequence processing runs in PostgreSQL. We will use Python after those CSV files exist, to draw the figures.
SELECT role, count(*) AS genomes,
min(bioseq.len(seq)) AS shortest,
max(bioseq.len(seq)) AS longest,
sum(bioseq.len(seq)) AS total_bases,
sum(length(seq::text) - length(replace(seq::text, 'N', ''))) AS unknown_n
FROM scoring.genomes GROUP BY role ORDER BY role;
SELECT accession, bioseq.len(seq) AS bases,
length(seq::text) - length(replace(seq::text, 'N', '')) AS unknown_n
FROM scoring.genomes WHERE role = 'sample'
ORDER BY accession LIMIT 5;
Result 1
Result rows (long values continue on indented lines):
Row 1
role / text: reference
genomes / bigint: 1
shortest / bigint: 16569
longest / bigint: 16569
total_bases / numeric: 16569
unknown_n / bigint: 1
Row 2
role / text: sample
genomes / bigint: 1000
shortest / bigint: 16561
longest / bigint: 16575
total_bases / numeric: 16568713
unknown_n / bigint: 8691
Result 2
accession / text | bases / bigint | unknown_n / integer
-----------------+----------------+--------------------
KF161060.1 | 16570 | 0
KF161061.1 | 16570 | 144
KF161062.1 | 16569 | 0
KF161063.1 | 16569 | 1
KF161065.1 | 16568 | 18
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
The website evaluates the SQL listings against a saved copy of these same input records. The downloadable workflow fetches them from NCBI. This keeps the displayed results tied to the complete cohort while letting us rebuild the page without requesting the same thousand sequences every time.
Measure N before excluding it
N records uncertainty about a nucleotide. Treating it as a fifth biological base would turn missing observations into apparent substitutions, so we count it separately and use A/C/G/T as the callable alphabet.
SELECT count(*) AS sample_genomes,
sum(bioseq.len(seq)) AS total_bases,
sum(length(seq::text) - length(replace(seq::text, 'N', ''))) AS unknown_n,
round(100.0 * sum(length(seq::text) - length(replace(seq::text, 'N', '')))
/ sum(bioseq.len(seq)), 4) AS percent_unknown
FROM scoring.genomes WHERE role = 'sample';
Result · 1 row
sample_genomes / bigint | total_bases / numeric | unknown_n / bigint | percent_unknown / numeric
------------------------+-----------------------+--------------------+--------------------------
1000 | 16568713 | 8691 | 0.0525
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
The 1,000 samples contain 16,568,713 bases in total. Their 8,691 Ns account for 0.0525% of those bases. The shortest record has 16,561 bases and the longest has 16,575.
!Sequence length and number of unresolved N bases in each of the 1,000 sample records.
Sequence length and number of unresolved N bases in each of the 1,000 sample records.
Each point is one accession. Sequence lengths help us spot a truncated input, while the second panel shows how much information each record will contribute after N is excluded. The reference also contains a historical N at position 3107. We keep that coordinate during alignment and exclude it when comparing nucleotide identities.
Align corresponding positions
An insertion changes the relationship between the two coordinate systems. If a sample has an extra base after position 309, comparing the raw strings at position 310, 311 and so on would shift every later comparison. We need the alignment path.
The query below asks BioseqDB for alignments and reads the coordinates and CIGAR attached to each match:
SELECT q.accession,
bioseq.target_start(t.seq)::integer AS target_start,
bioseq.target_end(t.seq)::integer AS target_end,
bioseq.query_start(t.seq)::integer AS query_start,
bioseq.query_end(t.seq)::integer AS query_end,
bioseq.is_reverse(t.seq) AS is_reverse,
bioseq.cigar(t.seq)::text AS cigar,
bioseq.edit(t.seq) AS edits,
bioseq.score(t.seq) AS score,
bioseq.is_primary(t.seq) AS is_primary
FROM scoring.genomes t
JOIN scoring.genomes q ON mmsearch(t.seq, q.seq, aligner_options(
kmer_size => 15, window_size => 10,
max_gap => 200, min_chain_score => 40))
WHERE t.role = 'reference' AND q.role = 'sample' AND q.accession = 'KF161060.1';
Result · 1 row
Result rows (long values continue on indented lines):
Row 1
accession / text: KF161060.1
target_start / integer: 0
target_end / integer: 16569
query_start / integer: 0
query_end / integer: 16570
is_reverse / boolean: f
cigar / text: 309M1I6M1I2791M1D13462M
edits / bigint: 17
score / integer: 33034
is_primary / boolean: t
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
[mmsearch](<bioseqdb.com>) finds a match between the reference sequence and a sample. The [aligner_options](<bioseqdb.com>) values set the [minimizer](<pubmed.ncbi.nlm.nih.gov>) and chaining parameters for these closely related genomes. Here the seeds are 15-base k-mers, with minimizers selected from windows of 10 consecutive k-mers. We request max_gap = 200. The default 500-base bandwidth raises the effective chaining distance to 500 bases, while alignment end extensions remain capped at 200. A minimum chain score of 40 filters weak candidates. [bioseq.target_start](<bioseqdb.com>) and [bioseq.target_end](<bioseqdb.com>) locate the reference interval, while the query accessors locate the corresponding sample interval. These intervals are zero-based and half-open.
[bioseq.cigar](<bioseqdb.com>) describes how the two intervals line up. An M consumes a base from both sequences and can contain either a match or a mismatch. An I consumes a sample base without advancing the reference, and a D advances the reference without consuming a sample base. The letters themselves determine which M positions differ.
The pipeline retains the best forward primary alignment for each accession and checks that it consumes the full sample and covers at least 99.9% of the reference. A failed check stops the run. In this cohort 999 alignments cover the whole reference. One leaves seven terminal reference positions uncovered, which we record as gaps. The [butterfly alignment tutorial](<bioseqdb.com>) explains these coordinate fields with a genomic ribbon plot, and the [final-alignment chapter](<bioseqdb.com>) describes how BioseqDB obtains the CIGAR.
Turn the CIGAR into aligned columns
sql/05_columns.sql splits each CIGAR into numbered runs. Two running sums track how far we have moved in the reference and sample. For a run consuming both sequences we copy the sample substring, for a deletion we write a gap, and for an insertion we retain its position in the CIGAR table but leave it out of the reference-column string.
This gives every sample a string with the reference's length. Character 500 now means the sample base aligned to reference position 500. We can group those strings by position and base:
SELECT p.pos::integer AS position1, p.base, count(*) AS samples
FROM scoring.projected_sequences s
CROSS JOIN LATERAL regexp_split_to_table(s.aligned, '')
WITH ORDINALITY AS p(base, pos)
GROUP BY p.pos, p.base
ORDER BY p.pos, p.base LIMIT 12;
SELECT base, sum(n) AS aligned_calls,
base IN ('A', 'C', 'G', 'T') AS included_in_matrix
FROM scoring.site_counts
GROUP BY base ORDER BY base;
Result 1
position1 / integer | base / text | samples / bigint
--------------------+-------------+-----------------
1 | - | 1
1 | G | 999
2 | - | 1
2 | A | 999
3 | - | 1
3 | T | 999
4 | C | 1000
5 | A | 1000
6 | C | 1000
7 | A | 1000
8 | G | 1000
9 | G | 1000
Result 2
base / text | aligned_calls / numeric | included_in_matrix / boolean
------------+-------------------------+-----------------------------
- | 1250 | f
A | 5114723 | t
C | 5179937 | t
G | 2174485 | t
N | 8616 | f
T | 4089989 | t
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
An inserted sample base has no reference column, so we report insertions separately. Of the 8,691 raw Ns, 75 occur in inserted sequence and 8,616 remain in the projected columns. The 1,250 gap calls consist of 1,243 deleted bases and the seven uncovered terminal reference positions. This reference-anchored arrangement makes the count at each position easy to inspect.
!Numbers of sample differences, N calls and gaps along the mitochondrial reference.
Numbers of sample differences, N calls and gaps along the mitochondrial reference.
The upper panel counts samples whose callable base differs from the reference. A tall line means many samples carry another base at that coordinate. The lower panel separates N calls from gaps, which helps distinguish a well-observed variant from a region with missing information.
Build the nucleotide matrix in SQL
At a column containing four As, one C and one G, there are six A/A pairs, four A/C pairs, four A/G pairs and one C/G pair. We can compute those counts from the letter totals without joining every sample to every other sample.
For the SQL implementation we count both orientations of each pair. Thus A/C and C/A each receive four observations, and A/A receives 4(4-1)=12. This ordered convention gives us a symmetric 4 × 4 table and the simple expected frequency p_a p_b.
SELECT a.base AS base_from,b.base AS base_to,
sum(CASE WHEN a.base=b.base THEN a.n*(a.n-1)
ELSE a.n*b.n END)::bigint AS pairs
FROM scoring.site_counts a JOIN scoring.site_counts b USING(position1)
JOIN scoring.reference_bases r USING(position1)
WHERE r.base IN ('A','C','G','T') AND a.base IN ('A','C','G','T') AND b.base IN ('A','C','G','T')
GROUP BY a.base,b.base
ORDER BY base_from, base_to;
Result · 16 rows
base_from / text | base_to / text | pairs / bigint
-----------------+----------------+---------------
A | A | 5099615260
A | C | 360392
A | G | 6862012
A | T | 59484
C | A | 360392
C | C | 5164850040
C | G | 58468
C | T | 6880251
G | A | 6862012
G | C | 58468
G | G | 2164099862
G | T | 47450
T | A | 59484
T | C | 6880251
T | G | 47450
T | T | 4076677940
Executed against Bioseq extension 1.0 on PostgreSQL 17.10
For a diagonal cell we use nₐ(nₐ − 1), which removes a sample's comparison with itself. For different letters we use n_a n_b. Summing those counts across co
Mentions
@fooliu @nilshomer (5) minibwa: a 2-3x improvement. @lh3lh3 has decided not to move forward with it, but it's there for anyone to use if they wish. minibwa#79
@PacBio Some of the optimizations went upstream too. The agents made hifiasm 2.9x and minimap2 1.5x faster on a whole genome, with identical output. Thanks to Haoyu Cheng and @lh3lh3 for building these incredible tools. hifiasm#956 minimap2#1364
PAM, BLOSUM and scoring matrices from 1,000 human mitochondrial genomes Two sequences can agree at most positions and still leave us with a few awkward decisions. Is a leucine opposite an isoleucine a reasonable match? Should an [alignment](<en.wikipedia.org>) insert a gap here, or accept three mismatches and keep going? A scoring scheme puts numbers on those choices, and changing the numbers can change the alignment you get. PAM and BLOSUM are families of substitution matrices for [proteins](<genome.gov>). They describe how amino-acid pairs occur in aligned related sequences, relative to the pairs we would expect from their individual abundance. We will work through those counts, then use BioseqDB to calculate [nucleotide](<genome.gov>) and protein comparisons from 1,000 human mitochondrial [genome](<genome.gov>) records. [SQL](<postgresql.org>) will do the counting and matrix calculations. Python will draw the results. A protein is a chain of [amino acids](<genome.gov>). A coding region of [DNA](<genome.gov>) supplies its instructions in three-letter groups called [codons](<genome.gov>), each specifying an amino acid or a stop signal. Mitochondria carry a small genome of their own, separate from the nuclear [chromosomes](<genome.gov>). That gives us complete mitochondrial records small enough to inspect and align in a tutorial, along with protein-coding [genes](<genome.gov>) that connect the DNA analysis to BLOSUM. Keep the two alphabets in mind: the A in a DNA sequence is adenine, while the A in BLOSUM is alanine, an amino acid. What a scoring matrix does An alignment places two sequences into corresponding columns. A column can contain the same letter twice, two different letters, or a letter opposite a gap. Those gaps let us represent insertions and deletions, usually shortened to indels. The scoring system has two parts. A substitution matrix assigns a score to each aligned letter pair, and a gap model charges for inserting or extending gaps. If we call the alignment 𝒜, its score is !S(\mathcal A)= \sum_{(a,b)\in\text{paired columns of }\mathcal A}s(a,b) -\sum_{g\in\text{gap runs of }\mathcal A}G(|g|). Here s(a,b) is the substitution score and G(|g|) is the cost of a gap run of length |g|. The algorithm searches for a high-scoring arrangement under that rule. A global alignment consumes both full sequences, whereas a local alignment can select matching regions within them. The [alignment foundations chapter](<bioseqdb.com>) gives the formal definitions and shows the different boundary conditions. Why use biological scores at all? Related proteins can accumulate substitutions while retaining a similar structure and function. A scoring matrix lets an alignment retain evidence from those changes instead of relying entirely on identical letters. We can then assess the resulting score against a chance model, which helps identify similarities consistent with common ancestry. The significance also depends on sequence length, composition and the size of the database being searched. Start with an identity matrix The smallest useful example rewards a match with +1 and charges -1 for a mismatch: !s(a,b)= \begin{cases} +1,&a=b,\ -1,&a\ne b. \end{cases} Some versions give mismatches zero instead. Either way, all mismatches receive the same treatment. For DNA this is a useful starting point, and simple match/mismatch schemes appear in many nucleotide aligners. It can also score proteins, although we lose the distinction between chemically similar and very different amino acids. Leucine and isoleucine, written L and I, have similar hydrophobic side chains, which tend to avoid water. Replacing one with the other can preserve the local environment of a protein. Glycine and tryptophan, G and W, are a different pairing: glycine is small and flexible, while tryptophan has a bulky aromatic side chain. Their exchange can demand a much larger structural adjustment. We call the first kind of replacement conservative and the second radical. The surrounding protein still matters. A substitution that is tolerated at one position can disrupt another, especially in a binding site or tightly packed core. A general substitution matrix averages observations across many positions and protein families. !An identity matrix gives all mismatches the same score, while BLOSUM62 distinguishes amino-acid pairs. An identity matrix gives all mismatches the same score, while BLOSUM62 distinguishes amino-acid pairs. The letters in the protein examples are amino-acid codes. BLOSUM62 rewards L/I with +2, while G/W receives -2. From observed pairs to log odds Suppose Q_ab is the probability of an ordered pair (a,b) in our aligned related sequences. Let pₐ and p_b be the background frequencies of those two residues. Under an independent chance model, the expected pair frequency is p_a p_b. The odds ratio compares the two: !O_{ab}=\frac{Q_{ab}}{p_ap_b}. An odds ratio above one means the pair appears more frequently in the related-sequence alignments than the background frequencies predict. Below one means it appears less frequently. With a model that combines evidence from individual columns, taking a logarithm turns multiplication into addition: !s_{ab}=\log_2\left(\frac{Q_{ab}}{p_ap_b}\right), \qquad \log_2\prod_r O_{a_rb_r}=\sum_r\log_2 O_{a_rb_r}. This is convenient for alignment algorithms: they can add a score as they advance through each paired column. Positive entries favor an alignment, negative entries reduce its score, and zero means that the observed and expected frequencies agree. BLOSUM62 stores integer scores in half-bit units. Before rounding, its scoring equation has the form used in the [ForensicMCQ BLOSUM guide](<forensicmcq.com>): !S_{ij}=2\log_2\left(\frac{\text{observed frequency}}{\text{expected frequency}}\right). For a hypothetical pair observed with frequency 0.04 and expected with frequency 0.01, the calculation is !2\log_2(0.04/0.01)=2\log_2(4)=2\times2=+4. Those numbers illustrate the arithmetic. Actual BLOSUM62 entries come from its measured pair frequencies. For example, the aligned protein fragments ALK and AID score !S(\mathrm{ALK},\mathrm{AID})= S(\mathrm A,\mathrm A)+S(\mathrm L,\mathrm I)+S(\mathrm K,\mathrm D) =4+2+(-1)=5. There are no gaps in that example, so we can simply add the three matrix entries. Longer alignments use the same operation repeatedly. Here are both calculations in SQL: SELECT round(2 * log(2::numeric, 0.04 / 0.01)) AS half_bit_score; SELECT sum(score) AS alignment_score FROM (VALUES ('A', 'A', 4), ('L', 'I', 2), ('K', 'D', -1)) AS aligned(first_residue, second_residue, score); Result 1 half_bit_score / numeric ------------------------ 4 Result 2 alignment_score / bigint ------------------------ 5 Executed against Bioseq extension 1.0 on PostgreSQL 17.10 PAM: estimating change over evolutionary distance Margaret Dayhoff and her colleagues developed the PAM approach during the 1960s and 1970s. Their original data came from 71 groups of closely related proteins, with alignments above roughly 85% identity. At that distance, ancestral reconstruction and counting changes along a tree are easier than they would be for highly divergent sequences. PAM stands for Point Accepted [Mutation](<genome.gov>). An accepted amino-acid replacement is one retained in an evolving lineage. It may preserve protein function well enough to persist, without giving the organism any particular advantage. One PAM unit represents an evolutionary interval with an expected one accepted replacement per 100 amino-acid positions. It is a measure of accumulated change. The amount of elapsed time needed to accumulate it depends on the proteins and lineages under study. The construction starts with global alignments of related proteins. From these alignments, a phylogenetic tree and likely ancestral sequences let us infer replacements along branches. Counts of those replacements, together with residue abundance, provide estimates of relative mutability, mₐ: how readily each amino acid changes in the fitted model. The probabilities are then normalized to the small amount of change represented by one PAM unit. This gives a PAM1 transition matrix, M(1). An entry describes the probability of moving from one amino acid to another over that unit of evolutionary distance. The diagonal contains the probabilities of retaining the same residue. To obtain a model for a longer distance, multiply the transition matrix by itself: !M(n)=M(1)^n. For PAM250, that is M(1)²⁵⁰. Each multiplication sums over possible intermediate residues, so a site can change several times and can return to its starting state. After calculating the transition probabilities at the chosen distance, we convert them into log-odds alignment scores using the background frequencies. The distinction between these two matrices matters: the matrix being raised to a power contains probabilities. The familiar printed PAM scoring table contains logarithms of odds. This also explains how sequences can be more than 100 PAM units apart. A site can undergo several replacements. PAM250 means 250 expected replacements per 100 positions, and the original model gives approximately 20% identity at that distance because repeated changes and reversions hide some of that history. Counting different letters at the endpoints only reveals the differences that remain. Lower-numbered matrices, such as PAM40, describe relatively little change. Higher-numbered matrices, such as PAM250, are useful when searching for more distant protein relationships. !PAM250 amino-acid scores, divided by three to display log odds in bits. PAM250 amino-acid scores, divided by three to display log odds in bits. The PAM250 plot puts the amino acids on both axes. Each cell scores one pair. Dividing its integer entries by three expresses them in bits, which lets us compare their colors with the BLOSUM62 plot below. What PAM assumes The original observations cover a small collection of closely related protein families. Extrapolating them over longer distances assumes that the same replacement process continues to apply. An error in a rare replacement estimate can therefore affect many later calculations. Protein positions also have different constraints. A residue buried in the core, one exposed to solvent and one involved in catalysis can tolerate different changes, even when they start with the same amino acid. A single general matrix averages across those environments. Rate differences between proteins and changes in constraints over time add further variation. These limitations [explain](<postgresql.org>) the interest in larger datasets and alternative models. They also matter when interpreting an alignment that extends beyond a credible homologous region: the matrix and gap costs may keep adding weakly favorable columns after the biological correspondence has ended. BLOSUM: count replacements in conserved blocks Steven and Jorja Henikoff introduced BLOSUM, short for BLOcks SUbstitution Matrix, in 1992. Their starting point was a collection of conserved, ungapped protein alignment blocks, drawn from hundreds of protein families. The [original paper](<pmc.ncbi.nlm.nih.gov>) describes about 2,000 blocks representing more than 500 groups. A block is a short region for which the correspondence between residues is sufficiently reliable to count pairs column by column. Conserved motifs help locate these regions. Selecting ungapped columns also keeps the substitution estimate separate from the problem of deciding where insertions and deletions belong. The family can contain both close and distant relatives. That makes it possible to observe replacement patterns over a broad range of divergence directly, instead of deriving every distance from a PAM1-style transition model. First reduce the weight of near duplicates Imagine a protein family represented by 90 nearly identical mammalian sequences and a handful of more distant sequences. Counting every pair equally would let those near duplicates dominate the result. BLOSUM groups sufficiently similar sequence segments into clusters and weights each cluster as one sequence when counting pairs between clusters. For BLOSUM80, the threshold is 80% identity within the block. BLOSUM62 uses 62%. Sequences at or above the threshold are grouped together, and their contributions are averaged. The number therefore describes how redundancy was handled during construction. A BLOSUM62 alignment can have much less than 62% identity. It also follows that lowering the threshold groups more similar sequences together and shifts the comparison toward more distant relationships. For the human cohort later in this tutorial, we will retain the sample counts because we want to describe variation within that cohort. The resulting log-odds matrix uses the same observed-versus-expected calculation. Its sampling rule and evolutionary scale will be visible alongside the results. Count the pairs in one column The [Columbia course presentation](<cs.columbia.edu>) walks through this construction with small examples. We can express the unweighted counting step directly in terms of residue counts, which is also convenient in SQL. Take one alignment column containing A A A A I L. There are four alanines, one isoleucine and one leucine. Every distinct pair of sequences contributes one pair of letters. The number of A/A pairs is 4 × 3 / 2 = 6. There are four A/I pairs, four A/L pairs and one I/L pair. The total is 15, exactly the number of ways to choose two of the six sequences. Grouping the letters first gives the counts directly: WITH residues(letter) AS ( VALUES ('A'), ('A'), ('A'), ('A'), ('I'), ('L') ), counts AS ( SELECT letter, count(*) AS n FROM residues GROUP BY letter ), pairs AS ( SELECT a.letter AS first_residue, b.letter AS second_residue, CASE WHEN a.letter = b.letter THEN a.n * (a.n - 1) / 2 ELSE a.n * b.n END AS observations FROM counts a JOIN counts b ON a.letter <= b.letter ) SELECT * FROM pairs WHERE observations > 0 ORDER BY first_residue, second_residue; Result · 4 rows first_residue / text | second_residue / text | observations / bigint ---------------------+-----------------------+---------------------- A | A | 6 A | I | 4 A | L | 4 I | L | 1 Executed against Bioseq extension 1.0 on PostgreSQL 17.10 !Six letters in one column produce 15 unordered sequence pairs, grouped by letter combination. Six letters in one column produce 15 unordered sequence pairs, grouped by letter combination. For residue a occurring n_(a,r) times in column r, the general formulas are !c^{(r)}_{aa}=\binom{n_{a,r}}{2} =\frac{n_{a,r}(n_{a,r}-1)}{2}, \qquad c^{(r)}_{ab}=n_{a,r}n_{b,r}\quad(a<b). Add the counts across all usable columns and normalize: !c_{ab}=\sum_r c^{(r)}_{ab},\qquad T=\sum_r\binom{n_r}{2},\qquad q_{ab}=\frac{c_{ab}}{T}, where nᵣ is the number of usable residues in column r. If every column contains the same n sequences and there are w columns, this becomes T=wn(n-1)/2. The column-specific version lets us exclude unresolved bases without pretending that every position has the same coverage. Keep ordered and unordered pairs straight We counted A/B and B/A together. For these unordered frequencies, the residue background and chance expectations are !p_a=q_{aa}+\frac12\sum_{b\ne a}q_{ab}, !e_{aa}=p_a^2,\qquad e_{ab}=2p_ap_b\quad(a\ne b). The factor of two accounts for the two possible orders of an unlike pair. With this convention, the score calculation is !s_{ab}=\log_2\frac{q_{ab}}{e_{ab}},\qquad S_{ab}=\operatorname{round}(2s_{ab}). The result agrees with the earlier ordered-pair equation: an off-diagonal ordered frequency is half its unordered counterpart. Keeping that factor consistent prevents a systematic error in every mismatch score. A pair that never occurs despite a positive expected frequency has an unsmoothed score of negative infinity. If a residue itself is absent, its background frequency is zero and the corresponding odds are undefined. Small datasets can have both kinds of cells. We will preserve their counts and mark unobserved cells in the plots, rather than displaying them as a neutral score. What the BLOSUM family gives us BLOSUM80 emphasizes closer relationships, while BLOSUM45 is commonly used for more distant ones. Each matrix is estimated from observations under its clustering threshold. BLOSUM62 became a widely used default for protein searches, including BLASTP. The larger and more varied training collection was an advantage over the original PAM dataset, and the Henikoffs reported improved alignment and search results in their comparisons. Performance still depends on the sequences and the scoring scheme used with the matrix. The blocks emphasize conserved regions. They provide reliable alignment columns, but they also select a particular part of protein sequence space. A matrix learned from those regions averages over many functions and structural environments. Short domains can present another difficulty: a permissive matrix aimed at distant relationships may accumulate too little information over a short alignment, or extend a good domain match into neighboring unrelated sequence. !BLOSUM62 amino-acid scores, divided by two to display log odds in bits. BLOSUM62 amino-acid scores, divided by two to display log odds in bits. The diagonal rewards identical residues. Off the diagonal, pairs such as L/I retain a positive score because they occur more frequently than their background abundances predict. Both reference matrix plots use the same bit scale. Comparing PAM and BLOSUM <table><tr><th>Property</th><th>PAM</th><th>BLOSUM</th></tr><tr><td>Starting observations</td><td>Closely related protein families in the original dataset</td><td>Conserved ungapped blocks from a broader set of families</td></tr><tr><td>Alignment context</td><td>Global alignments and reconstructed evolutionary changes</td><td>Local blocks and weighted residue-pair counts</td></tr><tr><td>Construction</td><td>An explicit transition model</td><td>Direct observed pair frequencies</td></tr><tr><td>Other family members</td><td>Extrapolate transition probabilities from PAM1</td><td>Repeat estimation with another clustering threshold</td></tr><tr><td>Closer relationships</td><td>Lower numbers, such as PAM40</td><td>Higher numbers, such as BLOSUM80</td></tr><tr><td>More distant relationships</td><td>Higher numbers, such as PAM250</td><td>Lower numbers, such as BLOSUM45</td></tr></table> One way to compare matrices is their average information per aligned pair, also called relative entropy: !H=\sum_{a\le b}q_{ab}\log_2\frac{q_{ab}}{e_{ab}}. Approximate comparisons based on information content put PAM120 near BLOSUM80, PAM160 near BLOSUM60-62, and PAM250 near BLOSUM45. These pairings describe a similar overall scale of discrimination. Individual entries and the alignments they favor can differ. Score units need attention too. Standard BLOSUM62 entries are half bits, while the standard PAM250 table is in thirds of a bit. Divide by two or three respectively when placing their values on a common color scale. Choose for the sequences and the alignment length A shallow matrix targets relatively similar sequences. Low PAM numbers and high BLOSUM numbers belong toward this end of their respective families. A deep matrix allows more divergence, as with PAM250 or BLOSUM45. For protein searches, shallow models can target roughly 50-90% identity, while deep models can target similarities around 20-30%, depending on the matrix and gap settings. If we expect close orthologs, a shallow matrix can focus the search on that range. Orthologs are genes whose relationship follows a speciation event. If we are looking for distant relatives across very different organisms, a deeper matrix can retain useful evidence from conservative replacements. Length changes the choice. A short domain has fewer positions over which to accumulate a score. [Pearson's discussion of matrix selection](<pmc.ncbi.nlm.nih.gov>) explains why shallow matrices can work better for short comparisons and how gap penalties affect the useful range of a matrix. An initial alignment and its identity can help establish the scale of the problem, especially when close relatives are expected. Inspect the aligned region as well: a high overall identity does not reveal whether an endpoint has extended into the wrong domain. BLOSUM62 is a practical starting point for many protein searches, with the matrix changed when the expected divergence or region length gives a reason to do so. Specialized matrices offer another option when the sequence environment differs from the mixture represented in general protein datasets. JTTtm was developed for transmembrane segments. [PHAT](<academic.oup.com>), short for predicted hydrophobic and transmembrane, uses observations from hydrophobic and transmembrane regions. Their relevance comes from the different residue composition and structural constraints of membrane-spanning sequence. Gap costs change what the matrix can do Suppose an alignment can avoid several unfavorable substitutions by inserting a gap. Whether that helps its total score depends on both the substitution matrix and the gap cost. Low costs make it cheap to add gaps. High costs can keep unlike residues opposite one another even when an insertion or deletion is a reasonable explanation. A constant penalty charges the same amount for every gap run: !G_{\mathrm{constant}}(L)=c. A linear penalty charges for every missing residue: !G_{\mathrm{linear}}(L)=eL. That is easy to compute, but it gives a five-residue deletion the same total cost as five separate one-residue gaps. An affine model distinguishes opening a gap from extending it: !G_{\mathrm{affine}}(L)=o+eL. Here o is the opening component and e the extension cost. With o=5 and e=1, a one-residue gap costs 6 and a five-residue gap costs 10. Five separate one-residue gaps cost 30. This favors a contiguous indel over several scattered gap openings, reflecting the possibility that a whole stretch was inserted or deleted in one event. Some software writes the same family of functions as o'+e(L-1). Its opening parameter includes the first extension cost, so o'=o+e makes the conventions agree. Check which convention an implementation uses before transferring parameter values. BioseqDB's [alignment scoring definition](<bioseqdb.com>) uses o+eL. For long indels, a two-piece affine model can use a cheaper extension slope once a gap grows beyond a crossover point: !G_{\mathrm{dual}}(L)=\min(o_1+e_1L,\ o_2+e_2L). !Constant, linear, affine and two-piece affine curves charge differently for gaps of increasing length. Constant, linear, affine and two-piece affine curves charge differently for gaps of increasing length. The second affine branch can charge more to open a gap but less for each additional base. It becomes the cheaper branch for long gaps. Affine scoring remains a common choice because [dynamic programming](<xlinux.nist.gov>) can evaluate it with a small set of states per matrix cell. Two-piece affine costs and other models extend that idea when the alignment problem calls for different treatment of long gaps. Where these scores appear in seed-and-extend aligners Modern nucleotide aligners frequently begin with exact matches. In [my explanation on Bioinformatics Stack Exchange](<bioinformatics.stackexchange.com>), I separated the process into finding short matching pieces, joining them into chains and optionally calculating a detailed alignment. It helps to keep the scores at those stages separate. First, an aligner extracts seeds, such as [k-mers](<sourmash.readthedocs.io>) or selected combinations of short substrings, and finds matching seeds in the reference. Each matching seed provides an anchor with coordinates in both sequences. The [Python read-mapper tutorial](<bioseqdb.com>) works through this construction with code. [Chaining](<lh3/minimap2>) then finds anchors that can plausibly belong to the same alignment. It rewards added matching coverage and penalizes awkward jumps between anchor coordinates. If the query advances 100 bases while the reference advances 120, that difference suggests a gap that the chaining score should account for. The [chaining chapter](<bioseqdb.com>) describes these scores and the algorithms used to evaluate them. !Seeds provide exact anchors, chaining joins their coordinates, and base alignment resolves the intervening letters. Seeds provide exact anchors, chaining joins their coordinates, and base alignment resolves the intervening letters. Anchor chaining can locate a corresponding region before a base-by-base alignment has been calculated. [Minimap2](<lh3/minimap2>) uses a concave gap cost in its chaining calculation. Its long-gap refinements can use the [minigraph chaining approach with range-minimum-query data structures](<bioseqdb.com>), which helps reconnect anchors across long indels. This is a geometric calculation on anchors. The two-piece affine cost described above belongs to its later base-level alignment. At the end of chaining, an aligner already has approximate mapping coordinates. Minimap2 can output those without calculating a [CIGAR string](<samtools.github.io>). Requesting base-level alignment makes it resolve the letters and gaps between anchors, using its [KSW2](<lh3/ksw2>) dynamic-programming implementation. The [minimap2 manual](<lh3.github.io>) distinguishes the chaining and alignment controls, and the [base-level alignment chapter](<bioseqdb.com>) covers the corresponding BioseqDB implementations. Several optimizations reduce that work. Banding restricts which part of the dynamic-programming grid is evaluated. X-drop and Z-drop stop an extension when its score falls sufficiently far below a previous best. Z-drop allows a larger drop when the path moves farther from that best-scoring cell's diagonal, which helps it continue across a gap. [Block Aligner](<pmc.ncbi.nlm.nih.gov>) evaluates selected regions of the grid as the alignment advances. Its [source code](<Daniel-Liu-c0deb0t/block-aligner>) shows how those blocks shift, grow and shrink. These nucleotide alignments use nucleotide scoring rules. Protein searches have a different alphabet and can use BLOSUM during both candidate selection and detailed alignment. For example, protein word-neighborhood seeding can accept a high-scoring word that differs from the query. The [MMseqs2 user guide](<mmseqs.com>) describes its protein-search stages. The [seed-and-extend](<pubmed.ncbi.nlm.nih.gov>) structure tells us how the search is organized, while the alphabet and scoring choices tell us what each stage considers similar. Count substitutions in 1,000 human mitochondrial genomes We can carry the same calculation through a real dataset. I chose 1,000 complete human mitochondrial genomes from the Danish cohort published by [Li and colleagues](<pmc.ncbi.nlm.nih.gov>). Mitochondria have their own small genome, about 16,600 bases long, which makes the download manageable while giving us enough samples to see repeated differences at individual positions. The study deposited 2,000 mitochondrial sequences. We take records in accession order and retain the first 1,000 whose [GenBank](<ncbi.nlm.nih.gov>) title explicitly says "complete genome". This excludes 11 partial records and gives a list beginning at KF161060.1 and ending at KF162070.1. Each has a different isolate identifier. The examples repository includes that list and its metadata, so the sample selection stays the same when somebody reruns the tutorial. These are mitochondrial genomes from one Danish study. Their frequencies describe that cohort, including its sampling and shared ancestry. A common variant carried by 400 people contributes 400 observations relative to the reference, even if those people inherited it from the same ancestor. That distinction will matter when we interpret the matrix. Our reference is the [revised Cambridge Reference Sequence, NC_012920.1](<ncbi.nlm.nih.gov>). We will align each sample to it, walk the CIGAR operations to recover corresponding positions, then count the letters in each column. The reference's gene annotations also tell us which groups of three bases form codons in its 13 protein-coding genes. <table><tr><th>Input</th><th>Role in the analysis</th></tr><tr><td>1,000 complete mitochondrial sample records</td><td>The observations used to estimate the cohort frequencies</td></tr><tr><td><code>NC_012920.1</code></td><td>A coordinate system for alignment and coding-region annotations</td></tr><tr><td>13 annotated coding regions</td><td>Corresponding codons for the amino-acid comparison</td></tr><tr><td>NCBI BLOSUM62 joint frequencies and scores</td><td>The protein reference matrix</td></tr><tr><td>NCBI PAM250 scores</td><td>A second historical matrix to inspect on the same bit scale</td></tr></table> The DNA calculation uses the four letters A, C, G and T. For the comparison with BLOSUM62 we translate the coding regions and use the 20 amino-acid alphabet. In BLOSUM, A means alanine, C cysteine, G glycine and T threonine. Keeping those two alphabets in separate tables makes the comparisons much easier to read. Install BioseqDB and get the example files Follow the [installation guide](<bioseqdb.com>) to create a [PostgreSQL](<postgresql.org>) database with BioseqDB. This tutorial uses its NCBI integration to download sequences through SQL. Warning Availability on 16 October 2026. The published 0.1.0-alpha.1 [Docker](<docs.docker.com>) images and native package builds omit the external-database integration used here (due to my oversight). It is planned for the next alpha release, but you can obviously compile source from Github directly. Currently using runtime SET, an API key, or a GRANT won't work as the code isn't present in the binary - when building one must use BIOSEQ_BUILD_EXTERNAL_API_CLIENTS=ON to enable this feature. The [source-build instructions in the butterfly example](<github.com/BioseqDB/examples/tree/main/butterfies…>) use the same external module. Check that it is present before starting the download: SELECT to_regnamespace('bioseq_external') IS NOT NULL AS external_schema, to_regclass('bioseq_external.sequences') IS NOT NULL AS remote_sequences; Result · 1 row external_schema / boolean | remote_sequences / boolean --------------------------+--------------------------- t | t Executed against Bioseq extension 1.0 on PostgreSQL 17.10 Both columns should be true. Your database role also needs USAGE on [bioseq_external](<bioseqdb.com>) and SELECT on [bioseq_external.sequences](<bioseqdb.com>), which the database owner can grant. Get the SQL files, accession metadata and plotting code from [BioseqDB/examples/scoring_matrices](<github.com/BioseqDB/examples/tree/main/scoring_ma…>), or download the [complete example bundle](<bioseqdb.com>): git clone https://github.com/BioseqDB/examples.git cd examples/scoring_matrices export DATABASE_URL=bioseq_examples Set DATABASE_URL to your database name. For a remote server, set PGHOST, PGPORT and PGUSER for your connection. Use a fresh database or one without a scoring schema. The pipeline keeps its intermediate tables in that schema so we can inspect each step. Download and inspect the sequences Our table keeps the accession, its role and the sequence itself: CREATE SCHEMA scoring; CREATE TABLE scoring.genomes ( accession text NOT NULL, role text NOT NULL, seq bioseq.seq('dna5', 64) NOT NULL ) USING bioseq; CREATE SCHEMA CREATE TABLE Executed against Bioseq extension 1.0 on PostgreSQL 17.10 USING [bioseq](<bioseqdb.com>) selects BioseqDB's table storage. [bioseq.seq](<bioseqdb.com>)('dna5', 64) stores DNA with the extra symbol N for an unresolved base. The 64 selects 64-bit packed storage units, and a sequence can occupy many of them. The [sequence representation chapter](<bioseqdb.com>) goes through the details. Here are the requests for the reference and one sample: INSERT INTO scoring.genomes (accession, role, seq) SELECT 'NC_012920.1', 'reference', sequence::bioseq.seq('dna5', 64) FROM bioseq_external.sequences WHERE provider = 'ncbi' AND accession = 'NC_012920.1'; INSERT INTO scoring.genomes (accession, role, seq) SELECT 'KF161060.1', 'sample', sequence::bioseq.seq('dna5', 64) FROM bioseq_external.sequences WHERE provider = 'ncbi' AND accession = 'KF161060.1'; SELECT accession, role, bioseq.len(seq) AS bases FROM scoring.genomes ORDER BY role, accession; INSERT 0 1 — 1 row affected INSERT 0 1 — 1 row affected accession / text | role / text | bases / bigint -----------------+-------------+--------------- NC_012920.1 | reference | 16569 KF161060.1 | sample | 16570 Executed against Bioseq extension 1.0 on PostgreSQL 17.10 The external table reads the NCBI response, and INSERT ... SELECT saves the sequence locally. The complete import repeats that request for the accession list and the reference. It pauses between requests and skips accessions already stored, so an interrupted download can resume. Run the full SQL workflow with: bash run.sh The script downloads the records, checks the sample count, aligns the genomes and writes the result tables to results/. All the sequence processing runs in PostgreSQL. We will use Python after those CSV files exist, to draw the figures. SELECT role, count(*) AS genomes, min(bioseq.len(seq)) AS shortest, max(bioseq.len(seq)) AS longest, sum(bioseq.len(seq)) AS total_bases, sum(length(seq::text) - length(replace(seq::text, 'N', ''))) AS unknown_n FROM scoring.genomes GROUP BY role ORDER BY role; SELECT accession, bioseq.len(seq) AS bases, length(seq::text) - length(replace(seq::text, 'N', '')) AS unknown_n FROM scoring.genomes WHERE role = 'sample' ORDER BY accession LIMIT 5; Result 1 Result rows (long values continue on indented lines): Row 1 role / text: reference genomes / bigint: 1 shortest / bigint: 16569 longest / bigint: 16569 total_bases / numeric: 16569 unknown_n / bigint: 1 Row 2 role / text: sample genomes / bigint: 1000 shortest / bigint: 16561 longest / bigint: 16575 total_bases / numeric: 16568713 unknown_n / bigint: 8691 Result 2 accession / text | bases / bigint | unknown_n / integer -----------------+----------------+-------------------- KF161060.1 | 16570 | 0 KF161061.1 | 16570 | 144 KF161062.1 | 16569 | 0 KF161063.1 | 16569 | 1 KF161065.1 | 16568 | 18 Executed against Bioseq extension 1.0 on PostgreSQL 17.10 The website evaluates the SQL listings against a saved copy of these same input records. The downloadable workflow fetches them from NCBI. This keeps the displayed results tied to the complete cohort while letting us rebuild the page without requesting the same thousand sequences every time. Measure N before excluding it N records uncertainty about a nucleotide. Treating it as a fifth biological base would turn missing observations into apparent substitutions, so we count it separately and use A/C/G/T as the callable alphabet. SELECT count(*) AS sample_genomes, sum(bioseq.len(seq)) AS total_bases, sum(length(seq::text) - length(replace(seq::text, 'N', ''))) AS unknown_n, round(100.0 * sum(length(seq::text) - length(replace(seq::text, 'N', ''))) / sum(bioseq.len(seq)), 4) AS percent_unknown FROM scoring.genomes WHERE role = 'sample'; Result · 1 row sample_genomes / bigint | total_bases / numeric | unknown_n / bigint | percent_unknown / numeric ------------------------+-----------------------+--------------------+-------------------------- 1000 | 16568713 | 8691 | 0.0525 Executed against Bioseq extension 1.0 on PostgreSQL 17.10 The 1,000 samples contain 16,568,713 bases in total. Their 8,691 Ns account for 0.0525% of those bases. The shortest record has 16,561 bases and the longest has 16,575. !Sequence length and number of unresolved N bases in each of the 1,000 sample records. Sequence length and number of unresolved N bases in each of the 1,000 sample records. Each point is one accession. Sequence lengths help us spot a truncated input, while the second panel shows how much information each record will contribute after N is excluded. The reference also contains a historical N at position 3107. We keep that coordinate during alignment and exclude it when comparing nucleotide identities. Align corresponding positions An insertion changes the relationship between the two coordinate systems. If a sample has an extra base after position 309, comparing the raw strings at position 310, 311 and so on would shift every later comparison. We need the alignment path. The query below asks BioseqDB for alignments and reads the coordinates and CIGAR attached to each match: SELECT q.accession, bioseq.target_start(t.seq)::integer AS target_start, bioseq.target_end(t.seq)::integer AS target_end, bioseq.query_start(t.seq)::integer AS query_start, bioseq.query_end(t.seq)::integer AS query_end, bioseq.is_reverse(t.seq) AS is_reverse, bioseq.cigar(t.seq)::text AS cigar, bioseq.edit(t.seq) AS edits, bioseq.score(t.seq) AS score, bioseq.is_primary(t.seq) AS is_primary FROM scoring.genomes t JOIN scoring.genomes q ON mmsearch(t.seq, q.seq, aligner_options( kmer_size => 15, window_size => 10, max_gap => 200, min_chain_score => 40)) WHERE t.role = 'reference' AND q.role = 'sample' AND q.accession = 'KF161060.1'; Result · 1 row Result rows (long values continue on indented lines): Row 1 accession / text: KF161060.1 target_start / integer: 0 target_end / integer: 16569 query_start / integer: 0 query_end / integer: 16570 is_reverse / boolean: f cigar / text: 309M1I6M1I2791M1D13462M edits / bigint: 17 score / integer: 33034 is_primary / boolean: t Executed against Bioseq extension 1.0 on PostgreSQL 17.10 [mmsearch](<bioseqdb.com>) finds a match between the reference sequence and a sample. The [aligner_options](<bioseqdb.com>) values set the [minimizer](<pubmed.ncbi.nlm.nih.gov>) and chaining parameters for these closely related genomes. Here the seeds are 15-base k-mers, with minimizers selected from windows of 10 consecutive k-mers. We request max_gap = 200. The default 500-base bandwidth raises the effective chaining distance to 500 bases, while alignment end extensions remain capped at 200. A minimum chain score of 40 filters weak candidates. [bioseq.target_start](<bioseqdb.com>) and [bioseq.target_end](<bioseqdb.com>) locate the reference interval, while the query accessors locate the corresponding sample interval. These intervals are zero-based and half-open. [bioseq.cigar](<bioseqdb.com>) describes how the two intervals line up. An M consumes a base from both sequences and can contain either a match or a mismatch. An I consumes a sample base without advancing the reference, and a D advances the reference without consuming a sample base. The letters themselves determine which M positions differ. The pipeline retains the best forward primary alignment for each accession and checks that it consumes the full sample and covers at least 99.9% of the reference. A failed check stops the run. In this cohort 999 alignments cover the whole reference. One leaves seven terminal reference positions uncovered, which we record as gaps. The [butterfly alignment tutorial](<bioseqdb.com>) explains these coordinate fields with a genomic ribbon plot, and the [final-alignment chapter](<bioseqdb.com>) describes how BioseqDB obtains the CIGAR. Turn the CIGAR into aligned columns sql/05_columns.sql splits each CIGAR into numbered runs. Two running sums track how far we have moved in the reference and sample. For a run consuming both sequences we copy the sample substring, for a deletion we write a gap, and for an insertion we retain its position in the CIGAR table but leave it out of the reference-column string. This gives every sample a string with the reference's length. Character 500 now means the sample base aligned to reference position 500. We can group those strings by position and base: SELECT p.pos::integer AS position1, p.base, count(*) AS samples FROM scoring.projected_sequences s CROSS JOIN LATERAL regexp_split_to_table(s.aligned, '') WITH ORDINALITY AS p(base, pos) GROUP BY p.pos, p.base ORDER BY p.pos, p.base LIMIT 12; SELECT base, sum(n) AS aligned_calls, base IN ('A', 'C', 'G', 'T') AS included_in_matrix FROM scoring.site_counts GROUP BY base ORDER BY base; Result 1 position1 / integer | base / text | samples / bigint --------------------+-------------+----------------- 1 | - | 1 1 | G | 999 2 | - | 1 2 | A | 999 3 | - | 1 3 | T | 999 4 | C | 1000 5 | A | 1000 6 | C | 1000 7 | A | 1000 8 | G | 1000 9 | G | 1000 Result 2 base / text | aligned_calls / numeric | included_in_matrix / boolean ------------+-------------------------+----------------------------- - | 1250 | f A | 5114723 | t C | 5179937 | t G | 2174485 | t N | 8616 | f T | 4089989 | t Executed against Bioseq extension 1.0 on PostgreSQL 17.10 An inserted sample base has no reference column, so we report insertions separately. Of the 8,691 raw Ns, 75 occur in inserted sequence and 8,616 remain in the projected columns. The 1,250 gap calls consist of 1,243 deleted bases and the seven uncovered terminal reference positions. This reference-anchored arrangement makes the count at each position easy to inspect. !Numbers of sample differences, N calls and gaps along the mitochondrial reference. Numbers of sample differences, N calls and gaps along the mitochondrial reference. The upper panel counts samples whose callable base differs from the reference. A tall line means many samples carry another base at that coordinate. The lower panel separates N calls from gaps, which helps distinguish a well-observed variant from a region with missing information. Build the nucleotide matrix in SQL At a column containing four As, one C and one G, there are six A/A pairs, four A/C pairs, four A/G pairs and one C/G pair. We can compute those counts from the letter totals without joining every sample to every other sample. For the SQL implementation we count both orientations of each pair. Thus A/C and C/A each receive four observations, and A/A receives 4(4-1)=12. This ordered convention gives us a symmetric 4 × 4 table and the simple expected frequency p_a p_b. SELECT a.base AS base_from,b.base AS base_to, sum(CASE WHEN a.base=b.base THEN a.n*(a.n-1) ELSE a.n*b.n END)::bigint AS pairs FROM scoring.site_counts a JOIN scoring.site_counts b USING(position1) JOIN scoring.reference_bases r USING(position1) WHERE r.base IN ('A','C','G','T') AND a.base IN ('A','C','G','T') AND b.base IN ('A','C','G','T') GROUP BY a.base,b.base ORDER BY base_from, base_to; Result · 16 rows base_from / text | base_to / text | pairs / bigint -----------------+----------------+--------------- A | A | 5099615260 A | C | 360392 A | G | 6862012 A | T | 59484 C | A | 360392 C | C | 5164850040 C | G | 58468 C | T | 6880251 G | A | 6862012 G | C | 58468 G | G | 2164099862 G | T | 47450 T | A | 59484 T | C | 6880251 T | G | 47450 T | T | 4076677940 Executed against Bioseq extension 1.0 on PostgreSQL 17.10 For a diagonal cell we use nₐ(nₐ − 1), which removes a sample's comparison with itself. For different letters we use n_a n_b. Summing those counts across co
Introducing BioseqDB BioseqDB is a PostgreSQL extension for storing, indexing and searching biological sequences with SQL. Combine genomic matches with relational data using compact sequence types and columnar storage! I wanted to keep sequences and their metadata together. If a sample has an organism, collection date or patient identifier, those values should be ordinary columns beside its DNA. Searching for a similar sequence should fit into the same query as filtering those columns. !Load FASTA or FASTQ and metadata into a BioseqDB table, then combine metadata filters and sequence alignment in one SQL query. How it started This project required several years of exploration and experiments. I began without a formal background in computational biology, thus learning the field was part of the research. - 2021–2022: Working with large GISAID datasets made FASTA files hard to manage. I proposed the database idea in May 2021 and directed a team prototype that connected BWA-MEM to PostgreSQL in May 2022. - 2023–2024: I wrote a Python aligner and then improved it, learning how the algorithms work along the way. - December 2025–2026: The current BioseqDB began with minimap2 and columnar storage. Custom seeding and chaining, indexes, sequence slices and end-to-end benchmarks followed. The current BioseqDB is a separate project, designed and written by me. It contains no code or intellectual property from the early BWA-MEM integration experiments. The introduction and full timeline have more of the story. Two queries to get a feel for it These examples use the three synthetic samples from the SQL walkthrough. The samples table has ordinary metadata columns and a bioseq.seq column. USING bioseq selects the sequence-aware storage. First, read metadata and a small piece of DNA. bioseq.len gives the sequence length; [:16] takes its first sixteen bases. SELECT patient_name, bioseq.len(seq) AS bases, seq[:16]::text AS prefix FROM samples WHERE age >= 30 ORDER BY age DESC; Result · 2 rows | patient_name | bases | prefix | | --- | --- | --- | | Ben | 512 | ACAAAAAGGGGCTTTC | | Anna | 4096 | AAAGAGTAGATCGTAG | Now search those samples for Ben's sequence, stored in row 2. mmsearch performs the alignment; bioseq.score and bioseq.cigar describe each match. SELECT t.patient_name, bioseq.target_start(t.seq) AS start, bioseq.score(t.seq) AS score, bioseq.cigar(t.seq) AS cigar FROM samples AS t, samples AS q WHERE t.age >= 30 AND q.id = 2 AND mmsearch(t.seq, q.seq::bioseq.seq('dna4', 64)) ORDER BY t.id, start; Result · 3 rows | patient_name | start | score | cigar | | --- | --- | --- | --- | | Anna | 256 | 1024 | 512M | | Anna | 2304 | 1024 | 512M | | Ben | 0 | 1024 | 512M | Anna has two occurrences and Ben has one. These are exact matches, so each CIGAR is 512M. Alignment can also handle differences between sequences: substitutions, insertions and deletions. The alignment chapter explains how seeding, chaining and base alignment fit together. Where it fits beside existing tools BioseqDB builds on ideas from existing aligners, especially minimap2. The difference I wanted was the workflow: keep the data in PostgreSQL and combine sequence search with SQL filters and joins. | Tool | Main focus | |:---|:---| | minimap2 | Mapping reads and aligning nucleotide sequences against a reference. | | BLAST | Finding local sequence similarity in sequence databases. | | Winnowmap | Long-read mapping, including repetitive regions. | | LexicMap | Nucleotide alignment against large collections of prokaryotic genomes. | | BioseqDB | Storing sequences beside relational metadata and running alignment inside SQL. | In my benchmarks, the preferred CPU configuration was typically 3–8× faster than minimap2 when a prepared reference index was available. Without that index, BioseqDB was usually slower than minimap2. These tests used 1,024 synthetic queries against one 2 or 4 Gbp reference, not a collection of real genomes. Several external tools also needed the reference split into smaller records, so the measurement details matter when comparing them. There is still work to do. The current GPU paths were slower than the CPU paths in these measurements, and indexing has known limitations. This first release is experimental software. Try the Docker image The image already includes PostgreSQL and BioseqDB. This is the Ubuntu 26.04 / PostgreSQL 19 prerelease image from the quick start: docker run -d --name bioseqdb-demo \ -e POSTGRES_PASSWORD=demo-change-me \ -p 127.0.0.1:55432:5432 \ styczynski/bioseqdb:0.1.0-alpha.1-pg19-ubuntu26.04 Wait until PostgreSQL is ready, then open its SQL prompt: until docker exec bioseqdb-demo pg_isready -h 127.0.0.1 -U postgres; do sleep 1; done docker exec -it -u postgres bioseqdb-demo psql -U postgres -d postgres Follow the walkthrough to create the sample table, load a FASTA file and try the queries above. I would like to hear how it works with your data: open an issue or contact me. --- Please check my OS project - BioseqDB. Read this post on BioseqDB · Read on Medium · Follow BioseqDB on LinkedIn.