Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • 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
    -@, --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)

  • #2
    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.

    Comment


    • #3
      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:

      Originally posted by zhanghao View Post
      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
      Code:
      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

      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

      Comment


      • #4
        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.

        Comment


        • #5
          @Brian: bs.sh is no longer in BBMap.

          Comment


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

            Comment


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

              Comment


              • #8
                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, 10:06 AM.

                Comment


                • #9
                  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.

                  Comment


                  • #10
                    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...

                    Comment


                    • #11
                      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-10-2016, 12:23 AM.

                      Comment

                      Latest Articles

                      Collapse

                      • seqadmin
                        Techniques and Challenges in Conservation Genomics
                        by seqadmin



                        The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                        Avian Conservation
                        Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                        03-08-2024, 10:41 AM
                      • seqadmin
                        The Impact of AI in Genomic Medicine
                        by seqadmin



                        Artificial intelligence (AI) has evolved from a futuristic vision to a mainstream technology, highlighted by the introduction of tools like OpenAI's ChatGPT and Google's Gemini. In recent years, AI has become increasingly integrated into the field of genomics. This integration has enabled new scientific discoveries while simultaneously raising important ethical questions1. Interviews with two researchers at the center of this intersection provide insightful perspectives into...
                        02-26-2024, 02:07 PM

                      ad_right_rmr

                      Collapse

                      News

                      Collapse

                      Topics Statistics Last Post
                      Started by seqadmin, 03-14-2024, 06:13 AM
                      0 responses
                      33 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, 03-08-2024, 08:03 AM
                      0 responses
                      72 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, 03-07-2024, 08:13 AM
                      0 responses
                      81 views
                      0 likes
                      Last Post seqadmin  
                      Started by seqadmin, 03-06-2024, 09:51 AM
                      0 responses
                      68 views
                      0 likes
                      Last Post seqadmin  
                      Working...
                      X