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_crlf )
 114 {
 115   fs::path example_path(EXAMPLE_DIR, fs::native);
 116   fs::path seq_path(example_path / "seq" / "mouse_mck_pro.fa");
 117   fs::path annot_path(example_path / "mm_mck3test.annot");
 118
 119   std::string mupa(
 120     "# hello\015\012"
 121     "ANA_NAME load_mupa_crlf\015\012");
 122   mupa += "SEQUENCE " + seq_path.native_file_string() + "\015\012";
 123   mupa += "ANNOTATION " + annot_path.native_file_string() + "\015\012";
 124
 125   istringstream mupa_stream(mupa);
 126   Mussa m;
 127   fs::path base;
 128   m.load_mupa_stream( mupa_stream, base );
 129   // Should run with no exceptions
 130 }
 131
 132 BOOST_AUTO_TEST_CASE( mussa_load_mupa_comment_character )
 133 {
 134   fs::path mupa_path(EXAMPLE_DIR, fs::native);
 135   fs::path seq_path = fs::initial_path() / "seq" / "mouse_mck_pro.fa";
 136   fs::path annot_path = fs::initial_path() / "mm_mck3test.annot";
 137
 138   std::string mupa(
 139     "# hello\015\012"
 140     "ANA_NAME load_mupa_crlf\015\012");
 141   mupa += "#SEQUENCE " + seq_path.native_file_string() + "\015\012";
 142   mupa += "#ANNOTATION " + annot_path.native_file_string() + "\015\012";
 143
 144   istringstream mupa_stream(mupa);
 145   Mussa m;
 146   fs::path base;
 147   m.load_mupa_stream( mupa_stream, base );
 148   // Should run with no exceptions
 149 }
 150
 151 BOOST_AUTO_TEST_CASE( mussa_load_mupa_exception )
 152 {
 153   std::string mupa(
 154     "# hello\015\012"
 155     "ANA_NAME load_mupa_crlf\015\012"
 156     "mwahhaha I broke you!\n"
 157   );
 158
 159   istringstream mupa_stream(mupa);
 160   Mussa m;
 161   fs::path base;
 162   BOOST_CHECK_THROW(m.load_mupa_stream( mupa_stream, base ), mussa_load_error);
 163 }
 164
 165 BOOST_AUTO_TEST_CASE( mussa_load_mupa )
 166 {
 167   fs::path mupa_path(EXAMPLE_DIR, fs::native);
 168   fs::path result_path = fs::initial_path() / "mck3test_w30_t20";
 169   mupa_path /= "mck3test.mupa";
 170
 171   Mussa m1;
 172   m1.load_mupa_file( mupa_path );
 173   m1.analyze();
 174   m1.save( result_path );
 175   BOOST_CHECK_EQUAL( m1.empty(), false);
 176   BOOST_CHECK_EQUAL( m1.get_name(), std::string("mck3test") );
 177   BOOST_CHECK( m1.size() > 0 );
 178   BOOST_CHECK_EQUAL( m1.get_analysis_path().string(), result_path.string());
 179
 180   Mussa m2;
 181   m2.load( result_path );
 182   BOOST_CHECK_EQUAL( m2.empty(), false);
 183   BOOST_CHECK_EQUAL( m2.get_name(), result_path.leaf() );
 184   BOOST_CHECK_EQUAL( m1.size(), m2.size() );
 185   BOOST_CHECK_EQUAL( result_path.string(), m2.get_analysis_path().string() );
 186
 187   // check clear a bit
 188   m2.clear();
 189   BOOST_CHECK_EQUAL( m2.empty(), true);
 190   BOOST_CHECK_EQUAL( m2.is_dirty(), false );
 191   BOOST_CHECK_EQUAL( m2.get_analysis_path().string(), fs::path().string());
 192 }
 193
 194 BOOST_AUTO_TEST_CASE( mussa_load_full_path )
 195 {
 196   Mussa m1;
 197   fs::path full_path(fs::path(EXAMPLE_DIR, fs::native) / "mck3test.mupa");
 198   m1.load_mupa_file( full_path );
 199   m1.analyze();
 200
 201   BOOST_CHECK( m1.size() > 0);
 202   BOOST_CHECK_EQUAL( m1.get_window(), 30 );
 203   BOOST_CHECK_EQUAL( m1.get_threshold(), 20);
 204   BOOST_CHECK_EQUAL( m1.is_dirty(), true);
 205   BOOST_CHECK_EQUAL( m1.get_analysis_path().string(), "");
 206 }
 207
 208 BOOST_AUTO_TEST_CASE( mussa_valid_motifs_in_new_analysis )
 209 {
 210   Mussa m1;
 211   fs::path full_path(fs::path(EXAMPLE_DIR, fs::native) / "mck3test.mupa");
 212   m1.load_mupa_file( full_path );
 213   m1.analyze();
 214   // check motifs
 215   BOOST_CHECK( m1.sequences().size() > 0 );
 216   BOOST_CHECK_EQUAL( m1.sequences()[0]->motifs().size(), 0 );
 217 }
 218
 219 // make sure we know that mupa files cannot be directories
 220 BOOST_AUTO_TEST_CASE( mussa_mupa_is_file_not_directory )
 221 {
 222   fs::path curdir(".");
 223   Mussa m1;
 224   BOOST_CHECK_THROW(m1.load_mupa_file( curdir ), mussa_load_error );
 225 }
 226
 227 // catch error if annotation isn't a file
 228 BOOST_AUTO_TEST_CASE( mussa_annotation_is_not_file )
 229 {
 230   Mussa m1;
 231   fs::path full_path(fs::path(EXAMPLE_DIR, fs::native) / "directory.mupa");
 232   BOOST_CHECK_THROW( m1.load_mupa_file( full_path ), mussa_load_error );
 233 }
 234
 235 BOOST_AUTO_TEST_CASE( mussa_load_analysis )
 236 {
 237   fs::path example_dir(EXAMPLE_DIR, fs::native);
 238   Mussa m1;
 239   m1.load_mupa_file( example_dir / "mck3test.mupa" );
 240   m1.analyze();
 241
 242   Mussa m2;
 243   fs::path analysis_path = fs::initial_path() / "mck3test_w30_t20";
 244   m2.load( analysis_path );
 245
 246   BOOST_CHECK_EQUAL( m1.size(), m2.size() );
 247   BOOST_CHECK_EQUAL( m1.get_window(), m2.get_window() );
 248   BOOST_CHECK_EQUAL( m1.get_threshold(), m2.get_threshold() );
 249   BOOST_CHECK_EQUAL( m2.get_analysis_path().string(), analysis_path.string());
 250 }
 251
 252 BOOST_AUTO_TEST_CASE( mussa_load_motif )
 253 {
 254   string data = "AAGG 1.0 1.0 0.0\n"
 255                 "GGTT 0.0 0.1 1.0 1.0\n";
 256
 257   istringstream test_istream(data);
 258
 259   Mussa m1;
 260   m1.append_sequence("AAAAGGGGTTTT");
 261   m1.append_sequence("GGGCCCCTTCCAATT");
 262   m1.load_motifs(test_istream);
 263
 264   BOOST_CHECK_EQUAL( m1.motifs().size(), 2);
 265   for (Mussa::vector_sequence_type::const_iterator seq_i = m1.sequences().begin();
 266        seq_i != m1.sequences().end();
 267        ++seq_i)
 268   {
 269     BOOST_CHECK( (*seq_i)->motifs().size() > 0 );
 270   }
 271 }
 272
 273 BOOST_AUTO_TEST_CASE( mussa_load_broken_motif )
 274 {
 275   string data = "AAGG 1.0 1.0 0.0\n"
 276                 "GGTT 0.0 0.1 1.0 1.0\n"
 277                 "ZZCTA 0.1 0.0 1.0\n";
 278
 279   istringstream test_istream(data);
 280
 281   Mussa m1;
 282   m1.append_sequence("AAAAGGGGTTTT");
 283   m1.append_sequence("GGGCCCCTTCCAATT");
 284   BOOST_CHECK_THROW(m1.load_motifs(test_istream), motif_load_error);
 285
 286   BOOST_CHECK_EQUAL( m1.motifs().size(), 0);
 287 }
 288
 289 BOOST_AUTO_TEST_CASE( mussa_named_motif )
 290 {
 291   string data = "CCAATT cat 0.1 0.2 0.3\n";
 292   istringstream test_istream(data);
 293
 294   Mussa m1;
 295   m1.append_sequence("AAAAGGGGTTTT");
 296   m1.append_sequence("GGGCCCCTTCCAATT");
 297   m1.load_motifs(test_istream);
 298
 299   std::set<Sequence> motifs = m1.motifs();
 300   BOOST_REQUIRE_EQUAL(motifs.size(), 1);
 301   BOOST_CHECK_EQUAL(motifs.begin()->get_name(), "cat");
 302 }
 303
 304 BOOST_AUTO_TEST_CASE( mussa_weirdly_spaced_named_motif )
 305 {
 306   string data = "CCAATT       cat_meow123     0.1    0.2 0.3\n";
 307   istringstream test_istream(data);
 308
 309   Mussa m1;
 310   m1.append_sequence("AAAAGGGGTTTT");
 311   m1.append_sequence("GGGCCCCTTCCAATT");
 312   m1.load_motifs(test_istream);
 313
 314   std::set<Sequence> motifs = m1.motifs();
 315   BOOST_REQUIRE_EQUAL(motifs.size(), 1);
 316   BOOST_CHECK_EQUAL(motifs.begin()->get_name(), "cat_meow123");
 317 }
 318
 319 BOOST_AUTO_TEST_CASE( mussa_name_quoted_motif )
 320 {
 321   string data = "CCAATT       \"cat meow 123\"     0.1    0.2 0.3\n";
 322   istringstream test_istream(data);
 323
 324   Mussa m1;
 325   m1.append_sequence("AAAAGGGGTTTT");
 326   m1.append_sequence("GGGCCCCTTCCAATT");
 327   m1.load_motifs(test_istream);
 328
 329   std::set<Sequence> motifs = m1.motifs();
 330   BOOST_REQUIRE_EQUAL(motifs.size(), 1);
 331   BOOST_CHECK_EQUAL(motifs.begin()->get_name(), "cat meow 123");
 332 }
 333
 334 BOOST_AUTO_TEST_CASE( mussa_name_embedded_quote_motif )
 335 {
 336   // pretty obviously this shouldn't work as " are our delimiter
 337   // and i'm too lazy to add support for \ in the parser
 338   string data = "ATA 0.5 0.5 0.5\n"
 339                 "CCAATT       \"cat \"meow 123\"     0.1    0.2 0.3\n";
 340   istringstream test_istream(data);
 341
 342   Mussa m1;
 343   m1.append_sequence("AAAAGGGGTTTT");
 344   m1.append_sequence("GGGCCCCTTCCAATT");
 345   BOOST_CHECK_THROW( m1.load_motifs(test_istream), motif_load_error);
 346
 347   std::set<Sequence> motifs = m1.motifs();
 348   BOOST_REQUIRE_EQUAL(motifs.size(), 0);
 349 }
 350
 351 BOOST_AUTO_TEST_CASE( mussa_save_motif )
 352 {
 353   string data = "ATA 1 1 1 1\n"
 354                 "CAT \"my name\" 1 0 0.5 0.5\n";
 355   istringstream data_istream(data);
 356
 357   Mussa m1;
 358   m1.append_sequence("AAAAGGGGTTTT");
 359   m1.append_sequence("GGGCCCCTTCCAATT");
 360   m1.load_motifs(data_istream);
 361
 362   string save;
 363   ostringstream save_ostream(save);
 364   m1.save_motifs(save_ostream);
 365
 366   istringstream reloaded_istream(save_ostream.str());
 367   Mussa m2;
 368   m2.append_sequence("AAAAGGGGTTTT");
 369   m2.append_sequence("GGGCCCCTTCCAATT");
 370   m2.load_motifs(reloaded_istream);
 371
 372   BOOST_REQUIRE_EQUAL(m1.motifs().size(), m2.motifs().size());
 373   Mussa::motif_set::const_iterator m1motif = m1.motifs().begin();
 374   Mussa::motif_set::const_iterator m2motif = m2.motifs().begin();
 375   for (;
 376        m1motif != m1.motifs().end() and m2motif != m2.motifs().end();
 377        ++m1motif, ++m2motif)
 378   {
 379     BOOST_CHECK_EQUAL(m1motif->get_sequence(), m2motif->get_sequence());
 380     BOOST_CHECK_EQUAL(m1motif->get_name(), m2motif->get_name());
 381     BOOST_CHECK_EQUAL(m1.colorMapper()->lookup("motif", m1motif->get_sequence()),
 382                       m2.colorMapper()->lookup("motif", m2motif->get_sequence()));
 383   }
 384 }
 385
 386 BOOST_AUTO_TEST_CASE( mussa_add_motif )
 387 {
 388   vector<Sequence> motifs;
 389   motifs.push_back("AAGG");
 390   vector<Color> colors;
 391   colors.push_back(Color(1.0, 0.0, 0.0));
 392
 393   Mussa m1;
 394   m1.append_sequence("AAAAGGGGTTTT");
 395   m1.append_sequence("GGGCCCCTTGGTT");
 396   m1.set_motifs(motifs, colors);
 397   int first_size = m1.motifs().size();
 398   BOOST_CHECK_EQUAL( first_size, 1 );
 399   BOOST_REQUIRE(first_size > 0);
 400   BOOST_CHECK_EQUAL(*(m1.motifs().begin()), motifs.front());
 401   // make sure that our sequences have the right number of motifs
 402   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 1);
 403   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1); // because of rc
 404
 405   // verify that setting the motif clears the arrays
 406   m1.set_motifs(motifs, colors);
 407   BOOST_CHECK_EQUAL( first_size, m1.motifs().size() );
 408   // make sure that our sequences have the right number of motifs
 409   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 1);
 410   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
 411
 412   // add a different motif
 413   motifs.clear();
 414   motifs.push_back("CCTTGG");
 415   BOOST_CHECK_EQUAL(motifs.size(), 1);
 416   m1.set_motifs(motifs, colors);
 417   BOOST_CHECK_EQUAL(m1.motifs().size(), 1);
 418   BOOST_REQUIRE(m1.motifs().size() > 0);
 419   BOOST_CHECK_EQUAL(*(m1.motifs().begin()), motifs.front());
 420   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 0);
 421   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
 422
 423   // try a motif that doesn't exist
 424   motifs.clear();
 425   motifs.push_back("CCTTGG");
 426   BOOST_CHECK_EQUAL(motifs.size(), 1);
 427   m1.set_motifs(motifs, colors);
 428   BOOST_CHECK_EQUAL(m1.motifs().size(), 1);
 429   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 0);
 430   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
 431
 432 }
 433
 434 static void
 435 two_way_local_align_test(const Mussa::vector_sequence_type &seqs,
 436                          const list<ConservedPath::path_type>& result,
 437                          const list<vector<bool> >& reversed)
 438 {
 439   map<char, vector <char> >  m;
 440   assign::insert(m)('A', assign::list_of('A')('T') )
 441                    ('T', assign::list_of('T')('A') )
 442                    ('G', assign::list_of('G')('C') )
 443                    ('C', assign::list_of('C')('G') );
 444   list<vector<bool> >::const_iterator rc_i = reversed.begin();
 445
 446   for(list<ConservedPath::path_type>::const_iterator base_i = result.begin();
 447       base_i != result.end();
 448       ++base_i, ++rc_i)
 449   {
 450     // since the reverse compliment flag is relative to the first sequence
 451     // the first one should always be false
 452     BOOST_CHECK_EQUAL( (*rc_i)[0], false );
 453     const int first_path_basepair_index = (*base_i)[0];
 454     const int second_path_basepair_index = (*base_i)[1];
 455     const char first_basepair = (*seqs[0])[first_path_basepair_index];
 456     const char second_basepair = (*seqs[1])[second_path_basepair_index];
 457     // get our index into our reverse compliment map m
 458     const int second_compliment_index = (*rc_i)[1];
 459     // lookup the forward or reverse compliment depending on our rc flag
 460     const char complimented_second = m[second_basepair][second_compliment_index];
 461
 462     BOOST_CHECK_EQUAL( first_basepair, complimented_second) ;
 463   }
 464 }
 465
 466 BOOST_AUTO_TEST_CASE( two_way_local_alignment )
 467 {
 468   string s0("GCGCATAT");
 469   string s1("AAAAAAAT");
 470   Sequence seq1(s1);
 471
 472   Mussa analysis;
 473   analysis.append_sequence(s0);
 474   analysis.append_sequence(s1);
 475   analysis.set_threshold(3);
 476   analysis.set_window(4);
 477   analysis.analyze();
 478   NwayPaths npath = analysis.paths();
 479   BOOST_REQUIRE_EQUAL( npath.pathz.size(), 2 );
 480
 481   list<ConservedPath::path_type> result;
 482   list<vector<bool> > reversed;
 483   list<ConservedPath>::iterator pathz_i = npath.pathz.begin();
 484
 485   list<ConservedPath> selected_paths;
 486   selected_paths.push_back(*pathz_i);
 487   analysis.createLocalAlignment(selected_paths.begin(),
 488                                 selected_paths.end(),
 489                                 result,
 490                                 reversed);
 491
 492   two_way_local_align_test(analysis.sequences(), result, reversed);
 493
 494   ++pathz_i;
 495   result.clear();
 496   reversed.clear();
 497   selected_paths.clear();
 498   selected_paths.push_back(*pathz_i);
 499   analysis.createLocalAlignment(selected_paths.begin(),
 500                                 selected_paths.end(),
 501                                 result,
 502                                 reversed);
 503   two_way_local_align_test(analysis.sequences(), result, reversed);
 504 }
 505
 506 BOOST_AUTO_TEST_CASE( three_way_local_alignment )
 507 {
 508   string s0("AGCAGGGAGGGTTTAAATGGCACCCAGCAGTTGGTGTGAGG");
 509   string s1("AGCGGGAAGGGTTTAAATGGCACCGGGCAGTTGGCGTGAGG");
 510   string s2("CAGCGCCGGGGTTTAAATGGCACCGAGCAGTTGGCGCAGGG");
 511
 512   Mussa analysis;
 513   analysis.append_sequence(s0);
 514   analysis.append_sequence(s1);
 515   analysis.append_sequence(s2);
 516   analysis.set_threshold(23);
 517   analysis.set_window(30);
 518   analysis.analyze();
 519   NwayPaths npath = analysis.paths();
 520   BOOST_CHECK_EQUAL( npath.refined_pathz.size(), 1 );
 521
 522   list<ConservedPath::path_type> result;
 523   list<vector<bool> > reversed;
 524   // grab 1 path (since there's only one)
 525   list<ConservedPath>::iterator pathz_i = npath.pathz.begin();
 526   list<ConservedPath> selected_paths;
 527   selected_paths.push_back(*pathz_i);
 528   analysis.createLocalAlignment(selected_paths.begin(),
 529                                 selected_paths.end(),
 530                                 result,
 531                                 reversed);
 532
 533   for(std::list<ConservedPath::path_type>::iterator result_i = result.begin();
 534       result_i != result.end();
 535       ++result_i)
 536   {
 537     ConservedPath::path_element first_element = *(result_i->begin());
 538     for (ConservedPath::path_type::iterator element_i = result_i->begin();
 539          element_i != result_i->end();
 540          ++element_i)
 541     {
 542       BOOST_CHECK_EQUAL( *element_i, first_element );
 543       BOOST_CHECK_EQUAL( s0[*element_i], s1[*element_i] );
 544       BOOST_CHECK_EQUAL( s1[*element_i], s2[*element_i] );
 545       BOOST_CHECK_EQUAL( s0[*element_i], s2[*element_i] );
 546     }
 547   }
 548 }
 549
 550 BOOST_AUTO_TEST_CASE( mussa_window_larger_than_sequence )
 551 {
 552   string s0("AGCAGGG");
 553   string s1("CAGCGGG");
 554
 555   Mussa analysis;
 556   analysis.append_sequence(s0);
 557   analysis.append_sequence(s1);
 558   analysis.set_threshold(23);
 559   analysis.set_window(30);
 560   BOOST_CHECK_THROW(analysis.analyze(), seqcomp_error);
 561 }
 562
 563 BOOST_AUTO_TEST_CASE( subanalysis )
 564 {
 565   Sequence s1("AATGAAGATTTTAATGCTTTAATTTTGTTTTGTAAACTTCGAATTTCCAAAATTTGAAA");
 566   Sequence s2("AGGAGCAAGTTCGCTTCATCGAGAATTTTTAATTTTTAGTCAAATTTTCCAATGTCTGA");
 567
 568   Mussa analysis;
 569   analysis.append_sequence(s1);
 570   analysis.append_sequence(s2);
 571   analysis.set_threshold(8);
 572   analysis.set_window(8);
 573   analysis.analyze();
 574
 575   NwayPaths perfect_path = analysis.paths();
 576   int perfect_match_count = perfect_path.pathz.size();
 577
 578   Sequence sub1 = s1.subseq(2, s1.size()-4);
 579   Sequence sub2 = s2.subseq(2, s2.size()-4);
 580   Mussa subanalysis;
 581   subanalysis.append_sequence(sub1);
 582   subanalysis.append_sequence(sub2);
 583   subanalysis.set_threshold(7);
 584   subanalysis.set_window(8);
 585   subanalysis.analyze();
 586   NwayPaths one_mismatch_path = subanalysis.paths();
 587   int one_mismatch_count = one_mismatch_path.pathz.size();
 588
 589   BOOST_CHECK( perfect_match_count < one_mismatch_count );
 590 }
 591
 592 BOOST_AUTO_TEST_CASE( dirty_flag )
 593 {
 594   Mussa m;
 595   BOOST_CHECK_EQUAL(m.is_dirty(), false);
 596   m.set_name("foo");
 597   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 598   m.clear();
 599   m.set_window(30);
 600   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 601   m.clear();
 602   m.set_threshold(1);
 603   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 604   m.clear();
 605   m.set_soft_threshold(1);
 606   BOOST_CHECK_EQUAL(m.is_dirty(), false);
 607   m.clear();
 608   m.append_sequence("AAGGCCTT");
 609   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 610 }
 611