1 //![main]
2 #include <iostream>
3 #include <seqan/align.h>
4 
5 using namespace seqan;
6 
main()7 int main()
8 {
9     typedef String<char> TSequence;
10     typedef Gaps<TSequence, ArrayGaps> TGaps;
11     typedef Iterator<TGaps>::Type TGapsIterator;
12     typedef Iterator<String<int> >::Type TIterator;
13 
14     TSequence text =    "MISSISSIPPIANDMISSOURI";
15     TSequence pattern = "SISSI";
16 
17     String<int> locations;
18     for (unsigned i = 0; i < length(text) - length(pattern); ++i)
19     {
20         // Compute the MyersBitVector in current window of text.
21         TSequence tmp = infix(text, i, i + length(pattern));
22 
23         // Report hits with at most 2 errors.
24         if (globalAlignmentScore(tmp, pattern, MyersBitVector()) >= -2)
25         {
26             appendValue(locations, i);
27         }
28     }
29 
30     TGaps gapsText;
31     TGaps gapsPattern;
32     assignSource(gapsPattern, pattern);
33     std::cout << "Text: " << text << "\tPattern: " << pattern << std::endl;
34     for (TIterator it = begin(locations); it != end(locations); ++it)
35     {
36         // Clear previously computed gaps.
37         clearGaps(gapsText);
38         clearGaps(gapsPattern);
39 
40         // Only recompute the area within the current window over the text.
41         TSequence textInfix = infix(text, *it, *it + length(pattern));
42         assignSource(gapsText, textInfix);
43 
44         // Use semi-global alignment since we do not want to track leading/trailing gaps in the pattern.
45         // Restirct search space using a band allowing at most 2 errors in vertical/horizontal direction.
46         int score = globalAlignment(gapsText, gapsPattern, Score<int>(0, -1, -1), AlignConfig<true, false, false, true>(), -2, 2);
47 
48         TGapsIterator itGapsPattern = begin(gapsPattern);
49         TGapsIterator itGapsEnd = end(gapsPattern);
50 
51         // Remove trailing gaps in pattern.
52         int count = 0;
53         while (isGap(--itGapsEnd))
54             ++count;
55         setClippedEndPosition(gapsPattern, length(gapsPattern) - count);
56 
57         // Remove leading gaps in pattern.
58         if (isGap(itGapsPattern))
59         {
60             setClippedBeginPosition(gapsPattern, countGaps(itGapsPattern));
61             setClippedBeginPosition(gapsText, countGaps(itGapsPattern));
62         }
63 
64         // Reinitilaize the iterators.
65         TGapsIterator itGapsText = begin(gapsText);
66         itGapsPattern = begin(gapsPattern);
67         itGapsEnd = end(gapsPattern);
68 
69         // Use a stringstream to construct the cigar string.
70         std::stringstream cigar;
71         while (itGapsPattern != itGapsEnd)
72         {
73             // Count insertions.
74             if (isGap(itGapsText))
75             {
76                 int numGaps = countGaps(itGapsText);
77                 cigar << numGaps << "I";
78                 itGapsText += numGaps;
79                 itGapsPattern += numGaps;
80                 continue;
81             }
82             // Count deletions.
83             if (isGap(itGapsPattern))
84             {
85                 int numGaps = countGaps(itGapsPattern);
86                 cigar << numGaps << "D";
87                 itGapsText += numGaps;
88                 itGapsPattern += numGaps;
89                 continue;
90             }
91             ++itGapsText;
92             ++itGapsPattern;
93         }
94         // Output the hit position in the text, the total number of edits and the corresponding cigar string.
95         std::cout << "Hit at position  " << *it << "\ttotal edits: " << abs(score) << std::endl;
96     }
97 
98     return 0;
99 }
100 //![main]
101