[Pedro]

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

[ttnguyen]

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

[seb567]

Ray 1.3.0 and de novo assembly of Illumina CEO genome in 11.5 h

Dear all,

Ray 1.3.0 is now available online.

http://sourceforge.net/projects/denovoassembler/files/Ray-1.3.0.tar.bz2

The most important change is the correction of a major bug that caused
parallel infinite loop on the human genome.

This, along concepts incorporated in Ray 1.2.4, allowed Ray to assemble
the genome of Illumina's CEO in 11.5 hours using 512 compute cores (see
below for the link).

What's new?

1.3.0

2011-03-22

* Vertices with less than 1 of coverage are ignored during the
computation of seeds and during the computation of extensions.
* Computation of library outer distances relies on the virtual
communicator.
* Expiry positions are used to toss away reads that are out-of-range
* When only one choice is given during the extension and some reads
are in-range, then the sole choice is picked up.
* Fixed a bug for empty reads.
* A read is not added in the active set if it is marked on a
repeated vertex and its mate was not encountered yet.
* Grouped messages in the extension of seeds.
* Reads marked on repeated vertices are cached during the extension.
* Paths are cached in the computation of fusions.
* Fixed an infinite loop in the extension of seeds.
* When fetching read markers for a vertex, send a list of mates to
meet if the vertex is repeated in order to reduce the communication.
* Updated the Instruction Manual
* Added a version of the logo without text.


I fixed a bug that caused an infinite loop. Now Ray can assemble large
genomes. See my blog post for more detail about that.

http://dskernel.blogspot.com/2011/03/de-novo-assembly-of-illumina-ceo-genome.html

Version 1.2.4 of Ray incorporated also new concepts that I will present
at RECOMB-Seq 2011.

The talk is available online:

http://boisvert.info/dropbox/recomb-seq-2011-talk.pdf

Sébastien Boisvert

[seb567]

What is the current status with color space data ?

It is mostly not tested, I would say.

[lcollado]

Hello,

I just began reading this whole thread (and skipped some of the longer posts :P) and I would like to know what is the current status with CS data. In some of the older posts you mention a few problems Ray has with CS data. Does it work correctly now?

Thank you and greetings
Leo

[seb567]

Ray 1.2.3: new heuristics for large distances (jumping libraries)

Hello,

Ray 1.2.3 is now available.

http://denovoassembler.sourceforge.net/

One of the new features is the capability of assembling genomes with
larger distances.

See this post for an example:

http://seqanswers.com/forums/showpost.php?p=34108&postcount=76

Other changes:

1.2.3 'cRAYfish'
2011-02-05

• Instruction manual
• Fixed a bug leading to a segmentation fault when providing an invalid file or no file at all. This bug was introduced in Ray 1.2.2.
• Fixed a bug leading to a hang when providing invalid arguments.
• Added reporting of memory usage after loading sequences and after distributing vertices.
• Removed read simulators because samtools include one called 'wgsim'.
• Changed the paired-end heuristics to accomodate larger distances.
• Removed RepeatedVertexWatchdog.
• Changed the RepeatedThreshold from 255 to 2*peakCoverage (or 255 if it is higher than 255)
• Now both versions of read pairs are utilized.
• Paired reads are utilized if the distance falls in mean-3 standard deviation;mean+3 standard deviation. (3 instead of 1 as described in the paper).
• If a pair of reads does not respect the constraints, the right read is removed from the used read set.
• Fusions are computed more cautiously.
• Added reporting of memory usage for sequence reads and for vertices. And for a few other places too.
• Improved memory usage for the extensions. Now calling destructors instead of clear() because the latter does not free memory.
• Messages are grouped in the computation of seeds, using VirtualCommunicator.
• Added a Doxygen configuration file.
• The paths in the distributed graph are now numbered with an integer of 64 bits. This fixes a segmentation fault occuring with large genomes -- when identifiers overflow on 32 bits. (need to verify this)
• The partition on sequence reads is outputted.
• Instead of having the Rank 0 (master) to request the contigs of each every other rank, now each rank appends its fusions to the contigs file and also to the AMOS file.
• Fixed a bug in the numbering of reads in the AMOS file.
• Fixed a memory usage problem in the memory consumption reducer
• Improved the running time for the extension of seeds by removing an inner loop.

