2
$\begingroup$

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?

$\endgroup$
  • 1
    $\begingroup$ Just as a general rule never do for var in $(cat file). Also known as bash pitfall number 1. $\endgroup$ – terdon Dec 6 '18 at 21:45
  • $\begingroup$ 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! $\endgroup$ – d_kennetz Dec 6 '18 at 21:51
  • $\begingroup$ Yeah, most of us use it very often until it bites us! :) $\endgroup$ – terdon Dec 6 '18 at 22:03
8
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – d_kennetz Dec 6 '18 at 20:18
  • $\begingroup$ just as a warning this can take much longer if the read name file is big. $\endgroup$ – Bioathlete Dec 6 '18 at 20:24
  • $\begingroup$ 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). $\endgroup$ – d_kennetz Dec 6 '18 at 20:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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