PDA

View Full Version : BWA sampe to bam?


dmacmillan
04-20-2012, 04:04 PM
So as everyone who uses bwa knows, the sampe function outputs a file in sam format. What I want to do is somehow convert that sam file to a bam file in some sort of pipe? It seems easy to implement, but I keep getting an error from samtools.

cat file.sam | samtools view -Sb

that does not work!

kopi-o
04-20-2012, 05:20 PM
Look at the samtools manualpage: http://samtools.sourceforge.net/samtools.shtml

You are looking for samtools view -bS or samtools view -bt

swbarnes2
04-21-2012, 10:23 PM
What you want is something like:

bwa sampe ref.fa r1.sai r2.sai r1.fq r2.fq | samtools view -bSho out.bam -;

xied75
04-22-2012, 11:43 AM
Dear all,

Just wondering, SAM is far bigger than BAM, and seems not much people will open the SAM and read it, if from BWA direct output BAM, it saves a lot effort and the disk I/O is faster due to smaller file size. Does this make sense or I forgot something?

Best,

dong

arvid
04-23-2012, 12:09 AM
Dear all,

Just wondering, SAM is far bigger than BAM, and seems not much people will open the SAM and read it, if from BWA direct output BAM, it saves a lot effort and the disk I/O is faster due to smaller file size. Does this make sense or I forgot something?

Best,

dong

Theoretically, when your server is more CPU-limited than I/O-limited and you only need to sequentially read the whole file, SAM will be faster than BAM (due to the compression overhead in BAM). I found that this is never the case for our applications and therefore pipe aligners directly into a samtools chain (with the -m option to samtools sort to fit most alignments in memory, thus avoiding temporary files to be written to disk), to directly get a sorted BAM on disk.

dmacmillan
04-26-2012, 10:01 AM
What you want is something like:

bwa sampe ref.fa r1.sai r2.sai r1.fq r2.fq | samtools view -bSho out.bam -;

I understand what you are doing here, but what is with the '-;' at the end (ignoring the single quotations)?

swbarnes2
04-26-2012, 10:52 AM
I understand what you are doing here, but what is with the '-;' at the end (ignoring the single quotations)?

the '-' means "the thing that's being piped". At least, that's how I understand it. That command works, I use it all the time just like I wrote it there, so would this:

bwa sampe ref.fa r1.sai r2.sai r1.fq r2.fq | samtools view -bSh - > out.bam;

sdriscoll
04-26-2012, 12:29 PM
I don't know if it's necessary from the BWA output or not but I like to use the -F option for output from bowtie to eliminate unaligned reads from making their way into the BAM file. Also the -h option isn't necessary in this example - the BAM header gets created appropriately..in fact I don't think samtools will allow you to create a BAM file from a SAM file without the SAM file already having the correct header information. I've only needed the -h option when I view BAM files. By default the header is left off when viewing a BAM file as SAM via Samtools.

So what I always use is this:

bwa sampe ref.fa r1.sai r2.sai r1.fq r2.fq | samtools view -bS -F 0x04 - > out.bam

sometimes followed by this:

samtools sort out.bam out-sorted

Bowtie doesn't properly sort its output and I don't remember if BWA does either. If you use the BAM file for any downstream analysis you usually need it to be sorted by chromosome and position.

dmacmillan
04-26-2012, 03:50 PM
Interesting tips, I will try both, thanks!

arvid
04-26-2012, 11:31 PM
To reduce the I/O load (and total CPU time as well) even further, this is my favourite:


bwa sampe ref.fa r1.sai r2.sai r1.fq r2.fq | samtools view -bSu -F 0x04 - | samtools sort -m 4294967296 - out.sorted
samtools index out.sorted.bam


Set -m as high as you can afford; in my hands samtools sort needs RAM up to 2x the value specified there in bytes (I set this to 16 GB when running on a server, which is enough for most BAMs to be sorted without writing temporary files to disk). -u removes the compression/decompression overhead in the pipe between view and sort.

swbarnes2
04-27-2012, 01:15 PM
piping into samtools sort works? I was afraid that that would get ugly.

How can I ask the server I'm on how much memory I can devote to sort?

nilshomer
04-27-2012, 07:48 PM
Use the "-m" option in samtools sort instead.