alg/test/test_mussa.cpp

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