Bwa Aln Output In .Bam

Question

[Dawnoflife]

How can I use the bwa aln function and get the output in .bam format.

I know bwasw gives the output in .sam so I just used SAMTOOLS to do the conversion but don't really know how to do it in this case as the output is in .sai.

What I have been trying to do is this:

1. Create the indexes for all the .FNA files.

2. Align them using aln command.

3. Take this output and convert it to .sam.

4. The sam files are very large to store. Hence convert them to .bam.

5. Data abstraction, like counting the ranges of values can't be performed on the binary (.bam) file, hence I converted that to a .txt file using SAMTOOLS.

6. Perform the desired data abstraction on the txt and get the results.

The tricky part is the placement of the files as I have to do this for about 1500 files. My code looks like this:

.././bwa index NC_008765.fna

.././bwa aln NC_008765.fna ../QUERY/updated_635_25bp.fasta > out.sai

../bwa samse NC_008765.fna out.sai ../QUERY/updated_635_25bp.fasta > aln.sam


../samtools-0.1.18/samtools-0.1.18/samtools view -bS -F 4 aln.sam > output.bam


../samtools-0.1.18/samtools-0.1.18/samtools view output2.bam > seemee.txt

awk '{ print $1 }' seemee.txt > comns.txt

sed 's/.\?.$//;s/$/00/;s/^00$/0/' comns.txt | sort | uniq -c | sort -nr > final.txt

I need to do this for all the 1000+ input files and 6 different query files. Does the general steps make sense? Can be it optimized? the code has been running on for 2 hours now and its only about 1/10th of its way through. Thanks for your help!

[toni]

Dear dawnoflife, you do not need to convert to BAM to, in the end, convert it to a text file. SAM is already a text file, so do what you want to do with your seemee.txt file directly with your aln.sam file.

[Dawnoflife]

SAM files are too large. they take up 36 mb as compared to < 1 mb for text files. unless I can perform the -F 4 operation when receiving the sam file itself.

[toni]

So just do not use "-bS -F 4" but "-S -F 4" , you're output will still be in SAM format (test), and then use it like your seemee.txt file

[toni]

you could then convert to BAM when you're done for persistent storage.

[Dawnoflife]

Is there any way I can pipe the sam file to bam directly and do the data analysis in b/w. I don't want to store the sam file.

Answer 1

[Vikas Bansal]

Use ->

bwa samse <in.db.fasta> <in.sai> <in.fq> > <out.sam> for single end reads.

bwa sampe <in.db.fasta> <in1.sai> <in2.sai> <in1.fq> <in2.fq> > <out.sam> for paired end

You can use some other parameters which are described here.

[Dawnoflife]

So we don't use the indexes created with bwa index -a

[Vikas Bansal]

Yes, the first step is to create index. and then you have 2 options -> aln and bwasw, depending on the type of reads. Since you used aln, it will give you binary file and then you have to use samse or sampe to get sam file.

[Dawnoflife]

Thanks for the help!

[Dawnoflife]

Actually, if I am trying to convert this bam file back to txt, it is giving me problems. The code below works when I convert the bam files obtained from

 bwasw .


for file in *.bam;  do /samtools-0.1.18/samtools-0.1.18/samtools view $file > ../$file.txt; done


The error I get is:
 EOF marker is absent. The input is probably truncated .


Will really appreciate the  help!

[Dawnoflife]

Edited the question, I am having troubles converting the bam to .txt. I need to get data from the columns and I was planning on using

awk for it. Thanks for your help!

[Vikas Bansal]

Have a look at Tony's answer. When you are using samtools view, you have to tell it that your file is in sam format (because by default it take input as bam) by adding -S argument.

[Dawnoflife]

I did use the -bS flag when using SamTools. Also how do I put code in the comments? It doesn't seem to work.

[toni]

Not -bS , just -S because -b would output in BAM !!! Please read the samtools documentation calmly :)

Answer 2

[toni]

'bwa aln' does not output in BAM format but in SAM which already is what you called a 'TXT' version of the BAM (binary format).

After that, [?]samtools[?] expects an input in BAM format so you have to specify the '-S' option to let it know that your input is in SAM format.

Synthesis :

bwa sampe <in.db.fasta> <in1.sai> <in2.sai> <in1.fq> <in2.fq> > <out.sam>

samtools view -S out.sam

You should be OK now.

[Dawnoflife]

Is <in.db.fasta> the indexes created with

bwa index ? When i was just using by .fna file for this it did not accept the input. 

[toni]

NO. As it is suggested by the name, it expects a FASTA file. This is the fasta sequence of your reference genome and this is already an input required by 'bwa index'.

[Vikas Bansal]

<in.db.fasta> is your reference genome (eg hg18, hg19 in fasta format). The first step is to create index. Use the same file in bwa samse. It should be fine if all other index files are in same directory.

[Dawnoflife]

Can I just use .fna since its a FASTA file as well. I'll try it again but the .fna didn't work for me.

[toni]

We cannot help you if you do not clarify what is not working.

[Dawnoflife]

Please see my edited question!

Answer 3

[Swbarnes2]

Are you doing something other than this to make out.bam?

bwa index ref.fa; 
bwa aln ref.fa read1.fq > r1.sai; 
bwa aln ref.fa read2.fq > r2.sai; 
bwa sampe ref.fa r1.sai r2.sai read1.fa read2.fq | samtools view -bSho out.bam -;

samtools view -S will do nothing. It will take a .sam file as input and output another .sam file.

(If you have a large index, like human, there's a slightly different index command)

[toni]

yes true, I wrote this without thinking but a simple 'cat' will then be enough

[Dawnoflife]

Kindly see the edited question!

[Swbarnes2]

Have you tried something like:

samtools view out.bam | awk '{print $1}' | sed 's/.\?.$//;s/$/00/;s/^00$/0/' - | sort | uniq -c | sort -nr > final.txt

And wouldn't it be easier to make an altered fastq and realign rather than altering the .bams?

[Swbarnes2]

Have you tried something like: samtools view out.bam | awk '{print $1}' | sed 's/.\?.$//;s/$/00/;s/^00$/0/' | sort | uniq -c | sort -nr > final.txt And wouldn't it be easier to make an altered fastq and realign rather than altering the .bams?