| Example Program Mummy Simple MUMmer clone. MUMmer is a tool to search for multiple exact matches (MUMs) between 2 given sequences. MUMs can be used as a starting point for a multiple genome alignment algorithm. This example shows how to implement a simple version of MUMer to find multiple exact matches of n sequences (n≥2) in SeqAn. File "index_mummy.cpp"
1 | #include <iostream>
2 | #include <fstream>
3 | #include <seqan/index.h>
4 |
5 | using namespace std;
6 | using namespace seqan;
7 |
8 |
9 | template <typename TIndex>
10 | void findMUMs(TIndex &esa, unsigned minLen)
11 | {
12 | typename Iterator<TIndex, Mums>::Type it(esa, minLen); // set minimum MUM length
13 | String< typename SAValue<TIndex>::Type > occs; // temp. string storing the hit positions
14 |
15 | cout << resetiosflags(ios::left);
16 | while (!atEnd(it))
17 | {
18 | occs = getOccurrences(it); // gives hit positions (seqNo,seqOfs)
19 | orderOccurrences(occs); // order them by seqNo
20 |
21 | for(unsigned i = 0; i < length(occs); ++i)
22 | cout << setw(8)
23 | << getValueI2(occs[i]) + 1 // output them in MUMmer's output format
24 | << " ";
25 |
26 | cout << setw(8)
27 | << repLength(it)
28 | << endl;
29 |
30 | ++it;
31 | }
32 | cout << setiosflags(ios::left) << endl;
33 | }
34 |
35 |
36 | template <typename TSpec>
37 | int runMummy(int argc, const char *argv[], unsigned seqCount, unsigned minLen)
38 | {
39 | typedef String<Dna5, TSpec> TString;
40 | typedef StringSet<TString> TStringSet;
41 | typedef Index<TStringSet> TIndex;
42 |
43 | TIndex index;
44 |
45 | // import sequences
46 | resize(indexText(index), seqCount);
47 | for(int arg = 1, seq = 0; arg < argc; ++arg)
48 | {
49 | // skip two argument options
50 | if (strcmp(argv[arg], "-p") == 0 || strcmp(argv[arg], "--profile") == 0 ||
51 | strcmp(argv[arg], "-l") == 0 || strcmp(argv[arg], "--minlen") == 0)
52 | {
53 | ++arg;
54 | continue;
55 | }
56 |
57 | if (argv[arg][0] != '-') {
58 | ifstream file;
59 | file.open(argv[arg], ios_base::in | ios_base::binary);
60 | if (!file.is_open()) {
61 | cerr << "Import of sequence " << argv[arg] << " failed." << endl;
62 | return 1;
63 | }
64 | read(file, indexText(index)[seq], Fasta());
65 | file.close();
66 | ++seq;
67 | }
68 | }
69 | cerr << lengthSum(indexText(index)) << " bps sequence imported." << endl;
70 |
71 | findMUMs(index, minLen);
72 |
73 | return 0;
74 | }
75 |
76 |
77 | void printHelp(int, const char *[], bool longHelp = false)
78 | {
79 | cerr << "***************************************" << endl;
80 | cerr << "*** Simple MUM finder ***" << endl;
81 | cerr << "*** written by David Weese (c) 2007 ***" << endl;
82 | cerr << "***************************************" << endl << endl;
83 | cerr << "Usage: mummy [OPTION]... <SEQUENCE FILE> ... <SEQUENCE FILE>" << endl;
84 | if (longHelp) {
85 | cerr << endl << "Options:" << endl;
86 | cerr << " -e, --extern \t" << "use external memory (for large datasets)" << endl;
87 | cerr << " -l, --minlen \t" << "set minimum MUM length" << endl;
88 | cerr << " \t" << "if not set, default value is 20" << endl;
89 | cerr << " -h, --help \t" << "print this help" << endl;
90 | } else {
91 | cerr << "Try 'mummy --help' for more information." << endl;
92 | }
93 | }
94 |
95 |
96 | int main(int argc, const char *argv[])
97 | {
98 | if (argc < 2) {
99 | printHelp(argc, argv);
100 | return 0;
101 | }
102 |
103 | unsigned optMinLen = 20;
104 | bool optExtern = false;
105 |
106 | unsigned seqCount = 0;
107 | for(int arg = 1; arg < argc; ++arg) {
108 | if (argv[arg][0] == '-') {
109 | // parse option
110 |
111 | if (strcmp(argv[arg], "-e") == 0 || strcmp(argv[arg], "--extern") == 0) {
112 | // use external memory algorithms
113 | optExtern = true;
114 | continue;
115 | }
116 | if (strcmp(argv[arg], "-l") == 0 || strcmp(argv[arg], "--minlen") == 0) {
117 | // minimum match length
118 | if (arg + 1 == argc) {
119 | printHelp(argc, argv);
120 | return 0;
121 | }
122 | ++arg;
123 | optMinLen = atoi(argv[arg]);
124 | continue;
125 | }
126 | if (strcmp(argv[arg], "-h") == 0 || strcmp(argv[arg], "--help") == 0) {
127 | // print help
128 | printHelp(argc, argv, true);
129 | return 0;
130 | }
131 | } else {
132 | // parse sequence file
133 | ++seqCount;
134 | }
135 | }
136 |
137 | if (optExtern)
138 | return runMummy<External<> >(argc, argv, seqCount, optMinLen);
139 | else
140 | return runMummy<Alloc<> >(argc, argv, seqCount, optMinLen);
141 | }
Output If you run the tool on 2 sequences it outputs exactly the same matches as MUMmer (called with -mum option), it only differs in the order of outputted matches. To get matches with increasing positions at the first sequence we piped the output to sort. As an example data set we used 3 strains of chlamydia bacterium (NC_002620.fna, NC_000117.fna, NC_007429.fna) and saved the Fasta files to the demos directory. weese@tanne:~/seqan$ cd demos weese@tanne:~/seqan/demos$ make index_mummy weese@tanne:~/seqan/demos$ ./index_mummy -h *************************************** *** Simple MUM finder *** *** written by David Weese (c) 2007 *** *************************************** Usage: mummy [OPTION]... <SEQUENCE FILE> ... <SEQUENCE FILE> Options: -e, --extern use external memory (for large datasets) -l, --minlen set minimum MUM length if not set, default value is 20 -h, --help print this help weese@tanne:~/seqan/demos$ ./index_mummy NC*.fna |sort > mums.txt 3159928 bps sequence imported. weese@tanne:~/seqan/demos$ head mums.txt 1565 323805 2159 48 1646 323886 2240 27 1722 323962 2316 37 1774 324014 2368 26 1941 324181 2535 23 2061 324301 2655 35 2102 324342 2696 29 2132 324372 2726 20 2183 324423 2777 24 weese@tanne:~/seqan/demos$ See Index, Mums Iterator, getOccurrences, representative, orderOccurrences, repLength, StringSet SeqAn - Sequence Analysis Library - www.seqan.de

Page built @2012/03/14 16:34:46