Biostar Beta. Not for public use.
Extract Alignment From Very Large Bam File
5
Entering edit mode
6.5 years ago
Plantae • 380

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?

bam • 21k views
ADD COMMENTlink
0
Entering edit mode

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

ADD REPLYlink
6
Entering edit mode
6.5 years ago
Bpow • 200
United States

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.

ADD COMMENTlink
0
Entering edit mode
ADD REPLYlink
0
Entering edit mode

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

ADD REPLYlink
0
Entering edit mode

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

ADD REPLYlink
0
Entering edit mode

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

ADD REPLYlink
5
Entering edit mode
16 months ago
Deutschland

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
ADD COMMENTlink
0
Entering edit mode

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

ADD REPLYlink
3
Entering edit mode
17 months ago
WCIP | Glasgow | UK

See also this question & answer where @matted and myself compared samtools view -L withsamtools 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
ADD COMMENTlink
1
Entering edit mode
21 months ago
Hypotheses • 70
Bangkok, Thailand

my two-cents

If this doesn't work samtools view -h reads.bam 1:1042000-1042010

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 
ADD COMMENTlink
0
Entering edit mode
19 months ago
Dublin / Ireland

use Bedtools, it accepts BAM files!

bamToBed -i reads.bam > reads.bed

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

ADD COMMENTlink

Login before adding your answer.

Similar Posts
Loading Similar Posts
Powered by the version 2.3.1