Piping bowtie2 output directly into BAM

Question

[hersh]

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.

[GenoMax]

A: Alignment with bowtie2 but empty .fastq files

[swbarnes2]

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.

Answer 1

[cmdcolin]

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.

[hersh]

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

[MBarcelo]

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!

[GenoMax]

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.

[MBarcelo]

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!!

[GenoMax]

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.

[MBarcelo]

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!

[GenoMax]

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

[MBarcelo]

I will do that!

Thanks a lot for your help!!