BWA mem and Samtools sort in Bash script
0
1
Entering edit mode
2.2 years ago

Hi, Biostars community.

I'm trying to use both BWA mem and Samtools sort at the same time via passing BWA mem output into samtools.

Command:

bwa mem reference.fa file.fq | samtools sort -o output.bam -


Well, it works fine while I'm using terminal, but the same command in bash script throws errors.

#!/bin/bash

bwa mem reference.fa file.fq | samtools sort -o output.bam -


Error message:

script.sh: line 3: 18: command not found
[main] unrecognized command 'sort'


Can someone show me the correct way to pass bwa output to samtools in bash script?

alignment bwa samtools bash • 1.8k views
0
Entering edit mode

do you have to load samtools first, e.g. via spack or the like? I.e. when you log into your server can you use samtools right away?

also, can you post the entire script?

0
Entering edit mode

I've posted all script. See "#!/bin/bash"

Well, yes, samtools is alredy installed and in works fine. I guess there is a different way to pass bwa output to samtools is bash script, not via vertical line | symbol.

0
Entering edit mode

The pipe | symbol should work. I can't currently replicate your issue. Can you provide more information about the environment you're working in (workstation/cluster/scheduler/HPC?), and the versions of bash, bwa, samtools you're using?

0
Entering edit mode

It works fine in terminal, but not when I try to call the script.

Ubuntu Server 18.04 bwa: 0.7.17-r1198-dirty samtools: 1.9

0
Entering edit mode

I think you need to also add -O BAM argument to write the output to file

0
Entering edit mode

No need for -O BAM, as recent versions of samtools will recognize the file format based on the extension passed with -o.

0
Entering edit mode

Thanks WoterDeCoster for this information

0
Entering edit mode

what is the output of

file script.sh


?

0
Entering edit mode

Do the commands work inside a script if you separate them in to different lines? i.e. output a file with bwa then read it separately with samtools sort.

There is no reason your pipe syntax shouldn't work, so there is something more 'sinister' going wrong behind the scenes I think.

0
Entering edit mode

but the same command in bash script throws errors.

How did you execute that script?

0
Entering edit mode

What makes this question interesting/amusing is that it takes a little bit of effort to write a code that deliberately produces the error. It is like a weird code golf: Write the shortest bash code that produces the error message shown above.

Here is my take:

bash script.sh


it prints:

script.sh: line 3: 18: command not found
[main] unrecognized command 'sort'


I feel I learned something nontrivial especially when it comes to reading error messages.

What is in my script? Is it possible to make it shorter?

wc script.sh


three lines, four words, 16 characters:

3       4      16 script.sh


As a response to the original poster: I believe that what ends up running is something completely different than what you think/show as being run.

0
Entering edit mode

To reduce the suspense:

cat script.sh


prints:

18 | bwa sort


Notably this example also helped us discover a bug in the markdown rendering for the javascript preview as well as the backend of the site.

The libraries we use won't properly render leading empty lines (two empty lines to push the error to line 3), hence I can't quite demonstrate my code correctly.