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