DEV Community: BioseqDB

PAM, BLOSUM and scoring matrices from 1,000 human mitochondrial genomes

BioseqDB — Thu, 22 Oct 2026 17:14:54 +0000

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 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. 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 and protein comparisons from 1,000 human mitochondrial genome records. SQL will do the counting and matrix calculations. Python will draw the results.

A protein is a chain of amino acids. A coding region of DNA supplies its instructions in three-letter groups called codons, each specifying an amino acid or a stop signal.

Mitochondria carry a small genome of their own, separate from the nuclear chromosomes. That gives us complete mitochondrial records small enough to inspect and align in a tutorial, along with protein-coding genes 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

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 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:

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.

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:

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:

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:

For a hypothetical pair observed with frequency 0.04 and expected with frequency 0.01, the calculation is

Those numbers illustrate the arithmetic. Actual BLOSUM62 entries come from its measured pair frequencies. For example, the aligned protein fragments ALK and AID score

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. 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:

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.

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 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 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 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.

For residue a occurring n_(a,r) times in column r, the general formulas are

Add the counts across all usable columns and normalize:

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

The factor of two accounts for the two possible orders of an unlike pair. With this convention, the score calculation is

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.

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

Property	PAM	BLOSUM
Starting observations	Closely related protein families in the original dataset	Conserved ungapped blocks from a broader set of families
Alignment context	Global alignments and reconstructed evolutionary changes	Local blocks and weighted residue-pair counts
Construction	An explicit transition model	Direct observed pair frequencies
Other family members	Extrapolate transition probabilities from PAM1	Repeat estimation with another clustering threshold
Closer relationships	Lower numbers, such as PAM40	Higher numbers, such as BLOSUM80
More distant relationships	Higher numbers, such as PAM250	Lower numbers, such as BLOSUM45

One way to compare matrices is their average information per aligned pair, also called relative entropy:

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 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, 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:

A linear penalty charges for every missing residue:

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:

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 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:

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 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, 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 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 works through this construction with code.

Chaining 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 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.

Anchor chaining can locate a corresponding region before a base-by-base alignment has been calculated.

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, 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. Requesting base-level alignment makes it resolve the letters and gaps between anchors, using its KSW2 dynamic-programming implementation. The minimap2 manual distinguishes the chaining and alignment controls, and the base-level alignment chapter 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 evaluates selected regions of the grid as the alignment advances. Its source code 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 describes its protein-search stages. The seed-and-extend 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. 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 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. 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.

Input	Role in the analysis
1,000 complete mitochondrial sample records	The observations used to estimate the cohort frequencies
`NC_012920.1`	A coordinate system for alignment and coding-region annotations
13 annotated coding regions	Corresponding codons for the amino-acid comparison
NCBI BLOSUM62 joint frequencies and scores	The protein reference matrix
NCBI PAM250 scores	A second historical matrix to inspect on the same bit scale

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 to create a PostgreSQL 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 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 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 and SELECT on bioseq_external.sequences, which the database owner can grant.

Get the SQL files, accession metadata and plotting code from BioseqDB/examples/scoring_matrices, or download the complete example bundle:

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

USINGbioseq selects BioseqDB's table storage. bioseq.seq('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 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

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.

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 finds a match between the reference sequence and a sample. The aligner_options values set the minimizer 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 and bioseq.target_end locate the reference interval, while the query accessors locate the corresponding sample interval. These intervals are zero-based and half-open.

bioseq.cigar 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 explains these coordinate fields with a genomic ribbon plot, and the final-alignment chapter 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.

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 columns produces the same totals as enumerating every distinct sample pair, while each column needs at most 16 combinations of letter counts.

With 1,000 samples there are (1,000 choose 2) = 499,500 unordered sample pairs at a fully callable position. That number falls wherever N or a gap removes a sample from the column. The denominator in the query comes from the retained counts themselves.

Now divide each pair count by the total, calculate the marginal frequencies and take the log odds:

WITH frequencies AS (
    SELECT *, pairs::numeric / sum(pairs) OVER () AS q
    FROM scoring.nucleotide_pairs
), marginals AS (
    SELECT base_from AS base, sum(q) AS p
    FROM frequencies GROUP BY base_from
)
SELECT f.base_from, f.base_to, pairs,
       round(q, 6) AS q,
       round(a.p * b.p, 6) AS expected,
       CASE WHEN q > 0 THEN
           round(2 * log(2::numeric, q / (a.p * b.p)), 2)
       END AS score_half_bits
FROM frequencies f
JOIN marginals a ON a.base = f.base_from
JOIN marginals b ON b.base = f.base_to
ORDER BY base_from, base_to;

Result · 16 rows

Result rows (long values continue on indented lines):

Row 1
base_from / text: A
base_to / text: A
pairs / bigint: 5099615260
q / numeric: 0.308436
expected / numeric: 0.095405
score_half_bits / numeric: 3.39

Row 2
base_from / text: A
base_to / text: C
pairs / bigint: 360392
q / numeric: 0.000022
expected / numeric: 0.096624
score_half_bits / numeric: -24.23

Row 3
base_from / text: A
base_to / text: G
pairs / bigint: 6862012
q / numeric: 0.000415
expected / numeric: 0.040559
score_half_bits / numeric: -13.22

Row 4
base_from / text: A
base_to / text: T
pairs / bigint: 59484
q / numeric: 0.000004
expected / numeric: 0.076289
score_half_bits / numeric: -28.74

Row 5
base_from / text: C
base_to / text: A
pairs / bigint: 360392
q / numeric: 0.000022
expected / numeric: 0.096624
score_half_bits / numeric: -24.23

Row 6
base_from / text: C
base_to / text: C
pairs / bigint: 5164850040
q / numeric: 0.312382
expected / numeric: 0.097858
score_half_bits / numeric: 3.35

Row 7
base_from / text: C
base_to / text: G
pairs / bigint: 58468
q / numeric: 0.000004
expected / numeric: 0.041077
score_half_bits / numeric: -27.01

Row 8
base_from / text: C
base_to / text: T
pairs / bigint: 6880251
q / numeric: 0.000416
expected / numeric: 0.077264
score_half_bits / numeric: -15.07

Row 9
base_from / text: G
base_to / text: A
pairs / bigint: 6862012
q / numeric: 0.000415
expected / numeric: 0.040559
score_half_bits / numeric: -13.22

Row 10
base_from / text: G
base_to / text: C
pairs / bigint: 58468
q / numeric: 0.000004
expected / numeric: 0.041077
score_half_bits / numeric: -27.01

Row 11
base_from / text: G
base_to / text: G
pairs / bigint: 2164099862
q / numeric: 0.130890
expected / numeric: 0.017243
score_half_bits / numeric: 5.85

Row 12
base_from / text: G
base_to / text: T
pairs / bigint: 47450
q / numeric: 0.000003
expected / numeric: 0.032432
score_half_bits / numeric: -26.93

Row 13
base_from / text: T
base_to / text: A
pairs / bigint: 59484
q / numeric: 0.000004
expected / numeric: 0.076289
score_half_bits / numeric: -28.74

Row 14
base_from / text: T
base_to / text: C
pairs / bigint: 6880251
q / numeric: 0.000416
expected / numeric: 0.077264
score_half_bits / numeric: -15.07

Row 15
base_from / text: T
base_to / text: G
pairs / bigint: 47450
q / numeric: 0.000003
expected / numeric: 0.032432
score_half_bits / numeric: -26.93

Row 16
base_from / text: T
base_to / text: T
pairs / bigint: 4076677940
q / numeric: 0.246567
expected / numeric: 0.061004
score_half_bits / numeric: 4.03

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

q is the observed ordered pair frequency. The two marginal frequencies give expected = p_a * p_b, and 2 * log(2, q / expected) expresses the score in half-bit units, the convention used by BLOSUM62. The chart divides that value by two so its color scale is in bits.

Four-by-four nucleotide log-odds matrix calculated from the 1,000 mitochondrial genomes.

All 16 nucleotide pairs occur in the data. The diagonal entries are positive, while off-diagonal scores are negative because differences between these closely related samples are much rarer than chance pairing of their base frequencies. Across the reference-aligned columns, 1,565 positions contain more than one observed canonical base.

These samples are closely related, so the diagonal contains most observations. The matrix describes their aligned base-pair frequencies, including unchanged positions. The numbers would change if we selected another population, another genomic region or a more divergent collection.

Separate the destination from the opportunity to differ

There are two useful ways to read an A-to-G count. We can divide by all callable sample observations at reference A positions, which tells us how often an A position differs as G. Or we can divide by all differences from reference A, which tells us how much of that difference spectrum goes to G rather than C or T.

Writing C_(a→b) for the reference-to-sample count, the first fraction is

For a destination among differences, with b ≠ a, we use

Both are useful, provided the denominator is clear. This query reports the counts and both percentages:

SELECT reference_base, sample_base, calls,
       round(100 * fraction_of_calls, 4) AS percent_of_callable_positions,
       round(100 * fraction_of_differences, 2) AS percent_of_differences
FROM scoring.reference_spectrum
WHERE reference_base <> sample_base
ORDER BY reference_base, sample_base;

Result · 12 rows

reference_base / text | sample_base / text | calls / bigint | percent_of_callable_positions / numeric | percent_of_differences / numeric
----------------------+--------------------+----------------+-----------------------------------------+---------------------------------
A                     | C                  | 95             | 0.0019                                  | 0.90
A                     | G                  | 10424          | 0.2035                                  | 98.70
A                     | T                  | 42             | 0.0008                                  | 0.40
C                     | A                  | 314            | 0.0061                                  | 8.43
C                     | G                  | 15             | 0.0003                                  | 0.40
C                     | T                  | 3397           | 0.0656                                  | 91.17
G                     | A                  | 3697           | 0.1705                                  | 98.48
G                     | C                  | 45             | 0.0021                                  | 1.20
G                     | T                  | 12             | 0.0006                                  | 0.32
T                     | A                  | 18             | 0.0004                                  | 0.35
T                     | C                  | 5111           | 0.1249                                  | 98.95
T                     | G                  | 36             | 0.0009                                  | 0.70

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

Four stacked bars showing destinations among differences from each reference nucleotide.

Each bar sums to 100% of the differences from its reference letter. A/G and C/T exchanges are called transitions, while exchanges between those two chemical groups are transversions. The chart lets us see those patterns without the much larger number of unchanged bases filling almost the whole plot.

For reference A, 98.70% of the observed differences are G. The corresponding transition fractions are 91.17% for C to T, 98.48% for G to A and 98.95% for T to C. There are 23,206 differences in 16,559,133 callable reference/sample comparisons, about 0.1401%.

The arrow here means reference base to sample base. Inferring which state was ancestral would require an evolutionary analysis, so the plotted quantities are frequencies of observed differences relative to rCRS.

We can also use the cohort matrix itself, treating each sample as a possible first partner rather than choosing rCRS. The same normalization gives the following proportions:

SELECT base_from, base_to,
       round(100 * fraction_of_differences, 2) AS percent_of_differences
FROM scoring.nucleotide_matrix
WHERE base_from <> base_to
ORDER BY base_from, base_to;

Result · 12 rows

base_from / text | base_to / text | percent_of_differences / numeric
-----------------+----------------+---------------------------------
A                | C              | 4.95
A                | G              | 94.23
A                | T              | 0.82
C                | A              | 4.94
C                | G              | 0.80
C                | T              | 94.26
G                | A              | 98.48
G                | C              | 0.84
G                | T              | 0.68
T                | A              | 0.85
T                | C              | 98.47
T                | G              | 0.68

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

Partner nucleotides among differing pairs across all 1,000 samples.

Here A/G accounts for 94.23% of differing pairs with A, and C/T accounts for 94.26% of those with C. The frequencies change because all sample pairs contribute. A variant that is common in the cohort but differs from rCRS receives different weight in these two summaries.

Translate the coding regions for the BLOSUM comparison

Our nucleotide matrix has four letters. BLOSUM62 has 20 amino acids, so the next table starts from the reference's annotated coding sequences rather than from arbitrary three-base chunks of the genome.

sql/07_proteins.sql stores each gene's interval, strand and annotated reference translation. Twelve coding regions lie on the forward strand. ND6 is on the reverse strand, so its projected bases are reverse-complemented before reading its codons. The translation table is NCBI genetic code 2, the vertebrate mitochondrial code.

That choice changes real amino-acid calls. For example, TGA encodes tryptophan in this code, while ATA encodes methionine. The script also handles the annotated initiation position and uses the reference protein length to leave incomplete terminal stop codons outside the amino-acid table.

Before translating a sample's coding region, we exclude that whole region if its alignment is incomplete or contains an insertion or deletion. A premature stop, an invalid initiation codon or an unresolved initiation codon also removes the region from the protein pair counts. The remaining codons give a conservative set of directly corresponding amino-acid positions. We count ambiguous codons separately, and keep the per-gene exclusions visible:

SELECT c.gene, count(*) AS sample_cds,
       count(*) FILTER (WHERE has_indel) AS indel_cds,
       count(*) FILTER (WHERE incompletely_aligned) AS incomplete_cds,
       count(*) FILTER (WHERE q.premature_stop) AS stop_cds,
       count(*) FILTER (WHERE q.invalid_start) AS invalid_start_cds,
       count(*) FILTER (WHERE q.ambiguous_start) AS ambiguous_start_cds,
       count(*) FILTER (WHERE NOT has_indel AND NOT incompletely_aligned
                       AND NOT q.premature_stop AND NOT q.invalid_start
                       AND NOT q.ambiguous_start) AS retained_cds
