alg/test/test_mussa.cpp

   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;
   8
   9 #include <string>
  10 #include <sstream>
  11 #include <vector>
  12
  13 #include "alg/mussa.hpp"
  14 #include "mussa_exceptions.hpp"
  15
  16 using namespace std;
  17
  18 //! can we initialize a mussa object?
  19 BOOST_AUTO_TEST_CASE( mussa_simple )
  20 {
  21   Mussa m;
  22   BOOST_CHECK_EQUAL(m.get_name(), "" );
  23   BOOST_CHECK_EQUAL(m.get_window(), 0);
  24   BOOST_CHECK_EQUAL(m.get_threshold(), 0);
  25   BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::TransitiveNway);
  26   m.set_name( "hello" );
  27   BOOST_CHECK_EQUAL(m.get_name(), "hello" );
  28   m.set_window(30);
  29   BOOST_CHECK_EQUAL(m.get_window(), 30);
  30   m.set_threshold(21);
  31   BOOST_CHECK_EQUAL(m.get_threshold(), 21);
  32   BOOST_CHECK_EQUAL(m.get_soft_threshold(), 21);
  33   m.set_soft_threshold(19);
  34   BOOST_CHECK_EQUAL(m.get_soft_threshold(), 21);
  35   m.set_soft_threshold(35);
  36   BOOST_CHECK_EQUAL(m.get_soft_threshold(), 30);
  37   m.set_soft_threshold(25);
  38   BOOST_CHECK_EQUAL(m.get_soft_threshold(), 25);
  39   m.set_analysis_mode(Mussa::RadialNway);
  40   BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::RadialNway);
  41
  42   m.clear();
  43   BOOST_CHECK_EQUAL(m.get_name(), "" );
  44   BOOST_CHECK_EQUAL(m.get_window(), 0);
  45   BOOST_CHECK_EQUAL(m.get_threshold(), 0);
  46   BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::TransitiveNway);
  47 }
  48
  49 BOOST_AUTO_TEST_CASE( mussa_analysis_name )
  50 {
  51   Mussa m;
  52   m.set_analysis_mode( Mussa::TransitiveNway );
  53   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Transitive" );
  54   m.set_analysis_mode( Mussa::RadialNway );
  55   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Radial" );
  56   m.set_analysis_mode( Mussa::EntropyNway );
  57   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Entropy" );
  58   m.set_analysis_mode( Mussa::RecursiveNway);
  59   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "[deprecated] Recursive" );
  60 }
  61
  62 BOOST_AUTO_TEST_CASE( mussa_sequences )
  63 {
  64   std::string s0("AAAANNNN");
  65   std::string s1("GGGGNNNN");
  66   std::string s2("TTTTNNNN");
  67
  68   Mussa analysis;
  69   analysis.append_sequence(s0);
  70   analysis.append_sequence(s1);
  71   analysis.append_sequence(s2);
  72
  73   BOOST_CHECK_EQUAL( analysis.sequences().size(), 3 );
  74   BOOST_CHECK_EQUAL( *(analysis.sequences()[0]), s0);
  75   BOOST_CHECK_EQUAL( *(analysis.sequences()[1]), s1);
  76   BOOST_CHECK_EQUAL( *(analysis.sequences()[2]), s2);
  77 }
  78
  79 // for some reason we can call nway once safely but it
  80 // somehow changed things and caused a segfault
  81 // fixed by adding a return statement in trans_path_search
  82 BOOST_AUTO_TEST_CASE ( empty_mussa_set_threshold )
  83 {
  84   Mussa m;
  85   m.set_soft_threshold(15);
  86   m.nway();
  87
  88   m.set_soft_threshold(25);
  89   m.nway();
  90 }
  91
  92 BOOST_AUTO_TEST_CASE( mussa_load_mupa )
  93 {
  94   fs::path mupa_path(EXAMPLE_DIR, fs::native);
  95   fs::path result_path = fs::initial_path() / "mck3test_w30_t20";
  96   mupa_path /= "mck3test.mupa";
  97
  98   Mussa m1;
  99   m1.load_mupa_file( mupa_path );
 100   m1.analyze();
 101   m1.save( result_path );
 102   BOOST_CHECK_EQUAL( m1.get_name(), std::string("mck3test") );
 103   BOOST_CHECK( m1.size() > 0 );
 104
 105   Mussa m2;
 106   m2.load( result_path );
 107   BOOST_CHECK_EQUAL( m2.get_name(), result_path.leaf() );
 108   BOOST_CHECK_EQUAL( m1.size(), m2.size() );
 109 }
 110
 111 BOOST_AUTO_TEST_CASE( mussa_load_full_path )
 112 {
 113   Mussa m1;
 114   fs::path full_path(fs::path(EXAMPLE_DIR, fs::native) / "mck3test.mupa");
 115   m1.load_mupa_file( full_path );
 116   m1.analyze();
 117
 118   BOOST_CHECK( m1.size() > 0);
 119   BOOST_CHECK_EQUAL( m1.get_window(), 30 );
 120   BOOST_CHECK_EQUAL( m1.get_threshold(), 20);
 121 }
 122
 123 BOOST_AUTO_TEST_CASE( mussa_mupa_directory )
 124 {
 125   fs::path curdir(".");
 126   Mussa m1;
 127   BOOST_CHECK_THROW(m1.load_mupa_file( curdir ), mussa_load_error );
 128 }
 129
 130 BOOST_AUTO_TEST_CASE( mussa_load_analysis )
 131 {
 132   fs::path example_dir(EXAMPLE_DIR, fs::native);
 133   Mussa m1;
 134   m1.load_mupa_file( example_dir / "mck3test.mupa" );
 135   m1.analyze();
 136
 137   Mussa m2;
 138   m2.load( fs::initial_path() / "mck3test_w30_t20");
 139
 140   BOOST_CHECK_EQUAL( m1.size(), m2.size() );
 141   BOOST_CHECK_EQUAL( m1.get_window(), m2.get_window() );
 142   BOOST_CHECK_EQUAL( m1.get_threshold(), m2.get_threshold() );
 143 }
 144
 145 BOOST_AUTO_TEST_CASE( mussa_load_motif )
 146 {
 147   string data = "AAGG 1.0 1.0 0.0\n"
 148                 "GGTT 0.0 0.1 1.0\n"
 149                 "ZXY 2 1.9 0\n";
 150
 151   istringstream test_istream(data);
 152
 153   Mussa m1;
 154   m1.append_sequence("AAAAGGGGTTTT");
 155   m1.append_sequence("GGGCCCCTTGGTT");
 156   m1.load_motifs(test_istream);
 157
 158   for (Mussa::vector_sequence_type::const_iterator seq_i = m1.sequences().begin();
 159        seq_i != m1.sequences().end();
 160        ++seq_i)
 161   {
 162     BOOST_CHECK( (*seq_i)->motifs().size() > 0 );
 163   }
 164 }
 165
 166 BOOST_AUTO_TEST_CASE( mussa_add_motif )
 167 {
 168   vector<Sequence> motifs;
 169   motifs.push_back("AAGG");
 170   vector<Color> colors;
 171   colors.push_back(Color(1.0, 0.0, 0.0));
 172
 173   Mussa m1;
 174   m1.append_sequence("AAAAGGGGTTTT");
 175   m1.append_sequence("GGGCCCCTTGGTT");
 176   m1.set_motifs(motifs, colors);
 177   int first_size = m1.motifs().size();
 178   BOOST_CHECK_EQUAL( first_size, 1 );
 179   BOOST_REQUIRE(first_size > 0);
 180   BOOST_CHECK_EQUAL(*(m1.motifs().begin()), motifs.front());
 181   // make sure that our sequences have the right number of motifs
 182   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 1);
 183   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1); // because of rc
 184
 185   // verify that setting the motif clears the arrays
 186   m1.set_motifs(motifs, colors);
 187   BOOST_CHECK_EQUAL( first_size, m1.motifs().size() );
 188   // make sure that our sequences have the right number of motifs
 189   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 1);
 190   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
 191
 192   // add a different motif
 193   motifs.clear();
 194   motifs.push_back("CCTTGG");
 195   BOOST_CHECK_EQUAL(motifs.size(), 1);
 196   m1.set_motifs(motifs, colors);
 197   BOOST_CHECK_EQUAL(m1.motifs().size(), 1);
 198   BOOST_REQUIRE(m1.motifs().size() > 0);
 199   BOOST_CHECK_EQUAL(*(m1.motifs().begin()), motifs.front());
 200   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 0);
 201   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
 202
 203   // try a motif that doesn't exist
 204   motifs.clear();
 205   motifs.push_back("CCTTGG");
 206   BOOST_CHECK_EQUAL(motifs.size(), 1);
 207   m1.set_motifs(motifs, colors);
 208   BOOST_CHECK_EQUAL(m1.motifs().size(), 1);
 209   BOOST_CHECK_EQUAL(m1.sequences()[0]->motifs().size(), 0);
 210   BOOST_CHECK_EQUAL(m1.sequences()[1]->motifs().size(), 1);
 211
 212 }
 213
 214 static void
 215 local_align_test(const Mussa::vector_sequence_type &seqs,
 216                  const list<ConservedPath::path_type>& result,
 217                  const list<vector<bool> >& reversed)
 218 {
 219   map<char, vector <char> >  m;
 220   assign::insert(m)('A', assign::list_of('A')('T') )
 221                    ('T', assign::list_of('T')('A') )
 222                    ('G', assign::list_of('G')('C') )
 223                    ('C', assign::list_of('C')('G') );
 224   list<vector<bool> >::const_iterator rc_i = reversed.begin();
 225
 226   for(list<ConservedPath::path_type>::const_iterator base_i = result.begin();
 227       base_i != result.end();
 228       ++base_i, ++rc_i)
 229   {
 230     // since the reverse compliment flag is relative to the first sequence
 231     // the first one should always be false
 232     BOOST_CHECK_EQUAL( (*rc_i)[0], false );
 233     const int first_path_basepair_index = (*base_i)[0];
 234     const int second_path_basepair_index = (*base_i)[1];
 235     const char first_basepair = (*seqs[0])[first_path_basepair_index];
 236     const char second_basepair = (*seqs[1])[second_path_basepair_index];
 237     // get our index into our reverse compliment map m
 238     const int second_compliment_index = (*rc_i)[1];
 239     // lookup the forward or reverse compliment depending on our rc flag
 240     const char complimented_second = m[second_basepair][second_compliment_index];
 241
 242     BOOST_CHECK_EQUAL( first_basepair, complimented_second) ;
 243   }
 244 }
 245
 246
 247 BOOST_AUTO_TEST_CASE( local_alignment )
 248 {
 249   string s0("GCGCATAT");
 250   string s1("AAAAAAAT");
 251   Sequence seq1(s1);
 252
 253   Mussa analysis;
 254   analysis.append_sequence(s0);
 255   analysis.append_sequence(s1);
 256   analysis.set_threshold(3);
 257   analysis.set_window(4);
 258   analysis.analyze();
 259   NwayPaths npath = analysis.paths();
 260   list<ConservedPath::path_type> result;
 261   list<vector<bool> > reversed;
 262   list<ConservedPath>::iterator pathz_i = npath.pathz.begin();
 263
 264   list<ConservedPath> selected_paths;
 265   selected_paths.push_back(*pathz_i);
 266   analysis.createLocalAlignment(selected_paths.begin(),
 267                                 selected_paths.end(),
 268                                 result,
 269                                 reversed);
 270
 271   local_align_test(analysis.sequences(), result, reversed);
 272
 273   ++pathz_i;
 274   result.clear();
 275   reversed.clear();
 276   selected_paths.clear();
 277   selected_paths.push_back(*pathz_i);
 278   analysis.createLocalAlignment(selected_paths.begin(),
 279                                 selected_paths.end(),
 280                                 result,
 281                                 reversed);
 282   local_align_test(analysis.sequences(), result, reversed);
 283
 284
 285 }
 286
 287 BOOST_AUTO_TEST_CASE( subanalysis )
 288 {
 289   Sequence s1("AATGAAGATTTTAATGCTTTAATTTTGTTTTGTAAACTTCGAATTTCCAAAATTTGAAA");
 290   Sequence s2("AGGAGCAAGTTCGCTTCATCGAGAATTTTTAATTTTTAGTCAAATTTTCCAATGTCTGA");
 291
 292   Mussa analysis;
 293   analysis.append_sequence(s1);
 294   analysis.append_sequence(s2);
 295   analysis.set_threshold(8);
 296   analysis.set_window(8);
 297   analysis.analyze();
 298
 299   NwayPaths perfect_path = analysis.paths();
 300   int perfect_match_count = perfect_path.pathz.size();
 301
 302   Sequence sub1 = s1.subseq(2, s1.size()-4);
 303   Sequence sub2 = s2.subseq(2, s2.size()-4);
 304   Mussa subanalysis;
 305   subanalysis.append_sequence(sub1);
 306   subanalysis.append_sequence(sub2);
 307   subanalysis.set_threshold(7);
 308   subanalysis.set_window(8);
 309   subanalysis.analyze();
 310   NwayPaths one_mismatch_path = subanalysis.paths();
 311   int one_mismatch_count = one_mismatch_path.pathz.size();
 312
 313   BOOST_CHECK( perfect_match_count < one_mismatch_count );
 314 }
 315
 316 BOOST_AUTO_TEST_CASE( dirty_flag )
 317 {
 318   Mussa m;
 319   BOOST_CHECK_EQUAL(m.is_dirty(), false);
 320   m.set_name("foo");
 321   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 322   m.clear();
 323   m.set_window(30);
 324   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 325   m.clear();
 326   m.set_threshold(1);
 327   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 328   m.clear();
 329   m.set_soft_threshold(1);
 330   BOOST_CHECK_EQUAL(m.is_dirty(), false);
 331   m.clear();
 332   m.append_sequence("AAGGCCTT");
 333   BOOST_CHECK_EQUAL(m.is_dirty(), true);
 334 }
 335