Re: [Seqan-dev] Bug in FragmentStore?

["Weese, David" <weese@campus.fu-berlin.de>]

Hi Mat,

the matePairId associates two reads to a mate-pair. Thus there are at most 2 reads with the same matePairId. The alignedReadStore stores matches of reads. Every read can have multiple matches. Each match has a unique id. Two matches of reads of the same mate-pair can form a mate-pair match. If so they have the same pairMatchId which is unique for these 2 matches.

Example:

X1 X2 Y1 Y2 Y3

2 reads A and B have matches at positions Xi and Yj. Theoretically 2*3 pairs of matches are possible. If only X1-Y1 and X2-Y3 form a mate pair match (due to their insert size) the alignedReadStore could look as follows:

read, position, pairMatchId

A, X1, 23

A, X2, 42

B, Y1, 23

B, Y3, 42

The pairMatchIds can be chosen arbitrarily (non-sequential) as long as there are exactly two matches with the same pairMatchId.

There is a function compactPairMatchIds that renumerates the pairMatchIds sequentially from 0 and returns the number of pair matches so that you can store additional information in an extra array using the pairMatchId.

If you need two matches of all pair matches you can call:

sortAlignedReads(me.alignedReadStore, SortPairMatchId());

to resort the alignedReadStore after the pairMatchId so that matches of a pair match are adjacent.

Cheers,

David

PS: Does the error still occur? If yes, please create a ticket with SAM and source files.

--
David Weese weese@inf.fu-berlin.de

Freie Universität Berlin http://www.inf.fu-berlin.de/

Institut für Informatik Phone: +49 30 838 75246

Takustraße 9 Algorithmic Bioinformatics

14195 Berlin Room 021 

Am 04.01.2011 um 16:47 schrieb Mat:

Hi Manuel!

Thanks for your reply. I can send you some source later, but i think its a more simple problem. I figured out that the pairMatchId of the FragmentStore's alignedReadStore is in general not the same as the matePairId from the matePairStore by looking into the source. What still is missing is: I would like to know on which contig the two reads of a mate-pair aligned. So far i got:

- First i iterate over the alignedReadStore to get a readId
- i get the matePairId via the readStore
- then i get the mate's id via the matePairStore

    - now i need a function to search for a readId in alignedReadStore efficiently - is there something like that?

The code looks now like this:

...
sortAlignedReads(fragStore.alignedReadStore, SortReadId());

    for (unsigned int n=0;n<nrAlignments;++n){

        //get the current readId
        TValue currentReadId=fragStore.alignedReadStore[n].readId;

        //get the correct mate-pair id
        TValue pair_id = fragStore.readStore[currentReadId].matePairId;

        if(pair_id!=TReadStoreElement::INVALID_ID){
            ++mateCounter;
            std::cout << "read-id: " << currentReadId << "n: " << n << std::endl;

            TValue r1 = fragStore.matePairStore[pair_id].readId[0];
            TValue r2 = fragStore.matePairStore[pair_id].readId[1];
            TValue other = 0;
           
            //now r1 or r2 should equal currentReadId
            if(r1==currentReadId){
                other = r2;
                std::cout << "other is r2:" << other << std::endl;
               
                //now get the alignId
                TValue alignId;
                ...how to do that?


So how do i get the alignId(s) of a mate's readId?

Thanks!

best,

mat


Am 1/4/11 4:16 PM, schrieb Holtgrewe, Manuel:

It turns out that the both pair_id's never match and the loop breaks with:

Trying to access an element behind the last one!

<<

This assertion is in the operator[]/value(str, pos) function for the String class. Did you try to run the program within a debugger to see where the problem occurs in your program?

Please provide a complete, compiling problem and a minimal SAM file (should be 2 lines after the header, right?) that allow us to reproduce the problem.

Bests,
Manuel
_______________________________________________
seqan-dev mailing list
seqan-dev@lists.fu-berlin.de
https://lists.fu-berlin.de/listinfo/seqan-dev

<ATT00001..txt>