ROC curve analysis of short read aligners

[Mchicken]

Hey guys,

i want to evaluate the performance of several aligners by applying them on a synthetic generated set of reads. Sadly i have not found a satisfying explanation for performing an analysis (including ROC curve analysis) on this kind of data and I am completely new to the ROC curve analysis and a little bit confused.

Here are some ideas:

Step1:

Generating 100.000 (or more) reads from my reference genome by usage of Mason or randomreads.sh (see BBTools) or an own script with illumina specific quality values and some sequencing error rate. What i then get are 100.000 reads of which i know the correct location in my reference genome.

Step2:

Run my set of aligners on this synthetic dataset. What i then have is for each aligner a sam-file containing the alignment locations of each read.


Step 3:

For my analysis i now have to define my set of TP,FP etc as well as terms like sensitivity, accuracy etc:

TP = correct aligned reads (if the alignment start is within 50bp (cause of clipping) of the correct position).

FP = incorrect aligned read (read is aligned at the wrong position)

FN = read was not aligned ?

TN = ???


Sensitivity = Number of correct aligned reads divided by the number of all reads (100.000) = TP/(TP+FP+FN)

Precision = Number of correct aligned reads divided by the number of all aligned reads = TP/(TP+FP)

FDR = Number of all incorrect aligned reads divided by the number of all aligned reads = FP (TP+FP)

Has someone other ideas?

Now i have the problem what should i plot? And how do i get my datapoints?

In the following lower figure http://www.ncbi.nlm.nih.gov/pmc/arti...-0065632-g002/ they plot FP vs TP and get the datapoints by calculating the cumulative number of correctly (TP) and incorrectly (FP) aligned reads from high to low mapping quality.

Here i have the problem that the aligners i use have varying mapping qualities for exactly the same alignments e.g. SSAHA2 has a value of 100 and BWA a value of 60 and bowtie2 a value of 255 in addition the number of different mapping quality values differs between the aligners leading to different numbers of datapoints. Don't the varying mapping qualities bring bias in my plot?

The mapping qualities achivied by bowtie2,BWA,SSAHA2 and NextGenMap and their frequencies in my set of 100.000 reads is provided here: mapping_qualities.pdf
As can be seen about 66.000 reads mapped by bowtie2 have a mapQ of 255, which should normally indicate (as can be seen in the SAM format specification) that there is no mapQ available. Nevertheless by looking at the same alignments in the sam files of the other aligners i see values of 60 for BWA and values about 100 for SSAHA2, such that i would say that a mapQ of 255 for bowtie2 tells me a good mapping quality.


What i would now do is that i make a figure with different cumulative bins for the mapping quality (0,>=10,>=20...,>=60) and use these bins on the y-axis and on the x-axis the sensitivity as given above (this is pretty similar to the strategy in the posted figure). Would this make sense for you?

Another figure might be that I again use the quality bins on the y-axis and on the x-axis the ratio of TP to FP. Does this also make sense for you?

Does anyone of you has other ideas for nice illustration of the mapping performance of aligners?

Thanks in advance for your help!

Mchicken

[Brian Bushnell]

One of the problems with ROC curves in mapping is that there are no 'true negatives' with normal data. You could throw in some completely random sequences, but those simply won't map. Also, ROC curves were designed to be for binary decisions: "Does the patient have the disease, or not?"

As a result, statisticians and electrical engineers may complain about details that seem very nit-picky, and some of the terms may not correlate perfectly to the same terms in binary-decision space. Also, in binary-decision space, a perfectly random decider would generate straight line of X=Y; while a perfectly random mapper would generate a cluster of points at X=100, Y=0.

I recommend that you do NOT do binning for the ROC curve - rather, each point should include data for all the reads with a mapping quality of up to some threshold. For example, if an aligner only produced MAPQ values of 2, 5, 9, and 255, the ROC curve should have 4 points:

1) (TP, FN) for all reads of MAPQ={2}
2) (TP, FN) for all reads of MAPQ={2, 5}
3) (TP, FN) for all reads of MAPQ={2, 5, 9}
4) (TP, FN) for all reads of MAPQ={2, 5, 9, 255}

As for Bowtie2, I don't know what the problem is; last time I tested it, the mapping qualities didn't look like that - there were no 255s. Edit: That was v2.1.0, so perhaps you should download an earlier version.

[Mchicken]

Probably I did not express myself clearly. With binning i ment a threshold like the one used in the following publications:

http://bioinformatics.oxfordjournals...7/20/2790.full - Figure 2a where it is called a quality threshold

or

http://www.plosone.org/article/info%...l.pone.0090581 - Figure1 where it is called a quality cutoff

[Brian Bushnell]

Oh, right, I misread you - I thought by binning you meant exclusive bins of 0-10, 11-20, etc.

[monia]

what's the commande line used for ROC construction?

[Brian Bushnell]

With BBTools you can generate ROC curves from any aligner like this:

(index the reference)
bbmap.sh ref=reference.fasta

(generate random reads)
randomreads.sh reads=100000 length=100 out=synth.fastq maxq=35 midq=25 minq=15

(map to produce a sam file)
...substitute this command with the appropriate one from your aligner of choice
bbmap.sh in=synth.fq out=mapped.sam

(generate ROC curve)
samtoroc.sh in=mapped.sam reads=100000


For randomreads.sh, run the script without any parameters for usage information. It allows options to add various different mutation types and rates in the reads.

[monia]

thank you it's just what i need.