Yesterday I was in the need of a tool for creating a REF_CACHE for some proprietary software package. It was documented, but links were dead, so I ended up reverse engineering the REF_CACHE script. A few hours later I found out that there is already a perl script in this repository that does the same thing because the REF_CACHE is a samtools thing.

However, this perl script is quite slow, so I decided to contribute back my script so others can benefit from it.

Benchmarks were done on a tmpfs in-memory filesystem for writing. The reference file was pulled through cat a few times so it was also in the memory cache.

For reference,

# Calculating the md5sum
$ /usr/bin/time md5sum ~/test/ucsc.hg19_nohap.fasta > /dev/null
4.09user 0.21system 0:04.30elapsed 99%CPU (0avgtext+0avgdata 2052maxresident)k
0inputs+0outputs (0major+108minor)pagefaults 0swaps

# Copying the whole file to the in-memory filesystem
$ /usr/bin/time cp ~/test/ucsc.hg19_nohap.fasta tmp/
0.01user 1.40system 0:01.41elapsed 100%CPU (0avgtext+0avgdata 2356maxresident)k
0inputs+0outputs (0major+163minor)pagefaults 0swaps

So in this setup colculating the md5sum and copying the file should take roughly 5.7 seconds. When the cache is already populated, no copying is required and the baseline is 4.3 seconds.

The perl script

# Running with no populated cache
$ /usr/bin/time ./seq_cache_populate.pl -root tmp/plcache ~/test/ucsc.hg19_nohap.fasta
36.43user 1.75system 0:38.19elapsed 99%CPU (0avgtext+0avgdata 272984maxresident)k
0inputs+0outputs (0major+1684minor)pagefaults 0swaps

# Running with already present cache
$ /usr/bin/time ./seq_cache_populate.pl -root tmp/plcache ~/test/ucsc.hg19_nohap.fasta
35.86user 0.52system 0:36.40elapsed 99%CPU (0avgtext+0avgdata 273024maxresident)k
0inputs+0outputs (0major+1686minor)pagefaults 0swaps

The perl script takes roughly 1.8 seconds for the copying, which is fair. The rest of the script takes quite long compared to a raw md5sum check.

The python script

# Running with no populated cache
$ /usr/bin/time ./seq_cache_populate.py -root tmp/pycache ~/test/ucsc.hg19_nohap.fasta
5.63user 1.53system 0:07.17elapsed 99%CPU (0avgtext+0avgdata 262480maxresident)k
0inputs+8outputs (0major+7520minor)pagefaults 0swaps

$ /usr/bin/time ./seq_cache_populate.py -root tmp/pycache ~/test/ucsc.hg19_nohap.fasta
# Running with already present cache
5.67user 0.27system 0:05.94elapsed 99%CPU (0avgtext+0avgdata 262556maxresident)k
0inputs+8outputs (0major+7518minor)pagefaults 0swaps

Some very slight overhead (1.7 seconds) compared to the raw operations. This is due to the newline removal and upper casing, which require some cost.

The python script is roughly 5 times faster for a new cache and roughly 6 times for an already populated cache. It uses the same CLI and provides the same functionality including the in-memory spooling.

I wouldn't have targeted seq_cache_populate.pl for improvement had I known of it before, but since I did most of the work already I might as well contribute it back and give other users an easier time generating their ref caches.