Background
A core challenge for bioinformatics in metagenomics is that many standard bioinformatics tools are not built to handle datasets of the size of our metagenomic libraries. At the same time, the continued rapid evolution of the field and development of new methods requires that we maintain compatibility with the standard bioinformatic working environment – in other words, if we want to keep up with the state of the field, we must be able to provide command-line access to standard file types (FASTA, GFF, GBK, etc) generated from our data. We also need to be able to subset our data quickly and easily. For example, while we may need or want to operate on the whole library when performing homology searches, once we’ve completed that search, we’ll often want to work with just a subset of the search results for further analysis. In other words, we need: our entire library accessible in a commonly acceptable format, the ability to rapidly subset our data for additional downstream processing, and we need this full set of capabilities available within something that looks like a standard Linux command-line environment. Finally, experience has shown that the closer we can hew to using standard Linux-based tools, the easier it is for our scientists and bioinformaticians to work with our data, the lower the support burden, and the faster we can iterate on methods.
Plainly stated, we need a fast and easy way to store and extract subsets of data from extremely large text files. Let’s investigate some options.
Tools
We’ll be relying on the time command to give us the time a program takes to run, top for memory, and pv to provide us the speed at which a file is being read. time and top are standard Linux tools, and pv is available via the pv package in most distributions.
Setup
The specific problem we’re tackling today is providing fast access and retrieval of gene annotations from a GFF file covering an entire metagenomic library. The GFF format was chosen because it’s simple, supports multiple sequences in one file, allows us to store and retrieve just the gene annotations, is widely supported by standard tooling, and can be converted into other formats as needed.
The CA29 Library
We’ll be working with one library for the duration of this post – the CA29 soil bacterial library. We’ve chosen this one because it’s one of our largest libraries – if we’ve got a solution that works for CA29, we’ve got a solution that works for everything. If you’d like more information about our libraries, our recent paper Precision Discovery of Novel Inhibitors of Cancer Target HsMetAP1 from Vast Metagenomic Diversity includes details on their construction and use.
Some quick stats on the library, for reference:
- Sequences: ~150M
- Annotations: ~1.6B
- Uncompressed Size: 235Gb
- Compressed Size: 24Gb
Note that depending on the context below I may use either GiB or Gb, depending on what the tool I’m using puts out. In both cases, the term refers to 2^30, or 1024^3, bytes.
Goals
Our goal here is to be able to easily retrieve the gene calls for an individual sequence from our combined library. We’d like to do this quickly enough for it to be practical as part of other pipelines, which we’ll define as “less than a minutes to retrieve the GFF for a sequence,” and we’d also like to keep the computational resources down to an arbitrarily-chosen limit of 2 CPUs and 4Gb of memory, as this will fit inside almost any process we run and stands a very low probability of bankrupting the company.
Approaches
Tabix
The standard tool for retrieval of subsets of data from GFF files is tabix, which is part of the htslib package. Tabix works by creating an index file that maps sequence IDs to their positions in the compressed GFF file, and then using this index to provide rapid lookup and retrieval by sequence ID or sequence ID and location. Tabix, however, suffers from two major limitations for our use case – first, as of the latest version, Tabix has a bug which causes it to crash if the total length of all of the sequence IDs in the GFF file it is indexing exceeds 2^31 bytes. We use UUIDs for sequence names, and so this imposes a limit of around 50 million sequences, which is on the low side for one of our libraries. Second, Tabix’s memory use is extremely high for large datasets – when searching a 12Gb compressed GFF file, Tabix requires 32Gb of memory.
Not all’s lost, though. Under the hood, Tabix uses another utility, BGZip, to facilitate its quick retrieval. BGZip turns out to be a real gem, both in its use and construction.
BGZip
BGZip makes use of a quirk in the widely-used and supported gzip file format – that you can concatenate multiple gzipped files together and the result is still a valid gzip archive. BGZip takes this to an extreme – when compressing a file, it splits the file into small chunks, gzips them independently, and then generates an index file mapping locations from the uncompressed file into the compressed file.
[ long_file.gff ]
↓
[ part 1 ][ part 2 ][ part 3 ][ part 4 ]
↓ Compressed and Indexed ↓
[a][b][c][d] + { part 1: a; part 2: b; ...}
This gives BGZip two neat tricks. One, which we’ll get back to later, is that it can compress and decompress multiple chunks at a time, since each chunk is fixed size and independent of its surroundings. The second is that it can retrieve subsections of the compressed file almost instantly:
# Relative size of the gzipped and uncompressed files:
$ ls -alh
total 259G
drwxrwxr-x 2 ubuntu ubuntu 4.0K Oct 4 21:31 .
drwxr-xr-x 9 ubuntu root 4.0K Oct 4 21:25 ..
-rw-rw-r-- 1 ubuntu ubuntu 235G Oct 4 21:35 ca29.gff
-rwxr-xr-x 1 ubuntu ubuntu 24G Oct 4 21:29 ca29.gff.gz
-rwxr-xr-x 1 ubuntu ubuntu 59M Oct 4 21:29 ca29.gff.gz.gzi
# ('.gzi' is bgzip's index - note it's only 60M for our 24Gb file)
# Retrieving a specific section of the uncompressed file, using `tail` and `head`:
$ time tail -c +251262670684 ca29.gff | head -c 60
##sequence-region fffff5f2-0bcb-4560-8399-c86c632969b8 1 125
real 0m0.001s
user 0m0.002s
sys 0m0.000s
# Retrieving a specific section of the compressed file using bgzip - note
# the start and size are the same as the uncompressed reads above:
$ time bgzip -b 251262670684 -s 60 ca29.gff.gz
#sequence-region fffff5f2-0bcb-4560-8399-c86c632969b8 1 1257
real 0m0.054s
user 0m0.025s
sys 0m0.029s
# Retrieval time does not scale with the size of the chunk retrieved:
$ time bgzip -b 251262670683 -s 62000000 ca29.gff.gz > /dev/null
real 0m0.054s
user 0m0.019s
sys 0m0.035s
So we’ve got the first part of our solution – nearly instantaneous retrieval of a subset of an excessively large gzipped file based on start location and size. By adding a sequence location index next to our GFF file with the sequence ID, start location, and size of the associated entries, we can get almost instant sequence retrieval:
# Our location index file: $ head ca29.gff.loc.tsv seq_id offset length 00000020-3ae5-4305-8de8-1bfceaa81f1c 16 1276 00000035-1d98-4718-960c-e1c0f9fda8c9 1292 937 00000057-0ac5-467f-a37f-3e64cf9928e8 2229 937 0000005c-da49-428c-95e5-342515731bc0 3166 937 00000076-c423-45c5-86b1-52fed94afbbf 4103 937 0000008a-7735-43f7-95fa-64ab84372ca1 5040 1617 000000f9-d5ca-4723-876d-b9942873d8cb 6657 602 0000015b-02d1-4c22-8b2f-8e723961f5fb 7259 939 00000197-45d6-4d40-8ab6-102b3e0000e8 8198 937 # Grab a random sequence from the end: $ tail -n50 ca29.gff.loc.tsv | head -n1 fffffb2f-878d-4b23-a82f-573fbc1fad5f 251262762041 935 # Retrieve the gene calls from the compressed GFF: $ time bgzip -b 251262762041 -s 935 ca29.gff.gz ##sequence-region fffffb2f-878d-4b23-a82f-573fbc1fad5f 1 1093 fffffb2f-878d-4b23-a82f-573fbc1fad5f annotation remark 1 1093 . . . molecule_type=DNA;topology=linear fffffb2f-878d-4b23-a82f-573fbc1fad5f feature region 1 1093 . . . ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f fffffb2f-878d-4b23-a82f-573fbc1fad5f feature CDS 1 87 . - 0 ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_number=1;name=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1 fffffb2f-878d-4b23-a82f-573fbc1fad5f feature exon 1 87 . - . gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_number=1 fffffb2f-878d-4b23-a82f-573fbc1fad5f feature CDS 212 1093 . - 0 ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_number=2;name=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2 fffffb2f-878d-4b23-a82f-573fbc1fad5f feature exon 212 1093 . - . gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_number=2 real 0m0.059s user 0m0.029s sys 0m0.025s
Using grep, we can retrieve the locations for an individual sequence id:
$ grep fffffb2f-878d-4b23-a82f-573fbc1fad5f ca29.gff.loc.tsv fffffb2f-878d-4b23-a82f-573fbc1fad5f 251262762041 935 $ bgzip -b 251262762041 -s 935 ca29.gff.gz ##sequence-region fffffb2f-878d-4b23-a82f-573fbc1fad5f 1 1093 fffffb2f-878d-4b23-a82f-573fbc1fad5f annotation remark 1 1093 . . . molecule_type=DNA;topology=linear fffffb2f-878d-4b23-a82f-573fbc1fad5f feature region 1 1093 . . . ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f fffffb2f-878d-4b23-a82f-573fbc1fad5f feature CDS 1 87 . - 0 ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_number=1;name=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1 fffffb2f-878d-4b23-a82f-573fbc1fad5f feature exon 1 87 . - . gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_number=1 fffffb2f-878d-4b23-a82f-573fbc1fad5f feature CDS 212 1093 . - 0 ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_number=2;name=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2 fffffb2f-878d-4b23-a82f-573fbc1fad5f feature exon 212 1093 . - . gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_number=2
We can even do this as a one-liner, if we’re clever:
# Isn't bash grand?
$ bgzip $( grep fffffb2f-878d-4b23-a82f-573fbc1fad5f ca29.gff.loc.tsv | awk '{ print "-b", $2, "-s", $3 }' ) ca29.gff.gz
##sequence-region fffffb2f-878d-4b23-a82f-573fbc1fad5f 1 1093
fffffb2f-878d-4b23-a82f-573fbc1fad5f annotation remark 1 1093 . . . molecule_type=DNA;topology=linear
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature region 1 1093 . . . ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature CDS 1 87 . - 0 ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_number=1;name=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature exon 1 87 . - . gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_number=1
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature CDS 212 1093 . - 0 ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_number=2;name=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature exon 212 1093 . - . gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_number=2
So, knowing nothing except the sequence ID we want to find and the library it’s in, we can quickly and efficiently grab all the gene calls for that sequence ID as a GFF file.
But, there’s a catch – our location file is Big! It contains one line for every sequence in the GFF file:
$ ls -alh ca29.gff.loc.tsv -rw-rw-r-- 1 ubuntu ubuntu 7.8G Oct 4 22:02 ca29.gff.loc.tsv
So let’s dig into some options for dealing with this file.
Compressing the Index File
We now have an index file formatted as follows:
seq_id offset length 00000020-3ae5-4305-8de8-1bfceaa81f1c 16 1276 00000035-1d98-4718-960c-e1c0f9fda8c9 1292 937
Note this file is extremely dense – we have the sequence ID, offset, and length, and very little else. Consequently, the file doesn’t compress particularly well:
# Uncompressed: $ du -hs ca29.gff.loc.tsv 7.8G ca29.gff.loc.tsv # Compressed: $ du -hs ca29.gff.loc.tsv.gz 3.7G ca29.gff.loc.tsv.gz
There’s just not a lot of slack to take out here. We could do something similar to BGZip and use a binary format, as opposed to just a standard ASCII file, and we might get some savings there, but then we’d need to build and distribute a toolchain for dealing with our proprietary index format, so that’s off the table for now.
Note also that the uncompressed index is comparable to the compressed GFF itself – 8Gb vs 24Gb, or about a third the size. We’d prefer to keep the compressed version. However, we quickly hit a problem: finding individual sequences from the gzipped file is a whole lot slower than the raw text format:
# Uncompressed: $ time cat ca29.gff.loc.tsv | grep fffffb2f-878d-4b23-a82f-573fbc1fad5f fffffb2f-878d-4b23-a82f-573fbc1fad5f 251262762041 935 real 0m4.004s user 0m2.328s sys 0m4.776s # Compressed: $ time zcat ca29.gff.loc.tsv.gz | grep fffffb2f-878d-4b23-a82f-573fbc1fad5f fffffb2f-878d-4b23-a82f-573fbc1fad5f 251262762041 935 real 1m0.754s user 1m0.947s sys 0m5.050s
Yowza! A minute to look up a single sequence in our compressed location file? What’s going on?
Let’s pull out pv and see where our bottleneck is. pv is a lovely little tool that works roughly like cat, but will show how fast it’s passing data along. Linux pipes affect our speed, so we’ll throw cat in between pv and grep just to make things fair:
# Uncompressed: $ time pv ca29.gff.loc.tsv | cat - | grep fffffb2f-878d-4b23-a82f-573fbc1fad5f 7.73GiB 0:00:03 [1.95GiB/s] [==========================================>] 100% fffffb2f-878d-4b23-a82f-573fbc1fad5f 251262762041 935 real 0m3.745s user 0m2.486s sys 0m2.580s # Compressed $ time pv ca29.gff.loc.tsv.gz | zcat - | grep fffffb2f-878d-4b23-a82f-573fbc1fad5f 3.62GiB 0:01:01 [60.3MiB/s] [==========================================>] 100% fffffb2f-878d-4b23-a82f-573fbc1fad5f 251262762041 935 real 1m1.473s user 1m1.319s sys 0m5.835s
Decompressing our compressed archive slows us down from nearly 2GiB/s to a piddling 60MiB/s. Our gzipped location file is half the size, but takes 30x as long to retrieve sequences from. That’s not great. The difference reduces a bit if we decide to search for multiple sequence IDs:
# 10 sequence IDs to search for $ cat seq_ids.txt fef21bda-479c-445b-affd-d04761795174 fdd86a15-1848-4a2b-be05-f9d62a0ee538 fcdd2245-17b1-4a7a-9f67-1a34802aefd3 ffc97fde-cf41-4440-b50b-86abad46f415 ff9fa0a4-7042-45bc-bf70-c502a1e49522 feed9c3a-9831-4c42-81a9-d072f728c9f2 fce7cd9b-8cf9-41f7-a4a2-04a7814e41b6 fd7fae42-f746-4934-b71e-d738d8da31bc fd55e51f-6773-4eb4-871d-4174e99547a9 fd581f3d-6ce0-4895-8b38-083c091a4729 # Uncompressed: $ time pv ca29.gff.loc.tsv | cat - | grep --file seq_ids.txt fcdd2245-17b1-4a7a-9f67-1a34802aefd3 248192112609 943 fce7cd9b-8cf9-41f7-a4a2-04a7814e41b6 248232660866 597 fd55e51f-6773-4eb4-871d-4174e99547a9 248655071744 597 fd581f3d-6ce0-4895-8b38-083c091a4729 248663484363 1276 fd7fae42-f746-4934-b71e-d738d8da31bc 248815519201 2291 fdd86a15-1848-4a2b-be05-f9d62a0ee538 249153912987 937 feed9c3a-9831-4c42-81a9-d072f728c9f2 250214540699 597 fef21bda-479c-445b-affd-d04761795174 250231619577 599 ff9fa0a4-7042-45bc-bf70-c502a1e49522 250893965638 941 ffc97fde-cf41-4440-b50b-86abad46f415 251054447078 937 7.73GiB 0:00:10 [ 781MiB/s] [==========================================>] 100% real 0m10.132s user 0m8.352s sys 0m6.810s # Compressed: $ time pv ca29.gff.loc.tsv.gz | zcat - | grep --file seq_ids.txt > /dev/null 3.62GiB 0:01:02 [58.8MiB/s] [==========================================>] 100% real 1m2.997s user 1m8.761s sys 0m5.691s
Not a big difference for our compressed reads – that’s obviously bottlenecked on decompression – but our uncompressed reads drop to around 800MiB/s, and it takes about 10s to get our results. This is probably as fast as we’re going to get with the index file approach, so let’s keep this as a benchmark.
First, though, let’s try a very different approach.
Using Sqlite3
What if we went with a database instead? Sqlite works fine as a command line utility, and it should certainly be faster than grepping through a text file – and indeed it is:
# The SQLite DB:
$ sqlite3 ca29.gff.sqlite3
SQLite version 3.38.2 2022-03-26 13:51:10
Enter ".help" for usage hints.
sqlite> .schema
CREATE TABLE location (
contig_id TEXT PRIMARY KEY,
start INTEGER NOT NULL,
length INTEGER NOT NULL
);
sqlite> select * from location limit 5;
seq_id|offset|length
00000020-3ae5-4305-8de8-1bfceaa81f1c|16|1276
00000035-1d98-4718-960c-e1c0f9fda8c9|1292|937
00000057-0ac5-467f-a37f-3e64cf9928e8|2229|937
0000005c-da49-428c-95e5-342515731bc0|3166|937
sqlite>
# Single sequence retrieval is almost instantaneous
$ time sqlite3 -tabs ca29.gff.sqlite3 'select * from location where contig_id = "fef21bda-479c-445b-affd-d04761795174"'
fef21bda-479c-445b-affd-d04761795174 250231619577 599
real 0m0.002s
user 0m0.002s
sys 0m0.000s
# Multiple sequence retrieval takes about the same time
$ time sqlite3 -tabs ca29.gff.sqlite3 "select * from location where contig_id in ('fef21bda-479c-445b-affd-d04761795174','fdd86a15-1848-4a2b-be05-f9d62a0ee538','fcdd2245-17b1-4a7a-9f67-1a34802aefd3','ffc97fde-cf41-4440-b50b-86abad46f415','ff9fa0a4-7042-45bc-bf70-c502a1e49522','feed9c3a-9831-4c42-81a9-d072f728c9f2','fce7cd9b-8cf9-41f7-a4a2-04a7814e41b6','fd7fae42-f746-4934-b71e-d738d8da31bc','fd55e51f-6773-4eb4-871d-4174e99547a9','fd581f3d-6ce0-4895-8b38-083c091a4729')"
fcdd2245-17b1-4a7a-9f67-1a34802aefd3 248192112609 943
fce7cd9b-8cf9-41f7-a4a2-04a7814e41b6 248232660866 597
fd55e51f-6773-4eb4-871d-4174e99547a9 248655071744 597
fd581f3d-6ce0-4895-8b38-083c091a4729 248663484363 1276
fd7fae42-f746-4934-b71e-d738d8da31bc 248815519201 2291
fdd86a15-1848-4a2b-be05-f9d62a0ee538 249153912987 937
feed9c3a-9831-4c42-81a9-d072f728c9f2 250214540699 597
fef21bda-479c-445b-affd-d04761795174 250231619577 599
ff9fa0a4-7042-45bc-bf70-c502a1e49522 250893965638 941
ffc97fde-cf41-4440-b50b-86abad46f415 251054447078 937
real 0m0.002s
user 0m0.002s
sys 0m0.000s
Great! The SQLite3 database lookup is almost immediate! So what’s the problem? Well:
$ du -hs ca29.gff.sqlite3 16G ca29.gff.sqlite3
The problem is the database is on the order of 2/3rds the size of the compressed file it’s searching. This might be a fine trade-off in some cases – after all, it’s providing nearly instantaneous lookups, which might be worth the storage space:
$ time bgzip $( sqlite3 -tabs ca29.gff.sqlite3 'select * from location where contig_id = "fffffb2f-878d-4b23-a82f-573fbc1fad5f"' | awk '{ print "-b", $2, "-s", $3 }' ) ca29.gff.gz
##sequence-region fffffb2f-878d-4b23-a82f-573fbc1fad5f 1 1093
fffffb2f-878d-4b23-a82f-573fbc1fad5f annotation remark 1 1093 . . . molecule_type=DNA;topology=linear
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature region 1 1093 . . . ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature CDS 1 87 . - 0 ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_number=1;name=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature exon 1 87 . - . gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_1;gene_number=1
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature CDS 212 1093 . - 0 ID=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_number=2;name=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2
fffffb2f-878d-4b23-a82f-573fbc1fad5f feature exon 212 1093 . - . gene_id=fffffb2f-878d-4b23-a82f-573fbc1fad5f_2;gene_number=2
real 0m0.058s
user 0m0.043s
sys 0m0.017s
But that’s an awful big trade, especially since our previous best of just grep-ing the raw text file for a single result was around 2 seconds. We’re not using this for an inline service, so the difference between 2 seconds and “immediate” really isn’t much for us. Let’s see if we can get any better from our index file.
Back to the Filesystem
So to recap, we’ve got an 8Gb uncompressed TSV file we can compress down to 4Gb, but then we see search times on the order of a minute or so for a set of sequence IDs.
Once again, taking a look at our timings, it’s clear the decompression time is our bottleneck. Taking grep out of the picture entirely makes this really clear:
# No pipes at all, just dump the file as fast as we can $ time pv ca29.gff.loc.tsv > /dev/null 7.73GiB 0:00:01 [7.61GiB/s] [==========================================>] 100% real 0m1.017s user 0m0.016s sys 0m1.000s # Pipe to cat - notice we cut our speed roughly in half right # off the bat if we want to consume our data: $ time pv ca29.gff.loc.tsv | cat - > /dev/null 7.73GiB 0:00:01 [4.38GiB/s] [==========================================>] 100% real 0m1.765s user 0m0.030s sys 0m2.838s # Decompressing the gzipped file is what's killing us, though: $ time pv ca29.gff.loc.tsv.gz | zcat - > /dev/null 3.62GiB 0:00:59 [62.4MiB/s] [==========================================>] 100% real 0m59.334s user 0m57.868s sys 0m1.973s
Let’s add grep back in, just so we’ve got a picture of where we’re trying to get:
# Grep for our set of 10 sequence IDs $ time pv ca29.gff.loc.tsv | grep --file seq_ids.txt > /dev/null 7.73GiB 0:00:09 [ 832MiB/s] [==========================================>] 100% real 0m9.511s user 0m8.164s sys 0m2.644s
So grep itself is limited to about 800Mb/s for multiple matches.
A Quick Diversion into Alternate Greps
There are a couple alternatives to grep that promise faster searches – let’s take a quick look at ag and ripgrep to see what they offer:
# RipGrep: $ time pv ca29.gff.loc.tsv | rg --file ./seq_ids.txt > /dev/null 7.73GiB 0:00:18 [ 422MiB/s] [==========================================>] 100% real 0m18.735s user 0m15.821s sys 0m10.066s # AG, the "silver searcher" - AG does not support a patterns file: $ time pv ca29.gff.loc.tsv | ag -F 'fef21bda-479c-445b-affd-d04761795174' fef21bda-479c-445b-affd-d04761795174 250231619577 599 7.73GiB 0:00:20 [ 380MiB/s] [==========================================>] 100% real 0m20.781s user 0m15.630s sys 0m19.217s
Both tools are half the speed of grep, and ag doesn’t support supplying multiple queries in a file. Generally, both of these tools are built to rapidly search a large tree of files, whereas we’re asking them to search one big file – as always, us humble bioinformaticians are left out in the cold. We’re not going to get saved by a different search utility – it looks like we’re stuck with grep.
Back to BGZip
So our bottleneck is decompression – remember earlier when I mentioned BGZip’s two big tricks? The multi-threaded decompression didn’t matter for us much before, since we were seeking specifically what we needed, but now we’re searching the whole file – let’s see if BGZip’s got something up its sleeve for us. As a reminder, here’s standard zcat:
# Standard zcat: $ time pv ca29.gff.loc.gz | zcat - > /dev/null 3.61GiB 0:01:03 [58.2MiB/s] [==========================================>] 100% real 1m3.615s user 1m2.497s sys 0m1.620s
And here’s bgzip with no threading, after we recompress the locations file using bgzip:
$ time pv ca29.gff.loc.tsv.gz | bgzip -c -d - > /dev/null 3.62GiB 0:00:24 [ 153MiB/s] [==========================================>] 100% real 0m24.102s user 0m22.340s sys 0m1.819s
Already we’re three times as fast – because BGZip can trade between reading and decompressing and make some better assumptions about the source data, we’re not bottlenecked quite as badly on decompression. But let’s see what happens when we actually use the parallel decompression:
# Two threads $ time pv ca29.gff.loc.tsv.gz | bgzip --threads 2 -c -d > /dev/null 3.62GiB 0:00:11 [ 320MiB/s] [==========================================>] 100% real 0m11.573s user 0m23.952s sys 0m2.653s # Four threads: $ time pv ca29.gff.loc.tsv.gz | bgzip --threads 4 -c -d > /dev/null 3.62GiB 0:00:05 [ 640MiB/s] [==========================================>] 100% real 0m5.788s user 0m23.714s sys 0m2.449s # Eight threads: $ time pv ca29.gff.loc.tsv.gz | bgzip --threads 8 -c -d > /dev/null 3.62GiB 0:00:02 [1.23GiB/s] [==========================================>] 100% real 0m2.953s user 0m24.297s sys 0m2.621s # 16 threads: $ time pv ca29.gff.loc.tsv.gz | bgzip --threads 16 -c -d > /dev/null 3.62GiB 0:00:02 [1.49GiB/s] [==========================================>] 100% real 0m2.439s user 0m29.208s sys 0m2.514s
These are enormous speedups in decompression, especially compared to the naive zcat approach – BGZip is doing some really fantastic work with threading here. You can see an almost linear increase in read speed right up through 8 threads, although we start to drop off above that – coordination takes its toll, eventually.
Let’s see what this looks like combined with grep. We’ll use 4 threads, and we’ll take pv out of the picture:
$ time bgzip --threads 4 -c -d ca29.gff.loc.tsv.gz | grep --file seq_ids.txt > /dev/null real 0m10.091s user 0m33.097s sys 0m5.670s
10 seconds – that’s roughly what we were getting with the uncompressed search. BGZip’s not as fast as just reading an uncompressed file, but it’s fast enough to keep up with grep – decompression is no longer our bottleneck. BGZip is fast enough that we can keep our compressed index file with no speed sacrifice.
A Couple More Grep Tricks
Two more notes on grep to close out our testing.
First, we’ve been searching for multiple results in our index file, since that’s a relatively normal use case for us, but that actually imposes a pretty steep penalty on grep:
$ time pv ca29.gff.loc.tsv | grep "fef21bda-479c-445b-affd-d04761795174" > /dev/null 7.73GiB 0:00:03 [1.96GiB/s] [==========================================>] 100% real 0m3.944s user 0m2.644s sys 0m2.616s
Fortunately, the real world effects here are pretty limited – BGZip is still fast enough that, roughly, this doesn’t matter:
$ time bgzip --threads 8 -c -d ca29.gff.loc.tsv.gz | grep --file just_one_seq_id.txt > /dev/null real 0m4.138s user 0m29.531s sys 0m5.396s
One change we can make, though, that does make a difference is telling grep when to stop. We’re looking for exactly one line per sequence, so we shouldn’t ever get more lines than we have in our input. The -m flag tells grep how many results to expect:
# Grab a sequence ID towards the front of our file: $ head -n100 ca29.gff.loc.tsv | tail -n1 00000c54-06f6-438b-9474-c56d4ff130b2 196289 932 # Grep without a hit of our count: (base) ubuntu@npd-edanielson-work:/data/ca29$ time cat ca29.gff.loc.tsv | grep "00000c54-06f6-438b-9474-c56d4ff130b2" > /dev/null real 0m3.134s user 0m0.026s sys 0m4.138s # Tell grep how many results to expect: (base) ubuntu@npd-edanielson-work:/data/ca29$ time cat ca29.gff.loc.tsv | grep -m1 "00000c54-06f6-438b-9474-c56d4ff130b2" > /dev/null real 0m0.002s user 0m0.000s sys 0m0.002s
We’ve been playing a bit of a trick with our timings so far, because we’ve been using the most pathological case we can – 10 sequence IDs all towards the end of the file, which means grep needs to search the entire file to find our results. In more realistic use cases, we’d have a more normally distributed set of sequence IDs, and adding the -m flag to our searches will make a much bigger difference. The takeaway here is that the 10 seconds we saw above could be considered our worst-case scenario – at worst, it will take 10 seconds to retrieve a set of results from our file. We can live with that.
Conclusion
Final Layout
Our final layout for our library includes four files:
- ca29.gff.gz (24Gb) – The full BGZip-compressed GFF file
- ca29.gff.gz.gzi (59Mb) – BGZip’s index file
- ca29.gff.loc.gz (3.7Gb) – The compressed sequence location TSV file
- ca29.gff.loc.gz.gzi (2.0Mb) – BGZip’s index file
The total space for all four files is around 28Gb, or around 12% the size of our original uncompressed GFF library. Retrieving a single sequence can be performed in one line of bash with only the sequence ID and the library in a matter of seconds:
$ time bgzip $( bgzip -c -d --threads 4 ./ca29.gff.loc.tsv.gz | grep -m1 -F "78ef59ad-94b9-4944-9a3e-e1c171fc0df5" | awk '{ print "-b", $2, "-s", $3 }' ) ca29.gff.gz > 78ef59ad-94b9-4944-9a3e-e1c171fc0df5.gff
real 0m2.803s
user 0m12.415s
sys 0m2.399s
$ cat 78ef59ad-94b9-4944-9a3e-e1c171fc0df5.gff
##sequence-region 78ef59ad-94b9-4944-9a3e-e1c171fc0df5 1 1620
78ef59ad-94b9-4944-9a3e-e1c171fc0df5 annotation remark 1 1620 . . . molecule_type=DNA;topology=linear
78ef59ad-94b9-4944-9a3e-e1c171fc0df5 feature region 1 1620 . . . ID=78ef59ad-94b9-4944-9a3e-e1c171fc0df5
78ef59ad-94b9-4944-9a3e-e1c171fc0df5 feature CDS 2 424 . - 0 ID=78ef59ad-94b9-4944-9a3e-e1c171fc0df5_1;gene_id=78ef59ad-94b9-4944-9a3e-e1c171fc0df5_1;gene_number=1;name=78ef59ad-94b9-4944-9a3e-e1c171fc0df5_1
78ef59ad-94b9-4944-9a3e-e1c171fc0df5 feature exon 2 424 . - . gene_id=78ef59ad-94b9-4944-9a3e-e1c171fc0df5_1;gene_number=1
78ef59ad-94b9-4944-9a3e-e1c171fc0df5 feature CDS 435 1619 . + 0 ID=78ef59ad-94b9-4944-9a3e-e1c171fc0df5_2;gene_id=78ef59ad-94b9-4944-9a3e-e1c171fc0df5_2;gene_number=2;name=78ef59ad-94b9-4944-9a3e-e1c171fc0df5_2
78ef59ad-94b9-4944-9a3e-e1c171fc0df5 feature exon 435 1619 . + . gene_id=78ef59ad-94b9-4944-9a3e-e1c171fc0df5_2;gene_number=2
Closing Remarks
Our original goal was to facilitate easy, fast, and cheap extraction of gene calls for a specific sequence from a library-level dataset. BGZip, part of the Samtools/HTSLib package, provides blazing fast indexed access into arbitrarily large gzip-compatible compressed files, with a very small space penalty for its index file. In addition, BGZip’s threaded compression & decompression, as well as some other optimizations, allow us to retrieve specific sequence locations from a compressed text location file as quickly as if the file were not compressed. The combination of these two features give us the ability to retrieve information for one or more sequences from our dataset in a matter of seconds.
Taken together, this allows us to keep our gene calls available on disk in a common bioinformatics format and to work with those files using standard *nix tools. Within our pipelines, we can go from a list of sequence IDs to the gene calls for those individual sequences in seconds. The approach generalizes to other file types as well – any file format in which the data for an individual item is written contiguously in the file and in which a single item can be read individually – such as FASTAs, JSONL, and others – is amenable to this approach.
In our next post, we will explore even faster ways to retrieve data from very large files!
(Feature photo by Joshua Sortino on Unsplash)
