Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • Compare mapped reads from different aligners

    I want to count the number of (single-end) reads that were mapped to approximately the same coordinates by different aligners.
    The problem is that the reads do not have identical IDs and may have shifted coordinates in a range of 1 bp (SOLiD mapped with BWA), for example:

    BWA:
    prefix_3_30_738 0 chr8 11162354 37 48M * 0 0 ...

    ABI BioScope:
    3_30_738 0 chr8 11162353 100 50M * 0 0 ...

    NovoalignCS:
    3_30_738_F3 0 chr8 11162353 150 50M * 0 0 ...

    Reads are in sorted, indexed BAM files. Of course I could change the read IDs and coordinates to find exact matches with Picard CompareSAMS, but I'd like to avoid redundance,
    reduce computational time and also output the matching reads. Besides, I'm interested in finding reads that may be aligned in a certain neighborhood.
    Has anyone already developed a tool that can handle such an issue? If not, what would be the most efficient strategy?

    Thank you for advice in advance!

    Barbara

  • #2
    Originally posted by epigen View Post
    I want to count the number of (single-end) reads that were mapped to approximately the same coordinates by different aligners.
    The problem is that the reads do not have identical IDs and may have shifted coordinates in a range of 1 bp (SOLiD mapped with BWA), for example:

    BWA:
    prefix_3_30_738 0 chr8 11162354 37 48M * 0 0 ...

    ABI BioScope:
    3_30_738 0 chr8 11162353 100 50M * 0 0 ...

    NovoalignCS:
    3_30_738_F3 0 chr8 11162353 150 50M * 0 0 ...

    Reads are in sorted, indexed BAM files. Of course I could change the read IDs and coordinates to find exact matches with Picard CompareSAMS, but I'd like to avoid redundance,
    reduce computational time and also output the matching reads. Besides, I'm interested in finding reads that may be aligned in a certain neighborhood.
    Has anyone already developed a tool that can handle such an issue? If not, what would be the most efficient strategy?

    Thank you for advice in advance!

    Barbara
    Yes, I developed the "dwgsim" toolset in DNAA (http://dnaa.sf.net). The "dwgsim" tool will create simulated reads, the "dwgsim_eval" function will give mapping sensitivity/accuracy statistics, and the "dwgsim_pileup_eval.pl" will give the sensitivity/accuracy of variant calling after samtools. Let me know if this works.

    Comment


    • #3
      storing read information: hash versus tree

      Thanks for the answer, Nils. I tried dwgsim_eval but it was not comfortable with the reads being not created by dwgsim and being single end:

      In function "process_bam": Warning[OutOfRange]. Variable/Value: 614 1806 266.
      Message: [dwgsim_eval] read was not generated by dwgsim?.
      ***** Warning *****
      ************************************************************
      ************************************************************
      In function "run": Fatal Error[OutOfRange]. Message: Found a read that was not paired.
      ***** Exiting due to errors *****
      ************************************************************

      By the way, some documentation on the DNAA tools would be really helpful ... As it seems, dwgsim_eval only makes statistics for one BAM or SAM file. What I want to do is to explicitely count how many reads were mapped to a similar coordinate by two different tools, and output these reads.
      For this purpose, I would obviously have to store the read information of one file. In a thread about duplicate removal someone recommended storing it in a trie structure instead of using one of these RAM-greedy Perl hashes. As far as I remember from my computer science lectures, a suffix tree (or eqivalent, a prefix tree) does not need less space than a hash. Searching in a tree is O(n) whereas in a hash it's O(1). Also, efficiency depends on how the tree structure is implemented.
      Now when it comes to putting that theoretical knowledge into code, read IDs are different from the strings always used to demonstrate how suffix trees work for finding matches. In my opinion, a giant Perl hash shouldn't be much of a problem with 32+ GB RAM available, ignoring the fact that one may call this approach "quick and dirty".

      What do the experts out there think?

      Barbara

      Comment


      • #4
        Originally posted by epigen View Post
        Thanks for the answer, Nils. I tried dwgsim_eval but it was not comfortable with the reads being not created by dwgsim and being single end:

        In function "process_bam": Warning[OutOfRange]. Variable/Value: 614 1806 266.
        Message: [dwgsim_eval] read was not generated by dwgsim?.
        ***** Warning *****
        ************************************************************
        ************************************************************
        In function "run": Fatal Error[OutOfRange]. Message: Found a read that was not paired.
        ***** Exiting due to errors *****
        ************************************************************

        By the way, some documentation on the DNAA tools would be really helpful ... As it seems, dwgsim_eval only makes statistics for one BAM or SAM file. What I want to do is to explicitely count how many reads were mapped to a similar coordinate by two different tools, and output these reads.
        For this purpose, I would obviously have to store the read information of one file. In a thread about duplicate removal someone recommended storing it in a trie structure instead of using one of these RAM-greedy Perl hashes. As far as I remember from my computer science lectures, a suffix tree (or eqivalent, a prefix tree) does not need less space than a hash. Searching in a tree is O(n) whereas in a hash it's O(1). Also, efficiency depends on how the tree structure is implemented.
        Now when it comes to putting that theoretical knowledge into code, read IDs are different from the strings always used to demonstrate how suffix trees work for finding matches. In my opinion, a giant Perl hash shouldn't be much of a problem with 32+ GB RAM available, ignoring the fact that one may call this approach "quick and dirty".

        What do the experts out there think?

        Barbara
        You have to use the read generator, otherwise it gets confused. The read name convention can be found at http://sourceforge.net/apps/mediawik...ome_Simulation. I am working on documentation, but you are welcome to help (go open source)!

        A perl hash would be easiest to implement, so try that.

        Comment

        Latest Articles

        Collapse

        • seqadmin
          Latest Developments in Precision Medicine
          by seqadmin



          Technological advances have led to drastic improvements in the field of precision medicine, enabling more personalized approaches to treatment. This article explores four leading groups that are overcoming many of the challenges of genomic profiling and precision medicine through their innovative platforms and technologies.

          Somatic Genomics
          “We have such a tremendous amount of genetic diversity that exists within each of us, and not just between us as individuals,”...
          05-24-2024, 01:16 PM
        • seqadmin
          Recent Advances in Sequencing Analysis Tools
          by seqadmin


          The sequencing world is rapidly changing due to declining costs, enhanced accuracies, and the advent of newer, cutting-edge instruments. Equally important to these developments are improvements in sequencing analysis, a process that converts vast amounts of raw data into a comprehensible and meaningful form. This complex task requires expertise and the right analysis tools. In this article, we highlight the progress and innovation in sequencing analysis by reviewing several of the...
          05-06-2024, 07:48 AM

        ad_right_rmr

        Collapse

        News

        Collapse

        Topics Statistics Last Post
        Started by seqadmin, 05-24-2024, 07:15 AM
        0 responses
        198 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 05-23-2024, 10:28 AM
        0 responses
        220 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 05-23-2024, 07:35 AM
        0 responses
        229 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 05-22-2024, 02:06 PM
        0 responses
        13 views
        0 likes
        Last Post seqadmin  
        Working...
        X