-seb

[corthay]

Hi seb567,

Thank you for the detailed answer.
I am looking forward to the new version
Also, it would make me happy if Ray can handle more large kmer
though it might not be important for Ray.

Thanks
Corthay

[seb567]

Hey,

I have uploaded the last version of Ray to can obtain an amos output format. It ran during ten hours but I obtain errors as:
MPI_Isend(145): MPI_Isend(buf=0x2b196b013ad0, count=1, MPI_UNSIGNED_LONG_LONG, dest=23612, tag=89, MPI_COMM_WORLD, request=0x7fff25927d24) failed
MPI_Isend(95).: Invalid rank has value 23612 but must be nonnegative and less than 16

Do you konw this problem (and have you a solution)?

Yes, there is bug in Ray 1.2.2 concerning the numbering of reads in the output of an AMOS file.

In Ray 1.2.2, I introduced the following feature:

Parallel reading of sequence files

But I forgot to modify the module that writes AMOS files.

In Ray, each sequence has a unique identifier, encoded in a unsigned integer of 64 bits.

Before Ray 1.2.2, the rank responsible for the sequence x was computed with x%NumberOfRanks. And the index of that sequence on the rank was computed with x/NumberOfRanks.

In Ray 1.2.2 and later, this changed.

Given a sequence numbered x, a number of ranks N and a number of sequences M, the rank of x is x/(M/N). If the computed rank exceeds the number of ranks minus 1, then it is decremented because the last rank will have a little more sequences.

The index of a sequence on a rank is computed by subtracting from x the total number of sequences stored by ranks having a lower rank number.

This issue is fixed in 1.2.3, which is set to be released in the next weeks, I guess.

Best regards

PS: on the same data with a default output format, it was ok...!

OK

[seb567]

Hi seb567,

I have two questions.

1) Can Ray handle sanger reads such as WGS, Fosmid end ?

Good question. Particularly relevant I would say !


I am currently modifying Ray heuristics to produce assemblies using paired sequences separated by very large physical distances.

I simulated 4 libraries from E. coli K-12 MG1655 genome sequences. These reads include 0.5% mismatch errors. ANd I used wgsim available in samtools.

Fragment standard deviation is 10% of the fragment average length.

These libraries are

200 +/- 20
1000 +/- 100
4000 +/- 400
10000 +/- 1000

/software/samtools-0.1.7/wgsim -e 0.005 -d 200 -s 20 -N 4000000 -1 30 -2 30 -r 0 -R 0 -X 0 Streptococcus-pneumoniae-R6.fasta 200-17_1.fastq 200-17_2.fastq 

/software/samtools-0.1.7/wgsim -e 0.005 -d 1000 -s 100 -N 4000000 -1 30 -2 30 -r 0 -R 0 -X 0 Ecoli-k12-mg1655.fasta 1000_1.fastq 1000_2.fastq 

/software/samtools-0.1.7/wgsim -e 0.005 -d 4000 -s 400 -N 4000000 -1 30 -2 30 -r 0 -R 0 -X 0 Ecoli-k12-mg1655.fasta 4000_1.fastq 4000_2.fastq 

/software/samtools-0.1.7/wgsim -e 0.005 -d 10000 -s 1000 -N 4000000 -1 30 -2 30 -r 0 -R 0 -X 0 Ecoli-k12-mg1655.fasta 10000_1.fastq 10000_2.fastq

Then, using Ray:

mpirun  -np 31 ~/Ray/trunk/code/Ray \
-p 200_1.fastq 200_2.fastq -p 1000_1.fastq 1000_2.fastq -p 4000_1.fastq 4000_2.fastq -p 10000_1.fastq 10000_2.fastq | tee log

Without scaffolding, Ray generates 65 contigs each having at least 500 nucleotides. There is a total of 4625792 assembled bases, all are A, T, C or G. The mean size is 71166, the N50 is 121756 and the largest is 371980.

This is without any scaffolding.

