Example Program
Interval Tree
Construction and usage demo for the interval tree.
An example for using interval trees.
1#include <iostream>
2#include <seqan/graph_align.h>
3
4using namespace seqan;
5
6int main()
7{
8
9    typedef CharString TCargo;  // id type
10    typedef int TValue;         // position type
11
12    typedef IntervalAndCargo<TValue, TCargo> TInterval;
13    typedef IntervalTree<TValue, TCargo> TIntervalTree;
14
15    String<TInterval> intervals;
16    resize(intervals, 5);
17
Store gene annotation in intervals.
18    intervals[0].i1 = 5;   intervals[0].i2 = 1000;
19    intervals[0].cargo = "gene";
20
21    intervals[1].i1 = 50;  intervals[1].i2 = 200;
22    intervals[1].cargo = "exon";
23
24    intervals[2].i1 = 600; intervals[2].i2 = 800;
25    intervals[2].cargo = "exon";
26
27    intervals[3].i1 = 100; intervals[3].i2 = 200;
28    intervals[3].cargo = "coding";
29
30    intervals[4].i1 = 600; intervals[4].i2 = 700;
31    intervals[4].cargo = "coding";
32
33    TIntervalTree tree(intervals);
34
Add another interval.
35    TInterval interval;
36    interval.i1 = 200; interval.i2 = 600;
37    interval.cargo = "intron";
38
39    addInterval(tree, interval);
40
Query a genomic region.
41    TValue delBegin = 300;
42    TValue delEnd   = 500;
43    String<TCargo> results;
44
45    findIntervals(tree, delBegin, delEnd, results);
46
47    std::cout << "Deletion " << delBegin << ".." << delEnd << " overlaps with ";
48    for (unsigned i = 0; i < length(results); ++i)
49        std::cout << results[i] << ",";
50    std::cout << std::endl;
51
Query a single position.
52    TValue snpPos = 150;
53    findIntervals(tree, snpPos, results);
54
55    std::cout << "SNP " << snpPos << " overlaps with ";
56    for (unsigned i = 0; i < length(results); ++i)
57        std::cout << results[i] << ",";
58    std::cout << std::endl;
59
60    CharString iCargo("exon");
61    bool res = removeInterval(tree, 50, 200, iCargo);
62    if (res)
63        std::cout << "Removed exon interval 50..200.\n";
64
Now, redo the query. This time one interval less should be returned.
65    String<TCargo> results2;
66    findIntervals(tree, snpPos, results2);
67
68    std::cout << "SNP " << snpPos << " overlaps with ";
69    for (unsigned i = 0; i < length(results2); ++i)
70        std::cout << results2[i] << ",";
71    std::cout << std::endl;
72
73    return 0;
74}
SeqAn - Sequence Analysis Library - www.seqan.de
 

Page built @2013/07/11 09:12:35