How to extract reads from bam file?

Question

[silas008]

Hi, guys!

I am trying to extract some reads that aligned in specific regions of the mitochondrial genome. I want to know what is the length of these reads, so I need to extract them before to measure, to avoid the measurement of other reads the aligned in other regions.

I did it by different ways, the last one was:

samtools view -bh original_sorted.bam chrM:1-55 > output.bam

(many times for different regions)

Then:

samtools merge merged.bam *.bam
samtools sort merged.bam > sorted.bam
samtools index sorted.bam

But I can't open the sorted.bam file on genome browser and when I try to transform it into a fastq I see this message:

[M::bam2fq_mainloop] processed 0 reads

The original_sorted.bam is ok. I can visualize the data and there are a lot of reads in the regions that I want to do the measurement.

Thanks, guys!!!

[Titus]

Hi, did you check if there is reads on your output.bam ?

[silas008]

Now i am using -L option and a bed file and the output.bam is ok. There are a lot of reads.

Thanks

[cpad0112]

samtools merge merged.bam *.bam

*.bam may pickup merged.bam as well. Try outputting some where outside the current directory @ silas008

[silas008]

I really didn't know why this was not working. I will try your suggestion.

Thanks

Answer 1

[finswimmer]

Hello silas008,

some questions and hints about this:

• What version of samtools are you using?
• If you want to extract multiple regions, you can define those regions in a bed file and use the -L option for samtools view. So there is no need to merge many file afterwards.
• How should your desired output look like? Depending on this there is a good chance that you don't have to make in intermediate bam file.
• Please use the formatting bar (especially the code option) to present your post better. I've done it for you this time.
  

Thank you!

fin swimmer

[silas008]

First of all, thank you for helpping.

I am sorry for the formatting, I am not so familiar with Biostars. I didn't know there is a code formatting bar. I will use it in next posts.

I am using Samtools version 1.4 and the option -L worked very well for me. I didn't know what I was doing wrong before but with -L option the output.bam is perfect.

Thank you very much

[finswimmer]

Hello again,

there's no need to say sorry. There are a lot of people who aren't aware of the formatting bar. That's why I'm pointing to it.

If you can do, I would recommend update samtools to the current version (v1.9).

fin swimmer

[silas008]

I will do that.

Thanks again