I then used MUMmer to validate contigs: 0 misassembly, 0 mismatch error and 0 small indel.

        %  & numberOfContigs & bases & meanSize  & n50  & max   & coverage   & misassembled & mismatches & indels
 Ray & 65 &  & 71166 & 121756 &  371980 &  0.9965 & 0 & 0 & 0 \\

99.65% of the genome is covered by Ray contigs.

Running time/31 processors:

[1,0]<stdout>: Beginning of computation: 2 seconds
[1,0]<stdout>: Distribution of sequence reads: 52 seconds
[1,0]<stdout>: Distribution of vertices & edges: 2 minutes, 56 seconds
[1,0]<stdout>: Calculation of coverage distribution: 0 seconds
[1,0]<stdout>: Indexing of sequence reads: 15 seconds
[1,0]<stdout>: Computation of seeds: 27 seconds
[1,0]<stdout>: Computation of library sizes: 4 minutes, 45 seconds
[1,0]<stdout>: Extension of seeds: 4 minutes, 15 seconds
[1,0]<stdout>: Computation of fusions: 37 seconds
[1,0]<stdout>: Collection of fusions: 0 seconds
[1,0]<stdout>: Completion of the assembly: 14 minutes, 9 seconds

The contig lengths:

388
641
1038
2164
2714
3438
4735
5845
7199
7939
9557
10364
11883
16323
16937
20213
24049
25788
26056
26433
27014
29449
30205
31159
33620
34875
36084
37445
38010
38550
39509
43086
43623
43841
46455
47970
58692
59083
59441
59889
64822
66096
66767
66783
67838
70713
82799
85321
94029
95691
97407
100961
102906
103167
121756
121858
138099
156439
161412
163121
169912
177402
213535
284796
318866
371980


Let me answer your question now. For seed extension, an algorithmic process that generates contigs, available versions of Ray will not accomodate very large fragment lengths, 400 at most I think.

It is because for larger distances, the right sequence in a pair most likely starts on a repeated k-mer. As such, constraints on the observed distances most be further enforced -- paired reads not satisfying the constraints must be set free in order to allow them to be used in the adequate time.

Ray 1.2.3 'cRAYfish' is set to be released soon. For this release, I want to give a running time on 512 processors for a human genome. But, for the time being, I am currently waiting for some compute time on a compute cluster.

1.2.3 fully utilizes large distances provided by jumping libraries (fosmid ends).

I also plan to incorporate a scaffolder, but that is not my priority now.

My priority is a memory utilization reducer that is utilized during the process of building the distributed de Bruijn graph.

2) Does Ray do scaffolding ? The results did not contain N even though
I used paired end as following.

Not yet, but it is a planned feature.

Note however that Ray make full use of paired reads to traverse repeated elements in the genome. And apparently it does that very well.


Meanwhile, I suggest you use SSPACE.

http://seqanswers.com/forums/showthread.php?t=8350

Bioinformatics paper http://bioinformatics.oxfordjournals...s.btq683.short

The header line of fasta file is modified by SOAPdenovo correction tool.
Is this problem ?

Thanks,
Corthay

no

Ray does not read the header in the fasta and fastq files containing sequences.

[corthay]

Sanger reads

Hi seb567,

I have two questions.

1) Can Ray handle sanger reads such as WGS, Fosmid end ?

2) Does Ray do scaffolding ? The results did not contain N even though
I used paired end as following.

Ray \
-s \
solexa/xxx.2.2.41.single.fasta \
-p \
solexa/xxx.2.2.41.1.fasta \
solexa/xxx.2.2.41.2.fasta \
-o \
ray/ray.kmer31.contig \
-k \
31

The header line of fasta file is modified by SOAPdenovo correction tool.
Is this problem ?

Thanks,
Corthay

[gstitan]

bug in Ray1.2.1?

Hey,

I have uploaded the last version of Ray to can obtain an amos output format. It ran during ten hours but I obtain errors as:
MPI_Isend(145): MPI_Isend(buf=0x2b196b013ad0, count=1, MPI_UNSIGNED_LONG_LONG, dest=23612, tag=89, MPI_COMM_WORLD, request=0x7fff25927d24) failed
MPI_Isend(95).: Invalid rank has value 23612 but must be nonnegative and less than 16

Do you konw this problem (and have you a solution)?

Best regards

PS: on the same data with a default output format, it was ok...!

