Extract Alignment From Very Large Bam File

Question

[Plantae]

I have bam files ~50Gb, which is read alignments onto a reference genome. now I want to extract some regions of alignment, but I found methods like samtools view -L region.bed all.bam

is too slow, usually taken more than 2 hours (even if very small region are submit)

Are there any other methods that can do this job faster?

[Arun]

Can you paste the header of your bam file (with samtools view -H file.bam)

Answer 1

[Bpow]

Extracting small regions from a bam file should be sub-seconds fast if the files are properly indexed. You don't say what version of samtools you are using (and I don't know if this is still an issue), but I recall that at least some point there was an issue where specifying the region on the command line

samtools view -h reads.bam 1:1042000-1042010

would use the index while specifying the -L option and providing a bed file of intervals would result in the entire file being read. You may try specifying an interval at the command line to see if that helps.

[Arun]

+1. Exactly what I thought as well. That's why I asked for the header to see if its coordinate sorted. Does samtools run on bam files that are not coordinate sorted? If so, probably the index is not of help?

[Plantae]

I indexed the bam file, then use your method. It is really very fast, much faster than bedtool

[Arun]

Plantae, in that case, do you mind accepting the answer? :)

[David Langenberger]

For only one position, it's the fastest way to do!

Answer 2

[David Langenberger]

BEDtools is a good tip, but I would recommend intersectbed here:

intersectBed -abam reads.bam -b regions.bed > reads.regions.bam

or, if you want output in BED format

intersectBed -abam reads.bam -b regions -bed > reads.regions.bed

[Chadi Saad]

intersectBed -abam reads.bam -b regions.bed > reads.regions.bam will report original reads right ? How can i clip this reads to get just the part that overlap my intervals in the bed file ? thanks

Answer 3

[dariober]

See also this question & answer where @matted and myself compared samtools view -L with samtools view <region>. The latter is much faster. If you have several regions in file references.txt use (credit @matted):

cat references.txt | xargs samtools view -b input.bam > output.bam

Answer 4

[Hypotheses]

My two cents:

samtools view -h reads.bam 1:1042000-1042010

If that doesn't work. try this:

samtools view -h reads.bam chr1:10420000-10421000

To output to another bam file format

samtools view -b reads.bam chr1:10420000-10421000 > subset.bam

Answer 5

[Raony Guimarães]

Use Bedtools, it accepts BAM files!

bamToBed -i reads.bam > reads.bed

http://code.google.com/p/bedtools/wiki/Usage

Answer 6

[Marco]

QuickBAM might be able to help, I think:

• https://academic.oup.com/bioinformatics/article/39/8/btad463/7232227
• https://gitlab.com/yiq/quickbam