Btukfyl

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?

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

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