[seb567]

Sounds interesting !

[talioto]

dnGASP

Seb,

If you haven't seen the posting already, I'd like to bring to your attention the de novo Genome Assembly Assessment Project (dnGASP). All details can be found at cnag.bsc.es, but briefly it's a project that solicits submissions of assemblies of a synthetic 1.8Gb diploid genome, which will be followed up with a workshop in April in Barcelona. Please take a look and see if you may be interested. We at the CNAG are looking at Ray as an option for assembling human genomes and other genomes that we have already sequenced here (but have been unable to assemble thus far). The project and associated workshop may be a good opportunity to gain exposure and to participate in a forum where issues regarding assembly and sequence data can be discussed. If you are interested, please sign up for the dngasp mailing list and register your team (Ray) on the cnag.bsc.es site as soon as possible in order to download the reads and submit assemblies.

Regards,
Tyler Alioto

[seb567]

Ray 1.2.1 'stringray'

Dear all,

Ray 1.2.1 is now available.

Source:

http://sourceforge.net/projects/denovoassembler/files/Ray-1.2.1.tar.bz2

This version fixes 2 critical flaws, can assemble polymorphic positions by
forcing bubble traversal, and adds an experimental feature: memory usage
reduction during the construction of the distributed graph.

A more detailed catalog of changes follows.

• SplayTreeIterator now iterates in preorder instead of inorder.
  
• Now the SplayTreeIterator is utilized to iterate over the vertices nstead ofstoring them in an array (takes too much memory).
  
• The forest of splay trees now contains 16384 trees instead of 4096. otethat as usual, each MPI rank has its own forest. Furthermore, the forest freezes when vertices distribution is properly finished -- which means no more splaying in the splay trees is to occur. The process ensures that vertices with low redundancy remain at leaves.
  
• Bubble traversal ensures no misassemblies as polymorphic positions (substitution, and indels) are assembled !!!
• Works with 454 data too (454 homopolymers error are interpreted as polymorphic positions).
• Added numeric indicator for files. (example: [1/9])
• Corrected a bug in library lengths messaging that leaded to hanging and/or Bus Error.
• Fixed a segmentation fault when the -a (output AMOS) is provided. Thanks to Daniel Brami from J. Craig Venter Institute (La Jolla, CA) for the timely report.
• Preliminary version of an algorithm to preemptively reduce the memory usage while building the distributed graph.
• Under GNU/Linux platforms, Ray outputs the virtual memory (VmData from /proc, that is the heap) utilized before exiting.

Thank you.

ps.

I am currently working on two human genome data sets: SRA000271
(African) and SRA010766 (Jay Flatley).

[seb567]

Dear caddymob,

I did some investigative work for you this morning.



You have 140 174 250 paired reads, each having a length of 105.

Would it be better to split these up into smaller chunks or does that matter?

Splitting the files would change nothing since Ray loads sequences in a lazy manner.

http://en.wikipedia.org/wiki/Lazy_loading

You use 256 MPI ranks, which are mapped on 8-core nodes, each having 16 GiB of memory.
You use 32 such nodes.

You have 512 GiB of distributed physical memory.



In modern operating systems, a program's addresses are virtual. If the computer has 2 GiB of random-access memory (RAM) and 100 GiB of swap memory,
then a program can allocate more than 2 GiB for its use. However, if the heap size (utilized memory) exceeds the random-access memory size, page faults occur.

"A page fault is a trap to the software raised by the hardware when a program accesses a page that is mapped in the virtual address space, but not loaded in physical memory."

http://en.wikipedia.org/wiki/Page_fault

It is followed by paging, 'one of the memory-management schemes by which a computer can store and retrieve data from secondary storage for use in main memory.'

http://en.wikipedia.org/wiki/Paging

Basically, if you exceed the physical memory, the job will just take forever because a lot of instructions will be only to swap pages between physical memory and swap memory.

Agreed, but this is something Ray did, I do not, to my knowledge have a way to specify how much swap to use. What can I do to prevent this?

Ray, presumably, requested memory, and page faults occured because physical memory was exausted.

You can't specify not to use swap because the programmer, and the program executing are unaware of which addresses are resident in physical memory and which are not.


