+
+static void
+two_way_local_align_test(const Mussa::vector_sequence_type &seqs,
+ const list<ConservedPath::path_type>& result,
+ const list<vector<bool> >& reversed)
+{
+ map<char, vector <char> > m;
+ assign::insert(m)('A', assign::list_of('A')('T') )
+ ('T', assign::list_of('T')('A') )
+ ('G', assign::list_of('G')('C') )
+ ('C', assign::list_of('C')('G') );
+ list<vector<bool> >::const_iterator rc_i = reversed.begin();
+
+ for(list<ConservedPath::path_type>::const_iterator base_i = result.begin();
+ base_i != result.end();
+ ++base_i, ++rc_i)
+ {
+ // since the reverse compliment flag is relative to the first sequence
+ // the first one should always be false
+ BOOST_CHECK_EQUAL( (*rc_i)[0], false );
+ const int first_path_basepair_index = (*base_i)[0];
+ const int second_path_basepair_index = (*base_i)[1];
+ const char first_basepair = (*seqs[0])[first_path_basepair_index];
+ const char second_basepair = (*seqs[1])[second_path_basepair_index];
+ // get our index into our reverse compliment map m
+ const int second_compliment_index = (*rc_i)[1];
+ // lookup the forward or reverse compliment depending on our rc flag
+ const char complimented_second = m[second_basepair][second_compliment_index];
+
+ BOOST_CHECK_EQUAL( first_basepair, complimented_second) ;
+ }
+}
+
+BOOST_AUTO_TEST_CASE( two_way_local_alignment )
+{
+ string s0("GCGCATAT");
+ string s1("AAAAAAAT");
+ Sequence seq1(s1);
+
+ Mussa analysis;
+ analysis.append_sequence(s0);
+ analysis.append_sequence(s1);
+ analysis.set_threshold(3);
+ analysis.set_window(4);
+ analysis.analyze();
+ NwayPaths npath = analysis.paths();
+ BOOST_REQUIRE_EQUAL( npath.pathz.size(), 2 );
+
+ list<ConservedPath::path_type> result;
+ list<vector<bool> > reversed;
+ list<ConservedPath>::iterator pathz_i = npath.pathz.begin();
+
+ list<ConservedPath> selected_paths;
+ selected_paths.push_back(*pathz_i);
+ analysis.createLocalAlignment(selected_paths.begin(),
+ selected_paths.end(),
+ result,
+ reversed);
+
+ two_way_local_align_test(analysis.sequences(), result, reversed);
+
+ ++pathz_i;
+ result.clear();
+ reversed.clear();
+ selected_paths.clear();
+ selected_paths.push_back(*pathz_i);
+ analysis.createLocalAlignment(selected_paths.begin(),
+ selected_paths.end(),
+ result,
+ reversed);
+ two_way_local_align_test(analysis.sequences(), result, reversed);
+}
+
+BOOST_AUTO_TEST_CASE( three_way_local_alignment )
+{
+ string s0("AGCAGGGAGGGTTTAAATGGCACCCAGCAGTTGGTGTGAGG");
+ string s1("AGCGGGAAGGGTTTAAATGGCACCGGGCAGTTGGCGTGAGG");
+ string s2("CAGCGCCGGGGTTTAAATGGCACCGAGCAGTTGGCGCAGGG");
+
+ Mussa analysis;
+ analysis.append_sequence(s0);
+ analysis.append_sequence(s1);
+ analysis.append_sequence(s2);
+ analysis.set_threshold(23);
+ analysis.set_window(30);
+ analysis.analyze();
+ NwayPaths npath = analysis.paths();
+ BOOST_CHECK_EQUAL( npath.refined_pathz.size(), 1 );
+
+ list<ConservedPath::path_type> result;
+ list<vector<bool> > reversed;
+ // grab 1 path (since there's only one)
+ list<ConservedPath>::iterator pathz_i = npath.pathz.begin();
+ list<ConservedPath> selected_paths;
+ selected_paths.push_back(*pathz_i);
+ analysis.createLocalAlignment(selected_paths.begin(),
+ selected_paths.end(),
+ result,
+ reversed);
+
+ for(std::list<ConservedPath::path_type>::iterator result_i = result.begin();
+ result_i != result.end();
+ ++result_i)
+ {
+ ConservedPath::path_element first_element = *(result_i->begin());
+ for (ConservedPath::path_type::iterator element_i = result_i->begin();
+ element_i != result_i->end();
+ ++element_i)
+ {
+ BOOST_CHECK_EQUAL( *element_i, first_element );
+ BOOST_CHECK_EQUAL( s0[*element_i], s1[*element_i] );
+ BOOST_CHECK_EQUAL( s1[*element_i], s2[*element_i] );
+ BOOST_CHECK_EQUAL( s0[*element_i], s2[*element_i] );
+ }
+ }
+}
+
+BOOST_AUTO_TEST_CASE( mussa_window_larger_than_sequence )
+{
+ string s0("AGCAGGG");
+ string s1("CAGCGGG");
+
+ Mussa analysis;
+ analysis.append_sequence(s0);
+ analysis.append_sequence(s1);
+ analysis.set_threshold(23);
+ analysis.set_window(30);
+ BOOST_CHECK_THROW(analysis.analyze(), seqcomp_error);
+}
+
+BOOST_AUTO_TEST_CASE( subanalysis )
+{
+ Sequence s1("AATGAAGATTTTAATGCTTTAATTTTGTTTTGTAAACTTCGAATTTCCAAAATTTGAAA");
+ Sequence s2("AGGAGCAAGTTCGCTTCATCGAGAATTTTTAATTTTTAGTCAAATTTTCCAATGTCTGA");
+
+ Mussa analysis;
+ analysis.append_sequence(s1);
+ analysis.append_sequence(s2);
+ analysis.set_threshold(8);
+ analysis.set_window(8);
+ analysis.analyze();
+
+ NwayPaths perfect_path = analysis.paths();
+ int perfect_match_count = perfect_path.pathz.size();
+
+ Sequence sub1 = s1.subseq(2, s1.size()-4);
+ Sequence sub2 = s2.subseq(2, s2.size()-4);
+ Mussa subanalysis;
+ subanalysis.append_sequence(sub1);
+ subanalysis.append_sequence(sub2);
+ subanalysis.set_threshold(7);
+ subanalysis.set_window(8);
+ subanalysis.analyze();
+ NwayPaths one_mismatch_path = subanalysis.paths();
+ int one_mismatch_count = one_mismatch_path.pathz.size();
+
+ BOOST_CHECK( perfect_match_count < one_mismatch_count );
+}
+
+BOOST_AUTO_TEST_CASE( dirty_flag )
+{
+ Mussa m;
+ BOOST_CHECK_EQUAL(m.is_dirty(), false);
+ m.set_name("foo");
+ BOOST_CHECK_EQUAL(m.is_dirty(), true);
+ m.clear();
+ m.set_window(30);
+ BOOST_CHECK_EQUAL(m.is_dirty(), true);
+ m.clear();
+ m.set_threshold(1);
+ BOOST_CHECK_EQUAL(m.is_dirty(), true);
+ m.clear();
+ m.set_soft_threshold(1);
+ BOOST_CHECK_EQUAL(m.is_dirty(), false);
+ m.clear();
+ m.append_sequence("AAGGCCTT");
+ BOOST_CHECK_EQUAL(m.is_dirty(), true);
+}
+