Seqanswers Leaderboard Ad

Collapse

Announcement

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

  • nilshomer
    replied
    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.

    Leave a comment:


  • epigen
    replied
    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

    Leave a comment:


  • nilshomer
    replied
    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.

    Leave a comment:


  • epigen
    started a topic Compare mapped reads from different aligners

    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

Latest Articles

Collapse

  • seqadmin
    Best Practices for Single-Cell Sequencing Analysis
    by seqadmin



    While isolating and preparing single cells for sequencing was historically the bottleneck, recent technological advancements have shifted the challenge to data analysis. This highlights the rapidly evolving nature of single-cell sequencing. The inherent complexity of single-cell analysis has intensified with the surge in data volume and the incorporation of diverse and more complex datasets. This article explores the challenges in analysis, examines common pitfalls, offers...
    06-06-2024, 07:15 AM
  • 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

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, 06-17-2024, 06:54 AM
0 responses
11 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-14-2024, 07:24 AM
0 responses
22 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-13-2024, 08:58 AM
0 responses
17 views
0 likes
Last Post seqadmin  
Started by seqadmin, 06-12-2024, 02:20 PM
0 responses
20 views
0 likes
Last Post seqadmin  
Working...
X