Re: [Seqan-dev] problem with interval tree - assertion failed

["Emde, Anne-Katrin" <Anne-Katrin.Emde@fu-berlin.de>]

Hi Mat,

 

there seems to be a problem with the addInterval function, when the interval tree is empty. I will fix that as soon as possible.

 

However, in your case it would make sense to create a string of intervals, and then construct the interval tree on the whole set of intervals at once. This will be more efficient and may construct a more balanced tree. addInterval was originally meant to add intervals later on to an existing tree (anyway, it should of course also work on an empty tree!).

You can see an example in demos/interval_tree.cpp. I’m copy&pasting the relevant part here:

 

 

        typedef CharString TCargo;      // id type

        typedef int TValue;             // position type

 

        typedef IntervalAndCargo<TValue, TCargo> TInterval;

        typedef IntervalTree<TValue, TCargo> TIntervalTree;

 

        String<TInterval> intervals;

        resize(intervals,3);

 

        // store gene annotation in intervals

        intervals[0].i1 = 5;   intervals[0].i2 = 1000;

        intervals[0].cargo = "gene";

 

        intervals[1].i1 = 50;  intervals[1].i2 = 200;

        intervals[1].cargo = "exon";

 

        intervals[2].i1 = 600; intervals[2].i2 = 800;

        intervals[2].cargo = "exon";

 

        TIntervalTree tree(intervals);

 

Best,

Anne-Katrin

 

 

 

 

---

From: Mat [mailto:matthias.dodt@mdc-berlin.de]
Sent: Freitag, 26. November 2010 13:52
To: seqan-dev@lists.fu-berlin.de
Subject: [Seqan-dev] problem with interval tree - assertion failed

 

Hey guys!

I got a problem with the interval tree. I want to load contigs and create intervals of the contigs sizes:

    //Interval tree construction
    typedef int TValue; //position type
    typedef gip::MappedFragment TCargo; // interval id type
    typedef IntervalAndCargo<TValue, TCargo*> TInterval;// interval type
    typedef IntervalTree<TValue, TCargo*> TIntervalTree;// interval tree type

    std::cout << "Loading contigs file..." << std::endl;
    if(!loadContigs(fragStore,genomeFileName)){
            std::cerr << "Error opening genome file!" << std::endl;
            return -1;
    }

    unsigned int numberOfContigs = length(fragStore.contigNameStore);
    std::cout << "Found " << numberOfContigs << " contigs in genome file." << std::endl;

    std::map<CharString, TIntervalTree*> contigMap;

    for(unsigned int n=0;n<numberOfContigs;++n){
        std::cout << "NAME: " << fragStore.contigNameStore[n] << std::endl;
        std::cout << "LENGTH: " << length(fragStore.contigStore[n].seq) << std::endl;
        TIntervalTree *tree = new TIntervalTree;
        TInterval in;
        in.i1 = 0;
        in.i2 = length(fragStore.contigStore[n].seq);
        in.cargo = NULL;

        addInterval(*tree, in);

     }

output:
=======
Loading contigs file...
Found 40 contigs in genome file.
/Users/mat/cpp_develop/seqan-trunk/projects/library/seqan/sequence/string_base.hNAME: 1
LENGTH: 6844
:219 Assertion failed : static_cast<TStringPos>(pos) < static_cast<TStringPos>(length(me)) was: 0 >= 0 (Trying to access an element behind the last one!)

if i remove the "addInterval" statement i get the output for each contig...

thanks!

best,

mat