[seb567]

I am not a Python guy but it looks to me like it would be fairly straightforward using Biopython's Bio/SeqIO/SffIO module. Read the file in, flip the value of 'key_sequence', write out a new file.

(Still waiting for Bioperl Bio::SeqIO::SFF )

Regardless, I guess it is correct to consider SFF files as containers, just like FASTA or FASTQ files.

Therefore, Ray will no longer try to match the key sequence. Instead, it will *simply* load all sequences in the SFF file and trim them using the clipping values therein.

See http://github.com/sebhtml/ray/commit/15826e290f1

+1

I know if you want something done you should probably take the initiative and contribute, but I have seen several posts where people have said they were working on this. So, like many people, I decided to wait, assuming it was in progress. (Sorry for taking things off track in the thread though )

You mean taking the initiative to write code changes to BioPython so that Bio/SeqIO/SffIO can change the key sequence, right ?


Is there a software tool from 454 that allows one to change header information in a SFF file ?


There is also this thing called flower (the code is pretty awesome by the way -- it is in Haskell)

Blog post: http://blog.malde.org/index.php/flower/
Source code: http://malde.org/~ketil/biohaskell/flower/


Also, the Ray git tree is now on github.

http://github.com/sebhtml/ray

Furthermore, Ray can now handle arbitrary large k-mers.

I am presently running some integration and unit tests on Ray v1.6.0-rc2.

You can download the latest development version of Ray with the following command *provided* that you have git.

git clone git://github.com/sebhtml/ray.git

To use large k-mers:

git clone git://github.com/sebhtml/ray.git
cd ray
make MAXKMERLENGTH=64 PREFIX=ray-git-master-kmax=64
make install
mpirun -np 128 ray-git-master-kmax=64/Ray -k 55 \
-p ABCD_1.fastq ABCD_2.fastq -o DeadlyBug,k=55

Enjoy !

[SES]

(Still waiting for Bioperl Bio::SeqIO::SFF )

+1

I know if you want something done you should probably take the initiative and contribute, but I have seen several posts where people have said they were working on this. So, like many people, I decided to wait, assuming it was in progress. (Sorry for taking things off track in the thread though )

[kmcarr]

Exactly my point.

Meanwhile, what do you think would be the best way to deal with these ill-encoded SFF files generated by sfffile <2.5 with the rapid library chemistry ?

I just don't see an easy way.

I am not a Python guy but it looks to me like it would be fairly straightforward using Biopython's Bio/SeqIO/SffIO module. Read the file in, flip the value of 'key_sequence', write out a new file.

(Still waiting for Bioperl Bio::SeqIO::SFF )

[seb567]

Sébastien, the 454 reads produced today will in many cases have the new key sequence GACT instead of TCAG. New library preparation kits (using the so-called 'Rapid Library' protocol) have this new key in the adaptors. It would be a great advantage if Ray could handle both key sequences!

Ray simply fetches the key sequence from the SFF header. Ray has no preference for GACT or TCAG.

I believe the problem is due to a bug in an earlier version of the Roche/454 software. With the switch to Rapid Library chemistry Roche switched the keytag to GACT and released new software (?2.3?). The gsRunProcessor produced properly formatted SFF files which reported GACT as the keytag in the common header section of the SFF. However if you used the program sfffile to manipulate those SFFs (e.g. decode MID tags, split files or merge files) the new common header would erroneously report TCAG as the keytag. This bug appears to have been corrected in the latest release (2.5) of sfffile.

What Sébastien seems to be saying is that Ray reads the common header of the SFF to determine what the keytag should be and in this case there is a mismatch between what the header reports the keytag to be and the keytag observed in the reads. It seems that figure002's SFF file(s) have fallen victim to this bug in sfffile.

Exactly my point.

Meanwhile, what do you think would be the best way to deal with these ill-encoded SFF files generated by sfffile <2.5 with the rapid library chemistry ?

I just don't see an easy way.

[kmcarr]

Sébastien, the 454 reads produced today will in many cases have the new key sequence GACT instead of TCAG. New library preparation kits (using the so-called 'Rapid Library' protocol) have this new key in the adaptors. It would be a great advantage if Ray could handle both key sequences!

I believe the problem is due to a bug in an earlier version of the Roche/454 software. With the switch to Rapid Library chemistry Roche switched the keytag to GACT and released new software (?2.3?). The gsRunProcessor produced properly formatted SFF files which reported GACT as the keytag in the common header section of the SFF. However if you used the program sfffile to manipulate those SFFs (e.g. decode MID tags, split files or merge files) the new common header would erroneously report TCAG as the keytag. This bug appears to have been corrected in the latest release (2.5) of sfffile.

What Sébastien seems to be saying is that Ray reads the common header of the SFF to determine what the keytag should be and in this case there is a mismatch between what the header reports the keytag to be and the keytag observed in the reads. It seems that figure002's SFF file(s) have fallen victim to this bug in sfffile.

[flxlex]