FROM scoring.cds_quality c
LEFT JOIN scoring.protein_cds_quality q USING (accession, gene)
GROUP BY c.gene ORDER BY c.gene;

SELECT count(*) AS eligible_codons,
       count(*) FILTER (WHERE aa IS NULL) AS ambiguous_codons,
       count(*) FILTER (WHERE aa IS NOT NULL AND aa <> '*') AS callable_codons,
       count(*) FILTER (WHERE aa = 'N') AS asparagine_codons
FROM scoring.amino_acids a
JOIN scoring.protein_cds_quality q USING (accession, gene)
WHERE NOT q.premature_stop AND NOT q.invalid_start AND NOT q.ambiguous_start;

Result 1

Result rows (long values continue on indented lines):

Row 1
gene / text: ATP6
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 1000

Row 2
gene / text: ATP8
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 1000

Row 3
gene / text: COX1
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 2
retained_cds / bigint: 998

Row 4
gene / text: COX2
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 1000

Row 5
gene / text: COX3
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 1000

Row 6
gene / text: CYTB
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 1000

Row 7
gene / text: ND1
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 1000

Row 8
gene / text: ND2
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 1000

Row 9
gene / text: ND3
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 2
retained_cds / bigint: 998

Row 10
gene / text: ND4
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 1000

Row 11
gene / text: ND4L
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 9
retained_cds / bigint: 991

Row 12
gene / text: ND5
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 1
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 999

Row 13
gene / text: ND6
sample_cds / bigint: 1000
indel_cds / bigint: 0
incomplete_cds / bigint: 0
stop_cds / bigint: 0
invalid_start_cds / bigint: 0
ambiguous_start_cds / bigint: 0
retained_cds / bigint: 1000

Result 2

eligible_codons / bigint | ambiguous_codons / bigint | callable_codons / bigint | asparagine_codons / bigint
-------------------------+---------------------------+--------------------------+---------------------------
3786259                  | 1540                      | 3784719                  | 163715

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

Retained and excluded sample coding regions for the 13 mitochondrial protein-coding genes.

The code checks its translation of the reference against all 13 annotated NCBI proteins. This catches strand, frame and genetic-code mistakes before they reach the matrix.

An N in a DNA sequence means an unresolved base. An N in an amino-acid sequence means asparagine, which is one of our 20 ordinary residues. We therefore exclude ambiguous DNA codons while retaining amino-acid N.

Count amino-acid pairs

The counting query follows the nucleotide version, grouping by gene and codon position first. This listing shows the leucine row; the script computes all 400 cells:

WITH observed AS (
    SELECT a.aa AS aa_from,b.aa AS aa_to,
       sum(CASE WHEN a.aa=b.aa THEN a.n*(a.n-1) ELSE a.n*b.n END)::bigint AS pairs
    FROM scoring.protein_site_counts a JOIN scoring.protein_site_counts b USING(gene,codon1)
    GROUP BY a.aa,b.aa
), alphabet AS (SELECT unnest(string_to_array('A R N D C Q E G H I L K M F P S T W Y V',' ')) AS aa)
SELECT a.aa AS aa_from,b.aa AS aa_to,coalesce(o.pairs,0)::bigint AS pairs
FROM alphabet a CROSS JOIN alphabet b
LEFT JOIN observed o ON o.aa_from=a.aa AND o.aa_to=b.aa
WHERE a.aa = 'L' ORDER BY b.aa;

Result · 20 rows

aa_from / text | aa_to / text | pairs / bigint
---------------+--------------+---------------
L              | A            | 6
L              | C            | 0
L              | D            | 0
L              | E            | 0
L              | F            | 161603
L              | G            | 0
L              | H            | 998
L              | I            | 185997
L              | K            | 0
L              | L            | 641458584
L              | M            | 997
L              | N            | 0
L              | P            | 6991
L              | Q            | 999
L              | R            | 0
L              | S            | 2985
L              | T            | 0
L              | V            | 1994
L              | W            | 0
L              | Y            | 0

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

Longer proteins contribute more columns because we pool the accepted codons across all 13 genes. Within each column, every retained sample has equal weight. The resulting cohort frequencies therefore reflect both the protein lengths and the number of callable samples at each position.

SELECT aa_from, aa_to, pairs,
       round(score_half_bits, 2) AS score_half_bits,
       round(100 * fraction_of_differences, 4) AS percent_of_differences
FROM scoring.protein_matrix
WHERE aa_from = 'L'
ORDER BY aa_to;

Result · 20 rows

aa_from / text | aa_to / text | pairs / bigint | score_half_bits / numeric | percent_of_differences / numeric
---------------+--------------+----------------+---------------------------+---------------------------------
L              | A            | 6              | -45.58                    | 0.0017
L              | C            | 0              | NULL                      | 0.0000
L              | D            | 0              | NULL                      | 0.0000
L              | E            | 0              | NULL                      | 0.0000
L              | F            | 161603         | -15.64                    | 44.5715
L              | G            | 0              | NULL                      | 0.0000
L              | H            | 998            | -28.02                    | 0.2753
L              | I            | 185997         | -16.37                    | 51.2996
L              | K            | 0              | NULL                      | 0.0000
L              | L            | 641458584      | 5.11                      | NULL
L              | M            | 997            | -30.22                    | 0.2750
L              | N            | 0              | NULL                      | 0.0000
L              | P            | 6991           | -24.75                    | 1.9282
L              | Q            | 999            | -27.80                    | 0.2755
L              | R            | 0              | NULL                      | 0.0000
L              | S            | 2985           | -27.85                    | 0.8233
L              | T            | 0              | NULL                      | 0.0000
L              | V            | 1994           | -27.58                    | 0.5500
L              | W            | 0              | NULL                      | 0.0000
L              | Y            | 0              | NULL                      | 0.0000

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

Twenty-by-twenty amino-acid log-odds matrix from corresponding mitochondrial coding positions.

The slashed cells mark pairs we never observed. For a zero observed frequency, the unsmoothed log odds is −∞. The SQL leaves those scores empty rather than replacing them with an arbitrary finite penalty. Smoothing would require an explicit prior and its weight, especially for rare pairs.

All 13 genes retain coding regions from the cohort. The start-codon filter removes 14 sample/gene combinations: one ND5 region starts with GTA, and 13 other regions have an unresolved first codon. Among the retained regions, 1,540 ambiguous codons are excluded, leaving 3,784,719 callable amino-acid observations. Amino-acid N remains present in 163,715 of those observations.

There are 317 positions with more than one observed amino acid. Of the 400 ordered matrix cells, 284 have no observed pair, which is why much of this heatmap is slashed. A thousand closely related samples still provide limited evidence about many possible amino-acid exchanges.

This calculation uses BLOSUM's log-odds idea with our cohort's frequencies. Actual BLOSUM62 construction includes clustering at its identity threshold to control redundant sequences. These human mitochondrial proteins are so similar that applying that threshold would merge much of our cohort, leaving very little independent between-cluster evidence. Equal sample weights instead answer the question we asked here: which aligned residue pairs occur in these thousand records?

Compare the differing amino-acid partners

The original BLOSUM62 joint probabilities are available in NCBI's BLAST source. They retain the frequency information that is lost when scores are rounded to integers. The example bundle includes those reference tables and their source information.

For each amino acid a, we remove its identical partner from the denominator and normalize the remaining ordered joint frequencies:

Apply that equation to both the cohort and BLOSUM62, and each row with observed differences becomes a distribution over differing partners. For example, the L row asks how much of the leucine pairing with a different residue involves isoleucine, valine and the other amino acids.

SELECT p.aa_from, p.aa_to,
       round(100 * p.fraction_of_differences, 4) AS cohort_percent,
       round(100 * b.fraction_of_differences, 4) AS blosum62_percent
FROM scoring.protein_matrix p
JOIN scoring.reference_protein_matrix b USING (aa_from, aa_to)
WHERE b.matrix = 'BLOSUM62'
  AND p.aa_from = 'L' AND p.aa_to <> 'L'
ORDER BY p.aa_to;

Result · 19 rows

aa_from / text | aa_to / text | cohort_percent / numeric | blosum62_percent / numeric
---------------+--------------+--------------------------+---------------------------
L              | A            | 0.0017                   | 7.1538
L              | C            | 0.0000                   | 2.5391
L              | D            | 0.0000                   | 2.4612
L              | E            | 0.0000                   | 3.2429
L              | F            | 44.5715                  | 8.7671
L              | G            | 0.0000                   | 3.3780
L              | H            | 0.2753                   | 1.5979
L              | I            | 51.2996                  | 18.4283
L              | K            | 0.0000                   | 3.9883
L              | M            | 0.2750                   | 7.9806
L              | N            | 0.0000                   | 2.2166
L              | P            | 1.9282                   | 2.2903
L              | Q            | 0.2755                   | 2.6186
L              | R            | 0.0000                   | 3.9170
L              | S            | 0.8233                   | 3.9305
L              | T            | 0.0000                   | 5.3812
L              | V            | 0.5500                   | 15.3459
L              | W            | 0.0000                   | 1.1852
L              | Y            | 0.0000                   | 3.5775

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

Conditional amino-acid partner frequencies calculated from NCBI's original BLOSUM62 probabilities.

Conditional amino-acid partner frequencies calculated from the mitochondrial cohort.

Both plots use the same alphabet and percentage scale. Open either matrix to inspect it at full size. Read across a row to compare the destinations. The underlying joint frequencies are symmetric, while these conditional rows can differ between directions because their denominators differ. In this cohort W, tryptophan, has no observed differing partner, so its denominator is zero and the whole row is undefined, shown by slashes.

The leucine row gives a visible example. In the cohort, isoleucine accounts for 51.30% of differing partners of leucine and phenylalanine for 44.57%. The corresponding BLOSUM62 proportions are 18.43% and 8.77%. Alanine pairs with threonine in 88.61% of its differing cohort pairs, compared with 7.05% in BLOSUM62.

The average log odds is also different: 3.985 bits per aligned pair for the cohort proteins, compared with 0.698 bits for BLOSUM62. That larger value reflects how concentrated the cohort observations are on identical and a few differing pairs. It describes these frequencies, rather than measuring how well the matrix would perform in a search for distant proteins.

BLOSUM62 pools conserved blocks from many protein families at much greater divergence. Our matrix pools a single cohort's mitochondrial proteins, most of which sit in membrane complexes. The two collections give different evidence about amino-acid pairing. The observed differences are a reason to examine the source of a scoring matrix before applying it to a particular alignment task.

Plot the SQL results with Python

At this point the CSV files contain the counts, marginal frequencies, log odds and conditional fractions. Python reads them and chooses how to display the axes and colors.

python3 -m venv .venv
.venv/bin/pip install matplotlib==3.10.7
.venv/bin/python render_charts.py --results results --output figures

from pathlib import Path
from matplotlib.colors import TwoSlopeNorm
import numpy as np
from render_charts import DNA, SCORE_MAP, figure, matrix, rows, save

records = rows(Path('results/nucleotide_matrix.csv'))
scores = matrix(records, DNA, 'base_from', 'base_to', 'score_half_bits')
fig, ax = figure(width=7, height=6)
edges = np.arange(len(DNA) + 1) - 0.5
image = ax.pcolormesh(edges, edges, scores, cmap=SCORE_MAP,
                      norm=TwoSlopeNorm(vcenter=0), shading='flat')
ax.invert_yaxis()
ax.set_xticks(range(4), DNA)
ax.set_yticks(range(4), DNA)
ax.set_xlabel('Partner nucleotide')
ax.set_ylabel('First nucleotide')
fig.colorbar(image, ax=ax, label='Log-odds score (half bits)')
fig.tight_layout()
save(fig, Path('nucleotide-matrix.svg'),
     'Human mitochondrial nucleotide scores',
     'Log-odds scores from aligned nucleotide pairs in the 1,000-genome cohort.')

The complete plotter writes the matrix heatmaps, nucleotide bars, quality plots and coordinate profiles shown here. Its SVG backgrounds are transparent, with separate text and accent colors for light and dark mode. Each plot keeps its numerical scale in its legend.

From pairwise differences to a distance table

A substitution score adds evidence across an alignment. To describe how much two of these samples differ, we can also calculate their mismatch fraction directly:

The pipeline calculates this for eight selected accessions, using the projected alignment columns and removing N or gaps separately for each sample pair:

WITH chosen AS MATERIALIZED (
    SELECT p.accession,r.pos::integer AS position1,r.base
    FROM scoring.projected_sequences p
    CROSS JOIN LATERAL regexp_split_to_table(p.aligned,'') WITH ORDINALITY r(base,pos)
    WHERE p.accession IN (SELECT 'KF'||(161060+142*i)||'.1' FROM generate_series(0,7) x(i))
), counts AS (
    SELECT a.accession AS sample_from,b.accession AS sample_to,
           count(*) AS compared_sites,count(*) FILTER(WHERE a.base<>b.base) AS differences
    FROM chosen a JOIN chosen b USING(position1)
    JOIN scoring.reference_bases r USING(position1)
    WHERE a.base IN ('A','C','G','T') AND b.base IN ('A','C','G','T')
      AND r.base IN ('A','C','G','T')
    GROUP BY a.accession,b.accession
)
SELECT *,round(differences::numeric/compared_sites, 6) AS p_distance FROM counts
ORDER BY sample_from, sample_to;

