BBmap can not produce sorted bam

[zhanghao]

I used BBmap to map RNAseq data with the command:
/home/bac/software/bbmap/bbmap.sh in1=CYP51C_2_1.fq.gz in2=CYP51C_2_2.fq.gz ref=genome.fasta nodisk out=cyp51c_bb.sam bamscript=bs.sh; sh bs.sh

samtools-1.3.1 is setted in PAYH. But it just produced sam file, seems samtools did not work. It showed:

bs.sh: 2: bs.sh: module: not found
bs.sh: 3: bs.sh: module: not found
Note: This script is designed to run with the amount of memory detected by BBMap.
If Samtools crashes, please ensure you are running on the same platform as BBMap,
or reduce Samtools' memory setting (the -m flag).
Note: Please ignore any warnings about 'EOF marker is absent'; this is a bug in samtools that occurs when using piped input.
[bam_sort] Use -T PREFIX / -o FILE to specify temporary and final output files
Usage: samtools sort [options...] [in.bam]
Options:
-l INT Set compression level, from 0 (uncompressed) to 9 (best)
-m INT Set maximum memory per thread; suffix K/M/G recognized [768M]
-n Sort by read name
-o FILE Write final output to FILE rather than standard output
-T PREFIX Write temporary files to PREFIX.nnnn.bam
-@, --threads INT
Set number of sorting and compression threads [1]
--input-fmt-option OPT[=VAL]
Specify a single input file format option in the form
of OPTION or OPTION=VALUE
-O, --output-fmt FORMAT[,OPT[=VAL]]...
Specify output format (SAM, BAM, CRAM)
--output-fmt-option OPT[=VAL]
Specify a single output file format option in the form
of OPTION or OPTION=VALUE
--reference FILE
Reference sequence FASTA FILE [null]
[E::hts_open_format] fail to open file 'cyp51c_bb_sorted.bam'
samtools index: failed to open "cyp51c_bb_sorted.bam": No such file or directory

I checked the bs.sh file:

#!/bin/bash
module unload samtools
module load samtools/0.1.19
echo "Note: This script is designed to run with the amount of memory detected by BBMap."
echo " If Samtools crashes, please ensure you are running on the same platform as BBMap,"
echo " or reduce Samtools' memory setting (the -m flag)."
echo "Note: Please ignore any warnings about 'EOF marker is absent'; this is a bug in samtools that occurs when using piped input."
samtools view -bSh1 cyp51c_bb.sam | samtools sort -m 13G -@ 3 - cyp51c_bb_sorted
samtools index cyp51c_bb_sorted.bam
bs.sh (END)

[zhanghao]

another question is:

I have 3 pairs of pair-end fastq data sets, how to set them as input? In hisat2, I can use "-1 read_A_1.fq,reads_B_1.fq,reads_C_1.fq -2 read_A_2.fq,reads_B_2.fq,reads_C_2.fq", in bbmap, it seems just in1 and in2 are supported.

[antoinefelden]

I have the same problem, looking at the SAMTools manual, I have the feeling there is a problem on the output option of this command:

samtools view -bSh1 cyp51c_bb.sam | samtools sort -m 13G -@ 3 - cyp51c_bb_sorted
samtools index cyp51c_bb_sorted.bam
bs.sh (END)

Could it be missing a "-o" before the filename? Have a look at the SAMTools page here : http://www.htslib.org/doc/samtools.html

Then it would be

samtools view -bSh1 cyp51c_bb.sam | samtools sort -m 13G -@ 3 -[B]o[/B] cyp51c_bb_sorted
samtools index cyp51c_bb_sorted.bam
bs.sh

However, I could not find any trace of this bs.sh file (I could only find it in my output), so I'm a bit puzzled about how to fix this.

Also, about your second question, you have to use bbwrap instead of bbmap, that will do the job I think

bbwrap.sh in1=read_A_1.fq,reads_B_1.fq,reads_C_1.fq in2=read_A_2.fq,reads_B_2.fq,reads_C_2.fq ref=genome.fasta nodisk out=cyp51c_bb.sam bamscript=bs.sh; sh bs.sh

[Brian Bushnell]

another question is:

I have 3 pairs of pair-end fastq data sets, how to set them as input? In hisat2, I can use "-1 read_A_1.fq,reads_B_1.fq,reads_C_1.fq -2 read_A_2.fq,reads_B_2.fq,reads_C_2.fq", in bbmap, it seems just in1 and in2 are supported.

Hi! Sorry, somehow I missed your post. BBMap only supports a single library per run; if you have different libraries with different insert sizes, I recommend NOT combining them, since BBMap's alignment score is affected by how close the insert size of a pair is to the average.

But you can use bbwrap.sh, as suggested by antoinefelden.

As for the failures of "bs.sh" - this is strange. I'll look into it tomorrow. Maybe it's a difference between Samtools v0.x versus Samtools 1.x - some of the behavior was altered; and bs.sh was designed for Samtools 0.1.19.

[GenoMax]

@Brian: bs.sh is no longer in BBMap.

[Brian Bushnell]

I just used it yesterday, and it worked fine for me... and it's documented on line 159 of bbmap.sh...?

[GenoMax]

Ah. It is inside bbmap.sh? I was expecting to see it at top level.

[Brian Bushnell]

Oh! Yes, technically bs.sh is not part of the BBMap package. Rather, it gets produced by BBMap if you add the flag "bs=bs.sh", because it is a custom script specific to the name of the sam/bam file BBMap produced. bs.sh simply calls Samtools with the correct (or, in the above case, apparently incorrect) commands to produce a sorted, indexed bam file.

Generally, you do this:

bbmap.sh in=reads.fq.gz out=mapped.bam bs=bs.sh; sh bs.sh

That produces sorted_mapped.bam and sorted_mapped.bam.bai.

[antoinefelden]

bs.sh simply calls Samtools with the correct (or, in the above case, apparently incorrect) commands to produce a sorted, indexed bam file.

Hi Brian,

Would you say that the problem comes from the missing output option?
If so, do you know how to fix that error?

If that helps, it may also be related to the samtools version, I'm using version 1.3 and when executing the script bs.sh, the module commands on line 2 and 3 of bs.sh did not work (command not found error, which I don't understand why so far).

I'm re-running the whole job with samtools 0.1.19, I'll post the outcome here.

[Brian Bushnell]

OK! I've tested this with samtools 0.1.19 and 1.3.1, and it fails with samtools 1.3.1. Apparently this is caused by some basic parameters that were changed when samtools migrated from 0.x to 1.x; this is something that I advocate against, but it it is sometimes unavoidable. Sorry about that! I'll release a new version soon that will work with both samtools 0.x and samtools 1.x, if that's possible. In the meantime, of course, you can use samtools manually...

[Brian Bushnell]

I'm re-running the whole job with samtools 0.1.19, I'll post the outcome here.

You don't need to rerun the whole job... the only problematic part is sorting and indexing the mapped reads. If you use samtools 0.x, it will work fine using the existing output of BBMap. But you can also use samtools 1.x if you look into the manual and find whatever arcane invocation is appropriate now. Hopefully there is a command that will work with all versions of samtools; I'll adjust the "bs" output in the next release.