BFast index creation & other SOLiD difficulties

[keebs42]

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.

[nilshomer]

Originally Posted by keebs42

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

[keebs42]

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

[nilshomer]

Originally Posted by keebs42

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.

[keebs42]

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!

[javijevi]

Originally Posted by nilshomer

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?

[drio]

Originally Posted by javijevi

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.

[javijevi]

Originally Posted by drio

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.

[nilshomer]

Originally Posted by javijevi

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.

[der_eiskern]

Originally Posted by nilshomer

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