I have a long list of read IDs of interest to me in a file called read_names.txt. it is simply in the format:

m54197_180831_211346/4981510/ccs
m54197_180831_211346/6226723/ccs
m54197_180831_211346/6619607/ccs
...

etc where these are the actual read names from a fastq file. I am then trying to find those in a bam file for which the fastq has been mapped, and I can accomplish this like:

for ID in `cat read_names.txt`
do
samtools view inbam.bam | grep $ID >> read_locs.sam
done

However, this method is obscenely slow because it is rerunning samtools view for every ID iteration (several hours now for 600 read IDs), and I was hoping to do this for several read_names.txt files.

I tried sort of flipping the script a bit and running samtools view first but it only returned the first read ID present in the file and stopped:

samtools view inbam.bam | for ID in `cat read_names.txt`; do grep $ID >> read_locs.txt; done

Am I missing something or is this the only way to accomplish this task?

  • 1
    Just as a general rule never do for var in $(cat file). Also known as bash pitfall number 1. – terdon Dec 6 at 21:45
  • thanks terdon, I will keep this in mind and find the better ways to handle this in the future. I use this method quite often so that is great for me to know! – d_kennetz Dec 6 at 21:51
  • Yeah, most of us use it very often until it bites us! :) – terdon Dec 6 at 22:03
up vote 5 down vote accepted

It is still slow but grep has a -f option to take in a file

samtools view inbam.bam | grep -f read_names.txt > read_locs.txt

  • This took 30 seconds as opposed to several hours, so I'd say the overhead was definitely in samtools view. I feel silly because of how simple this was. Thank you! – d_kennetz Dec 6 at 20:18
  • just as a warning this can take much longer if the read name file is big. – Bioathlete Dec 6 at 20:24
  • I appreciate the heads up! Comparatively, I do not think it will ever take longer than the method I was using so it is still a great answer. (This is serving as a caveat to your warning). – d_kennetz Dec 6 at 20:26

The reason your approach is slow is because you need to process the entire bam file from start to finish for every ID you are searching for. Here are some other solutions and a time comparison:

$ time bash -c 'for id in `cat names`; do 
    samtools view foo.bam | grep  \"$id\" >> out; 
done'

real    0m14.310s
user    0m14.170s
sys     0m2.540s

A very simple improvement, stop reading the file when you find the first match, will already speed things up a little:

$ time bash -c 'for id in `cat names`; do 
    samtools view foo.bam | grep -m1 \"$id\" >> out; 
done'

real    0m14.091s
user    0m13.927s
sys     0m2.506s

So will a safer shell loop (you should never do for var in `command`; do) although, in any case, using shell loops for text parsing is inherently a bad idea):

$ time bash -c "while read -r id; do 
    samtools view foo.bam | grep -m1  '$id'; 
done < names > out"

real    0m0.302s
user    0m0.355s
sys     0m0.126s

I don't quite understand why that one is so much faster, but presumably the shel is doing some sort of optimization.

Now, Bioathlete's answer:

$ time bash -c "samtools view foo.bam | grep -f names > out"

real    0m0.176s
user    0m0.210s
sys     0m0.015s

Or, on a larger (336M) file so the times are worth measuring:

$ time bash -c "samtools view bar.bam  | grep -F names > out"

real    0m11.137s
user    0m10.653s
sys     0m1.972s

And we can probably improve on that slightly if we add -w so it will only match entire words (which you want anyway since one read's name might be a substring of another's) and -F so that grep doesn't try to interpret the name as regular expressions:

$ time bash -c "samtools view bar.bam  | grep -m 1 -wFf names > out"

real    0m4.972s
user    0m6.527s
sys     0m0.595s
  • oh good call on the -w forgot about that – Bioathlete Dec 6 at 22:43
  • @Bioathlete it's the -m that really speeds it up though. But yeah, the -w is needed for safety. – terdon Dec 6 at 22:49

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.