SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Smatools mpileup thinks my sorted bam files are not sorted plumb_r Bioinformatics 6 04-28-2014 08:22 AM
what's the difference between a bam and a sorted bam KnowNothing2 Bioinformatics 7 10-17-2013 04:00 PM
convert sorted bam to sorted sam for htseq-count ugolino Bioinformatics 3 10-19-2012 06:30 AM
sorted bam larger than unsorted bam ians Bioinformatics 4 05-30-2012 05:51 AM

Reply
 
Thread Tools
Old 10-25-2016, 05:19 PM   #1
zhanghao
Member
 
Location: BEIJING

Join Date: Nov 2015
Posts: 17
Default BBmap can not produce sorted bam

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
[email protected], --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 [email protected] 3 - cyp51c_bb_sorted
samtools index cyp51c_bb_sorted.bam
bs.sh (END)
zhanghao is offline   Reply With Quote
Old 10-25-2016, 05:26 PM   #2
zhanghao
Member
 
Location: BEIJING

Join Date: Nov 2015
Posts: 17
Default

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.
zhanghao is offline   Reply With Quote
Old 12-08-2016, 03:26 PM   #3
antoinefelden
Junior Member
 
Location: Wellington, NZ

Join Date: Dec 2016
Posts: 3
Default

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:

Quote:
Originally Posted by zhanghao View Post
samtools view -bSh1 cyp51c_bb.sam | samtools sort -m 13G [email protected] 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
Code:
samtools view -bSh1 cyp51c_bb.sam | samtools sort -m 13G [email protected] 3 -o 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

Code:
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
antoinefelden is offline   Reply With Quote
Old 12-08-2016, 05:05 PM   #4
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,695
Default

Quote:
Originally Posted by zhanghao View Post
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.
Brian Bushnell is offline   Reply With Quote
Old 12-09-2016, 03:22 AM   #5
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,492
Default

@Brian: bs.sh is no longer in BBMap.
GenoMax is offline   Reply With Quote
Old 12-09-2016, 05:32 AM   #6
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,695
Default

I just used it yesterday, and it worked fine for me... and it's documented on line 159 of bbmap.sh...?
Brian Bushnell is offline   Reply With Quote
Old 12-09-2016, 06:24 AM   #7
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,492
Default

Ah. It is inside bbmap.sh? I was expecting to see it at top level.
GenoMax is offline   Reply With Quote
Old 12-09-2016, 06:27 AM   #8
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,695
Default

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:

Code:
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.

Last edited by Brian Bushnell; 12-09-2016 at 09:06 AM.
Brian Bushnell is offline   Reply With Quote
Old 12-09-2016, 02:03 PM   #9
antoinefelden
Junior Member
 
Location: Wellington, NZ

Join Date: Dec 2016
Posts: 3
Default

Quote:
Originally Posted by Brian Bushnell View Post
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.
antoinefelden is offline   Reply With Quote
Old 12-09-2016, 02:38 PM   #10
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,695
Default

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 is offline   Reply With Quote
Old 12-09-2016, 02:57 PM   #11
Brian Bushnell
Super Moderator
 
Location: Walnut Creek, CA

Join Date: Jan 2014
Posts: 2,695
Default

Quote:
Originally Posted by antoinefelden View Post
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.

Last edited by Brian Bushnell; 12-09-2016 at 11:23 PM.
Brian Bushnell is offline   Reply With Quote
Reply

Tags
bbtools, samtools

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 02:37 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2017, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO