Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • #31
    Each line consists of mapping quality threshold, #mapped reads with mapQ no less than the 1st column and #mismapped reads. It does not show #reads with mapQ=0. If we include mapQ=0 mappings, the sensitivity of bwa is also good for simulated data, but on single-end real data, the low-quality tail on reads makes bwa much worse. This is what Steven and Ben have observed. This is also why it is recommended to enable trimming when using bwa.

    BWA always gives mapQ 0 to repetitive hits, but other mappers (gsnap, bowtie2 and novoalign) may give mapQ<=3 to repetitive hits. This is theoretically correct. I may further set a mapQ threshold 1-4 when plotting.

    Comment


    • #32
      Originally posted by cjp View Post
      @nickloman

      The output of the wgsim_eval.pl program looks a bit like the data below - bowtie 1 always gives a mapping score of 255 (column1). I'm guessing that bowtie 2 has many FP's at a mapping score of 1 (column3 if column1 == 1), but cumulatively finds more TP's with all mapping scores (column2 if column1 == 1). But I was also wondering the exact meaning from the output of the wgsim_eval.pl script.

      % tail *.roc
      ==> bowtie2.roc <==
      14 172922 11
      13 172925 12
      12 177943 27
      11 177945 28
      10 179990 37
      9 179995 40
      4 180250 40
      3 187273 578
      2 187324 580
      1 199331 5877

      ==> bowtie.roc <==
      255 86206 1740

      ==> bwa.roc <==
      10 192354 72
      9 192560 107
      8 192595 107
      7 192628 110
      6 192652 115
      5 192669 116
      4 192681 117
      3 192731 117
      2 192741 118
      1 192762 119

      Chris
      I forked Heng's code a while back into the dwgsim project (links below). I also added user documentation:

      Download dnaa for free. DNAA is the DNA analysis package, for analyzing next-generation post-alignment whole genome resequencing data. Specifically, DNAA is able to find structural variation, SNP and indel variants, as well as evaluating the mapping and data quality.

      Comment


      • #33
        Originally posted by rskr View Post
        I disagree. If you look at hash based aligners there are certain patterns of indels, mismatches and errors, where they won't find the right result even if it is unique. For example if the word size is 15, and there are are two mismatches 10 bases apart in a 50mer, the hash won't return the region at all. Likewise for longer reads the number of mismatches is likely to be higher and the Suffix Array search will terminate before finding the ideal match.
        That's a ridiculous statement! Most hashed aligners using a 15-mer hash would need 3 equally spaced mismatches in a 50mer to miss an alignment. But there are some hashed based aligners that are can find a 50-mer alignment with 5 or 6 mismatches even if equally spaced. Novoalign can do this.

        Comment


        • #34
          Looking at those ROC curves, it appears to me that Novoalign is the best mapper in the specified simulation that was run with respect to sensitivity and specificity. Is this a correct interpretation?

          Comment


          • #35
            Originally posted by Heisman View Post
            Looking at those ROC curves, it appears to me that Novoalign is the best mapper in the specified simulation that was run with respect to sensitivity and specificity. Is this a correct interpretation?
            From this comparison yes that would be the case, with smalt showing good performance as well.What would be interesting would be to do this with genericforms' suggestion of a 30x coverage genome.

            Comment


            • #36
              Originally posted by Heisman View Post
              Looking at those ROC curves, it appears to me that Novoalign is the best mapper in the specified simulation that was run with respect to sensitivity and specificity. Is this a correct interpretation?
              For sure, we believe it's still the best

              Colin

              PS. Some mapper comparisons have shown different results but this can be the result of targeting the simulation at the mapper the developer is promoting. One recent example used high simulated indel rates and then didn't adjust the other mappers gap open & extend penalties to suit. Their mapper, with default low gap penalties, came out as the clear winner.
              It can be a problem to optimise the parameters for every mapper when doing these comparisons and the tendency is to use defaults which is probably reasonable.

              I think Heng Li is doing an honest and unbiased comparison at mutation rates you'd expect in a resequencing project.

              Comment


              • #37
                Heng, thanks for making the ROC plots available. I think they're pretty interesting.

                About GSNAP having only a single point on the plots: actually, GSNAP does calculate mapping quality, but my understanding was that the quality should translate into the probability that the given read alignment is the correct one. So GSNAP does a Bayesian normalization of all of the raw mapping qualities, and then reports the normalized mapping quality. This tends to produce a dichotomous result, where if one alignment is much better than the others, it gets a very high mapping quality of 40 or so, but if there are two or more roughly similar multimappers, they all get a mapping quality of 3 or so (where 3 corresponds to the log probability of 0.5). Perhaps I have the wrong understanding of mapping quality here (and maybe someone can correct me), but I am told that GATK has a similar expectation.

                To get around this, I have added a new field in the SAM output of GSNAP called XQ, which gives the non-normalized mapping quality, although it is still scaled so the best alignment gets a mapping quality of 40. (There are certain reasons why the scaling is important, mostly having to do with GSTRUCT needing to know this information.)

                Regarding the comment someone made about multimappers, GSNAP is designed to report all multimappers, but I think it has a different meaning than other programs of what a multimapper is. Some programs expect to give the program a specific search parameter, like 5 mismatches or less, and then expect to get all mappings that satisfy that parameter. On the other hand, GSNAP interprets the search parameter to be the maximum perimeter of the search space, so has a hard limitation of looking for 6 mismatches or more. However, if GSNAP finds an exact match, it will also place a soft limitation at that point and just report all multimappers that are also exact matches, and not report the 1-mismatch through 5-mismatch answers. The exception is that if GSNAP is given a value for suboptimal mismatches, and then it will proceed to search at that many mismatches past its optimal finding. For example, if GSNAP cannot find an exact match, but finds a 1-mismatch, and is given a suboptimal-mismatch value of 1, then it will report all 1-mismatch and 2-mismatch alignments, but will still not go past that to report 3-mismatch through 5-mismatch answers.
                Last edited by twu; 11-18-2011, 02:41 PM.

                Comment


                • #38
                  Heng, in cases where the distribution of positive vs negative examples is very skewed (such as variant calling), the ROC curve can also be misleading. The ROC curve typically only looks at positive examples (false positive rate vs true positive rate), but one should also look at the corresponding ROC curve for negative examples (TNR vs FNR), or, look at precision-recall curves.

                  Comment


                  • #39
                    @kopi-o As Steve has argued, the experiment I am doing now is not a classical binary classification. There are only "wrongly mapped", but actually no "false positives" in the strict sense.

                    I guess what you mean is to look at reads generated from regions with simulated polymorphisms. This is also Gerton insists doing and I agree it is a good thing to do. The current ROC does not tell you if a mapper uniformly misses a hit or consistently misses a hit in some regions. Similarly, the ROC does not tell you whether a mapper tends to produce consistent errors or random errors. All these matter in variant calling.

                    If you are interested in variant calling, the right evaluation is to plot the ROC for variant calls. This is a standard binary classification and is more telling.

                    Comment


                    • #40
                      @Thomas I quite like the GSNAP algorithm as well as its implementation and I have recommended it to others already. It is one of the top NGS mappers nowadays. My opinion is for GSNAP the only thing might be improved is a more useful mapping quality. I know GSNAP gives mapQ, but for "unique" hits, the vast majority get a mapQ 40 and the few hits with higher mapQ are actually not more accurate. Perhaps having higher mapQ may not help the standard SNP calling too much, but there are areas where extremely high mapping accurate is preferred.

                      I am not sure how to improve mapQ for single-end mapping, but I kinda think you should be able to derive better mapQ for paired-end mapping. It seems to me that GSNAP will visit more suboptimal hits in the PE mode. By seeing more hits and using the pairing information, you can know some hits can barely wrong.

                      Comment


                      • #41
                        BTW, here is an interesting observation on speed. I simulated 100k single-end (SE) reads and 100k pairs of paired-end (PE) reads (200k reads). One would think a program should run about twice as slow in the PE mode simply because there are twice as many reads. This is true for bowtie2/bwa/bwa-sw. Nonetheless, gsnap is much slower in the PE mode. My guess is in the PE mode, gsnap visits more suboptimal hits to get more reads paired. It is slower, but the accuracy is also higher. On the other hand, both novoalign and smalt are faster in the PE mode, but at a cost the false positive rate also goes slightly higher. My explanation is that they do not map each end separately and then pair them (bwa/bwa-sw does this), but rather map the pair as a whole.

                        Comment


                        • #42
                          Heng, I have a quick question. In trying to use simulated data to recall mutations, would you base recalibrate BWA-mapped reads? With real human data you can use dbSNP, however with simulated data what would you use?

                          Comment


                          • #43
                            BWA does not use base quality during alignment except for trimming. I have not been convinced that the difference between using base quality or not has a significant effect on downstream data analyses.

                            Comment


                            • #44
                              I have been wondering the same thing. So if I was to compare the recall of mutations from BWA mapped reads to a mapper that does recalibrate base qualities, you do think it would matter if I use GATK to first recall rate the reads mapped by BWA? You I think did this previously in a paper with Nils, right? How did you do the comparison?

                              Sorry for all the questions!

                              Comment


                              • #45
                                Updated to bowtie2-beta4. On accuracy, bowtie2-beta4 is similar to bwa-sw overall. I have also done the comparison on real data following the way I used in the bwa-sw paper. Out of 138k 454 reads with average read length 355bp, bwa-sw misses 1094+58 good alignments (~90% shorter than 100bp) and gives 31 questionable alignments, while bowtie2-beta4 misses 13+91 good alignments and gives 65 questionable alignments. The accuracy is largely indistinguishable for practical applications. On speed, Bowtie2 is about 20% faster and uses less memory.

                                In conclusion, bowtie2-beta4 has similar accuracy to bwa-sw for both 100bp simulated data and 350bp real 454 data. It is one of the best (accuracy+speed) mappers for hiseq and 454 reads. I will start to recommend it to others along with smalt/novoalign/gsnap. I think a missing feature in bowtie2 is to properly report chimeric alignments, which is essential to mapping even longer sequences. This should be fairly easy to implement.
                                Last edited by lh3; 12-06-2011, 12:24 PM. Reason: typo

                                Comment

                                Latest Articles

                                Collapse

                                • seqadmin
                                  Techniques and Challenges in Conservation Genomics
                                  by seqadmin



                                  The field of conservation genomics centers on applying genomics technologies in support of conservation efforts and the preservation of biodiversity. This article features interviews with two researchers who showcase their innovative work and highlight the current state and future of conservation genomics.

                                  Avian Conservation
                                  Matthew DeSaix, a recent doctoral graduate from Kristen Ruegg’s lab at The University of Colorado, shared that most of his research...
                                  03-08-2024, 10:41 AM
                                • seqadmin
                                  The Impact of AI in Genomic Medicine
                                  by seqadmin



                                  Artificial intelligence (AI) has evolved from a futuristic vision to a mainstream technology, highlighted by the introduction of tools like OpenAI's ChatGPT and Google's Gemini. In recent years, AI has become increasingly integrated into the field of genomics. This integration has enabled new scientific discoveries while simultaneously raising important ethical questions1. Interviews with two researchers at the center of this intersection provide insightful perspectives into...
                                  02-26-2024, 02:07 PM

                                ad_right_rmr

                                Collapse

                                News

                                Collapse

                                Topics Statistics Last Post
                                Started by seqadmin, 03-14-2024, 06:13 AM
                                0 responses
                                34 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 03-08-2024, 08:03 AM
                                0 responses
                                72 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 03-07-2024, 08:13 AM
                                0 responses
                                81 views
                                0 likes
                                Last Post seqadmin  
                                Started by seqadmin, 03-06-2024, 09:51 AM
                                0 responses
                                68 views
                                0 likes
                                Last Post seqadmin  
                                Working...
                                X