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
          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
        • seqadmin
          Essential Discoveries and Tools in Epitranscriptomics
          by seqadmin




          The field of epigenetics has traditionally concentrated more on DNA and how changes like methylation and phosphorylation of histones impact gene expression and regulation. However, our increased understanding of RNA modifications and their importance in cellular processes has led to a rise in epitranscriptomics research. “Epitranscriptomics brings together the concepts of epigenetics and gene expression,” explained Adrien Leger, PhD, Principal Research Scientist...
          04-22-2024, 07:01 AM

        ad_right_rmr

        Collapse

        News

        Collapse

        Topics Statistics Last Post
        Started by seqadmin, 05-14-2024, 07:03 AM
        0 responses
        20 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 05-10-2024, 06:35 AM
        0 responses
        44 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 05-09-2024, 02:46 PM
        0 responses
        54 views
        0 likes
        Last Post seqadmin  
        Started by seqadmin, 05-07-2024, 06:57 AM
        0 responses
        42 views
        0 likes
        Last Post seqadmin  
        Working...
        X