[Seqan-dev] Bug in FragmentStore?

[Mat <matthias.dodt@mdc-berlin.de>]

Hi guys!

I think there is a problem with the fragment store reading paired-end reads from SAM. The pairMatchIDs seem to be assigned incorrectly. As an input i used a sam file with the following statistics:
87323 in total
0 QC failure
0 duplicates
86119 mapped (98.62%)
87323 paired in sequencing
43662 read1
43661 read2
86965 properly paired (99.59%)
86002 with itself and mate mapped
117 singletons (0.13%)
3341 with mate mapped to a different chr
3242 with mate mapped to a different chr (mapQ>=5)


The code...:

...
unsigned int nrAlignments = length(fragStore.alignedReadStore),mateCounter=0,spanMateCnt=0;
    sortAlignedReads(fragStore.alignedReadStore, SortReadId());

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

        TValue pair_id = fragStore.alignedReadStore[n].pairMatchId;
        TValue currentReadId=fragStore.alignedReadStore[n].readId;

        if(pair_id!=TReadStoreElement::INVALID_ID){

            std::cout << "read-id: " << currentReadId << std::endl;
            std::cout << "pair-id: " << pair_id << std::endl;

            //get mate read-id
            TValue r1 = fragStore.matePairStore[pair_id].readId[0];
            TValue r2 = fragStore.matePairStore[pair_id].readId[1];

            //now r1 or r2 should equal currentReadId
            if((r1!=currentReadId)&&(r2!=currentReadId)){
                std::cout << "Wrong mate association" << std::endl;
            }
            else {
                ++spanMateCnt;
            }

}
}
}
...

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!
<<

Any suggestions?

thanks!

best,

mat