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