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
| 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
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 % 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.
Top comments (0)