Result · 64 rows

sample_from / text | sample_to / text | compared_sites / bigint | differences / bigint | p_distance / numeric
-------------------+------------------+-------------------------+----------------------+---------------------
KF161060.1         | KF161060.1       | 16568                   | 0                    | 0.000000
KF161060.1         | KF161202.1       | 16567                   | 29                   | 0.001750
KF161060.1         | KF161344.1       | 16565                   | 9                    | 0.000543
KF161060.1         | KF161486.1       | 16568                   | 15                   | 0.000905
KF161060.1         | KF161628.1       | 16568                   | 11                   | 0.000664
KF161060.1         | KF161770.1       | 16567                   | 10                   | 0.000604
KF161060.1         | KF161912.1       | 16566                   | 12                   | 0.000724
KF161060.1         | KF162054.1       | 16526                   | 13                   | 0.000787
KF161202.1         | KF161060.1       | 16567                   | 29                   | 0.001750
KF161202.1         | KF161202.1       | 16567                   | 0                    | 0.000000
KF161202.1         | KF161344.1       | 16564                   | 25                   | 0.001509
KF161202.1         | KF161486.1       | 16567                   | 31                   | 0.001871
KF161202.1         | KF161628.1       | 16567                   | 27                   | 0.001630
KF161202.1         | KF161770.1       | 16566                   | 24                   | 0.001449
KF161202.1         | KF161912.1       | 16565                   | 28                   | 0.001690
KF161202.1         | KF162054.1       | 16525                   | 29                   | 0.001755
KF161344.1         | KF161060.1       | 16565                   | 9                    | 0.000543
KF161344.1         | KF161202.1       | 16564                   | 25                   | 0.001509
KF161344.1         | KF161344.1       | 16565                   | 0                    | 0.000000
KF161344.1         | KF161486.1       | 16565                   | 8                    | 0.000483
KF161344.1         | KF161628.1       | 16565                   | 6                    | 0.000362
KF161344.1         | KF161770.1       | 16564                   | 5                    | 0.000302
KF161344.1         | KF161912.1       | 16563                   | 7                    | 0.000423
KF161344.1         | KF162054.1       | 16523                   | 8                    | 0.000484
KF161486.1         | KF161060.1       | 16568                   | 15                   | 0.000905
KF161486.1         | KF161202.1       | 16567                   | 31                   | 0.001871
KF161486.1         | KF161344.1       | 16565                   | 8                    | 0.000483
KF161486.1         | KF161486.1       | 16568                   | 0                    | 0.000000
KF161486.1         | KF161628.1       | 16568                   | 10                   | 0.000604
KF161486.1         | KF161770.1       | 16567                   | 11                   | 0.000664
KF161486.1         | KF161912.1       | 16566                   | 13                   | 0.000785
KF161486.1         | KF162054.1       | 16526                   | 14                   | 0.000847
KF161628.1         | KF161060.1       | 16568                   | 11                   | 0.000664
KF161628.1         | KF161202.1       | 16567                   | 27                   | 0.001630
KF161628.1         | KF161344.1       | 16565                   | 6                    | 0.000362
KF161628.1         | KF161486.1       | 16568                   | 10                   | 0.000604
KF161628.1         | KF161628.1       | 16568                   | 0                    | 0.000000
KF161628.1         | KF161770.1       | 16567                   | 7                    | 0.000423
KF161628.1         | KF161912.1       | 16566                   | 9                    | 0.000543
KF161628.1         | KF162054.1       | 16526                   | 10                   | 0.000605
KF161770.1         | KF161060.1       | 16567                   | 10                   | 0.000604
KF161770.1         | KF161202.1       | 16566                   | 24                   | 0.001449
KF161770.1         | KF161344.1       | 16564                   | 5                    | 0.000302
KF161770.1         | KF161486.1       | 16567                   | 11                   | 0.000664
KF161770.1         | KF161628.1       | 16567                   | 7                    | 0.000423
KF161770.1         | KF161770.1       | 16567                   | 0                    | 0.000000
KF161770.1         | KF161912.1       | 16565                   | 8                    | 0.000483
KF161770.1         | KF162054.1       | 16525                   | 9                    | 0.000545
KF161912.1         | KF161060.1       | 16566                   | 12                   | 0.000724
KF161912.1         | KF161202.1       | 16565                   | 28                   | 0.001690
KF161912.1         | KF161344.1       | 16563                   | 7                    | 0.000423
KF161912.1         | KF161486.1       | 16566                   | 13                   | 0.000785
KF161912.1         | KF161628.1       | 16566                   | 9                    | 0.000543
KF161912.1         | KF161770.1       | 16565                   | 8                    | 0.000483
KF161912.1         | KF161912.1       | 16566                   | 0                    | 0.000000
KF161912.1         | KF162054.1       | 16524                   | 11                   | 0.000666
KF162054.1         | KF161060.1       | 16526                   | 13                   | 0.000787
KF162054.1         | KF161202.1       | 16525                   | 29                   | 0.001755
KF162054.1         | KF161344.1       | 16523                   | 8                    | 0.000484
KF162054.1         | KF161486.1       | 16526                   | 14                   | 0.000847
KF162054.1         | KF161628.1       | 16526                   | 10                   | 0.000605
KF162054.1         | KF161770.1       | 16525                   | 9                    | 0.000545
KF162054.1         | KF161912.1       | 16524                   | 11                   | 0.000666
KF162054.1         | KF162054.1       | 16526                   | 0                    | 0.000000

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

Pairwise callable-site mismatch percentages for eight sample mitochondrial genomes.

This distance table is the kind of input a distance-based tree method can use. A phylogenetic analysis can then choose a substitution model, estimate branch lengths and test how strongly the data support the relationships. Models that correct for repeated changes become more relevant as divergence grows. A raw BLOSUM alignment score alone does not specify that evolutionary distance: alignment length, composition and the scoring settings all affect it.

What to try with the tables

The 1,000 complete mitochondrial records give us a nucleotide matrix with a strong transition pattern and a protein matrix concentrated on a small set of observed exchanges. We have also kept the missing observations visible: 8,691 raw nucleotide Ns, alignment gaps and the coding-region exclusions all have their own counts.

PAM and BLOSUM both turn substitution evidence into additive alignment scores, with different choices about their input alignments and how they represent divergence. Gap penalties then determine how those substitution scores compete with opening or extending an insertion or deletion. That combination decides which alignment path wins.

In our database the counts remain beside the accession list, coding annotations and alignment paths. You can select one mitochondrial gene, change which sample records contribute, or inspect the coordinates behind a particular matrix cell with another SQL query. Recalculate its marginal frequencies as well as its pair counts, then rerun the plotter to see how the result changes.

References and further reading

The scoring-matrix explanations follow the topics in Scoring Matrices in Bioinformatics: PAM, BLOSUM, Log-Odds, Gap Penalties. The ForensicMCQ guide supplies the small log-odds and alignment-score examples, and the Columbia course presentation is useful for working through the column counts by hand.

For the original construction, read Henikoff and Henikoff, 1992. Pearson, 2013 discusses matrix selection and alignment length. Trivedi and Nagarajaram, 2020 covers general and specialized protein matrices, and Mount, 2008 discusses using BLOSUM in sequence alignments.

The cohort records come from Li and colleagues' study of 2,000 Danish mitochondrial genomes. NCBI provides the reference sequence and its annotations, the vertebrate mitochondrial genetic code, and the reference matrix frequencies.

Please check my OS project - BioseqDB

Read this tutorial on Medium. Follow BioseqDB on LinkedIn.

Working with sequence files in BioseqDB and SeqKit

BioseqDB — Mon, 19 Oct 2026 17:16:18 +0000

A FASTA header can carry quite a lot of information. The sequence has a name, somebody adds an organism, then an assembly and perhaps some coordinates, and eventually you need a little regular expression just to recover something that could have been a column in a table.

A FASTA record starts with a header line beginning with >, followed by the sequence letters on one or more lines.

I like having those columns. This is one reason I wanted sequence operations inside PostgreSQL when I started working on BioseqDB. You can keep the DNA next to its metadata and write a query that uses both, without exporting another intermediate file every time you change a filter.

SeqKit's tutorial has plenty of useful examples of the work people do with sequence files. I picked several of them and some tasks from the FAQ to try in SQL: checking lengths, counting bases, finding motifs, selecting records and dealing with names.

We will use the same sequences in both tools. I checked the examples with SeqKit 2.14.0 and a BioseqDB development build with external database support. You can find the SQL, shell commands and chart code in BioseqDB/examples/seqkit_workflows.

1. Get a small dataset from UCSC

For this exercise I chose three mitochondrial chromosomes: human, mouse and zebrafish. Mitochondria have their own DNA, separate from the much larger nuclear chromosomes, so we get complete chromosome records small enough to inspect without waiting for a large download. These three records total 49,464 bases.

Our source is the UCSC Genome Browser API. BioseqDB's external tables can fetch sequence intervals from UCSC directly through SQL. We will use assembly names hg38, mm39 and danRer11, each with chromosome chrM.

Warning

Start with the installation instructions. The source-build procedure in the butterfly example also applies here. In an already configured source checkout, enable the external module before building and installing:

cmake -S . -B build/seqkit -DBIOSEQ_BUILD_EXTERNAL_API_CLIENTS=ON
cmake --build build/seqkit --target postgresql_bioseq postgresql_bioseq_cache

Connect to a fresh database where that build is installed, run CREATE EXTENSION bioseq;, then check for the external sequence table:

SELECT to_regclass('bioseq_external.sequences') IS NOT NULL
       AS external_sequences_available;

Result · 1 row

external_sequences_available / boolean
--------------------------------------
t

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

The result must be true. If the table exists but access is denied, your database owner needs to grant your role USAGE on bioseq_external and SELECT on bioseq_external.sequences.

Create a schema for this exercise and a table for the chromosome records:

