[ShaunMahony]

Thanks for the info Ih3. I agree that it is very useful to report the locations of repetitive / non-uniquely mappable reads. However, can MAQ be set to report ALL of the repetitive locations rather than just a single random one? I know that this is off-topic; apologies.

[lh3]

Q0;Q10;Q20 is threshold on the alignment quality score assigned by MAQ.

MAQ is initially designed for resequencing and keeping these repetitive reads is quite useful for the subsequent SNP calling with MAQ. This also helps CNV calling. I could understand that a lot of people do not want to see all these repetitive reads, but putting a threshold on mapping quality is very easy anyway. In addition, different people may want to set different different threshold.

As for the calculation of mapping quality, it just follows a very simple Bayesian procedure. You can calculate p(z|x,u) of read z mapped to u on the reference x. With Bayesian formula you get p(u|x,z). The mapping quality is -10log10(1-p(u|x,z)).

[ECO]

I'm definitely out of my league in this discussion, but if anyone needs hosting for some of these sample datasets, let me know!

[ShaunMahony]

I'm also interested in how MAQ assigns quality scores. Can you confirm what you meant by "Q0;Q10;Q20" in the MAQ tests? Is this a threshold on the quality score that MAQ gives the alignment (as opposed to the quality score of the read in the Fastq file)? If a read maps to multiple locations, MAQ reports one location at random and assigns a quality score of 0. Therefore the Q0 accuracy should be much less than if you had excluded these alignments. I think that this behavior is a bit strange; it would be less confusing if MAQ didn't report any matches for the non-uniquely mapping reads and instead reported the number of places that the read maps (the whole read, not just the first 25 bases).

[zee]

How does maq calculate mapping quality? I looked in the documentation and couldn't find anything. Is it based on the ASCII quality code in the query sequence?

[zee]

Thanks, I just got the maq package and I'd like to simulate data according to the same rules you used.
Is it correct to assume that in cases where I use eland that I won't be able to make use of the fastq format and instead I need to convert to fasta using a script?
The simulation steps seem quite straightforward.

[lh3]

The simulated data are free to use, of course, but unfortunately I could not find a FTP site to upload the data. You may try maq to simulate the reads by yourself at the moment. The read positions are coded in read names and so you can write your own script to evaluate the accuracy.

[zee]

Dataset

Nice work Ih3,

I think this is an invaluable exercise amidst all the emerging tools. Is the benchmark dataset available anywhere because I'd like to test similar metrics on tools that we have in our lab?

Thanks for your help.

[brudno]

Hi --

Any interest on trying this with SHRiMP? I am guessing it'll be much slower than the other tools, though we are working to improve this. While the main use for shrimp is probably going to be SOLiD, we do support solexa as well.

Cheers, -Mike

[Chipper]

I was not aware that there is more than one clock to watch... Anyway, the main advantage for SOAP over Eland is that you can actually download and use it, plus the ability to use longer reads, or trim them if no match is found. Which happens quite frequently...

[lh3]

I just tried SOAP with and without -p 4 on a quadcore, it was about 3 times faster with all four cores active (2.5 M reads, chr1).

You must be couting the wall-clock time or the time soap was printing out. Although wall-clock time should be considered, it is more important to evaluate on processor time. If you split your read file into 4 chunks and align them with 4 eland, I am sure eland will be still several times faster.

Here is the follow-up. When '-p 4 -g 0 -c 0 -r 0' is applied:

real 6m12.443s
user 23m37.292s (=1417.292 sec)
sys 0m3.799s

As we would expect, it takes more user time (vs. 1228.12 sec). Anyway, multithreading is still useful if you do not have many reads and want to get the result quickly. This is definitely the gain of multithreading. SOAP is also a great software. Thank you for pointing this out.

[Chipper]

I just tried SOAP with and without -p 4 on a quadcore, it was about 3 times faster with all four cores active (2.5 M reads, chr1).

[lh3]

Eland supports multithreading in its source codes, but apparently it is not activated by default (if I am right). In fact, as Eland is fast and small, we can invoke several eland on a node at the same time. Parallelization on a large cluster can be easily managed by LSF or SGE.

On a node, SOAP can be parallelized with -p (the CPU time should be similar with or without -p) and on a cluster you can also use LSF/SGE. However, as soap may require huge memory, you will have to design a clever strategy, based on the hardware configurations, to run it efficiently without breaking your clusters. In addition, LSF/SGE usually prefers single-thread jobs. Multi-thread jobs may reduce the overall efficiency of a cluster unless all the nodes in the cluster has dozens of processors.

In my view, indexing the reference helps to get faster speed but tends to be memory demanding for human alignment. Indexing reads is more scalable but may sacrify some speed. Eland is still faster firstly because Tony Cox at Illumina is one of the best programmers in this field and secondly because soap has to trade speed for memory. Using long seed for human alignment is impractical.

Anyway, all the software in this benchmark are great. Sometimes which to use is just subjected to your appetite.

[Malarky67]

I tried SOAP some more and it seems to adjust seed size to available memory, which is likely to affect the run time. As I understand it Eland places reads not reference in memory.

Yes. That is what I have heard. Does anyone understand how these algorithms are parallelised across multiple processors or even nodes of a cluster? (especially in the case of nodes are reference tables built for each node?)

[Chipper]

Hi
SOAP does seem slow. You haven't set it to work on more than one processor. Was this simulation done on a single processor machine?

e.g. from SOAP usage

-p <int> number of processors to use, default=1


Does Eland automatically have the same behaviour?

I tried SOAP some more and it seems to adjust seed size to available memory, which is likely to affect the run time. As I understand it Eland places reads not reference in memory.