Faster tetranucleotide (k-mer) frequencies!
Tags: perl, bioinformatics, profiling, softwaredevelopment
Published: 2024-03-15 08:31:22

I saw [Jennifer's post](https://dev.to/jmeneghin/calculating-tetranucleotide-k-mer-frequencies-969) about re-writing her [perl scripts](https://github.com/jmeneghin/perl-for-reysenbach-lab/blob/master/get_kmer_frequencies.pl) in python and how she saw a 2.5 times improvement.

How could this be?  My favourite language can't be that slow.
It must be programmer error.

I have an interest in Perl and Science, so time to roll up sleeves and learn me some profiling/benchmarking. What follows is my internal monologue and the notes I scribbled down during the learning process. For those that want to follow along, I've [created a small repo](https://github.com/duffee/faster-perl-for-reysenbach/tree/main/get_kmer_frequencies).

## Getting started

Voltaire said that Hell is other people's code. My first step was to [re-write it](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies1.pl) into [Modern Perl](http://modernperlbooks.com/) and in the process, understand what each line does. When it's written idiomatically, it's easier to refactor and I should be able to make some minor performance improvements along the way.

Assume that the original script has been tested enough.  For me to be correct, I've got to produce the exact same output.  I got close, except for the header line.
[line 99](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies_original.pl#L99) `print OUT "\t$prefix_$j";` becomes
[line 89](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies1.pl#L89) `print $out_fh "\t$j";`  Yes, that's a bug because `$prefix_` doesn't exist.

Search "benchmarking tools for linux" and decide that [hyperfine](https://github.com/sharkdp/hyperfine) is good for what I'm doing. Run Jennifer's new python script against my refactored perl and find that the python is 1.26 times faster for k=3 and 1.47 times faster for k=4. For the Covid-19 sequence, these are both on the order of hundreds of milliseconds.
```sh
hyperfine --warmup 3 'perl/get_kmer_frequencies.pl Covid-19_seq.fasta 3 boyd1' 'python/get_kmer_frequencies.py -i Covid-19_seq.fasta -k 3 -p boyd2'
```
Ok, not bad. Better than 2.9 times faster, but that's probably down to the way that **hyperfine** warms the cache and separates out User time from System time.

Oh, I should just check how much I improved when I refactored.  Run it against Jennifer's original perl script and ... hers was 1.1 times faster. Well, that was a bit embarrassing.

_ahem_ I was ... aiming at improving readability, .. maintainability, y'know best practice and all that. That's my story and I'm sticking to it. ;)

## For sanity's sake

Check that the output of the new file is the same as the original, otherwise you've messed up the refactoring. I started using this test script with [prove](https://metacpan.org/dist/Test-Harness/view/bin/prove) to make it quick and easy.
Saved as **i.t**, I run it with `prove i.t` for the lols.
It gets noisy when there's a problem, so I go back to running it by hand.

```perl
use Test2::V0;

my $standard = 'get_kmer_frequencies.pl';
my @files = sort { -M $a <=> -M $b } glob 'get_kmer*';
ok my $latest = shift @files;
isnt $latest, $standard, 'Files to compare';

my @args = qw'Covid-19.fasta 3';
ok system('perl', $standard, @args, 'A') == 0, 'Make A_kmers.txt';
ok system('perl', $latest, @args, 'B') == 0, 'Make B_kmers.txt';

is `diff A_kmers.txt B_kmers.txt`, q{}, 'No differences in output';

done_testing();
```
Clever people do this from the start.
I did this after a bug I introduced messed up the output and I hadn't immediately noticed. What it was is that I changed the key separator to a character that was found in some of those keys and it then split those keys. Oops.

## NYTProf time

When you get serious about optimizing programs, trying to enhance performance, you reach for [profiling tools](https://en.wikipedia.org/wiki/Profiling_(computer_programming)) that can analyze your code's memory or time complexity. In Perl, [Devel::NYTProf](https://metacpan.org/pod/Devel::NYTProf) comes highly recommended. I use it to collect data on the number of times each statement is called and how long it spends executing it. That way I can work out where to invest the effort making the script faster, what gives the most bang for the buck.

Grab the profiler and run
```
perl -d:NYTProf get_kmer_frequencies.pl Covid-19_seq.fasta 3 boyd1
```
and open up the **nytprof/index.html** using `nytprofhtml --open` to see
```
Calls   P   F   Exclusive Time  Inclusive Time  Subroutine
9653    1   1   31.6ms  31.6ms  main::rc_seq
25      2   1   28.1ms  59.7ms  main::process_it
82498   7   1   9.46ms  9.46ms  main::CORE:print (opcode)
3175    4   1   7.89ms  7.89ms  main::CORE:sort (opcode)
```

### Sorting out the sort

Obviously, the `rc_seq` is the big sub that needs attention, but what about that `sort`? Quickly looking at the sort on [Line 78](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies1.pl#L78) `for my $i (keys %knucs)` I see that there's no reason to sort those keys. Saved one sort and the script runs about the same. There's another sort [inside a loop](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies1.pl#L96) which can be extracted out of the loop. [Extracting that](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies2.pl#L84) made it run 1.15 times faster!

Changing the header line ([line 87](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies1.pl#L87)) from a for loop to a [join over a list](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies3.pl#L87) is 1 or 2 percent faster.

### How do I print thee? Let me count the ways.

Messing about with printing in the inner loop didn't gain much, but changing the key separator from a tab `"\t"` (interpolated string) to an underscore `'_'` (a string literal) made a 10% improvement. (it also introduced the bug noted above because the keys used the underscore. changed it to a colon - bug gone)

`say` is marginally slower than `print` so use [print inside the loop](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies5.pl#L93) that gets called a lot to save maybe 10% on that call. From 32ms to 28ms is a small, but nice gain for a one line change.

### rc_seq - transforming the sequence

The `rc_seq` sub is an [if-elsif block](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies4.pl#L151) that splits a string into individual characters, translates **ACGT** into their complement (**TGCA**), reverses the array and joins it back into a string.

Being Perl, we can [manipulate and reverse the string](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies5.pl#L144) in-place. The change makes it shorter and more obvious (sometimes it runs faster). Actually, I ran this through the profiler and the sub now runs 5 times faster.

### process_it - collecting the frequencies

This sub does the work of splitting the sequence into **kmers** and counting them. The longest time spent here is [incrementing](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies5.pl#L136) the `%knucs` hash.

The second longest time is spent turning the sequence into an array of letters to create all the kmer substrings. Splitting isn't bad, but [joining sets of letters](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies5.pl#L129) together is. Use the [string function,](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies6.pl#L130) `substr` instead and speed that line up by 5 times.

Now marginally faster the python script in speed. Over 20% faster for k=3, and 5% (+/- 5%) faster for k=4. That's a decent improvement.

### Like the end of a great song ... a Key change!

[line 91](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies6.pl#L91)
```
my @items = map { my $key = join ':', $_, $i; $knucs{$key} // 0 } @record_keys;
```
The script spends most of its time (55ms!!!), longer than _anything else_, on this line.

Assume that the problem isn't the `map` but constructing the key for the lookup. Change the key to a [2 dimensional lookup](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies7.pl#L92) and see if that improves things.

**WHEN** you finally get it right (and remember the correct order that you construct the keys in), [line 92](https://github.com/duffee/faster-perl-for-reysenbach/blob/main/get_kmer_frequencies/get_kmer_frequencies7.pl#L92) is now 2.5 times faster than before and the perl script is now 40% faster than the python script.

Keys are constructed/used on lines 79, 91, 135

## STOP!!!!

Know when to stop.

There are no more obvious or easy gains here. Any more work is likely to yield small returns. Go outside, have a life or at the least consult [the relevant chart](https://xkcd.com/1205/).

Well, after thinking a while, maybe constructing the output could be improved, but I'm moving on. I've exceeded my goal of making the perl script as fast as the python script and learned more about refactoring and profiling. A bit like how audiophiles use your music to listen to their equipment, I've used Jennifer's science to better understand my Perl and had fun doing it.

There's a niggling thought at the back of my mind, now that I feel I better understand the purpose of the script, whether [BioPerl](https://bioperl.org) can do this even faster. I will leave that for another day. Oh, look _glycine_ has already done [most of the hard work](https://blogs.perl.org/users/glycine/2024/03/reading-sequences-from-fasta-foramt-alignment-by-bioperl.html) for me. Many thanks!

## Lessons learned

In summary, these are reflections on the changes that I made in chronological order. This may be someone's first time considering performance, so I include basic rules of thumb I used along with the things I did not know before.

* Modern Perl style adds a small amount of overhead, but the sanity it brings is a price worth paying.
* Streamline a method of checking the output hasn't changed
* Don't `sort` when order is not important
* Calculate constant values _outside_ of loops
* Use built-in list functions over loops (`join` instead of `for`)
* Interpolated strings are slower than string literals (prefer single quotes over double quotes)
* `say` is slightly slower than `print`.  Avoid it in heavy loops.
* `substr` is faster than `split`ing and `join`ing
* Creating a single hash key is slower than using a 2 level hash

The Human Genome is _way_ too large. Grab the protein sequence for Caenorhabditis elegans. It takes about 5 minutes to run.

Run your frequent tests with the Covid-19 sequence. Repeated runs with anything larger take too long for rapid turnaround.

**WARNING:** hyperfine will run each program 20 to 40 times to get decent statistics. You won't want to wait around for a file that takes 5 minutes to process a single run.

I'll leave you with a couple of related references for further reading, _chrisarg_'s work on [parsing FastQ files fast](https://blogs.perl.org/users/chrisarg/2023/09/of-go-c-perl-and-fastq-file-conversion-vol-i-intro.html)
and a marketsplash tutorial on [Perl code profiling tools](https://marketsplash.com/tutorials/perl/perl-code-profiling-tools/)

## In Conclusion

My corollary to Cunningham's Law:
Don't ask people how to make your code run faster;
[Tell them their language is slow](https://xkcd.com/386/)

It's taken a lot longer to reply Jennifer's post than I'd anticipated, but right now I have the warm glow that comes from being able to say, (until someone iterates on the above corollary)

... **Python is SLOW!**