This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Faster extraction of reads mapped to a region

So I think I've just read every single post about this on the internet.

I simply want to extract reads mapped to a certain region. This is the pipeline I'm currently using:

samtools view sorted.bam "chr3:1000-5000" | cut -f1 > ids
LC_ALL=C grep -wFf ids < original.sam > subset.sam  # this is the bottleneck
echo -e "$(samtools view -H original.sam)\n$(cat subset.sam)" > subset.sam
samtools view -b subset.sam | samtools fastq -1 reads1.fq -2 reads2.fq -

This is not too bad for a single run, but I actually have to do this for 1000+ regions. Is there possibly a more efficient way to extract the reads, or is this as fast as it gets?

Edit: added updating-the-header step
Edit: Help

samtools alignment
• 2,825 views • link
written by novice novice

Could you clarify what the desired output is? Do you want fastq files, read IDs, a subset bam or something else?

• 0 views • link
written by WouterDeCoster WouterDeCoster

The fastq files (paired end).

• 0 views • link
written by novice novice

1 answer

I would do something like this

cat <(samtools view -H merged.mdup.bam) <(samtools view merged.mdup.bam "chr3:20000-200000") | samtools view -bS - > subset.bam

You can loop the above command on every region giving region as name to the output file instead of subset.bam.

BTW, are you trying to generate fastq from your subset.bam ? Because I did not find fastq argument in my version (samtools-1.2), instead I found bam2fq.

• 0 views • link
written by venu venu

samtools fastq is in samtools v.1.3.1

• 0 views • link
written by GenoMax GenoMax

fastq is the new bam2fq.

The problem with this solution is that, if I pipe it into fastq, I get discordant fastq files. I need to make sure to extract the pairs of reads as well.

• 0 views • link
written by novice novice

Log in to answer this question.

More posts like this