Piping bowtie2 output directly into BAM
1
6
Entering edit mode
5.1 years ago
hersh ▴ 60

I am using bowtie2, and due to the size of the SAM output, I would like to entirely avoid the SAM file, and pipe the results into samtools to convert the output to a BAM file. I am getting a little bit hung up on how to accomplish this.

I have tried a few different options:

The approach here: http://slhogle.github.io/2014/bowtie-and-samtools/

and here: http://www.metagenomics.wiki/tools/samtools/breadth-of-coverage

are both giving me an empty BAM file.

I have also tried:

bowtie2 -x [index] -U [reads] -S ${TMPDIR}/temp.sam | samtools view -bS > out.bam

To see if I would get anything different using the temp directory to store the SAM temporarily.

Edit in response to a comment to add the commands I have tried:

bowtie2 -x [index] -U [reads] -S - | samtools view -bS > out.bam

And

bowtie2 -x [index] -U [reads] -S - | samtools view -bS -  > out.bam

I'm sure this has been done before, is there an obvious way to accomplish this task?

Thanks.

alignment • 22k views
ADD COMMENT
0
Entering edit mode
ADD REPLY
0
Entering edit mode

Why don't you write all the command lines you tried. It seems to me that the one you did write won't work, because it's not going to send output to samtools, it's writing it to temp.sam. Your first link is to bowtie, which might not work with bowtie2.

ADD REPLY
13
Entering edit mode
5.1 years ago
cmdcolin ★ 3.8k

You shouldn't need to use the -S command for the output, just omit this arguement

bowtie2 -x genomeindex -U out.fq|samtools view -bS - > out.bam

Versions

% bowtie2 --version
/usr/bin/bowtie2-align-s version 2.3.4.2
64-bit
% samtools --version
samtools 1.7
Using htslib 1.7-2
Copyright (C) 2018 Genome Research Ltd.
ADD COMMENT
0
Entering edit mode

Thanks! I tried this string and my BAM was correct. Was the issue just that I was using the -S flag?

ADD REPLY
0
Entering edit mode

Hey!

I did that and worked great. This is the code I used:

bowtie2 -p 10 -q --no-unal -k 20 -x bowtieIndex -1 out_1.fq  -2 out_2.fq | samtools view -@ 10 -Sb -o bowtie2.bam 1>1.txt 2>2.txt &

The problem is that I was running it on the background and now I only have the BAM file but it did not save the alignment stats in either the stdout or stderr... Does anyone know how to get the alignment stats (the log output from bowtie2) from a BAM file?

Thanks a lot!

ADD REPLY
1
Entering edit mode

Did you look in the two files above to see what was captured? You can get alignment stats using samtools idxstats. You will need to sort and index the BAM file.

ADD REPLY
0
Entering edit mode

The .txt files were both empty... I think because of the pipe, I was not getting the output from the bowtie2 but from the samtools (which has no stdout or stderr).

I tried also with samtools flagstats and bamtools stats but since for the bowtie2 run I dropped the unaligned reads with the option --no-unal I don't get the stats for the alignment of all reads, only for the reads that are contained in the bam file...

I will try with samtools idxstats and see what I can get

Thanks a lot!!

ADD REPLY
0
Entering edit mode

Since you used --no-unal idxstats output is only going to show you the number of reads that are aligned. You could simply subtract the total aligned from total reads to get the number that were excluded.

ADD REPLY
0
Entering edit mode

I was just doing that and I got extremely confused since the number of aligned reads is 2,340M and the number of total reads I had as input is 596M... I think is because it mapped each read to up to 20 possible sites with the option -k 20 I will investigate a way to get info of one mapping site per read to assess the number of mapped reads. Any ideas no how to do that are welcome!

Thanks a lot!

ADD REPLY
1
Entering edit mode

Correct. You could filter out secondary alignments: Regarding filter primary alignment

ADD REPLY
0
Entering edit mode

I will do that!

Thanks a lot for your help!!

ADD REPLY

Login before adding your answer.

Traffic: 3296 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6