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