Question: Extract Alignment From Very Large Bam File
5
gravatar for Plantae
8.3 years ago by
Plantae380
Plantae380 wrote:

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 • 24k views
ADD COMMENTlink modified 6.0 years ago by dariober11k • written 8.3 years ago by Plantae380

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

ADD REPLYlink written 8.3 years ago by Arun2.3k
6
gravatar for Bpow
8.3 years ago by
Bpow240
United States
Bpow240 wrote:

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 written 8.3 years ago by Bpow240

+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?

ADD REPLYlink written 8.3 years ago by Arun2.3k

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

ADD REPLYlink modified 8.3 years ago • written 8.3 years ago by Plantae380

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

ADD REPLYlink written 8.3 years ago by Arun2.3k

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

ADD REPLYlink written 8.3 years ago by David Langenberger9.5k
5
gravatar for David Langenberger
8.3 years ago by
Deutschland
David Langenberger9.5k wrote:

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 written 8.3 years ago by David Langenberger9.5k

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 written 3.2 years ago by Chadi Saad80
4
gravatar for dariober
6.0 years ago by
dariober11k
WCIP | Glasgow | UK
dariober11k wrote:

See also this A: Create smaller bam file based on the list of sequences 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 A: Create smaller bam file based on the list of sequences):

cat references.txt | xargs samtools view -b input.bam > output.bam
ADD COMMENTlink written 6.0 years ago by dariober11k
3
gravatar for Hypotheses
6.0 years ago by
Hypotheses90
Bangkok, Thailand
Hypotheses90 wrote:

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 written 6.0 years ago by Hypotheses90
0
gravatar for Raony Guimarães
8.3 years ago by
Dublin / Ireland
Raony Guimarães1.1k wrote:

use Bedtools, it accepts BAM files!

bamToBed -i reads.bam > reads.bed

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

ADD COMMENTlink written 8.3 years ago by Raony Guimarães1.1k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1375 users visited in the last hour