Comparing butterfly genomes with BioseqDB, SQLAlchemy and Plotly

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

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

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:

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; 

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

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:

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

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