Re: [Seqan-dev] SAM issue?

["Holtgrewe, Manuel" <manuel.holtgrewe@fu-berlin.de>]

Hi Mat,

I think I know understand your problem. The begin and end position the aligned read store are in gap-space, that is the position in the global reads-to-reference alignment. If you want to get the position in the reference without gaps, you have to use the positionGapToSeq function.

typedef typename TFragmentStore::TContigStore                        TContigStore;
typedef typename Value<TContigStore>::Type                           TContig;
typedef typename TFragmentStore::TContigSeq                          TContigSeq;
typedef Gaps<TContigSeq, AnchorGaps<typename TContig::TGapAnchors> > TContigGaps;

typedef typename TFragmentStore::TAlignedReadStore                   TAlignedReadStore;
typedef typename Value<TAlignedReadStore>::Type                      TAlignedRead;
typedef typename TAlignedRead::TPos                                  TAlignedReadPos;

TContigGaps contigGaps(contigStore[contigId].seq, contigStore[contigId].gaps);
TAlignedRead const & alignedRead = alignedReadStore[idx];
// Translate end position from aligned read record to sequence space in reference.
TAlignedReadPos endPos = positionGapToSeq(contigGaps, alignedRead.endPos);


In the build_gold_standard.h file of app rabema, you can see this in action (from line 470).

Bests,
Manuel

Am 21.03.2011 um 11:11 schrieb Mat:

> Hi Manuel!
> 
> I tried this but the generated SAM looks fine. Still i get the wrong coordinates displayed- I simplified the loop to a basic output (cpp attached).
> 
> So i tried to increase the filesize to see what happens:
> 
> head -n 70000 test.sam > bug2.sam
> illumina_80bp_3kb.000000571
> 19626,19706
> 
> head -n 80000 test.sam > bug2.sam
> illumina_80bp_3kb.000000571
> 19628,19708
> 
> head -n 90000 test.sam > bug2.sam
> illumina_80bp_3kb.000000571
> 19629,19709
> 
> Pretty weired! Looks to me like an overflow problem? I am running a Mac with OSX...
> 
> Could you check if you run into the same problems on windows?
> 
> I don't want to spam the list with parts the SAM file but you can download the entire sam file from here:
> 
> http://rapidshare.com/files/453617711/test.sam
> 
> Thanks!
> 
> best,
> 
> mat
>  
> 
> 
> 
> Am 3/18/11 5:06 PM, schrieb Manuel Holtgrewe:
>> Another thing you could try is to read the SAM file, and write it out immediately after and see whether there are any problems with the generated SAM. 
>> 
>> *m 
>> 
>> Am 18.03.2011 um 17:05 schrieb Mat: 
>> 
>>> Hi Manuel&thanks for the reply! 
>>> 
>>> I only send you the head of the SAM file i used- i don't get the issue with that one either - so now i have to figure out what's wrong there... 
>>> 
>>> Ill text you guys if i find an issue - thanks a lot! 
>>> 
>>> cheers 
>>> 
>>> mat 
>>> 
>>> Am 3/18/11 4:46 PM, schrieb Manuel Holtgrewe: 
>>>> 
>>>> Hi, Mat. 
>>>> 
>>>> Here is what I get: 
>>>> 
>>>> $ ./demos/SAMtest -s bug.sam | grep -A 1 571 
>>>> illumina_80bp_3kb.000000571 
>>>> illumina_80bp_3kb.000000571 
>>>> 19625,19705 and 22480,22400 
>>>> 
>>>> Which is OK since indices are 1-based in SAM and 0-based in SeqAn. 
>>>> 
>>>> I just saw that there were uninitialized variables in your code, try to compile it with all warnings enabled and in optimized mode (the compiler only detects certain things when           you tell him to think harder about the code). Me getting the right values (probably my OS initializes data differently from yours), and you wrong indicates problems with initializations. 
>>>> 
>>>> Normally, I would also point you at valgrind, but it did not find the uninitialized memory on the stack :) 
>>>> 
>>>> Below is the C++ output: 
>>>> 
>>>> Building CXX object demos/CMakeFiles/SAMtest.dir/Users/manuel/Development/seqan-trunk/projects/library/demos/SAMtest.cpp.o
>>>> /Users/manuel/Development/seqan-trunk/projects/library/demos/SAMtest.cpp: In function ‘int main(int, const char**)’: 
>>>> /Users/manuel/Development/seqan-trunk/projects/library/demos/SAMtest.cpp:57: warning: ‘cId1’ may be used uninitialized in this function 
>>>> /Users/manuel/Development/seqan-trunk/projects/library/demos/SAMtest.cpp:57: warning: ‘id1’ may be used uninitialized in this function 
>>>> 
>>>> 
>>>> 
>>>> HTH 
>>>> 
>>>> Am 18.03.2011 um 15:52 schrieb Mat: 
>>>> 
>>>>> Hey guys! 
>>>>> 
>>>>> I have a strange issue loading a SAM-File into a FragmentStore. The SAM file has the following paired-read alignment: 
>>>>> ... 
>>>>> illumina_80bp_3kb.000000571     99      contig00002     19626   60      80M     =       22401   2855 
>>>>> illumina_80bp_3kb.000000571     147     contig00002     22401   60      80M     =       19626   -285 
>>>>> 
>>>>> 
>>>>> I load this SAM file into a FragmentStore and sort it by PairMatchId. Then i extract all the pairs that map on the same contig: 
>>>>> 
>>>>> illumina_80bp_3kb.000000571 
>>>>> illumina_80bp_3kb.000000571 
>>>>> 19638,19718 and 22495,22415 
>>>>> 
>>>>> I don't really get why these coordinates differ?! 
>>>>> 
>>>>> Thanks a lot!! 
>>>>> 
>>>>> (see sourcecode attached)... 
>>>>> <bug.sam><SAMtest.cpp><ATT00001..txt> 
>>>> 
>>> <ATT00001..txt> 
>> 
> <SAMtest.cpp><ATT00001..txt>