Compare mapped reads from different aligners

[epigen]

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

[nilshomer]

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.

[epigen]

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

[nilshomer]

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.