Working with sequence files in BioseqDB and SeqKit

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

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.

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

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

INSERT 0 1 — 1 row affected

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 

49 

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 Source：https://dev.to/bioseqdb/working-with-sequence-files-in-bioseqdb-and-seqkit-1df2