Further speedup in iVar Trim primer finding

[niemasd]

Is your feature request related to a problem? Please describe. I noticed that iVar Trim was recently optimized quite a bit (awesome stuff!), and I noticed that (assuming I'm understanding the release descriptions properly) one of the optimizations is that a faster sorting implementation was incorporated to sort primers, and if I recall correctly, the primers need to be sorted to enable binary search to find the end of a primer that spans a certain position (for primer trimming). However, you can actually avoid binary search entirely, which could dramatically speed things up

Describe the solution you'd like Because viral genomes are so short, at the beginning of iVar Trim's execution, you can precompute an array where index i contains the end position of the primer covering position i of the reference genome. Then, while trimming a given read that mapped to position i, you can perform a constant-time lookup (rather than binary search) to find the primer end position. A general summary of the preprocessing algorithm is as follows:

primer_end = array of length genome size
primers = sorted list of (start, end) tuples
primers_ind = 0
curr_primers = empty queue
for i from 0 to genome size:
    while curr_primers.front()[1] < i:
        curr_primers.dequeue()
    while primers[primers_ind][0] >= i:
        curr_primers.enqueue(primers[primers_ind])
        primers_ind = primers_ind + 1
    primer_end[i] = curr_primers.back()[1]

This slightly overly complicated approach allows overlapping primers, but if you can assume primers do not overlap:

primer_end = array of length genome size
primers = sorted list of (start, end) tuples
for start, end in primers:
    for i from start to end:
        primer_end[i] = end

Then, in your code, if a read maps to position i, you primer trim until the base that mapped to position primer_end[i]

Describe alternatives you've considered N/A

Additional context @kga1978 told me to make a GitHub Issue; happy to chat with you folks about this approach if you're interested 😄

[cmaceves]

hi @niemasd! thank you for the suggestion, it seems pretty reasonable for future development.

[gkarthik]

Hey @niemasd, this will certainly reduce the lookup time, thanks! A few additional things Chrissy and I discussed yesterday,

• Using a hashmap instead of an array since the array will be pretty sparse.
• We'll need to preprocess this for forward and reverse reads (if paired).
• We'll need to take the offset into account while creating the hashmap.

[niemasd]

Agreed re: 2nd and 3rd bullets, but re: first bullet, I think the memory overhead a hash map has with respect to load factor (for constant-time average-case time complexity) + buckets will result in similar memory as using an array, and you might not get the benefits of memory locality. But regardless, I think hash map or just a basic preallocated array will have similar runtime + memory (both should be extremely fast and low memory), so whichever is easier for you to implement is probably the better choice (I think the difference will be negligible compared to the rest of the runtime/memory of the tool)

But yeah, glad to hear you'll look into it! Feel free to ping me if you ever want to chat about it