alg/test/test_mussa.cpp

   1 #define BOOST_AUTO_TEST_MAIN
   2 #include <boost/test/auto_unit_test.hpp>
   3 #include <boost/filesystem/path.hpp>
   4 #include <boost/filesystem/operations.hpp>
   5 namespace fs = boost::filesystem;
   6 #include <boost/assign/list_of.hpp>
   7 #include <boost/assign/list_inserter.hpp>
   8 #include <boost/assign.hpp>
   9 namespace assign = boost::assign;
  10
  11 #include <string>
  12 #include <sstream>
  13 #include <vector>
  14
  15 #include "alg/mussa.hpp"
  16 #include "mussa_exceptions.hpp"
  17
  18 using namespace std;
  19
  20 //! can we initialize a mussa object?
  21 BOOST_AUTO_TEST_CASE( mussa_simple )
  22 {
  23   Mussa m;
  24   BOOST_CHECK_EQUAL(m.empty(), true);
  25   BOOST_CHECK_EQUAL(m.get_name(), "" );
  26   BOOST_CHECK_EQUAL(m.get_window(), 0);
  27   BOOST_CHECK_EQUAL(m.get_threshold(), 0);
  28   BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::TransitiveNway);
  29   m.set_name( "hello" );
  30   BOOST_CHECK_EQUAL(m.get_name(), "hello" );
  31   m.set_window(30);
  32   BOOST_CHECK_EQUAL(m.get_window(), 30);
  33   m.set_threshold(21);
  34   BOOST_CHECK_EQUAL(m.get_threshold(), 21);
  35   BOOST_CHECK_EQUAL(m.get_soft_threshold(), 21);
  36   m.set_soft_threshold(19);
  37   BOOST_CHECK_EQUAL(m.get_soft_threshold(), 21);
  38   m.set_soft_threshold(35);
  39   BOOST_CHECK_EQUAL(m.get_soft_threshold(), 30);
  40   m.set_soft_threshold(25);
  41   BOOST_CHECK_EQUAL(m.get_soft_threshold(), 25);
  42   m.set_analysis_mode(Mussa::RadialNway);
  43   BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::RadialNway);
  44   // make sure our path is empty
  45   BOOST_CHECK_EQUAL(m.get_analysis_path().string(), fs::path().string() );
  46
  47   m.clear();
  48   BOOST_CHECK_EQUAL(m.get_name(), "" );
  49   BOOST_CHECK_EQUAL(m.get_window(), 0);
  50   BOOST_CHECK_EQUAL(m.get_threshold(), 0);
  51   BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::TransitiveNway);
  52 }
  53
  54 BOOST_AUTO_TEST_CASE ( mussa_title )
  55 {
  56   Mussa m;
  57
  58   BOOST_CHECK_EQUAL( m.get_title(), "Unnamed");
  59   string foo("foo");
  60   m.set_name(foo);
  61   BOOST_CHECK_EQUAL( m.get_title(), foo);
  62   string foopath_name("/my/silly/path");
  63   fs::path foopath(foopath_name);
  64   m.set_analysis_path(foopath);
  65   BOOST_CHECK_EQUAL( m.get_title().size(), 14);
  66 }
  67
  68 BOOST_AUTO_TEST_CASE( mussa_analysis_name )
  69 {
  70   Mussa m;
  71   m.set_analysis_mode( Mussa::TransitiveNway );
  72   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Transitive" );
  73   m.set_analysis_mode( Mussa::RadialNway );
  74   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Radial" );
  75   m.set_analysis_mode( Mussa::EntropyNway );
  76   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Entropy" );
  77   m.set_analysis_mode( Mussa::RecursiveNway);
  78   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "[deprecated] Recursive" );
  79 }
  80
  81 BOOST_AUTO_TEST_CASE( mussa_sequences )
  82 {
  83   std::string s0("AAAANNNN");
  84   std::string s1("GGGGNNNN");
  85   std::string s2("TTTTNNNN");
  86
  87   Mussa analysis;
  88   BOOST_CHECK_EQUAL(analysis.empty(), true);
  89   analysis.append_sequence(s0);
  90   analysis.append_sequence(s1);
  91   analysis.append_sequence(s2);
  92
  93   BOOST_CHECK_EQUAL( analysis.empty(), false);
  94   BOOST_CHECK_EQUAL( analysis.sequences().size(), 3 );
  95   BOOST_CHECK_EQUAL( *(analysis.sequences()[0]), s0);
  96   BOOST_CHECK_EQUAL( *(analysis.sequences()[1]), s1);
  97   BOOST_CHECK_EQUAL( *(analysis.sequences()[2]), s2);
  98 }
  99
 100 // for some reason we can call nway once safely but it
 101 // somehow changed things and caused a segfault
 102 // fixed by adding a return statement in trans_path_search
 103 BOOST_AUTO_TEST_CASE ( empty_mussa_set_threshold )
 104 {
 105   Mussa m;
 106   m.set_soft_threshold(15);
 107   m.nway();
 108
 109   m.set_soft_threshold(25);
 110   m.nway();
 111 }
 112
 113 BOOST_AUTO_TEST_CASE( mussa_load_mupa )
 114 {
 115   fs::path mupa_path(EXAMPLE_DIR, fs::native);
 116   fs::path result_path = fs::initial_path() / "mck3test_w30_t20";
 117   mupa_path /= "mck3test.mupa";
 118
 119   Mussa m1;
 120   m1.load_mupa_file( mupa_path );
 121   m1.analyze();
 122   m1.save( result_path );
 123   BOOST_CHECK_EQUAL( m1.empty(), false);
 124   BOOST_CHECK_EQUAL( m1.get_name(), std::string("mck3test") );
 125   BOOST_CHECK( m1.size() > 0 );
 126   BOOST_CHECK_EQUAL( m1.get_analysis_path().string(), result_path.string());
 127
 128   Mussa m2;
 129   m2.load( result_path );
 130   BOOST_CHECK_EQUAL( m2.empty(), false);
 131   BOOST_CHECK_EQUAL( m2.get_name(), result_path.leaf() );
 132   BOOST_CHECK_EQUAL( m1.size(), m2.size() );
 133   BOOST_CHECK_EQUAL( result_path.string(), m2.get_analysis_path().string() );
 134
 135   // check clear a bit
 136   m2.clear();
 137   BOOST_CHECK_EQUAL( m2.empty(), true);
 138   BOOST_CHECK_EQUAL( m2.is_dirty(), false );
 139   BOOST_CHECK_EQUAL( m2.get_analysis_path().string(), fs::path().string());
 140 }
 141
 142 BOOST_AUTO_TEST_CASE( mussa_load_full_path )
 143 {
 144   Mussa m1;
 145   fs::path full_path(fs::path(EXAMPLE_DIR, fs::native) / "mck3test.mupa");
 146   m1.load_mupa_file( full_path );
 147   m1.analyze();
 148
 149   BOOST_CHECK( m1.size() > 0);
 150   BOOST_CHECK_EQUAL( m1.get_window(), 30 );
 151   BOOST_CHECK_EQUAL( m1.get_threshold(), 20);
 152   BOOST_CHECK_EQUAL( m1.is_dirty(), true);
 153   BOOST_CHECK_EQUAL( m1.get_analysis_path().string(), "");
 154 }
 155
 156 BOOST_AUTO_TEST_CASE( mussa_valid_motifs_in_new_analysis )
 157 {
 158   Mussa m1;
 159   fs::path full_path(fs::path(EXAMPLE_DIR, fs::native) / "mck3test.mupa");
 160   m1.load_mupa_file( full_path );
 161   m1.analyze();
 162   // check motifs
 163   BOOST_CHECK( m1.sequences().size() > 0 );
 164   BOOST_CHECK_EQUAL( m1.sequences()[0]->motifs().size(), 0 );
 165 }
 166
 167 // make sure we know that mupa files cannot be directories
 168 BOOST_AUTO_TEST_CASE( mussa_mupa_is_file_not_directory )
 169 {
 170   fs::path curdir(".");
 171   Mussa m1;
 172   BOOST_CHECK_THROW(m1.load_mupa_file( curdir ), mussa_load_error );
 173 }
 174
 175 // catch error if annotation isn't a file
 176 BOOST_AUTO_TEST_CASE( mussa_annotation_is_not_file )
 177 {
 178   Mussa m1;
 179   fs::path full_path(fs::path(EXAMPLE_DIR, fs::native) / "directory.mupa");
 180   BOOST_CHECK_THROW( m1.load_mupa_file( full_path ), mussa_load_error );
 181 }
 182
 183 BOOST_AUTO_TEST_CASE( mussa_load_analysis )
 184 {
 185   fs::path example_dir(EXAMPLE_DIR, fs::native);
 186   Mussa m1;
 187   m1.load_mupa_file( example_dir / "mck3test.mupa" );
 188   m1.analyze();
 189
 190   Mussa m2;
 191   fs::path analysis_path = fs::initial_path() / "mck3test_w30_t20";
 192   m2.load( analysis_path );
 193
 194   BOOST_CHECK_EQUAL( m1.size(), m2.size() );
 195   BOOST_CHECK_EQUAL( m1.get_window(), m2.get_window() );
 196   BOOST_CHECK_EQUAL( m1.get_threshold(), m2.get_threshold() );
 197   BOOST_CHECK_EQUAL( m2.get_analysis_path().string(), analysis_path.string());
 198 }
 199
 200 BOOST_AUTO_TEST_CASE( mussa_load_motif )
 201 {
 202   string data = "AAGG 1.0 1.0 0.0\n"
 203                 "GGTT 0.0 0.1 1.0 1.0\n";
 204
 205   istringstream test_istream(data);
 206
 207   Mussa m1;
 208   m1.append_sequence("AAAAGGGGTTTT");
 209   m1.append_sequence("GGGCCCCTTCCAATT");
 210   m1.load_motifs(test_istream);
 211
 212   BOOST_CHECK_EQUAL( m1.motifs().size(), 2);
 213   for (Mussa::vector_sequence_type::const_iterator seq_i = m1.sequences().begin();
 214        seq_i != m1.sequences().end();
 215        ++seq_i)
 216   {
 217     BOOST_CHECK( (*seq_i)->motifs().size() > 0 );
 218   }
 219 }
 220
 221 BOOST_AUTO_TEST_CASE( mussa_load_broken_motif )
 222 {
 223   string data = "AAGG 1.0 1.0 0.0\n"
 224                 "GGTT 0.0 0.1 1.0 1.0\n"
 225                 "ZZCTA 0.1 0.0 1.0\n";
 226
 227   istringstream test_istream(data);
 228
 229   Mussa m1;
 230   m1.append_sequence("AAAAGGGGTTTT");
 231   m1.append_sequence("GGGCCCCTTCCAATT");
 232   BOOST_CHECK_THROW(m1.load_motifs(test_istream), motif_load_error);
 233
 234   BOOST_CHECK_EQUAL( m1.motifs().size(), 0);
 235 }
 236
 237 BOOST_AUTO_TEST_CASE( mussa_named_motif )
 238 {
 239   string data = "CCAATT cat 0.1 0.2 0.3\n";
 240   istringstream test_istream(data);
 241
 242   Mussa m1;
 243   m1.append_sequence("AAAAGGGGTTTT");
 244   m1.append_sequence("GGGCCCCTTCCAATT");
 245   m1.load_motifs(test_istream);
 246
 247   std::set<Sequence> motifs = m1.motifs();
 248   BOOST_REQUIRE_EQUAL(motifs.size(), 1);
 249   BOOST_CHECK_EQUAL(motifs.begin()->get_name(), "cat");
 250 }
 251
 252 BOOST_AUTO_TEST_CASE( mussa_weirdly_spaced_named_motif )
 253 {
 254   string data = "CCAATT       cat_meow123     0.1    0.2 0.3\n";
 255   istringstream test_istream(data);
 256
 257   Mussa m1;
 258   m1.append_sequence("AAAAGGGGTTTT");
 259   m1.append_sequence("GGGCCCCTTCCAATT");
 260   m1.load_motifs(test_istream);
 261
 262   std::set<Sequence> motifs = m1.motifs();
 263   BOOST_REQUIRE_EQUAL(motifs.size(), 1);
 264   BOOST_CHECK_EQUAL(motifs.begin()->get_name(), "cat_meow123");
 265 }
 266
 267 BOOST_AUTO_TEST_CASE( mussa_name_quoted_motif )
 268 {
 269   string data = "CCAATT       \"cat meow 123\"     0.1    0.2 0.3\n";
 270   istringstream test_istream(data);
 271
 272   Mussa m1;
 273   m1.append_sequence("AAAAGGGGTTTT");
 274   m1.append_sequence("GGGCCCCTTCCAATT");
 275   m1.load_motifs(test_istream);
 276
 277   std::set<Sequence> motifs = m1.motifs();
 278   BOOST_REQUIRE_EQUAL(motifs.size(), 1);
 279   BOOST_CHECK_EQUAL(motifs.begin()->get_name(), "cat meow 123");
 280 }
 281
 282 BOOST_AUTO_TEST_CASE( mussa_name_embedded_quote_motif )
 283 {
 284   // pretty obviously this shouldn't work as " are our delimiter
 285   // and i'm too lazy to add support for \ in the parser
 286   string data = "ATA 0.5 0.5 0.5\n"
 287                 "CCAATT       \"cat \"meow 123\"     0.1    0.2 0.3\n";
 288   istringstream test_istream(data);
 289
 290   Mussa m1;
 291   m1.append_sequence("AAAAGGGGTTTT");
 292   m1.append_sequence("GGGCCCCTTCCAATT");
 293   BOOST_CHECK_THROW( m1.load_motifs(test_istream), motif_load_error);
 294
 295   std::set<Sequence> motifs = m1.motifs();
 296   BOOST_REQUIRE_EQUAL(motifs.size(), 0);
 297 }
 298
 299 BOOST_AUTO_TEST_CASE( mussa_save_motif )
 300 {
 301   string data = "ATA 1 1 1 1\n"
 302                 "CAT \"my name\" 1 0 0.5 0.5\n";
 303   istringstream data_istream(data);
 304
 305   Mussa m1;
 306   m1.append_sequence("AAAAGGGGTTTT");
 307   m1.append_sequence("GGGCCCCTTCCAATT");
 308   m1.load_motifs(data_istream);
 309
 310   string save;
 311   ostringstream save_ostream(save);
 312   m1.save_motifs(save_ostream);
 313
 314   istringstream reloaded_istream(save_ostream.str());
 315   Mussa m2;
 316   m2.append_sequence("AAAAGGGGTTTT");
 317   m2.append_sequence("GGGCCCCTTCCAATT");
 318   m2.load_motifs(reloaded_istream);
 319
 320   BOOST_REQUIRE_EQUAL(m1.motifs().size(), m2.motifs().size());
 321   Mussa::motif_set::const_iterator m1motif = m1.motifs().begin();
 322   Mussa::motif_set::const_iterator m2motif = m2.motifs().begin();
 323   for (;
 324        m1motif != m1.motifs().end() and m2motif != m2.motifs().end();
 325        ++m1motif, ++m2motif)
 326   {
 327     BOOST_CHECK_EQUAL(m1motif->get_sequence(), m2motif->get_sequence());
 328     BOOST_CHECK_EQUAL(m1motif->get_name(), m2motif->get_name());
 329     BOOST_CHECK_EQUAL(m1.colorMapper()->lookup("motif", m1motif->get_sequence()),
 330                       m2.colorMapper()->lookup("motif", m2motif->get_sequence()));
 331   }
 332 }
 333
 334 BOOST_AUTO_TEST_CASE( mussa_add_motif )
 335 {
 336   vector<Sequence> motifs;
 337   motifs.push_back("AAGG");
 338   vector<Color> colors;
 339   colors.push_back(Color(1.0, 0.0, 0.0));
 340
 341   Mussa m1;
 342   m1.append_sequence("AAAAGGGGTTTT");
 343   m1.append_sequence("GGGCCCCTTGGTT");
 344   m1.set_motifs(motifs, colors);
 345   int first_size = m1.motifs().size();
 346   BOOST_CHECK_EQUAL( first_size, 1 );
 347   BOOST_REQUIRE(first_size > 0);
 348   BOOST_CHECK_EQUAL(*(m1.motifs().begin()), motifs.front());
 349   // make sure that our sequences have the right number of motifs
 350   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 1);
 351   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1); // because of rc
 352
 353   // verify that setting the motif clears the arrays
 354   m1.set_motifs(motifs, colors);
 355   BOOST_CHECK_EQUAL( first_size, m1.motifs().size() );
 356   // make sure that our sequences have the right number of motifs
 357   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 1);
 358   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
 359
 360   // add a different motif
 361   motifs.clear();
 362   motifs.push_back("CCTTGG");
 363   BOOST_CHECK_EQUAL(motifs.size(), 1);
 364   m1.set_motifs(motifs, colors);
 365   BOOST_CHECK_EQUAL(m1.motifs().size(), 1);
 366   BOOST_REQUIRE(m1.motifs().size() > 0);
 367   BOOST_CHECK_EQUAL(*(m1.motifs().begin()), motifs.front());
 368   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 0);
 369   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
 370
 371   // try a motif that doesn't exist
 372   motifs.clear();
 373   motifs.push_back("CCTTGG");
 374   BOOST_CHECK_EQUAL(motifs.size(), 1);
 375   m1.set_motifs(motifs, colors);
 376   BOOST_CHECK_EQUAL(m1.motifs().size(), 1);
 377   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 0);
 378   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
 379
 380 }
 381
 382 static void
 383 two_way_local_align_test(const Mussa::vector_sequence_type &seqs,
 384                          const list<ConservedPath::path_type>& result,
 385                          const list<vector<bool> >& reversed)
 386 {
 387   map<char, vector <char> >  m;
 388   assign::insert(m)('A', assign::list_of('A')('T') )
 389                    ('T', assign::list_of('T')('A') )
 390                    ('G', assign::list_of('G')('C') )
 391                    ('C', assign::list_of('C')('G') );
 392   list<vector<bool> >::const_iterator rc_i = reversed.begin();
 393
 394   for(list<ConservedPath::path_type>::const_iterator base_i = result.begin();
 395       base_i != result.end();
 396       ++base_i, ++rc_i)
 397   {
 398     // since the reverse compliment flag is relative to the first sequence
 399     // the first one should always be false
 400     BOOST_CHECK_EQUAL( (*rc_i)[0], false );
 401     const int first_path_basepair_index = (*base_i)[0];
 402     const int second_path_basepair_index = (*base_i)[1];
 403     const char first_basepair = (*seqs[0])[first_path_basepair_index];
 404     const char second_basepair = (*seqs[1])[second_path_basepair_index];
 405     // get our index into our reverse compliment map m
 406     const int second_compliment_index = (*rc_i)[1];
 407     // lookup the forward or reverse compliment depending on our rc flag
 408     const char complimented_second = m[second_basepair][second_compliment_index];
 409
 410     BOOST_CHECK_EQUAL( first_basepair, complimented_second) ;
 411   }
 412 }
 413
 414 BOOST_AUTO_TEST_CASE( two_way_local_alignment )
 415 {
 416   string s0("GCGCATAT");
 417   string s1("AAAAAAAT");
 418   Sequence seq1(s1);
 419
 420   Mussa analysis;
 421   analysis.append_sequence(s0);
 422   analysis.append_sequence(s1);
 423   analysis.set_threshold(3);
 424   analysis.set_window(4);
 425   analysis.analyze();
 426   NwayPaths npath = analysis.paths();
 427   BOOST_REQUIRE_EQUAL( npath.pathz.size(), 2 );
 428
 429   list<ConservedPath::path_type> result;
 430   list<vector<bool> > reversed;
 431   list<ConservedPath>::iterator pathz_i = npath.pathz.begin();
 432
 433   list<ConservedPath> selected_paths;
 434   selected_paths.push_back(*pathz_i);
 435   analysis.createLocalAlignment(selected_paths.begin(),
 436                                 selected_paths.end(),
 437                                 result,
 438                                 reversed);
 439
 440   two_way_local_align_test(analysis.sequences(), result, reversed);
 441
 442   ++pathz_i;
 443   result.clear();
 444   reversed.clear();
 445   selected_paths.clear();
 446   selected_paths.push_back(*pathz_i);
 447   analysis.createLocalAlignment(selected_paths.begin(),
 448                                 selected_paths.end(),
 449                                 result,
 450                                 reversed);
 451   two_way_local_align_test(analysis.sequences(), result, reversed);
 452 }
 453
 454 BOOST_AUTO_TEST_CASE( three_way_local_alignment )
 455 {
 456   string s0("AGCAGGGAGGGTTTAAATGGCACCCAGCAGTTGGTGTGAGG");
 457   string s1("AGCGGGAAGGGTTTAAATGGCACCGGGCAGTTGGCGTGAGG");
 458   string s2("CAGCGCCGGGGTTTAAATGGCACCGAGCAGTTGGCGCAGGG");
 459
 460   Mussa analysis;
 461   analysis.append_sequence(s0);
 462   analysis.append_sequence(s1);
 463   analysis.append_sequence(s2);
 464   analysis.set_threshold(23);
 465   analysis.set_window(30);
 466   analysis.analyze();
 467   NwayPaths npath = analysis.paths();
 468   BOOST_CHECK_EQUAL( npath.refined_pathz.size(), 1 );
 469
 470   list<ConservedPath::path_type> result;
 471   list<vector<bool> > reversed;
 472   // grab 1 path (since there's only one)
 473   list<ConservedPath>::iterator pathz_i = npath.pathz.begin();
 474   list<ConservedPath> selected_paths;
 475   selected_paths.push_back(*pathz_i);
 476   analysis.createLocalAlignment(selected_paths.begin(),
 477                                 selected_paths.end(),
 478                                 result,
 479                                 reversed);
 480
 481   for(std::list<ConservedPath::path_type>::iterator result_i = result.begin();
 482       result_i != result.end();
 483       ++result_i)
 484   {
 485     ConservedPath::path_element first_element = *(result_i->begin());
 486     for (ConservedPath::path_type::iterator element_i = result_i->begin();
 487          element_i != result_i->end();
 488          ++element_i)
 489     {
 490       BOOST_CHECK_EQUAL( *element_i, first_element );
 491       BOOST_CHECK_EQUAL( s0[*element_i], s1[*element_i] );
 492       BOOST_CHECK_EQUAL( s1[*element_i], s2[*element_i] );
 493       BOOST_CHECK_EQUAL( s0[*element_i], s2[*element_i] );
 494     }
 495   }
 496 }
 497
 498 BOOST_AUTO_TEST_CASE( mussa_window_larger_than_sequence )
 499 {
 500   string s0("AGCAGGG");
 501   string s1("CAGCGGG");
 502
 503   Mussa analysis;
 504   analysis.append_sequence(s0);
 505   analysis.append_sequence(s1);
 506   analysis.set_threshold(23);
 507   analysis.set_window(30);
 508   BOOST_CHECK_THROW(analysis.analyze(), seqcomp_error);
 509 }
 510
 511 BOOST_AUTO_TEST_CASE( subanalysis )
 512 {
 513   Sequence s1("AATGAAGATTTTAATGCTTTAATTTTGTTTTGTAAACTTCGAATTTCCAAAATTTGAAA");
 514   Sequence s2("AGGAGCAAGTTCGCTTCATCGAGAATTTTTAATTTTTAGTCAAATTTTCCAATGTCTGA");
 515
 516   Mussa analysis;
 517   analysis.append_sequence(s1);
 518   analysis.append_sequence(s2);
 519   analysis.set_threshold(8);
 520   analysis.set_window(8);
 521   analysis.analyze();
 522
 523   NwayPaths perfect_path = analysis.paths();
 524   int perfect_match_count = perfect_path.pathz.size();
 525
 526   Sequence sub1 = s1.subseq(2, s1.size()-4);
 527   Sequence sub2 = s2.subseq(2, s2.size()-4);
 528   Mussa subanalysis;
 529   subanalysis.append_sequence(sub1);
 530   subanalysis.append_sequence(sub2);
 531   subanalysis.set_threshold(7);
 532   subanalysis.set_window(8);
 533   subanalysis.analyze();
 534   NwayPaths one_mismatch_path = subanalysis.paths();
 535   int one_mismatch_count = one_mismatch_path.pathz.size();
 536
 537   BOOST_CHECK( perfect_match_count < one_mismatch_count );
 538 }
 539
 540 BOOST_AUTO_TEST_CASE( dirty_flag )
 541 {
 542   Mussa m;
 543   BOOST_CHECK_EQUAL(m.is_dirty(), false);
 544   m.set_name("foo");
 545   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 546   m.clear();
 547   m.set_window(30);
 548   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 549   m.clear();
 550   m.set_threshold(1);
 551   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 552   m.clear();
 553   m.set_soft_threshold(1);
 554   BOOST_CHECK_EQUAL(m.is_dirty(), false);
 555   m.clear();
 556   m.append_sequence("AAGGCCTT");
 557   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 558 }
 559