Luckily, it seems you can limit the physical and virtual memory available to any process using a properly configured PBS scheduler: http://wiki.hpc.ufl.edu/index.php/PBS_Directives#Memory



You use PBS to launch the job.

According to the post numbered 69, your script is the following:

http://seqanswers.com/forums/showpost.php?p=31722&postcount=69

###Parameterized PBS Script ####
#PBS -S /bin/bash
#PBS -N LUNDE.Ray4
#PBS -l nodes=256
#PBS -l walltime=25:00:00
#PBS -q normal
#PBS -j oe
#PBS -o LUNDE.Ray4.o
#PBS -M ---redacted---
#PBS -m abe

wd=/scratch/myfiles/LUNDE_ASSEMBLE/
cd $wd

use intel-openmpi-1.4.2
use Ray-1.2.0

mpirun -np 256 Ray -p $wd\Lunde_1.fastq $wd\Lunde_2.fastq -o Lunde-contigs

To properly understand the parameters of the PBS scheduler, I have read http://wiki.hpc.ufl.edu/index.php/PBS_Directives
Your setup seems OK to me although I never worked with PBS.



In the post numbered 61 http://seqanswers.com/forums/showpos...8&postcount=61,

you wrote this:

Req[0] TaskCount: 256 Partition: anon
Utilized Resources Per Task: PROCS: 120.18 MEM: 2596M SWAP: 881G
Avg Util Resources Per Task: PROCS: 120.18
Max Util Resources Per Task: PROCS: 237.11 MEM: 2596M SWAP: 881G
Average Utilized Memory: 1641.30 MB
Average Utilized Procs: 48473.59
NodeSet=ONEOF:FEATURE:awesometown
NodeAccess: SINGLEJOB
NodeCount: 32

As I understand it, you have 256 tasks, and a task utilizes 120 processors, 2596 MiB of physical memory, and 881 GiB of swap space.

For all that matters, the node set is named 'awesometown'.

That does not make sense.

Too, 'Average Utilized Procs: 48473.59' does not make sense neither.



Your PBS script might be wrongly written or PBS might be misconfigured.

To ascertain the wholly truth, get support from your compute department.



As I promised in the post numbered 66 http://seqanswers.com/forums/showpos...4&postcount=66, I tested Ray with the Intel compiler.

I compiled Ray 1.2.1 with Open-MPI 1.4.3 and the Intel compiler, version 11.1.059.

I then launched a job, using mpirun (from Open-MPI 1.4.3) compiled with the Intel compiler. The script for Sun Grid Engine follows.

#!/bin/bash
#$ -N Ray1.2.1
#$ -P nne-790-aa
#$ -l h_rt=0:20:00
#$ -pe node 64
#$ -M sebastien.boisvert.3@<removed>
#$ -R y
#$ -m bea
module load compilers/intel/11.1.059 mpi/openmpi/1.4.3_intel
/software/MPI/openmpi-1.4.3_gcc/bin/mpirun /home/sboisver12/Ray/trunk/code/Ray  \
-p /home/sboisver12/nne-790-aa/SRA001125/SRR001665_1.fastq /home/sboisver12/nne-790-aa/SRA001125/SRR001665_2.fastq \
-p /home/sboisver12/nne-790-aa/SRA001125/SRR001666_1.fastq /home/sboisver12/nne-790-aa/SRA001125/SRR001666_2.fastq \
-o  Intel

The output in stdout is http://pastebin.com/yprcXwDe

Each MPI rank needed, on average, 284 MiB. The utilized distributed memory was 17 GiB.

I then used MUMmer (via the script named print-latex.sh in the scripts directory of the Ray distribution) to assess the quality of the assembly.

[sboisver12@colosse1 ~]$ print-latex.sh nne-790-aa/nuccore/Ecoli-k12-mg1655.fasta Intel.fasta  Ray-1.2.1
        %  & numberOfContigs & bases & meanSize  & n50  & max   & coverage   & misassembled & mismatches & indels
 Ray-1.2.1 & 123 & 4616336 & 37531 & 72499 &  176360 &  0.9819 & 0 & 2 & 4 \\

Everything's OK with the Intel compiler.


Cheers.

-seb