Header Leaderboard Ad

Collapse

BFast index creation & other SOLiD difficulties

Collapse

Announcement

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

  • BFast index creation & other SOLiD difficulties

    Hi all,

    I've been having lots of trouble getting 36 million SOLiD reads aligned to hg18. I started using bwa, but less than 40% of my reads align, which seems pretty low. Posts on here suggesting there are bugs in the SOLiD module of BWA have made me look towards other solutions. The fastq file for the bwa input was created using the csfastaToFastq script.

    Here's the bwa aln command Ive been using:
    bwa aln -c hg18.cs.fa reads.part1.fastq


    I've also tried bowtie, and can get over 54% of the reads to align, however when using samtools to import the resulting sam file I'm getting this:

    [samopen] SAM header is present: 25 sequences.
    Parse error at line 28: sequence and quality are inconsistent

    Here's my bowtie command and the samtools import command:
    bowtie_0.12.2/bowtie -t -C -S --nomaqround --best --snpfrac .001 -n 2 --chunkmbs 256 -e 90 -l 28 -p 6 hg18_c bfast_built_validation_reads_91.fastq bowtie_run1.sam

    samtools import /home/jekeeble/storage/hg18/hg18.fa.fai ./bowtie_run1.sam bowtie_run1.bam

    As you might guess, the fastq file used for the bowtie input was built using the bfast solid2fastq program.

    Finally, i've also been trying to use BFast based on lh3's mention that many people are converging to that or bioscope for solid data. I can't seem to get beyond getting the bfast index built.

    running:
    bfast index -f hg18.fa -w 14 -n 5 -A 1 -m 1011..111 -i 2

    produces this output:
    ************************************************************
    Checking input parameters supplied by the user ...
    Validating fastaFileName hg18.fa.
    Validating tmpDir path ./.
    Input arguments look good!
    ************************************************************
    ************************************************************
    Printing Program Parameters:
    programMode: [ExecuteProgram]
    fastaFileName: hg18.fa
    space: [Color Space]
    mask: 10111111011001100011111000111111
    depth: 0
    hashWidth: 14
    indexNumber: 2
    repeatMasker: [Not Using]
    startContig: 0
    startPos: 0
    endContig: 2147483647
    endPos: 2147483647
    exonsFileName: [Not Using]
    numThreads: 5
    tmpDir: ./
    timing: [Not Using]
    ************************************************************
    ************************************************************
    Reading in reference genome from hg18.fa.cs.brg.
    In total read 25 contigs for a total of 3080436051 bases
    ************************************************************
    Creating the index...
    ************************************************************
    Warning: startContig was less than zero.
    Defaulting to contig=1 and position=1.
    ************************************************************
    ************************************************************
    Warning: endContig was greater than the number of contigs in the reference genome.
    Defaulting to reference genome's end contig=25 and position=16571.
    ************************************************************
    Currently on [contig,pos]:
    [------10,---68000000]
    [------25,------16571]
    Sorting by thread...
    bfast: RGIndex.c:743: RGIndexSort: Assertion `IsAPowerOfTwo(numThreads)==1' failed.
    Aborted

    Anyone have suggestions to where I should go from here? It's been a pretty frustrating day . I"m also running shrimp, but it's taking its sweet time.

    Thanks,
    Jonathan K.

  • #2
    Originally posted by keebs42 View Post
    Finally, i've also been trying to use BFast based on lh3's mention that many people are converging to that or bioscope for solid data. I can't seem to get beyond getting the bfast index built.

    running:
    bfast index -f hg18.fa -w 14 -n 5 -A 1 -m 1011..111 -i 2

    produces this output:
    ************************************************************
    Checking input parameters supplied by the user ...
    Validating fastaFileName hg18.fa.
    Validating tmpDir path ./.
    Input arguments look good!
    ************************************************************
    ************************************************************
    Printing Program Parameters:
    programMode: [ExecuteProgram]
    fastaFileName: hg18.fa
    space: [Color Space]
    mask: 10111111011001100011111000111111
    depth: 0
    hashWidth: 14
    indexNumber: 2
    repeatMasker: [Not Using]
    startContig: 0
    startPos: 0
    endContig: 2147483647
    endPos: 2147483647
    exonsFileName: [Not Using]
    numThreads: 5
    tmpDir: ./
    timing: [Not Using]
    ************************************************************
    ************************************************************
    Reading in reference genome from hg18.fa.cs.brg.
    In total read 25 contigs for a total of 3080436051 bases
    ************************************************************
    Creating the index...
    ************************************************************
    Warning: startContig was less than zero.
    Defaulting to contig=1 and position=1.
    ************************************************************
    ************************************************************
    Warning: endContig was greater than the number of contigs in the reference genome.
    Defaulting to reference genome's end contig=25 and position=16571.
    ************************************************************
    Currently on [contig,pos]:
    [------10,---68000000]
    [------25,------16571]
    Sorting by thread...
    bfast: RGIndex.c:743: RGIndexSort: Assertion `IsAPowerOfTwo(numThreads)==1' failed.
    Aborted

    Anyone have suggestions to where I should go from here? It's been a pretty frustrating day . I"m also running shrimp, but it's taking its sweet time.

    Thanks,
    Jonathan K.
    As the error indicates from BFAST, the # of threads must be a power of two. Can you try with 4 or 8 threads (match it to the number of cores)? Remember, you only need to index once

    Nils
    Last edited by nilshomer; 02-07-2010, 05:15 PM. Reason: I can't spell

    Comment


    • #3
      Yeah I thought that might be the case, but apparently my head was so far into the problem I didn't notice I was requesting 5 threads. Bumping it up to 8 should take care of it. Thanks!

      Jonathan

      Comment


      • #4
        Originally posted by keebs42 View Post
        Yeah I thought that might be the case, but apparently my head was so far into the problem I didn't notice I was requesting 5 threads. Bumping it up to 8 should take care of it. Thanks!

        Jonathan
        I should probably code it to "round down" to the nearest power of two. It's there for a parallel index sort.

        Comment


        • #5
          That would work, or just make a note of it in the Bfast book in section 3.3.1 or 3.3.2. Also the tutorial in 7.7.1 and 7.7.2 don't have any mention of an even thread count. If -n isn't specified does it default to 2?

          Otherwise I've found your documentation pretty helpful. Great work!

          Comment


          • #6
            Originally posted by nilshomer View Post
            As the error indicates from BFAST, the # of threads must be a power of two. Can you try with 4 or 8 threads (match it to the number of cores)? Remember, you only need to index once

            Nils
            I used -n 6 (my box has 8 cores, but I wanted to have 2 of them free), and also get an error indicating that number of thread must be a power of two. Why do you propose 4 or 8 threads in the previous message? Isn't 6 cores a valid power of two for bfast?

            Comment


            • #7
              Originally posted by javijevi View Post
              I used -n 6 (my box has 8 cores, but I wanted to have 2 of them free), and also get an error indicating that number of thread must be a power of two. Why do you propose 4 or 8 threads in the previous message? Isn't 6 cores a valid power of two for bfast?
              6 is not a power of two.
              -drd

              Comment


              • #8
                Originally posted by drio View Post
                6 is not a power of two.
                Oops... I feel my ears growing like a donkey... My poor english... I misunderstood 'power of two'. I thought it was 'multiple of two'. I strongly apologize.

                Comment


                • #9
                  Originally posted by javijevi View Post
                  Oops... I feel my ears growing like a donkey... My poor english... I misunderstood 'power of two'. I thought it was 'multiple of two'. I strongly apologize.
                  No problem. It only needs to be a power of two for the index creation step.

                  Comment


                  • #10
                    Originally posted by nilshomer View Post
                    As the error indicates from BFAST, the # of threads must be a power of two. Can you try with 4 or 8 threads (match it to the number of cores)? Remember, you only need to index once

                    Nils
                    perhaps in the next release the error report could look a little bit more like this

                    Comment

                    Working...
                    X