CREATE SCHEMA seqkit_lab;
CREATE TABLE seqkit_lab.genomes (
    assembly text NOT NULL,
    species 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

Column	What we put there
`assembly`	UCSC assembly identifier, such as `hg38`
`species`	A short label for grouping and plotting
`seq`	The DNA sequence, stored as a BioseqDB value

USING bioseq selects BioseqDB's table storage. The dna5 alphabet accepts A, C, G, T and N, where N means an unresolved base. The 64 selects 64-bit storage blocks, and a sequence can occupy many such blocks. The sequence representation documentation explains the storage details.

Now fetch the three records. Each INSERT ... SELECT reads a remote interval and saves it locally. Subsequent analysis will use the saved rows.

INSERT INTO seqkit_lab.genomes (assembly, species, seq)
SELECT 'hg38', 'human', sequence::bioseq.seq('dna5', 64)
FROM bioseq_external.sequences
WHERE provider = 'ucsc' AND assembly_database = 'hg38'
  AND chromosome = 'chrM' AND start0 = 0 AND end0 = 16569;

INSERT INTO seqkit_lab.genomes (assembly, species, seq)
SELECT 'mm39', 'mouse', sequence::bioseq.seq('dna5', 64)
FROM bioseq_external.sequences
WHERE provider = 'ucsc' AND assembly_database = 'mm39'
  AND chromosome = 'chrM' AND start0 = 0 AND end0 = 16299;

INSERT INTO seqkit_lab.genomes (assembly, species, seq)
SELECT 'danRer11', 'zebrafish', sequence::bioseq.seq('dna5', 64)
FROM bioseq_external.sequences
WHERE provider = 'ucsc' AND assembly_database = 'danRer11'
  AND chromosome = 'chrM' AND start0 = 0 AND end0 = 16596;

SELECT assembly, species, bioseq.len(seq) AS bp
FROM seqkit_lab.genomes ORDER BY species;

INSERT 0 1 — 1 row affected

assembly / text | species / text | bp / bigint
----------------+----------------+------------
hg38            | human          | 16569
mm39            | mouse          | 16299
danRer11        | zebrafish      | 16596

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

UCSC's start0 is zero-based and end0 is exclusive. The human interval [0, 16569) therefore contains 16,569 bases, including the first and last base of that record. These coordinates refer to each assembly separately.

Lengths of the three mitochondrial chromosome records retrieved from UCSC.

The records differ by only a few hundred bases. Their lengths match the intervals requested from UCSC, which gives us a quick check that all three imports returned complete records.

2. Make records we can filter and export

Three rows would make a rather short filtering exercise. We will split each chromosome into non-overlapping windows of 1,000 bases and keep the remaining shorter piece at the end. This gives us 51 records to filter, with each record's position on its original chromosome saved alongside it.

CREATE TABLE seqkit_lab.windows (
    id text NOT NULL,
    species text NOT NULL,
    assembly text NOT NULL,
    start0 integer NOT NULL,
    end0 integer NOT NULL,
    seq bioseq.seq('dna5', 64) NOT NULL
) USING bioseq;
INSERT INTO seqkit_lab.windows
SELECT species || '_' || lpad((start0 / 1000 + 1)::text, 2, '0'),
       species, assembly, start0,
       least(start0 + 1000, bioseq.len(seq)),
       seq[start0:least(start0 + 1000, bioseq.len(seq))]
FROM seqkit_lab.genomes
CROSS JOIN LATERAL generate_series(0, bioseq.len(seq) - 1, 1000) AS p(start0);

SELECT species, count(*) AS windows FROM seqkit_lab.windows
GROUP BY species ORDER BY species;

CREATE TABLE

INSERT 0 51 — 51 rows affected

species / text | windows / bigint
---------------+-----------------
human          | 17
mouse          | 17
zebrafish      | 17

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

generate_series supplies the start positions. The slice seq[start0:end0] uses the same zero-based, half-open convention as the UCSC request, and least stops the final window at the end of the chromosome. We retain the coordinates because a motif position inside human_05 is otherwise easy to confuse with a position on the original chromosome.

Each row has an id, species, assembly, start0, end0 and seq. Here are the first six:

SELECT id, assembly, start0, end0, bioseq.len(seq) AS bp,
       seq[:12]::text AS prefix
FROM seqkit_lab.windows
ORDER BY id
LIMIT 6;

Result · 6 rows

Result rows (long values continue on indented lines):

Row 1
id / text: human_01
assembly / text: hg38
start0 / integer: 0
end0 / integer: 1000
bp / bigint: 1000
prefix / text: GATCACAGGTCT

Row 2
id / text: human_02
assembly / text: hg38
start0 / integer: 1000
end0 / integer: 2000
bp / bigint: 1000
prefix / text: CCAGTTGACACA

Row 3
id / text: human_03
assembly / text: hg38
start0 / integer: 2000
end0 / integer: 3000
bp / bigint: 1000
prefix / text: CGAGCCTGGTGA

Row 4
id / text: human_04
assembly / text: hg38
start0 / integer: 3000
end0 / integer: 4000
bp / bigint: 1000
prefix / text: GGACATCCCGAT

Row 5
id / text: human_05
assembly / text: hg38
start0 / integer: 4000
end0 / integer: 5000
bp / bigint: 1000
prefix / text: TTATAATAAACA

Row 6
id / text: human_06
assembly / text: hg38
start0 / integer: 5000
end0 / integer: 6000
bp / bigint: 1000
prefix / text: ATCTTAGCATAC

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

For the command-line comparisons, install SeqKit and export these same rows. Run the following in psql. The \o command saves the query output to a file in your terminal's current directory.

\pset format unaligned
\pset tuples_only on
\o genomes.fa
SELECT '>' || assembly || E'\n' || seq::text
FROM seqkit_lab.genomes ORDER BY assembly;
\o windows.fa
SELECT '>' || id || ' species=' || species || ' assembly=' || assembly
       || ' start0=' || start0 || ' end0=' || end0 || E'\n' || seq::text
FROM seqkit_lab.windows ORDER BY id;
\o

For example, a header will start with >human_01 species=human assembly=hg38 start0=0 end0=1000. SeqKit uses the first whitespace-delimited part as the ID by default. The SQL table already has each field separately.

The example repository includes these export commands in 02_export.sql. Run them after 01_import.sql to create genomes.fa and windows.fa for the SeqKit examples below.

3. Check lengths before doing anything clever

A missing record or a truncated import can ruin all the later analysis. Start by counting the rows and their bases.

BioseqDB

SELECT count(*) AS sequences, sum(bioseq.len(seq)) AS total_bp,
       min(bioseq.len(seq)) AS min_bp,
       round(avg(bioseq.len(seq)), 1) AS mean_bp,
       max(bioseq.len(seq)) AS max_bp
FROM seqkit_lab.windows;

Result · 1 row

sequences / bigint | total_bp / numeric | min_bp / bigint | mean_bp / numeric | max_bp / bigint
-------------------+--------------------+-----------------+-------------------+----------------
51                 | 49464              | 299             | 969.9             | 1000

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit stats -T windows.fa

file    format  type    num_seqs    sum_len min_len avg_len max_len
windows.fa  FASTA   DNA 51  49464   299 969.9   1000

Both tools report 51 records and 49,464 bases. The shortest record is 299 bases, which is the final mouse window. The mean is about 970 bases because all three final windows are shorter than 1,000.

To find those short records, use the sequence length as a filter:

BioseqDB

SELECT id, bioseq.len(seq) AS bp
FROM seqkit_lab.windows
WHERE bioseq.len(seq) < 1000
ORDER BY id;

Result · 3 rows

id / text    | bp / bigint
-------------+------------
human_17     | 569
mouse_17     | 299
zebrafish_17 | 596

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit seq --max-len 999 windows.fa | seqkit fx2tab -n -i -l

human_17    569
mouse_17    299
zebrafish_17    596

bioseq.len counts sequence symbols. SeqKit's --max-len (-M) accepts an inclusive maximum length, so -M 999 corresponds to our SQL condition < 1000.

4. Count bases and calculate GC content

GC content is the percentage of known DNA bases that are G or C. A sequence with 30 Gs and 20 Cs among 100 known bases has 50% GC content. We can calculate it for a whole chromosome or for smaller windows to see where its base composition changes.

First count the letters in each complete chromosome:

SELECT species,
       bioseq.symbol_count(seq, 'A') AS a,
       bioseq.symbol_count(seq, 'C') AS c,
       bioseq.symbol_count(seq, 'G') AS g,
       bioseq.symbol_count(seq, 'T') AS t,
       bioseq.symbol_count(seq, 'N') AS n
FROM seqkit_lab.genomes ORDER BY species;

Result · 3 rows

Result rows (long values continue on indented lines):

Row 1
species / text: human
a / bigint: 5124
c / bigint: 5181
g / bigint: 2169
t / bigint: 4094
n / bigint: 1

Row 2
species / text: mouse
a / bigint: 5629
c / bigint: 3976
g / bigint: 2013
t / bigint: 4681
n / bigint: 0

Row 3
species / text: zebrafish
a / bigint: 5301
c / bigint: 3969
g / bigint: 2658
t / bigint: 4668
n / bigint: 0

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

Counts of A, C, G and T in the three mitochondrial records, with the unresolved-base count shown separately.

There is one N in the human record and none in the other two. Dropping it would change the coordinate system, so we leave it in the stored sequence. For GC content, however, we exclude N from the denominator, matching SeqKit's calculation on these dna5 records.

BioseqDB

SELECT species, bioseq.len(seq) AS bp,
       round(100.0 * (bioseq.symbol_count(seq, 'G') + bioseq.symbol_count(seq, 'C'))
             / nullif(bioseq.len(seq) - bioseq.symbol_count(seq, 'N'), 0), 2) AS gc_percent,
       bioseq.symbol_count(seq, 'N') AS unknown_bases
FROM seqkit_lab.genomes ORDER BY species;

Result · 3 rows

species / text | bp / bigint | gc_percent / numeric | unknown_bases / bigint
---------------+-------------+----------------------+-----------------------
human          | 16569       | 44.36                | 1
mouse          | 16299       | 36.74                | 0
zebrafish      | 16596       | 39.93                | 0

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit fx2tab -n -i -l -g -C N genomes.fa

danRer11    16596   39.93   0
hg38    16569   44.36   1
mm39    16299   36.74   0

bioseq.symbol_count counts one letter. The 100.0 makes this a percentage calculation with a decimal result, and nullif(..., 0) avoids dividing by zero if a future record contains only Ns. With alphabets containing other ambiguity symbols, the denominator should explicitly count A, C, G and T instead.

The whole human mitochondrial record has 44.36% GC, the mouse record 36.74%, and zebrafish 39.93%. Now calculate the same statistic per window:

BioseqDB

SELECT id,
       round(100.0 * (bioseq.symbol_count(seq, 'G') + bioseq.symbol_count(seq, 'C'))
             / nullif(bioseq.len(seq) - bioseq.symbol_count(seq, 'N'), 0), 2) AS gc_percent
FROM seqkit_lab.windows ORDER BY id LIMIT 6;

Result · 6 rows

id / text | gc_percent / numeric
----------+---------------------
human_01  | 45.60
human_02  | 44.00
human_03  | 42.70
human_04  | 47.75
human_05  | 42.60
human_06  | 43.90

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit head -n 6 windows.fa | seqkit fx2tab -n -i -g

human_01    45.60
human_02    44.00
human_03    42.70
human_04    47.75
human_05    42.60
human_06    43.90

GC percentage along each mitochondrial chromosome, measured in the tutorial's non-overlapping windows.

Each point sits at the midpoint of a window. The final window is shorter, and the lines connect neighbouring measurements so we can follow how GC content changes along each chromosome. The horizontal position is an offset within each species' own chromosome. To find corresponding regions between species, we would align the sequences, as in the butterfly tutorial.

5. Select records by ID and by their contents

The SeqKit FAQ starts several workflows with a list of wanted IDs. SQL's IN expression is enough for a short list:

BioseqDB

SELECT id, bioseq.len(seq) AS bp
FROM seqkit_lab.windows
WHERE id IN ('human_01', 'mouse_02')
ORDER BY id;

Result · 2 rows

id / text | bp / bigint
----------+------------
human_01  | 1000
mouse_02  | 1000

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

printf 'human_01\nmouse_02\n' > wanted.txt
seqkit grep -f wanted.txt windows.fa | seqkit fx2tab -n -i -l

human_01    1000
mouse_02    1000

If you have thousands of IDs, put them in a table and join on id. The SeqKit equivalent accepts a file with grep -f ids.txt. To exclude the two IDs instead, reverse the condition:

BioseqDB

SELECT count(*) AS retained
FROM seqkit_lab.windows
WHERE id NOT IN ('human_01', 'mouse_02');

Result · 1 row

retained / bigint
-----------------
49

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

printf 'human_01\nmouse_02\n' > wanted.txt
seqkit grep -v -f wanted.txt -C windows.fa

Both retain 49 records. The SeqKit tutorial uses this filtering step to remove contaminated reads. Once you have identified those reads, put their IDs in the exclusion list and use the same query.

We can also find the window containing the unresolved base:

BioseqDB

SELECT id, bioseq.symbol_count(seq, 'N') AS unknown_bases
FROM seqkit_lab.windows
WHERE bioseq.symbol_count(seq, 'N') > 0
ORDER BY id;

Result · 1 row

id / text | unknown_bases / bigint
----------+-----------------------
human_04  | 1

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit grep -s -P -r -p '[^ACGT]' windows.fa | seqkit fx2tab -n -i -C N

human_04    1

The result is human_04. On the SeqKit side, -s searches the sequence rather than its ID, -r enables a regular expression and -P searches only the positive strand. [^ACGT] means any character outside those four letters. Our SQL test is equivalent for this dna5 dataset, whose only possible extra symbol is N.

A motif is a short sequence pattern we want to find inside a longer sequence. For a literal motif, try GAATTC:

BioseqDB

SELECT id
FROM seqkit_lab.windows
WHERE seq::text LIKE '%GAATTC%'
ORDER BY id;

Result · 8 rows

id / text
------------
human_05
human_06
human_13
mouse_02
mouse_04
mouse_05
zebrafish_09
zebrafish_15

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit grep -s -P -p GAATTC windows.fa | seqkit seq -n -i

human_05
human_06
human_13
mouse_02
mouse_04
mouse_05
zebrafish_09
zebrafish_15

Eight windows contain it. SQL's % matches any text before or after the motif, while SeqKit's grep -s performs the corresponding sequence search. We keep -P because an ordinary SQL text search examines the stored strand. SeqKit also checks the reverse complement by default.

A variable position in a motif can be written as a regular expression. TATA[AT]AA allows A or T in the fifth position:

BioseqDB

SELECT id FROM seqkit_lab.windows
WHERE seq::text ~ 'TATA[AT]AA'
ORDER BY id;

Result · 10 rows

id / text
------------
human_07
human_13
mouse_01
mouse_02
mouse_03
mouse_05
mouse_08
mouse_11
mouse_12
zebrafish_01

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit grep -s -P -r -p 'TATA[AT]AA' windows.fa | seqkit seq -n -i

human_07
human_13
mouse_01
mouse_02
mouse_03
mouse_05
mouse_08
mouse_11
mouse_12
zebrafish_01

This matches ten windows. The FAQ also demonstrates IUPAC degenerate symbols with SeqKit's -d option. Here we use an explicit character class so the two listings have the same meaning.

These searches find motifs contained within individual windows. To include matches spanning two windows, search the complete chromosomes in seqkit_lab.genomes. Mitochondrial DNA is circular, so a search across the chromosome's origin also needs a sequence slice joining its last bases to its first bases.

6. Locate a motif and take a slice

To extract the bases around a motif, we need its position inside the record:

BioseqDB

SELECT id, strpos(seq::text, 'GAATTC') - 1 AS first_start0
FROM seqkit_lab.windows
WHERE seq::text LIKE '%GAATTC%'
ORDER BY id;

Result · 8 rows

id / text    | first_start0 / integer
-------------+-----------------------
human_05     | 120
human_06     | 273
human_13     | 639
mouse_02     | 749
mouse_04     | 802
mouse_05     | 12
zebrafish_09 | 777
zebrafish_15 | 448

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit locate -P -p GAATTC windows.fa

seqID   patternName pattern strand  start   end matched
human_05    GAATTC  GAATTC  +   121 126 GAATTC
human_06    GAATTC  GAATTC  +   274 279 GAATTC
human_13    GAATTC  GAATTC  +   640 645 GAATTC
mouse_02    GAATTC  GAATTC  +   750 755 GAATTC
mouse_04    GAATTC  GAATTC  +   803 808 GAATTC
mouse_05    GAATTC  GAATTC  +   13  18  GAATTC
zebrafish_09    GAATTC  GAATTC  +   778 783 GAATTC
zebrafish_15    GAATTC  GAATTC  +   449 454 GAATTC

strpos returns the first occurrence using PostgreSQL's one-based string positions. Subtracting one converts it to a zero-based offset inside the window. SeqKit reports one-based inclusive start and end coordinates, so a SQL offset of 120 corresponds to SeqKit start 121. Add the window's start0 to get the chromosome offset.

There is one GAATTC occurrence in each of these eight matching windows, so both listings return eight positions. With repeated occurrences in a window, strpos returns the first one and seqkit locate returns all of them.

Taking the first twelve bases is shorter to write:

BioseqDB

SELECT id, seq[:12]::text AS first_12
FROM seqkit_lab.windows ORDER BY id LIMIT 3;

Result · 3 rows

id / text | first_12 / text
----------+----------------
human_01  | GATCACAGGTCT
human_02  | CCAGTTGACACA
human_03  | CGAGCCTGGTGA

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit head -n 3 windows.fa | seqkit subseq -r 1:12 | seqkit fx2tab -Q -i

human_01    GATCACAGGTCT
human_02    CCAGTTGACACA
human_03    CGAGCCTGGTGA

BioseqDB's seq[:12] selects offsets 0 through 11. SeqKit's 1:12 names the same bases using inclusive positions. Similarly, seq[-10:] takes the last ten bases, while seq[9:20] takes bases 10 through 20.

To inspect the opposite strand, reverse the sequence and complement A with T and C with G:

BioseqDB

SELECT id, seq[:20]::text AS forward,
       bioseq.reverse_complement(seq[:20])::text AS reverse_complement
FROM seqkit_lab.windows
WHERE id IN ('human_01', 'mouse_01', 'zebrafish_01')
ORDER BY id;

Result · 3 rows

id / text    | forward / text       | reverse_complement / text
-------------+----------------------+--------------------------
human_01     | GATCACAGGTCTATCACCCT | AGGGTGATAGACCTGTGATC
mouse_01     | GTTAATGTAGCTTAATAACA | TGTTATTAAGCTACATTAAC
zebrafish_01 | ACGGCCGGCGACAATTTATA | TATAAATTGTCGCCGGCCGT

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit grep -r -p '_01$' windows.fa | seqkit subseq -r 1:20 | seqkit seq -t DNA -r -p | seqkit fx2tab -Q -i

human_01    AGGGTGATAGACCTGTGATC
mouse_01    TGTTATTAAGCTACATTAAC
zebrafish_01    TATAAATTGTCGCCGGCCGT

The query uses only the first twenty bases of the first window from each species so that the result fits on screen. SeqKit's -r -p options apply the same reverse-complement operation. Reversing alone changes their order without complementing A to T or C to G.

7. Group by species, then choose an order

The FAQ's header-splitting example extracts a group name from each FASTA header and writes separate files. We already have a species column, so we can inspect the groups directly:

BioseqDB

SELECT species, count(*) AS records
FROM seqkit_lab.windows
GROUP BY species ORDER BY species;

Result · 3 rows

species / text | records / bigint
---------------+-----------------
human          | 17
mouse          | 17
zebrafish      | 17

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit split --by-id --id-regexp '^([^_]+)_' --out-dir by-species --force windows.fa
seqkit stats -T by-species/*.fa

file    format  type    num_seqs    sum_len min_len avg_len max_len
by-species/windows.part_human.fa    FASTA   DNA 17  16569   569 974.6   1000
by-species/windows.part_mouse.fa    FASTA   DNA 17  16299   299 958.8   1000
by-species/windows.part_zebrafish.fa    FASTA   DNA 17  16596   596 976.2   1000

There are seventeen windows per species. To work with the human group, filter the column. The SeqKit command below writes three files by extracting the prefix before the underscore in each ID.

BioseqDB

SELECT id, bioseq.len(seq) AS bp
FROM seqkit_lab.windows WHERE species = 'human'
ORDER BY id LIMIT 3;

Result · 3 rows

id / text | bp / bigint
----------+------------
human_01  | 1000
human_02  | 1000
human_03  | 1000

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit split --by-id --id-regexp '^([^_]+)_' --out-dir by-species --force windows.fa
seqkit head -n 3 by-species/windows.part_human.fa | seqkit fx2tab -n -i -l

human_01    1000
human_02    1000
human_03    1000

The SQL listing previews three human rows. You can apply WHERE species = 'human' to any subsequent query to work with that group. To save it as a FASTA file for another program, add the same condition to the export query from section 2.

For length sorting, the three short tail windows make the result easy to check:

BioseqDB

SELECT id, bioseq.len(seq) AS bp
FROM seqkit_lab.windows
ORDER BY bioseq.len(seq), id
LIMIT 3;

Result · 3 rows

id / text    | bp / bigint
-------------+------------
mouse_17     | 299
human_17     | 569
zebrafish_17 | 596

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit sort -l windows.fa | seqkit head -n 3 | seqkit fx2tab -n -i -l

mouse_17    299
human_17    569
zebrafish_17    596

SQL adds id as a tie-breaker because the other 48 records all have length 1,000. ORDER BYbioseq.len(seq), id gives these equal-length records a consistent order too. In SQL, use ORDER BY whenever the order of the results matters.

That last detail matters when selecting a range of records:

BioseqDB

SELECT id, bioseq.len(seq) AS bp
FROM seqkit_lab.windows
ORDER BY id
LIMIT 5 OFFSET 3;

Result · 5 rows

id / text | bp / bigint
----------+------------
human_04  | 1000
human_05  | 1000
human_06  | 1000
human_07  | 1000
human_08  | 1000

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit range -r 4:8 windows.fa | seqkit fx2tab -n -i -l

human_04    1000
human_05    1000
human_06    1000
human_07    1000
human_08    1000

These are records four through eight in ID order. OFFSET 3 skips the first three and LIMIT 5 takes the next five. SeqKit's range selects whole records in file order, which agrees here because we exported with ORDER BY id.

8. Rename, edit and join small sequence fragments

Adding a prefix to a name is a string operation:

BioseqDB

SELECT id, 'mito_' || id AS renamed
FROM seqkit_lab.windows ORDER BY id LIMIT 3;

Result · 3 rows

id / text | renamed / text
----------+---------------
human_01  | mito_human_01
human_02  | mito_human_02
human_03  | mito_human_03

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit head -n 3 windows.fa | seqkit replace -p '^' -r 'mito_' | seqkit seq -n -i

mito_human_01
mito_human_02
mito_human_03

|| concatenates ordinary PostgreSQL text here. For more complicated header changes, PostgreSQL has regexp_replace. For a lookup table of old and new names, use a join. SeqKit's replace also supports regular expressions and key-value replacement files.

We can make a copy with its fifth base set to T:

BioseqDB

SELECT id, seq[:20]::text AS before,
       overlay(seq[:20]::text placing 'T' from 5 for 1) AS after
FROM seqkit_lab.windows ORDER BY id LIMIT 3;

Result · 3 rows

id / text | before / text        | after / text
----------+----------------------+---------------------
human_01  | GATCACAGGTCTATCACCCT | GATCTCAGGTCTATCACCCT
human_02  | CCAGTTGACACAAAATAGAC | CCAGTTGACACAAAATAGAC
human_03  | CGAGCCTGGTGATAGCTGGT | CGAGTCTGGTGATAGCTGGT

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit head -n 3 windows.fa | seqkit subseq -r 1:20 | seqkit mutate -p 5:T | seqkit fx2tab -Q -i

human_01    GATCTCAGGTCTATCACCCT
human_02    CCAGTTGACACAAAATAGAC
human_03    CGAGTCTGGTGATAGCTGGT

PostgreSQL's text function overlay uses one-based positions, so from 5 for 1 agrees with seqkit mutate -p 5:T. This SELECT returns the edited twenty-base text beside the original while leaving the stored sequence unchanged. Cast a complete edited string back to bioseq.seq('dna5', 64) when you want to store it as DNA.

SeqKit's concat joins sequences with matching IDs across files. We can do the same with a SQL join:

BioseqDB

WITH left_part AS (
  SELECT id, seq[:10] AS seq
  FROM seqkit_lab.windows ORDER BY id LIMIT 3
), right_part AS (
  SELECT id, seq[-10:] AS seq
  FROM seqkit_lab.windows ORDER BY id LIMIT 3
)
SELECT l.id,
       (l.seq OPERATOR(bioseq.||) r.seq)::text AS joined
FROM left_part l JOIN right_part r USING (id)
ORDER BY l.id;

Result · 3 rows

id / text | joined / text
----------+---------------------
human_01  | GATCACAGGTGTAAAAAACT
human_02  | CCAGTTGACAACAAACCTAC
human_03  | CGAGCCTGGTTGTTGGATCA

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit head -n 3 windows.fa > trio.fa
seqkit subseq -r 1:10 trio.fa > left.fa
seqkit subseq -r -10:-1 trio.fa > right.fa
seqkit concat left.fa right.fa | seqkit sort | seqkit fx2tab -Q -i

human_01    GATCACAGGTGTAAAAAACT
human_02    CCAGTTGACAACAAACCTAC
human_03    CGAGCCTGGTTGTTGGATCA

Here left.fa contains the first ten bases of the first three human windows, and right.fa contains their last ten bases with the same IDs. The repository's seqkit.sh script prepares those files. The SQL common table expressions prepare the same two inputs, then JOIN ... USING (id) pairs them. OPERATOR(bioseq.||) concatenates the BioseqDB sequence values. The result is a twenty-base sequence containing the two ends of each window.

seqkit rename adds suffixes to repeated IDs, while seqkit rmdup -s compares sequence contents. To try both, seqkit.sh creates three records in duplicates.fa: two records share the ID copy, the second has an edited final base, and the third has ID other but repeats the first sequence exactly.

For that input, seqkit rename duplicates.fa returns IDs copy, copy_2, other. In SQL, number the rows within each ID using row_number() OVER (PARTITION BY id ORDER BY ordinal) and append the number after the first occurrence. Keep an input-order column such as ordinal if the identity of the first record matters.

To remove repeated sequence contents, seqkit rmdup -s -P -i duplicates.fa retains the two distinct sequences. PostgreSQL's DISTINCT ON (seq) does the same forward-strand grouping, with ORDER BYseq, ordinal choosing the first copy. -P restricts SeqKit's comparison to the forward strand, matching the SQL grouping. The complete runnable queries for both operations are in 03_queries.sql.

9. Shuffle and take a reproducible subset

For an exploratory subset in SQL, you can write ORDER BY random() LIMIT 5. Re-running it changes the sample. Here we use a fixed ordering derived from each ID and the label '2026:', so you can repeat the query and get the same five records:

BioseqDB

SELECT id
FROM seqkit_lab.windows
ORDER BY md5('2026:' || id)
LIMIT 5;

Result · 5 rows

id / text
------------
mouse_09
zebrafish_04
human_06
human_13
human_09

Executed against Bioseq extension 1.0 on PostgreSQL 17.10

SeqKit

seqkit sample -p 0.1 -s 11 windows.fa | seqkit seq -n -i

human_01
human_11
mouse_10
zebrafish_08
zebrafish_09
zebrafish_16

Changing '2026:' changes the SQL ordering. Remove LIMIT 5 to reorder all rows using the same method. SeqKit has shuffle -s 23 windows.fa for a seeded shuffle and sample -p 0.1 -s 11 windows.fa for a seeded fraction sample.

The SQL example returns exactly five rows. SeqKit's -p 0.1 gives each record a 10% chance of inclusion and returned six rows on this input with seed 11. The chosen IDs differ because these are two different ways of selecting a subset. For sampling a fixed number from large files, read the SeqKit sample documentation and its notes on the count option.

A short command map

Here are the commands and expressions from the examples above, gathered in one table for reference.

Task	BioseqDB / PostgreSQL expression	SeqKit command
Count records and inspect lengths	`count(*)`, `min(bioseq.len(seq))`, `max(bioseq.len(seq))`	stats
Reverse-complement DNA or filter length	`bioseq.reverse_complement(seq)`, `WHEREbioseq.len(seq) < 1000`	seq
Keep IDs or sequences matching a pattern	`WHERE id IN (...)`, `WHEREseq::text LIKE '%GAATTC%'`	grep
Select records by position in an ordered input	`ORDER BY id LIMIT 5 OFFSET 3`	range
Order records by sequence length	`ORDER BYbioseq.len(seq), id`	sort
Change record order	`ORDER BY random()`	shuffle
Select a subset	`ORDER BY random() LIMIT 5` for five rows	sample
Separate groups encoded in headers	Filter a metadata column, export groups if files are needed	split
Change names with a pattern or mapping	Text concatenation, `regexp_replace`, or a metadata join	replace
Make repeated IDs distinct	`row_number() OVER (PARTITION BY id ORDER BY ordinal)` plus a suffix	rename
Edit a sequence copy	`overlay(seq::text placing 'T' from 5 for 1)`	mutate
Join fragments sharing an ID	`JOIN ... USING (id)` and `OPERATOR(bioseq.\	\
Keep one copy of each stored sequence	{% raw %}`DISTINCT ON (seq)` with an explicit choice of retained row	rmdup

For me the useful difference appears when the next question involves metadata. Once the sequences are in seqkit_lab.windows, adding a species condition or joining a list of interesting regions is an ordinary SQL query. The data stays beside the information explaining where it came from. SeqKit remains convenient when I need to inspect or transform a file from the terminal, and the FASTA export means I can use it on exactly the same records.

For a comparison of longer regions that allows substitutions, insertions and deletions, try the BioseqDB alignment tutorial. The queries there return matching intervals you can plot, along with the species and coordinate columns used throughout the analysis.

Sources

This tutorial was inspired by the awesome work of Wei Shen on SeqKit. His FAQ and command reference explain the SeqKit commands used here and many more.

The sequences come from the UCSC Genome Browser API. You can find the import queries, SeqKit commands and Python code for all three charts in BioseqDB/examples/seqkit_workflows.

Please check my OS project - BioseqDB.

You can also follow BioseqDB on LinkedIn.

Read on Medium.

Download the SQL, SeqKit commands and chart scripts.

Comparing butterfly genomes with BioseqDB, SQLAlchemy and Plotly

BioseqDB — Sun, 18 Oct 2026 16:20:44 +0000

Two butterflies can look remarkably similar without belonging to the same species. For Heliconius, part of the story why that happens, involves optix, a gene whose activity helps determine where red patterns appear on the wings. It encodes a transcription factor - so a protein which helps control the activity of other genes, that in very simplified terms works like sompphisticated switch, because changes in regulatory regions can affect not only if, but where and when a gene is used. Reed and colleagues studied this connection across Heliconius species.

I decided to write this tutorial to showcase very simple example of analysis that one can do for single selected organism. As similar topics were always an interest of mine, I dedicated few last years to create a software that combines power of usual relational databases and software used commonly in compuatational biology to perform alignments (one can view them as problem of searching text in the database and understanding text similarity in specific biological way).
In this tutorial we will align a piece of chromosome or scaffold from two Heliconius species that includes optix. We will fetch the sequences from NCBI, keep them in PostgreSQL, ask BioseqDB (the database software I created and do little self-promotion :)) to find matching regions, and use SQLAlchemy and Plotly to inspect the results in Python. We won't need any separate alignment executable or tools.

I inspired myself by excellent blog post: butterfly genome-alignment tutorial using Minimap2 and R. I wanted to work through the same kind of comparison with BioseqDB. This version uses newer NCBI assemblies and their own gene annotations, so the coordinates and results below will be different obviously!

You can find all the full code, SQL, dependencies and Python notebook in BioseqDB/examples/butterfies_alignments. You can also find this tutorial’s code snapshot on the website, but please note that neither provide FASTA files as I rely on NCBI datasets.

1. Very short biological introduction on what we will be doing

A DNA sequence is commonly written with four letters: A, C, G and T, representing adenine, cytosine, guanine and thymine. These are the bases of its nucleotides. A sequence file may also contain N, meaning the base is unresolved. Removing those letters would move every coordinate after them, so we keep them for simplicity.

DNA has two complementary strands - A pairs with T, and C with G and the strands also run in opposite directions. To read the same region from the other strand, we need to reverse the sequence and complement each letter. For example, the reverse complement of AACGT is ACGTT. I mentioning this, because this detail will matter quite a lot in our butterfly plots.

A genome is an organism's complete set of DNA, and an assembly is a reconstruction of it from sequencing data collected from some instruments (this process goes into "wet lab" so you ask your fellow collegues to run Illumina machine, which yields many many parts of the sequence you are interested in along with some noise. I won't go into details now how exactly this is done, but we have our sequence already in the dataset!).
A chromosome record can represent a reconstructed chromosome, while a scaffold is an assembled sequence whose placement may still be incomplete. We will be usiong record for two species - our H. melpomene rosina record is chromosome 18, while the H. erato demophoon record is an unplaced scaffold named ntLink_18

Now the core of any genomic/metagenomic analysis pipeline is the problem of alignment. An alignment is a process that places two sequences in columns so that corresponding bases (nucleotites A,C,T,G) can be compared. Equal letters are matches, different letters are substitutions, and gaps mean insertions or deletions relative to the other sequence. This allows us to understand for example where specific nucleotides were changed during evolutionary processes (how mutations arise, where etc.)

We cannot sensibly compare every position in a million-base sequence with every position in another million-base sequence (that is a lot of text, so our alignment algorithm would be ridiciously slow!). Thankfully in BioseqDB's I provided search functionmmsearch. Internally it starts with short shared fragments, connects promising matches into chains, then aligns selected regions in more detail. This approach is called the seed-and-extend idea and it allows more efficient processing so you can align ridiciously long sequences with no issue. If you are interested in writing your own aligner, my Python read-mapper tutorial goes through the algorithms specifically written in very simple Python.

Here we start analysis with a pairwise comparison across two species, caluclating alignemnt of fragments per every pair of sequences. It is a useful first exercise in comparative genomics.

2. Setting things up!

We will import two one-million-base sequence fragments, inspect their lengths and nucleotide counts, calculate GC content, align them, and progressively draw one match, all retained matches, and the neighbourhood of optix.

PostgreSQL via BioseqDB extension stores the sequences and BioseqDB performs the alignment. We will write Python code using SQLAlchemy, that will send SQL and receive rows, which later will be used by Plotly to constrtuct some figures. You can do that same thing with plain SQL client like psql as every step of data processing will happen inside the database itself.

3. Install BioseqDB and check external access

Warning
Availability on 16 October 2026. The published 0.1.0-alpha.1 Docker 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 installation page covers the regular Docker and native packages. For this tutorial, use the external-enabled development build described in the example's setup instructions:

cmake -S . -B build/butterflies -DBIOSEQ_BUILD_EXTERNAL_API_CLIENTS=ON
cmake --build build/butterflies --target postgresql_bioseq postgresql_bioseq_cache

These lines modify an already configured source build to enable the external clients module, which is responsible for directly interfacing with NCBI database. Obviously a workaround at this time, would be directly downloading datasets by yourself and loading FASTA using bioseq.load function.
After installing that build and enabling the extension in a fresh database, check its SQL objects:

SELECT to_regnamespace('bioseq_external') IS NOT NULL AS external_schema,
       to_regclass('bioseq_external.catalog') IS NOT NULL AS catalog,
       to_regclass('bioseq_external.sequences') IS NOT NULL AS sequences;

Result · 1 row

external_schema	catalog	sequences
t	t	t

All three checks should be true. If the tables exist but a query reports permission denied, ask the database owner to grant access to your tutorial role:

GRANT USAGE ON SCHEMA bioseq_external TO butterfly_reader;
GRANT SELECT ON bioseq_external.catalog, bioseq_external.sequences
TO butterfly_reader;

Replace butterfly_reader with your existing database role.

4. Search NCBI via SQL

We will start with some metadata to make sure we refer to correct organism etc. The following query will use NCBI provider in BioseqDB and send requests to puiblic NCBI API to search organisms for you (all in SQL - how amazing it is instead using hard to understand CLIs?)

SELECT accession, scientific_name, assembly_level, assembly_status
FROM bioseq_external.catalog
WHERE provider = 'ncbi'
  AND scientific_name = 'Heliconius erato demophoon'
ORDER BY accession;

Result · 2 rows

accession	scientific_name	assembly_level	assembly_status
GCA_018249695.1	Heliconius erato demophoon	Contig	current
GCA_056362945.2	Heliconius erato demophoon	Scaffold	current

The provider conditionis important as it selects NCBI and BioseqDB supports manyh other public databases of sequences. Field scientific_name selects the organism, in this case our butterfly species and naturally all seasrch fields are compatible what official NCBI API provides. Field accession is the assembly identifier and it's common name for most of public data sources. The last important one - assembly_level - is a field that distinguishes such things as scaffold and chromosome assemblies (there may be different kind of objects in the database so we would like to see what's there!)

When you run this query, you can probably get different results as the database is updated. For this comparison, I selected these two assemblies and sequence windows:

Organism	Assembly	Sequence record	Imported interval, zero-based
H. erato demophoon	GCA_056362945.2	JBTWUK020000023.1, unplaced scaffold	`[3700000, 4700000)`
H. melpomene rosina	GCF_056956725.1	NC_140406.1, chromosome 18	`[600000, 1600000)`

Each interval contains one million bases, or 1 Mb (in this notation kb is one thousand bases and 1Mb means sequence is one million A/C/T/G/maybe others characters-long).

We use zero-based, half-open coordinates throughout: [start, end) includes start and excludes end, so its length is end - start. NCBI GFF annotations use one-based inclusive coordinates; convert them with start0 = start1 - 1, leaving the end unchanged. NCBI documents the GFF convention here.

5. Now the storage!

Here is our own table that we can create as in ordinary SQL you should be familiar with. USING bioseq gives it BioseqDB's sequence storage and search support. The other fields are tehre just to save the information needed to interpret a coordinates later.

CREATE SCHEMA IF NOT EXISTS butterflies;

CREATE TABLE IF NOT EXISTS butterflies.regions (
    key text NOT NULL,
    accession text NOT NULL,
    organism text NOT NULL,
    assembly text NOT NULL,
    start0 bigint NOT NULL,
    end0 bigint NOT NULL,
    seq bioseq.seq('dna5', 64) NOT NULL
) USING bioseq;

Field	Meaning
`key`	Short name used in this example: `erato` or `melpomene`.
`accession`	Versioned NCBI sequence record, not the assembly identifier.
`organism`	Full scientific name, including the subspecies.
`assembly`	Assembly accession used to select the matching annotations.
`start0`, `end0`	Boundaries of the imported window on the original record.
`seq`	The actual DNA, stored as `bioseq.seq('dna5', 64)`

Oner important gotcha: you might be familiar with SQL VARCHAR where the number signifies how long the text can be. Here it's not the maximum bound on sequence length, but because of my decisions, this it how the sequence is encoded. It can be encoded on units of 8-bits, 16-bits, 32-bits or 64-bits and by default I recommend using the 64. If you try to plug here some other number you will gen an error.

Now we import the two intervals. The external start0 and end0 filters request the selected part of each accession sequence. Beatuful thing is that every step of the process so far happens inside the database, so you don't have to import anything in Python or use BioPython at all. Th SQL is self-explanatory (at least this is what I think - we just do previous query but along with INSERT INTO statement so it's persisted in the table we created):

INSERT INTO butterflies.regions
    (key, accession, organism, assembly, start0, end0, seq)
SELECT 'erato', 'JBTWUK020000023.1', 'Heliconius erato demophoon', 'GCA_056362945.2', 3700000, 4700000,
       sequence::bioseq.seq('dna5', 64)
FROM bioseq_external.sequences
WHERE provider = 'ncbi'
  AND accession = 'JBTWUK020000023.1'
  AND start0 = 3700000
  AND end0 = 4700000;

INSERT INTO butterflies.regions
    (key, accession, organism, assembly, start0, end0, seq)
SELECT 'melpomene', 'NC_140406.1', 'Heliconius melpomene rosina', 'GCF_056956725.1', 600000, 1600000,
       sequence::bioseq.seq('dna5', 64)
FROM bioseq_external.sequences
WHERE provider = 'ncbi'
  AND accession = 'NC_140406.1'
  AND start0 = 600000
  AND end0 = 1600000;

INSERT ... SELECT copies the returned sequences into our table. Once imported, the composition and alignment queries use local data.
Let's inspect a short prefix to see what's there out of curiousity:

SELECT key, accession, bioseq.len(seq) AS length,
       seq[0:20]::text AS first_20_bases
FROM butterflies.regions
ORDER BY key;

Result · 2 rows

key	accession	length	first_20_bases
erato	JBTWUK020000023.1	1000000	TCGCGTGAATGCAAATGCGG
melpomene	NC_140406.1	1000000	ATATCAAATACCAATGTTAG

bioseq.len(seq) gives the number of bases. seq[0:20] takes the first twenty, using the same half-open convention as our genomic windows. A coordinate returned by an alignment will be relative to this imported sequence. To get back to the NCBI record, add its stored start0.

6. Let's do some exploration and count bases!

Two records can have similar overall composition and still contain very different sequences. I don't want to read entire million of bases as this would crowd up our screen, so lets stick to some meta properties, like the sequence lengths and counts for every nucleotide:

WITH counts AS (
    SELECT key, accession, organism, start0 AS interval_start,
           bioseq.len(seq) AS length,
           bioseq.symbol_count(seq, 'A') AS a,
           bioseq.symbol_count(seq, 'C') AS c,
           bioseq.symbol_count(seq, 'G') AS g,
           bioseq.symbol_count(seq, 'T') AS t,
           bioseq.symbol_count(seq, 'N') AS n
    FROM butterflies.regions
)
SELECT *, round(100.0 * (g + c) / NULLIF(a + c + g + t, 0), 4)
          AS gc_percent
FROM counts
ORDER BY key;

Result · 2 rows

Row 1
key: erato
accession: JBTWUK020000023.1
organism: Heliconius erato demophoon
interval_start: 3700000
length: 1000000
a: 329927
c: 172322
g: 170873
t: 326878
n: 0
gc_percent: 34.3195

Row 2
key: melpomene
accession: NC_140406.1
organism: Heliconius melpomene rosina
interval_start: 600000
length: 1000000
a: 327164
c: 171599
g: 171833
t: 329404
n: 0
gc_percent: 34.3432

bioseq.symbol_count counts a selected letter and GC content is the percentage of called bases which are G or C:

GC% = 100 \frac{G + C}{A + C + G + T} .

We exclude N from the denominator as N represents an unresolved base and in practice it can be anything of A/C/T/G, but we never know for sure. Usually Ns are small part of the sequence overall and we are not interested in them so much. NULLIF here avoids dividing by zero if a sequence contains no called bases. That case is unknown GC content, not zero GC. Lucily for us, both windows we saved and used here happen to contain no Ns (hooray!)

Both analysis windows are exactly 1 Mb. The source erato scaffold is 5,549,739 bases, the source melpomene chromosome is 16,989,447 bases.

A and T are more frequent than C and G in both windows. Erato contains 329,927 A, 172,322 C, 170,873 G and 326,878 T bases. Melpomene contains 327,164 A, 171,599 C, 171,833 G and 329,404 T bases.

The GC values are close: 34.3195% and 34.3432%, but note that this isn't reflecting their similarity as we don't measure here any position-based properties.

More interesting is measuring GC content in windows of specified length, because they we have more than only one single data point per entire sequence of interest (they are long so any mutations will average to very similar total GC content percentage. Using smaller windows on sequences allow us to see some regions where "something happens" and the structure of sequence changes - all this without actually comparing all nucleotides themselves on our screen - in practice this can be used as very primitive tool to quickly search for regions of interest)

WITH windows AS (
    SELECT r.accession, w.start,
           least(w.start + 10000, bioseq.len(r.seq)) AS "end",
           r.seq[w.start:least(w.start + 10000, bioseq.len(r.seq))]
               AS piece
    FROM butterflies.regions AS r
    CROSS JOIN LATERAL generate_series(
        0, bioseq.len(r.seq) - 1, 10000) AS w(start)
), counts AS (
    SELECT accession, start, "end", bioseq.len(piece) AS length,
           bioseq.symbol_count(piece, 'G') + bioseq.symbol_count(piece, 'C') AS gc,
           bioseq.symbol_count(piece, 'N') AS n
    FROM windows
)
SELECT accession, start, "end",
       round(100.0 * gc / NULLIF(length - n, 0), 4) AS gc_percent,
       1.0 * n / length AS n_fraction
FROM counts
ORDER BY accession, start
LIMIT 6;

Result · 6 rows

accession	start	end	gc_percent
JBTWUK020000023.1	0	10000	33.6800
JBTWUK020000023.1	10000	20000	34.1000
JBTWUK020000023.1	20000	30000	35.7900
JBTWUK020000023.1	30000	40000	35.7600
JBTWUK020000023.1	40000	50000	33.2900
JBTWUK020000023.1	50000	60000	35.6600

These curves show where the composition changes within each imported region - peaks/valleyhs or other interesting shapes are not, by themselves, genes, enhancers or anything else and we need the appropriate annotations or experiments in "wet labs" to check, but nevertheless it's some information we can use in our analysis.

7. Align the regions inside PostgreSQL

We will use melpomene as the target and erato as the query. Both come from the same local table, selected by their key fields. This listing shows the first eight result (06_align.sql in the repository omits the final LIMIT and returns every block. I used LIMIT here just for simpler illustrative purposes)

SELECT t.accession AS target_accession,
       q.accession AS query_accession,
       bioseq.target_start(t.seq) AS target_start,
       bioseq.target_end(t.seq) AS target_end,
       bioseq.query_start(t.seq) AS query_start,
       bioseq.query_end(t.seq) AS query_end,
       bioseq.is_reverse(t.seq) AS is_reverse,
       bioseq.score(t.seq) AS score,
       bioseq.cigar(t.seq) AS cigar,
       bioseq.edit(t.seq) AS edits,
       bioseq.block_length(t.seq) AS block_length,
       bioseq.mapq(t.seq) AS mapq,
       bioseq.is_primary(t.seq) AS is_primary
FROM butterflies.regions AS t
JOIN butterflies.regions AS q
  ON mmsearch(t.seq, q.seq, aligner_options(
      kmer_size => 15, window_size => 10,
      max_gap => 200, max_gap_ref => 200, min_chain_score => 40))
WHERE t.key = 'melpomene' AND q.key = 'erato'
ORDER BY score DESC, target_start, query_start
LIMIT 8;

Result · 8 rows

Row 1
target_accession: NC_140406.1
query_accession: JBTWUK020000023.1
target_start: 639721
target_end: 644608
query_start: 332091
query_end: 337300
is_reverse: t
score: 5874
cigar: 27M4D14M2I9M2I39M1I29M3I15M1I3M1D8M3I53M4D3M1D12M1D24M7D2M1D8M1I3M4I15M1D14M8I89M2D78M16D5M5D3M4D4M6D1M4D6M1I13M2D1M2D1M1D32M4D15M1I33M2I6M8I5M3I1M1I17M3D17M12I4M7I5M1I10M7I2M4I6M2D5M1I6M5I12M5I5M5I4M4I10M6I5M11I2M8I1M4I5M10I4M12I3M41I3M1I2M5I4M1I3M1I1M1I114M6I41M3I13M5I89M1I29M1I9M3I9M5I6M3I11M1D1M2D5M1D4M1I20M2I2M1I21M2I3M1D47M4D71M3I3M2I46M2D101M3D67M6D65M2I103M12I7M1I4M1I65M60I5M19I6M32I3M37I5M3I1M3I26M8I4M1I21M1I29M1I2M2D36M2I37M1I687M3I2217M1D1M1D3M1D15M1D1M2I12M1D15M
edits: 894
block_length: 5307
mapq: 60
is_primary: t

Row 2
target_accession: NC_140406.1
query_accession: JBTWUK020000023.1
target_start: 243317
target_end: 247028
query_start: 793769
query_end: 797544
is_reverse: t
score: 5190
cigar: 32M3I11M1I2738M3I21M23D3M2D4M18D1M5D1M6D5M6D5M15D13M6D11M2I11M3I3M8I25M2I125M4D8M1I17M1I4M1D27M3I14M22I3M18I12M5I3M5I4M1D11M3I1M11I3M5I6M4I2M3I1M10I4M1D1M3I5M1I1M1I2M4I5M6I1M6I10M1I2M1I9M15I2M10I22M1D3M3I13M1I1M1I4M1D6M1D59M1I8M5D12M2D34M2D7M1I5M2I36M6D228M
edits: 484
block_length: 3881
mapq: 60
is_primary: t

Row 3
target_accession: NC_140406.1
query_accession: JBTWUK020000023.1
target_start: 436485
target_end: 440921
query_start: 582051
query_end: 586426
is_reverse: t
score: 4902
cigar: 55M2I4M2D399M1I2M1D977M1D36M4I6M1D2M3D34M2D2M5D6M2D4M3D21M4I32M2D1M2D13M1I3M1I30M4I44M2I3M2D13M1D77M1I4M2D6M1D2M3D29M1D5M1D6M2D10M2I15M1I1M2I6M2D11M1I9M2D20M2I2M3I5M1D4M1D2M1I1M1I3M1D3M2I257M2I9M3I37M1D27M3D89M3D1M7D13M5D13M12D3M1I14M1D42M2I102M2D4M2D21M1D39M6I18M6D5M2D2M1D5M3D6M18D87M1D5M2D5M11I11M12I1191M3D50M4I9M7I5M2I24M5I1M3I72M2I2M2I5M3I11M2I2M1I1M1I12M2I4M4D4M2D3M1D4M3D2M2D5M1D20M6D2M4D3M15D2M4D5M8D107M
edits: 728
block_length: 4542
mapq: 60
is_primary: t

Row 4
target_accession: NC_140406.1
query_accession: JBTWUK020000023.1
target_start: 735657
target_end: 740842
query_start: 214190
query_end: 219875
is_reverse: t
score: 4868
cigar: 5M2D10M4I194M3I4M10I6M1I30M10I3M9I7M36I5M62I4M55I4M8I2M8I56M7I8M1I58M2I71M4I28M8I20M1D12M8D15M2I1M1D19M1I16M1D49M17D3M5D3M9D6M5D8M9D2M5D2M1D62M6I5M1I11M3I2M1I2M7I2M7I6M9I11M2I2M1I5M4I17M5I3M7I4M2I4M1I2M2I9M1I7M2D3M3D11M5I4M1I1M2I7M7I4M3I5M12I4M1D8M21I115M1I1M3I94M1I235M6D44M2I16M18I2M1I2M24I62M1I5M5D4M1D4M5I2M1I25M1D34M1I48M2I8M2I5M1D49M1D4M3I36M1I114M1D91M1I10M1D210M3I38M1I16M1D30M2I5M4D10M4D49M1I252M3I31M2I150M1I14M2I138M3I5M1I151M3I9M1I8M12I34M10I1M2I22M2I93M2D1M1D22M1I35M4I27M2I25M1I140M1D33M2I12M1D73M1D63M2I51M3I23M9I2M1I4M2I44M1D102M5I25M4I2M5I4M10I55M14I10M6I4M3I40M2D5M1D6M4D3M1I5M3I23M1D9M2D6M1I5M1I3M3I8M1I3M2D4M2D4M5I5M2I2M1I2M1I3M1I7M6I9M5I1M5I4M3I13M4D14M2D4M3I13M2D5M1D88M6I44M1I1M2D31M8D132M1I94M6I90M17D4M1D41M4I3M4I12M7I4M7I63M3D8M1I2M1I16M1D21M3I116M2D58M1I14M11I2M1I7M1D3M21I24M1I3M1I9M14I16M3D24M
edits: 1297
block_length: 5849
mapq: 60
is_primary: t

Row 5
target_accession: NC_140406.1
query_accession: JBTWUK020000023.1
target_start: 782603
target_end: 787530
query_start: 156816
query_end: 162219
is_reverse: t
score: 4858
cigar: 20M3I41M7I8M4I26M1D3M1I10M1I2M3D38M1D222M1I134M1D16M1I98M3D2M2D5M2D9M2D81M2I8M1D1M1D9M1I2M2I7M1I8M4I35M11I5M2I2M3I3M5D5M3D7M10I15M1I2M2I9M2I5M10I5M3I15M1D2M4I79M7I215M1I46M9I19M1I19M1I44M1D61M5D145M3I57M1I6M1D3M1I6M1D2M1I50M18I9M1I1M1I24M12I4M7I2M2I4M4I4M1I65M8I6M2I7M2D5M1I3M1I69M2D1M2D148M6I165M10I38M1D9M3I1M2I33M1D52M2D33M4I22M1I14M1D16M4D109M7I1M1I3M1I12M8I1M1I16M4I32M1D434M2I80M1I12M3I23M1D70M1I1M1I37M2D60M2D278M1I73M16I99M1D49M1D7M1D11M3I4M1I3M1I2M1I4M1I2M1D5M6I6M1D6M6I3M1I11M3I9M2I4M3I5M5I1M4I7M3I8M3D6M3D2M1D10M2D4M3I17M6I10M2I2M1D3M5I2M2I1M2I8M2I21M1I2M1I1M4I6M2I1M3I12M1I7M4I5M1I7M1I6M5I5M1I3M1I2M11I1M1I6M9I154M14I5M2D3M4D56M1D14M3I2M3I7M2I2M2I2M2I3M11I3M5I3M13I13M7I1M6I13M4I4M2I10M1I11M2I2M10I11M3I2M1I4M18I7M6I2M7I2M1I4M5I2M3I4M3I2M4I1M1I2M2I10M8I9M9I1M2I4M26I4M6I2M1I1M1I3M1D2M7I11M1I25M1I1M1I11M1D29M4D257M1I34M7I2M1I43M6I59M
edits: 1109
block_length: 5486
mapq: 60
is_primary: t

Row 6
target_accession: NC_140406.1
query_accession: JBTWUK020000023.1
target_start: 123169
target_end: 127102
query_start: 931912
query_end: 935729
is_reverse: t
score: 4824
cigar: 23M4D1M1D3M2D14M1D21M1D7M3D45M1I15M3D7M1D19M3D2M13D2M4D47M1I13M3D56M2I23M1D43M4I9M9I90M1I86M1I179M1I2M1I14M7D2M1D678M1I1M1I32M1D14M3I4M1D9M1I7M1I76M7I21M3I29M1I62M3I7M4I2M3I30M4D21M1I66M16D5M7D8M8D6M40D4M19D347M1D120M1I1M1I36M2I7M2D5M4D9M1D7M2D2M6D61M8D61M2D15M1D204M5I5M2I64M3D1M1D726M2I15M1D94M1I34M1I5M1D1M1D11M1I2M1I2M2I85M1I4M3D2M4D2M1D14M
edits: 567
block_length: 4003
mapq: 60
is_primary: t

Row 7
target_accession: NC_140406.1
query_accession: JBTWUK020000023.1
target_start: 421610
target_end: 427235
query_start: 593427
query_end: 599306
is_reverse: t
score: 4738
cigar: 35M2I9M2D36M1D4M1I1M1I76M1I10M9I99M1D30M2I18M5D562M1D95M4D54M9D8M1D1M1D10M2I4M3I19M3D3M1D3M1D46M1D10M14D129M1I3M1D25M8D3M1D475M1D21M1D70M1D4M2I2M2I16M3I81M1I38M9I9M1D4M2I39M4I166M3D1M5D33M2D46M2D6M6D1M1D7M1I4M1D2M1D12M1D10M5I6M5D12M1D3M4D3M3I10M4D3M1D2M4I11M1I3M4D6M1I3M2D3M1D6M2D2M11D72M1I10M3I21M2I78M4D62M3D19M3I8M9I8M14I2M3I62M33I4M24I7M36I5M3I1M4I3M19I3M22I4M32I2M1I2M10I7M4I2M13I3M11I7M20I5M38I18M1I4M2I992M1D7M3I8M1I21M2D6M11I3M3I31M1D1M1D4M2D12M2D1M1I2M2D50M10D14M7D19M20D7M5D1M2D4M7D3M2D4M4D115M1I1M3I20M2I4M2I76M2D1M1D1M1D13M4I9M6D1M5D14M8D4M15D5M1D3M2D3M3I6M3I8M2D6M1I2M2D5M3D6M1I2M7D13M1I56M11I10M6I88M8I13M5I9M4I9M3I6M4I5M1I3M1I7M3I5M7I5M9I3M11I5M22I9M3I9M8I4M6I4M13I7M1D3M3I16M1I3M3I2M3I2M2I1M4I3M4I3M8I4M6I2M1I2M4I8M3I3M1I16M1I8M1I1M1I2M1I3M2I5M1I4M2I11M21I3M4I3M1I3M2I3M14I2M2I6M1D5M1D2M1I4M11I238M2D6M16D4M23D1M6D1M5D3M1D3M19D6M5D9M5D3M10D1M2D5M18D9M1D2M2D3M4D5M9D6M5D4M4D4M1D195M
edits: 1509
block_length: 6266
mapq: 60
is_primary: t

Row 8
target_accession: NC_140406.1
query_accession: JBTWUK020000023.1
target_start: 409490
target_end: 413107
query_start: 608855
query_end: 612268
is_reverse: t
score: 4586
cigar: 68M30I1M1I2M1I7M1I24M3D8M1D772M7D15M1I3M4D11M6D8M1I10M7D92M2D13M6D23M4I10M3I63M3D8M2I69M1I1M1I5M4D24M2D8M16D227M2I2M4I947M19D3M1D6M2D4M1D42M4D2M4D2M45D6M49D4M37D3M3D1M1D3M10D6M23D55M1D17M1I2M1I5M2D54M1I39M4I679M
edits: 539
block_length: 3676
mapq: 60
is_primary: t

mmsearch(t.seq, q.seq, ...) is SQL function you can use to produce local sequence matches (perform sequence alignment). Here kmer_size => 15 uses 15-base seeds and window_size => 10 controls their minimizer windows (I don't want to go too deep into details, but a k-mer is a run of k consecutive bases and a minimizer is one representative k-mer selected from a window by a fixed ordering, which reduces the number of seeds we keep. The longer k-mers can increase performance by reducing sensitivity, but here we use BioseqDB's default k-mer size 15 and window size 10, so we will stick to them. They are presented here to let you know what you can change. These two values can be omitted when the session defaults are unchanged, but keep the gap options below to reproduce these results).
The gap and minimum chain-score options affect which candidate chains are explored and can do the same - improve performance at the cost of sensitivity. Here max_gap => 200 differs from the default 5000; max_gap_ref is accepted but is not yet used by the current chaining implementation.

The bioseq.* projections describe the match we get from aligning oft.seq in the current result row:

Projection	What we use it for
`target_start`, `target_end`	Half-open interval on the imported melpomene sequence.
`query_start`, `query_end`	Half-open interval on the imported erato sequence, in its deposited orientation.
`is_reverse`	Whether the query aligns as its reverse complement.
`cigar`	Base-level alignment operations: for example, `M` consumes a base from each sequence, `I` from the query, and `D` from the target.
`edit`	Number of substitutions and inserted/deleted bases in that alignment.
`block_length`	Alignment-column count, including gap columns.
`score`	Alignment score under the scoring scheme; not a percentage.
`mapq`, `is_primary`	Estimated placement confidence and primary/secondary information used when inspecting competing matches.

Mapping quality, or MAPQ, concerns confidence in a placement among alternatives. It is different from the fraction of identical bases in an alignment, it says us how many different regions where there for some part of the query, so a high MAPQ does not mean the two sequences are identical, but means that for each region we were aligning there weren't many regions that can match as well so result is likely around correct place.

Requesting bioseq.cigar asks BioseqDB for the base-level alignment rather than only an estimates of positions (other base-level projections, including bioseq.edit, also request this step; without them the start/end of query/target might change! ASking for cigar improves the accuracy to the positions of individual bases). In a CIGAR string, M consumes a base from each sequence and can include both matches and mismatches; I marks inserted query bases and D marks deleted target bases. The numbers are run lengths. Use bioseq.edit() to retrieve the number of substitutions and inserted/deleted bases.

The validated run returned 661 blocks, of which 592 had reverse query orientation (the aligners consider reverse strand which means that we not only search for some sequence X, but also it's reverse-complement. In very simple terms this is because of how DNA work and the reverse complement of ACTG is CAGT). For the plots we keep blocks with at least 300 alignment columns, mapping quality at least 20 and identity at least 70%. That leaves 355 blocks, including 342 reverse matches.

8. Set up Python with Poetry

Run these commands in a terminal with Python 3.11–3.13 and Poetry installed:

git clone https://github.com/BioseqDB/examples.git
cd examples/butterfies_alignments
poetry install
export DATABASE_URL='postgresql+psycopg://postgres:bioseq@localhost:55432/butterflies'
poetry run butterfly-tutorial

Use the connection details of your external-enabled database. The password shown here is only a local example. poetry install installs SQLAlchemy, psycopg, Plotly, Kaleido and JupyterLab. If starting a separate project instead, the corresponding dependency command is:

poetry add sqlalchemy 'psycopg[binary]' plotly kaleido jupyterlab nbclient

Static SVG export uses Plotly with Kaleido and requires Chrome Plotly explains that dependency here. The database analysis can run without image export using --no-figures. After the first successful import, you can use --reuse to repeat the analysis without fetching those sequences again.

For the notebook interface:

poetry run jupyter lab notebooks/butterfly_alignments.ipynb

9. Mastering SQL now

You don't have to manually enter SQL into the database, you can do that the same thing from Python via the industry-standardised SQLAlchemy library:

import os
from pathlib import Path
from sqlalchemy import create_engine, text

root = Path.cwd()
if not (root / "sql").is_dir():
    root = root.parent
sql_dir = root / "sql"
engine = create_engine(
    os.environ["DATABASE_URL"],
    connect_args={"prepare_threshold": None},
)

statement = text((sql_dir / "04_inspect.sql").read_text())
with engine.connect() as connection:
    preview = connection.execute(statement).mappings().all()
for row in preview:
    print(dict(row))

I use create_engineto tell SQLAlchemy how and where to connect and later the with engine.connect() block acquires a connection and returns it to the pool when finished. .mappings() gives rows with column names, which is convenient for tables and figures. We disable psycopg's automatic prepared statements here because the current external-query implementation needs the accession and range values available in its custom query plan.

The first two returned rows are:

{'key': 'erato', 'accession': 'JBTWUK020000023.1', 'length': 1000000, 'first_20_bases': 'TCGCGTGAATGCAAATGCGG'}
{'key': 'melpomene', 'accession': 'NC_140406.1', 'length': 1000000, 'first_20_bases': 'ATATCAAATACCAATGTTAG'}

For a changing value, bind a parameter instead of building SQL with string interpolation:

statement = text((sql_dir / "01_discover.sql").read_text())
with engine.connect() as connection:
    assemblies = connection.execute(
        statement, {"organism": "Heliconius erato demophoon"}
    ).mappings().all()
print([dict(row) for row in assemblies])

If you want to read more about SQLAlchemy transactions and what is engine.connect() you are welcome to read more here: SQLAlchemy's tutorial describes these transaction patterns.

To repeat the imports from Python, use the transactional helper:

import json
from butterfly_tutorial.workflow import import_regions

sources = json.loads((root / "sources.json").read_text())
import_regions(engine, sources["records"])

The alignment is SQL, so it's another query through the same connection (nothing special here really!)

statement = text((sql_dir / "06_align.sql").read_text())
with engine.connect() as connection:
    matches = [dict(row) for row in connection.execute(statement).mappings()]
print(f"{len(matches)} alignment blocks")

The matches list contains the same 661 rows as the SQL query. The full runner in the GitHub repo also validates lengths, nucleotide counts, sequence hashes, CIGAR operations and the optix coding sequences before, but I'm not going to present it all here as this is just some additional boilerplate not very relevant to the analysis itself.

10. First drawing!

There's a think called ribbon diagram, where a sequence track is just a rectangle whose x coordinates are sequence positions and its height and y position are display choices:

An alignment connects two intervals with a four-corner polygon. For a forward match, the target's left endpoint corresponds to the query's left endpoint. For a reverse match, it corresponds to the query's right endpoint. No rocket science here, right? Only butterfly science! Huh.

Here is the small plotting operation behind a ribbon:

import plotly.graph_objects as go

hit = matches[0]
left, right = hit["query_start"], hit["query_end"]
if hit["is_reverse"]:
    left, right = right, left
figure = go.Figure(go.Scatter(
    x=[hit["target_start"], hit["target_end"], right, left],
    y=[1, 1, 0, 0], fill="toself", mode="lines",
    line={"color": "#65b9dd"},
    fillcolor="rgba(101,185,221,0.25)",
))
figure.update_layout(xaxis_title="Position in imported window (bp)",
                     yaxis={"tickvals": [0, 1], "ticktext": ["Erato", "Melpomene"]})
figure.show()

The figures on this page are static images, but the notebook and generated HTML versions also let you zoom and hover over a match to inspect its coordinates and it feel much cooler to see them update after your code changes. You can give it a try.

Its target interval is [639721, 644608) and its query interval is [332091, 337300), with reverse orientation. The spans need not have equal lengths because the alignment includes insertions and deletions. This block has 5,307 alignment columns and 894 edits, giving about 83.15% identity.

11. We need more plots! Mooore!

One plot won't satisfy my needs for colourful graphs. We can use a dot plot. A dot plot is a useful first view of many matches - each line segment plots target position against query position and its slope shows orientation.

The broad reverse trend agrees with the opposite orientations of the annotated optix regions in records, but it does not by itself mean a biological inversion between the species and somehow they reversed everyrthing just by randon accident - an assembly record can be deposited in either orientation to NCBI so it's likely that.

Now draw all retained ribbons in the original coordinate systems:

Overlapping blocks are possible, especially around repeated sequence, so summing their lengths would overstate coverage. After merging overlapping intervals, the displayed blocks cover 374,441 target bases and 387,528 query bases. These are the spans touched by the selected alignments.

For a closer look at optix, we centre each gene on zero and reverse the erato display axis. In local window coordinates, the gene midpoints are 493,223.5 for melpomene and 495,840 for erato. The display transformations are therefore x = target_position - 493223.5 and x = -(query_position - 495840).

We need to make a choice how to center the genome and it's compltely arbitrary - their annotated transcript spans differ, so corresponding coding bases need not fall at exactly the same x position. This transformation makes the comparison easier to read. In the interactive notebook and HTML versions, the hover text retains the original local coordinates. When selecting a zoom, we keep blocks that overlap it, including blocks crossing its boundary, requiring both endpoints to lie inside the zoom would silently discard useful matches. If you have better idea how to draw this, go ahead!

12. Add gene models

The erato feature AC2CFC_001707 is annotated as homeobox protein SIX6; I treat it as the optix homolog here because its coding sequence shares 778/804 bases and its protein 266/267 amino acids with melpomene Optix.

The annotation loader on GitHub downloads the exact GFF files listed in sources.json, selects the appropriate accession, converts coordinates and clips features to our imported windows. There are 35 gene or pseudogene models overlapping the erato window and 38 overlapping the melpomene window, including one pseudogene in each. A pseudogene is a gene-like sequence annotated as no longer encoding the usual functional product. A gene model is an annotation of a gene's structure.

The notebook assembles these views and make_figures returns a dictionary of Plotly figures.

from butterfly_tutorial.workflow import analyze
from butterfly_tutorial.annotations import load_annotations, load_optix_exons
from butterfly_tutorial.plots import make_figures

analysis = analyze(engine, sources["records"])
annotation_cache = root / "output" / "downloads"
genes = load_annotations(sources["records"], annotation_cache)
exons = load_optix_exons(sources["records"], annotation_cache)
figures = make_figures(
    analysis["sequences"], analysis["plotted_matches"],
    analysis["gc_windows"], genes, optix_exons=exons,
)
print(json.dumps(analysis["summary"], indent=2))

The exon panel is a closer view of the annotations. Exons are the parts retained in a mature transcript, The annotated erato transcript has two exons and the melpomene transcript has one in these files and I put everything into one simple graphic. The gene here tracks place the matching segments in their annotated neighbourhood. The ribbons show sequence resemblance under the alignment and the last one - GC tracks add a simple measurement of GC content made directly from the imported DNA.

13. What did we learn?

The two regions have almost the same overall GC content, yet their alignment is a collection of local blocks with gaps and competing matches. The optix coding sequences are especially similar, while the surrounding DNA provides a larger and less uniform comparison. None of these plots alone identifies the regulatory change responsible for a wing pattern, but hey, we got some cool plots!

What I find useful about this workflow is that the sequences, their accession identifiers and the analysis stay all withing SQL. Selecting another region means changing explicit coordinates in your code without the need to look for CLI docs in the internet or using different tools blues together via bash scripts. You can even search NCBI via SQL layer that directly uses their API. I didn't show it here but you can export the table into usual FASTA as well.

To continue, try changing the display thresholds and compare the dot plot before interpreting an apparent break. Then select another gene neighbourhood, preserving the accession versions and coordinate conventions. Adding more genomes can turn this small pairwise exercise into a broader comparative analysis. with even more cool charts!

Please check my OS project - BioseqDB.

Read this post on BioseqDB · Read on Medium · Follow BioseqDB on LinkedIn.

Introducing BioseqDB

BioseqDB — Sun, 18 Oct 2026 16:17:07 +0000

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.

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.