1 #include <boost/test/auto_unit_test.hpp>
2 #include <boost/filesystem/operations.hpp>
3 namespace fs = boost::filesystem;
4 #include <boost/assign/list_of.hpp>
5 #include <boost/assign/list_inserter.hpp>
6 #include <boost/assign.hpp>
7 namespace assign = boost::assign;
13 #include "alg/mussa.hpp"
17 //! can we initialize a mussa object?
18 BOOST_AUTO_TEST_CASE( mussa_simple )
21 BOOST_CHECK_EQUAL(m.get_name(), "" );
22 BOOST_CHECK_EQUAL(m.get_window(), 0);
23 BOOST_CHECK_EQUAL(m.get_threshold(), 0);
24 BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::TransitiveNway);
25 m.set_name( "hello" );
26 BOOST_CHECK_EQUAL(m.get_name(), "hello" );
28 BOOST_CHECK_EQUAL(m.get_window(), 30);
30 BOOST_CHECK_EQUAL(m.get_threshold(), 21);
31 BOOST_CHECK_EQUAL(m.get_soft_threshold(), 21);
32 m.set_soft_threshold(19);
33 BOOST_CHECK_EQUAL(m.get_soft_threshold(), 21);
34 m.set_soft_threshold(35);
35 BOOST_CHECK_EQUAL(m.get_soft_threshold(), 30);
36 m.set_soft_threshold(25);
37 BOOST_CHECK_EQUAL(m.get_soft_threshold(), 25);
38 m.set_analysis_mode(Mussa::RadialNway);
39 BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::RadialNway);
42 BOOST_CHECK_EQUAL(m.get_name(), "" );
43 BOOST_CHECK_EQUAL(m.get_window(), 0);
44 BOOST_CHECK_EQUAL(m.get_threshold(), 0);
45 BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::TransitiveNway);
48 BOOST_AUTO_TEST_CASE( mussa_analysis_name )
51 m.set_analysis_mode( Mussa::TransitiveNway );
52 BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Transitive" );
53 m.set_analysis_mode( Mussa::RadialNway );
54 BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Radial" );
55 m.set_analysis_mode( Mussa::EntropyNway );
56 BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Entropy" );
57 m.set_analysis_mode( Mussa::RecursiveNway);
58 BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "[deprecated] Recursive" );
61 BOOST_AUTO_TEST_CASE( mussa_sequences )
63 std::string s0("AAAANNNN");
64 std::string s1("GGGGNNNN");
65 std::string s2("TTTTNNNN");
68 analysis.append_sequence(s0);
69 analysis.append_sequence(s1);
70 analysis.append_sequence(s2);
72 BOOST_CHECK_EQUAL( analysis.sequences().size(), 3 );
73 BOOST_CHECK_EQUAL( *(analysis.sequences()[0]), s0);
74 BOOST_CHECK_EQUAL( *(analysis.sequences()[1]), s1);
75 BOOST_CHECK_EQUAL( *(analysis.sequences()[2]), s2);
78 // for some reason we can call nway once safely but it
79 // somehow changed things and caused a segfault
80 // fixed by adding a return statement in trans_path_search
81 BOOST_AUTO_TEST_CASE ( empty_mussa_set_threshold )
84 m.set_soft_threshold(15);
87 m.set_soft_threshold(25);
91 BOOST_AUTO_TEST_CASE( mussa_load_mupa )
93 fs::path mupa_path(EXAMPLE_DIR, fs::native);
94 fs::path result_path = fs::initial_path() / "mck3test_w30_t20";
95 mupa_path /= "mck3test.mupa";
98 m1.load_mupa_file( mupa_path );
100 m1.save( result_path );
101 BOOST_CHECK_EQUAL( m1.get_name(), std::string("mck3test") );
102 BOOST_CHECK( m1.size() > 0 );
105 m2.load( result_path );
106 BOOST_CHECK_EQUAL( m2.get_name(), result_path.leaf() );
107 BOOST_CHECK_EQUAL( m1.size(), m2.size() );
110 BOOST_AUTO_TEST_CASE( mussa_load_full_path )
113 fs::path full_path(fs::path(EXAMPLE_DIR, fs::native) / "mck3test.mupa");
114 m1.load_mupa_file( full_path );
117 BOOST_CHECK( m1.size() > 0);
118 BOOST_CHECK_EQUAL( m1.get_window(), 30 );
119 BOOST_CHECK_EQUAL( m1.get_threshold(), 20);
122 BOOST_AUTO_TEST_CASE( mussa_load_analysis )
124 fs::path example_dir(EXAMPLE_DIR, fs::native);
126 m1.load_mupa_file( example_dir / "mck3test.mupa" );
130 m2.load( fs::initial_path() / "mck3test_w30_t20");
132 BOOST_CHECK_EQUAL( m1.size(), m2.size() );
133 BOOST_CHECK_EQUAL( m1.get_window(), m2.get_window() );
134 BOOST_CHECK_EQUAL( m1.get_threshold(), m2.get_threshold() );
137 BOOST_AUTO_TEST_CASE( mussa_load_motif )
139 string data = "AAGG 1.0 1.0 0.0\n"
143 istringstream test_istream(data);
146 m1.append_sequence("AAAAGGGGTTTT");
147 m1.append_sequence("GGGCCCCTTGGTT");
148 m1.load_motifs(test_istream);
150 for (Mussa::vector_sequence_type::const_iterator seq_i = m1.sequences().begin();
151 seq_i != m1.sequences().end();
154 BOOST_CHECK( (*seq_i)->motifs().size() > 0 );
158 BOOST_AUTO_TEST_CASE( mussa_add_motif )
160 vector<Sequence> motifs;
161 motifs.push_back("AAGG");
162 vector<Color> colors;
163 colors.push_back(Color(1.0, 0.0, 0.0));
166 m1.append_sequence("AAAAGGGGTTTT");
167 m1.append_sequence("GGGCCCCTTGGTT");
168 m1.set_motifs(motifs, colors);
169 int first_size = m1.motifs().size();
170 BOOST_CHECK_EQUAL( first_size, 1 );
171 BOOST_REQUIRE(first_size > 0);
172 BOOST_CHECK_EQUAL(*(m1.motifs().begin()), motifs.front());
173 // make sure that our sequences have the right number of motifs
174 BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 1);
175 BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1); // because of rc
177 // verify that setting the motif clears the arrays
178 m1.set_motifs(motifs, colors);
179 BOOST_CHECK_EQUAL( first_size, m1.motifs().size() );
180 // make sure that our sequences have the right number of motifs
181 BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 1);
182 BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
184 // add a different motif
186 motifs.push_back("CCTTGG");
187 BOOST_CHECK_EQUAL(motifs.size(), 1);
188 m1.set_motifs(motifs, colors);
189 BOOST_CHECK_EQUAL(m1.motifs().size(), 1);
190 BOOST_REQUIRE(m1.motifs().size() > 0);
191 BOOST_CHECK_EQUAL(*(m1.motifs().begin()), motifs.front());
192 BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 0);
193 BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
195 // try a motif that doesn't exist
197 motifs.push_back("CCTTGG");
198 BOOST_CHECK_EQUAL(motifs.size(), 1);
199 m1.set_motifs(motifs, colors);
200 BOOST_CHECK_EQUAL(m1.motifs().size(), 1);
201 BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 0);
202 BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
207 local_align_test(const Mussa::vector_sequence_type &seqs,
208 const list<ConservedPath::path_type>& result,
209 const list<vector<bool> >& reversed)
211 map<char, vector <char> > m;
212 assign::insert(m)('A', assign::list_of('A')('T') )
213 ('T', assign::list_of('T')('A') )
214 ('G', assign::list_of('G')('C') )
215 ('C', assign::list_of('C')('G') );
216 list<vector<bool> >::const_iterator rc_i = reversed.begin();
218 for(list<ConservedPath::path_type>::const_iterator base_i = result.begin();
219 base_i != result.end();
222 // since the reverse compliment flag is relative to the first sequence
223 // the first one should always be false
224 BOOST_CHECK_EQUAL( (*rc_i)[0], false );
225 const int first_path_basepair_index = (*base_i)[0];
226 const int second_path_basepair_index = (*base_i)[1];
227 const char first_basepair = (*seqs[0])[first_path_basepair_index];
228 const char second_basepair = (*seqs[1])[second_path_basepair_index];
229 // get our index into our reverse compliment map m
230 const int second_compliment_index = (*rc_i)[1];
231 // lookup the forward or reverse compliment depending on our rc flag
232 const char complimented_second = m[second_basepair][second_compliment_index];
234 BOOST_CHECK_EQUAL( first_basepair, complimented_second) ;
239 BOOST_AUTO_TEST_CASE( local_alignment )
241 string s0("GCGCATAT");
242 string s1("AAAAAAAT");
246 analysis.append_sequence(s0);
247 analysis.append_sequence(s1);
248 analysis.set_threshold(3);
249 analysis.set_window(4);
251 NwayPaths npath = analysis.paths();
252 list<ConservedPath::path_type> result;
253 list<vector<bool> > reversed;
254 list<ConservedPath>::iterator pathz_i = npath.pathz.begin();
256 list<ConservedPath> selected_paths;
257 selected_paths.push_back(*pathz_i);
258 analysis.createLocalAlignment(selected_paths.begin(),
259 selected_paths.end(),
263 local_align_test(analysis.sequences(), result, reversed);
268 selected_paths.clear();
269 selected_paths.push_back(*pathz_i);
270 analysis.createLocalAlignment(selected_paths.begin(),
271 selected_paths.end(),
274 local_align_test(analysis.sequences(), result, reversed);
279 BOOST_AUTO_TEST_CASE( subanalysis )
281 Sequence s1("AATGAAGATTTTAATGCTTTAATTTTGTTTTGTAAACTTCGAATTTCCAAAATTTGAAA");
282 Sequence s2("AGGAGCAAGTTCGCTTCATCGAGAATTTTTAATTTTTAGTCAAATTTTCCAATGTCTGA");
285 analysis.append_sequence(s1);
286 analysis.append_sequence(s2);
287 analysis.set_threshold(8);
288 analysis.set_window(8);
291 NwayPaths perfect_path = analysis.paths();
292 int perfect_match_count = perfect_path.pathz.size();
294 Sequence sub1 = s1.subseq(2, s1.size()-4);
295 Sequence sub2 = s2.subseq(2, s2.size()-4);
297 subanalysis.append_sequence(sub1);
298 subanalysis.append_sequence(sub2);
299 subanalysis.set_threshold(7);
300 subanalysis.set_window(8);
301 subanalysis.analyze();
302 NwayPaths one_mismatch_path = subanalysis.paths();
303 int one_mismatch_count = one_mismatch_path.pathz.size();
305 BOOST_CHECK( perfect_match_count < one_mismatch_count );
308 BOOST_AUTO_TEST_CASE( dirty_flag )
311 BOOST_CHECK_EQUAL(m.is_dirty(), false);
313 BOOST_CHECK_EQUAL(m.is_dirty(), true);
316 BOOST_CHECK_EQUAL(m.is_dirty(), true);
319 BOOST_CHECK_EQUAL(m.is_dirty(), true);
321 m.set_soft_threshold(1);
322 BOOST_CHECK_EQUAL(m.is_dirty(), false);
324 m.append_sequence("AAGGCCTT");
325 BOOST_CHECK_EQUAL(m.is_dirty(), true);