Maybe they changed the SFF standard or you are using multiplex identifiers. In either case, I suggest you convert your SFF files to FASTA (or FASTQ) and supply the resulting files to Ray instead of the SFF files.

Sébastien, the 454 reads produced today will in many cases have the new key sequence GACT instead of TCAG. New library preparation kits (using the so-called 'Rapid Library' protocol) have this new key in the adaptors. It would be a great advantage if Ray could handle both key sequences!

[seb567]

Dear Sébastien,

I recently started testing Ray on read data of a ~450Mb genome. I can load FASTQ files without a problem, but I can't figure out how to load SFF files. The Instruction Manual only has examples for loading FASTQ files.

I just started a job with the following command,

mpirun -np 16 time Ray \
-s /home/sp/data/454/shotgun/F0A0H9G01.sff \
-i /home/sp/data/454/pairedend_20k/FPFSKVK01.sff \
-i /home/sp/data/454/pairedend_3k/FO2K76101.sff \
-k 17 -o melon_454_small_test_20110604

But the output contains a lot of the following,

$ tail mpirun.o22658
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG

So I'm probably doing something wrong. Could you please explain how I should do this?

In the SFF specification, it is said that the header of the file contains the prefix of sample sequences.

The message you encountered means that your SFF file contains sequence reads with a sequence key that does not match with the one listed in the header.

Maybe they changed the SFF standard or you are using multiplex identifiers. In either case, I suggest you convert your SFF files to FASTA (or FASTQ) and supply the resulting files to Ray instead of the SFF files.

Sébastien

[figure002]

How to load SFF files?

Dear Sébastien,

I recently started testing Ray on read data of a ~450Mb genome. I can load FASTQ files without a problem, but I can't figure out how to load SFF files. The Instruction Manual only has examples for loading FASTQ files.

I just started a job with the following command,

mpirun -np 16 time Ray \
-s /home/sp/data/454/shotgun/F0A0H9G01.sff \
-i /home/sp/data/454/pairedend_20k/FPFSKVK01.sff \
-i /home/sp/data/454/pairedend_3k/FO2K76101.sff \
-k 17 -o melon_454_small_test_20110604

But the output contains a lot of the following,

$ tail mpirun.o22658
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG
Not KEY, was GACT expected TCAG

So I'm probably doing something wrong. Could you please explain how I should do this?

[seb567]

Ray 1.4.0: built-in scaffolder &amp; more

Dear Ray users,


Ray 1.4.0 is now available.


The most significant change is the built-in scaffolder.



The second most significant change is the new algorithm that finds
assembly seeds.

Also, I added a lot of output files in Ray.

They are listed here:

http://denovoassembler.sourceforge.net/manual.html#SECTION00050000000000000000

Finally, our new website http://denovoassembler.sf.net is hopefully
easier to browse.

On the website, there is a manual for Ray.

http://denovoassembler.sourceforge.net/manual.html

Sébastien

1.4.0
2011-05-30

* A built-in scaffolder is now available -- Thanks to Dr.
Jean-Francois Pombert (University of British Columbia) for the
suggestion.
* The maximum number of libraries is now 499 instead of 250.
* The number of seeds is now divided by 2 to speed up their
extension.
* Fixed a bug in the depth first search that leaded to vertices
having no coverage values.
* Removed the configure script, now Ray must be compiled with the
provided Makefile.
* Added a switch to enable the profiler: -run-profiler
* Added a switch to debug seed generation: -debug-seeds
* Added a switch to debug bubble detection: -debug-bubbles
* Added a switch to show memory usage: -show-memory-usage
* Added a switch to show the ending context of extensions:
-show-ending-context
* Devised a new algorithm that finds the peak coverage, minimum
coverage and repeat coverage in distributions.
* Ray now writes the peak, minimum and repeat coverages to a file.
* Ray now writes the statistics for libraries to a file.
* Fixed a bug that disallowed mixing manual and automatic
detection of outer distances.
* Ray now writes the statistics for seed lengths to a file.
* Devised a new algorithm that computes longer seeds to bootstrap
assemblies.
* Slave modes, master modes and MPI tags are generated with macros
for method prototypes, enumerations and assignments in arrays.
* Added some changes for Microsoft Windows compatibility. Thanks
to Hannes Pouseele (Applied Maths, Inc.) for some suggestions.
* Added instructions regarding mpic++ and CXX environment
variable. Thanks to Dr. Harry Mangalam from UC Irvice for
pointing that out.
* Changed the merger behavior for ends of contigs.
* Added a script to validate scaffolds.

[seb567]

Hi seq567 do you think it's a wonderful or terrible idea to try Ray on RNASeq data since coverage is not uniform?

Worth a try although I believe the coverage distribution may have more than one peak because, as you said, non-uniformity is present.

You don't do bubble bursting, but when choosing paths could non-uniform coverage create mis-assemblies?

Choosing the next vertex to visit is not performed with the coverage. Instead, pairs of sequences are utilised for that purpose.

Bubbles are not bursted, but when Ray encounters one (in general, an heterozygous site), it will attempt to pass through it by selecting one child.

In my opinion, bubble bursting and bubble merging are bad approaches because you either lose information (bursting) or create misassemblies that are otherwise not so hard to avoid (merging).

(I could read all the source, but I thought it would be faster this way :-) )

However, the bits related to message passing may interfere with your reading.

Things are easier to grasp in natural languages than in computer languages !

***
Sébastien Boisvert

http://boisvert.info

http://twitter.com/sebhtml

http://denovoassembler.sf.net

http://genome.ulaval.ca/corbeillab

[lletourn]

Hi seq567 do you think it's a wonderful or terrible idea to try Ray on RNASeq data since coverage is not uniform?

You don't do bubble bursting, but when choosing paths could non-uniform coverage create mis-assemblies?

(I could read all the source, but I thought it would be faster this way :-) )

[seb567]

Cool, I was looking forward for a MPI-based assembler like this. Will have a try now. Thanks a lot!

I hope you enjoy Ray.

Hi Sebastien,

thank you very much for this tool. I've been using several assemblers to process solexa and roche data on 3 bacterial genomes and your tool is one of the best performers. In fact, with my solexa datasets is the best . Anyway, I wonder if you could include more output formats in a next release. More precisely, I really wanted to have FASTQ output files.
Indeed, this would not be needed if there was an (almost) universal sequence/assembly format conversion tool. There are several around but are often very targeted.

cheers

Ray is very good when errors are randomly occuring. 454's homopolymers are not random.


I guess you would like the nucleotides to be accompanied by quality values in the output. Is that so ?

Hi seb567, glad to have met you in Quebec at the Rendez-vous. I'm playing with Ray on some of my "troublesome" species and I'm having very good results.

Like corthay, I was wondering if you had any plans to support bigger kmers than 31 either through a compiler switch or something else.

One could argue that past a point (kmer > 100) an overlap based assembler might do a very good job, but where speed is concerned graph based assemblers are hard to beat.

Also, with illumina that now has high 3' quality long reads, maybe bigger kmer can be justified.

Thanks

I agree with you. I will check how easy (or hard) it is for me to change the upper bound for k-mer length in Ray.

seb, Ray looks like a great assembly tool but I have a question -- could I expect Ray to do a better job than Newbler if the input reads are 454 only (single and paired end)? Or is the advantage mainly seen when mixed platform reads are input (e.g. 454 and Illumina)?

With only 454 reads, Ray will not produce a very good assembly because the errors are not randomly occuring in reads. Ray will go through some of these pesky homopolymer bubbles in the graph.

However, should you mix 454 with another sequencing technology with random errors, you will obtain a better that one with only 454 alone.

I would also like to put in a request for longer kmers. Assembling paired illumina reads of 100 bp (forward) and 83 bp (reverse, quality-trimmed) and scaffolding with SSPACE, I get 537 scaffolds with an N50 of 19163 when I use the default kmer setting, but 254 scaffolds with an N50 of 52083 when I use -k 31. Velvet Optimiser using the same data picks a kmer length of 67. I can only guess how awesome an assembly Ray could give me with the ability to use longer kmers...

I understand, I will see how easy (or hard) it is to change the maximum k-mer length in Ray.




Thank you for your interest in Ray !

p.s.: I am presently optimising Ray by studying its communicational behavior.

http://i.imgur.com/g9NFV.png

x axis: time; blue: number of iterations; red: number of sent messages; green: number of received messages; data are collected 10 times per second

I am also in the process of devising a better way to compute seeds -- with are relatively similar to the concept of unipaths observed in classic assemblers.

[Adjuvant]

Like corthay, I was wondering if you had any plans to support bigger kmers than 31 either through a compiler switch or something else.

One could argue that past a point (kmer > 100) an overlap based assembler might do a very good job, but where speed is concerned graph based assemblers are hard to beat.

I would also like to put in a request for longer kmers. Assembling paired illumina reads of 100 bp (forward) and 83 bp (reverse, quality-trimmed) and scaffolding with SSPACE, I get 537 scaffolds with an N50 of 19163 when I use the default kmer setting, but 254 scaffolds with an N50 of 52083 when I use -k 31. Velvet Optimiser using the same data picks a kmer length of 67. I can only guess how awesome an assembly Ray could give me with the ability to use longer kmers...

[ssully]

seb, Ray looks like a great assembly tool but I have a question -- could I expect Ray to do a better job than Newbler if the input reads are 454 only (single and paired end)? Or is the advantage mainly seen when mixed platform reads are input (e.g. 454 and Illumina)?

[lletourn]

Hi seb567, glad to have met you in Quebec at the Rendez-vous. I'm playing with Ray on some of my "troublesome" species and I'm having very good results.

Like corthay, I was wondering if you had any plans to support bigger kmers than 31 either through a compiler switch or something else.

One could argue that past a point (kmer > 100) an overlap based assembler might do a very good job, but where speed is concerned graph based assemblers are hard to beat.

Also, with illumina that now has high 3' quality long reads, maybe bigger kmer can be justified.

Thanks