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   m.set_analysis_mode(Mussa::RadialNway);
  32   BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::RadialNway);
  33
  34   m.clear();
  35   BOOST_CHECK_EQUAL(m.get_name(), "" );
  36   BOOST_CHECK_EQUAL(m.get_window(), 0);
  37   BOOST_CHECK_EQUAL(m.get_threshold(), 0);
  38   BOOST_CHECK_EQUAL(m.get_analysis_mode(), Mussa::TransitiveNway);
  39 }
  40
  41 BOOST_AUTO_TEST_CASE( mussa_analysis_name )
  42 {
  43   Mussa m;
  44   m.set_analysis_mode( Mussa::TransitiveNway );
  45   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Transitive" );
  46   m.set_analysis_mode( Mussa::RadialNway );
  47   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Radial" );
  48   m.set_analysis_mode( Mussa::EntropyNway );
  49   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "Entropy" );
  50   m.set_analysis_mode( Mussa::RecursiveNway);
  51   BOOST_CHECK_EQUAL( m.get_analysis_mode_name(), "[deprecated] Recursive" );
  52 }
  53
  54 BOOST_AUTO_TEST_CASE( mussa_sequences )
  55 {
  56   std::string s0("AAAANNNN");
  57   std::string s1("GGGGNNNN");
  58   std::string s2("TTTTNNNN");
  59
  60   Mussa analysis;
  61   analysis.add_a_seq(s0);
  62   analysis.add_a_seq(s1);
  63   analysis.add_a_seq(s2);
  64
  65   BOOST_CHECK_EQUAL( analysis.sequences().size(), 3 );
  66   BOOST_CHECK_EQUAL( analysis.sequences()[0].get_seq(), s0);
  67   BOOST_CHECK_EQUAL( analysis.sequences()[1].get_seq(), s1);
  68   BOOST_CHECK_EQUAL( analysis.sequences()[2].get_seq(), s2);
  69 }
  70
  71 // for some reason we can call nway once safely but it
  72 // somehow changed things and caused a segfault
  73 // fixed by adding a return statement in trans_path_search
  74 BOOST_AUTO_TEST_CASE ( empty_mussa_set_threshold )
  75 {
  76   Mussa m;
  77   m.set_soft_thres(15);
  78   m.nway();
  79
  80   m.set_soft_thres(25);
  81   m.nway();
  82 }
  83
  84 BOOST_AUTO_TEST_CASE( mussa_load_mupa )
  85 {
  86   fs::path mupa_path(EXAMPLE_DIR);
  87   mupa_path /= "mck3test.mupa";
  88
  89   Mussa m1;
  90   m1.load_mupa_file( mupa_path );
  91   m1.analyze(0, 0);
  92   BOOST_CHECK_EQUAL( m1.get_name(), std::string("mck3test") );
  93   BOOST_CHECK( m1.size() > 0 );
  94
  95   Mussa m2;
  96   fs::path result_path = fs::initial_path() / "mck3test_w30_t20";
  97   m2.load( result_path );
  98   BOOST_CHECK_EQUAL( m2.get_name(), result_path.leaf() );
  99   BOOST_CHECK_EQUAL( m1.size(), m2.size() );
 100
 101 }
 102
 103 BOOST_AUTO_TEST_CASE( mussa_load_full_path )
 104 {
 105   Mussa m1;
 106   fs::path full_path(fs::path(EXAMPLE_DIR) / "mck3test.mupa");
 107   m1.load_mupa_file( full_path );
 108   m1.analyze(0, 0);
 109
 110   BOOST_CHECK( m1.size() > 0);
 111   BOOST_CHECK_EQUAL( m1.get_window(), 30 );
 112   BOOST_CHECK_EQUAL( m1.get_threshold(), 20);
 113 }
 114
 115 BOOST_AUTO_TEST_CASE( mussa_load_analysis )
 116 {
 117   fs::path example_dir(EXAMPLE_DIR);
 118   Mussa m1;
 119   m1.load_mupa_file( example_dir / "mck3test.mupa" );
 120   m1.analyze(0, 0);
 121
 122   Mussa m2;
 123   m2.load( fs::initial_path() / "mck3test_w30_t20");
 124
 125   BOOST_CHECK_EQUAL( m1.size(), m2.size() );
 126   BOOST_CHECK_EQUAL( m1.get_window(), m2.get_window() );
 127   BOOST_CHECK_EQUAL( m1.get_threshold(), m2.get_threshold() );
 128 }
 129
 130 BOOST_AUTO_TEST_CASE( mussa_load_motif )
 131 {
 132   string data = "AAGG 1.0 1.0 0.0\n"
 133                 "GGTT 0.0 0.1 1.0\n"
 134                 "ZXY 2 1.9 0\n";
 135
 136   istringstream test_istream(data);
 137
 138   Mussa m1;
 139   m1.add_a_seq("AAAAGGGGTTTT");
 140   m1.add_a_seq("GGGCCCCTTGGTT");
 141   m1.load_motifs(test_istream);
 142
 143   for (vector<Sequence>::const_iterator seq_i = m1.sequences().begin();
 144        seq_i != m1.sequences().end();
 145        ++seq_i)
 146   {
 147     BOOST_CHECK( seq_i->motifs().size() > 0 );
 148   }
 149 }
 150
 151 BOOST_AUTO_TEST_CASE( mussa_add_motif )
 152 {
 153   vector<string> motifs;
 154   motifs.push_back("AAGG");
 155   vector<Color> colors;
 156   colors.push_back(Color(1.0, 0.0, 0.0));
 157
 158   Mussa m1;
 159   m1.add_a_seq("AAAAGGGGTTTT");
 160   m1.add_a_seq("GGGCCCCTTGGTT");
 161   m1.add_motifs(motifs, colors);
 162   int first_size = m1.motifs().size();
 163   BOOST_CHECK_EQUAL( first_size, 1 );
 164   m1.add_motifs(motifs, colors);
 165   BOOST_CHECK_EQUAL( first_size, m1.motifs().size() );
 166 }
 167
 168 static void
 169 local_align_test(const vector<Sequence> &seqs,
 170                  const list<ConservedPath::path_type>& result,
 171                  const list<vector<bool> >& reversed)
 172 {
 173   map<char, vector <char> >  m;
 174   assign::insert(m)('A', assign::list_of('A')('T') )
 175                    ('T', assign::list_of('T')('A') )
 176                    ('G', assign::list_of('G')('C') )
 177                    ('C', assign::list_of('C')('G') );
 178   list<vector<bool> >::const_iterator rc_i = reversed.begin();
 179
 180   for(list<ConservedPath::path_type>::const_iterator base_i = result.begin();
 181       base_i != result.end();
 182       ++base_i, ++rc_i)
 183   {
 184     // since the reverse compliment flag is relative to the first sequence
 185     // the first one should always be false
 186     BOOST_CHECK_EQUAL( (*rc_i)[0], false );
 187     const int first_path_basepair_index = (*base_i)[0];
 188     const int second_path_basepair_index = (*base_i)[1];
 189     const char first_basepair = seqs[0][first_path_basepair_index];
 190     const char second_basepair = seqs[1][second_path_basepair_index];
 191     // get our index into our reverse compliment map m
 192     const int second_compliment_index = (*rc_i)[1];
 193     // lookup the forward or reverse compliment depending on our rc flag
 194     const char complimented_second = m[second_basepair][second_compliment_index];
 195
 196     BOOST_CHECK_EQUAL( first_basepair, complimented_second) ;
 197   }
 198 }
 199
 200
 201 BOOST_AUTO_TEST_CASE( local_alignment )
 202 {
 203   string s0("GCGCATAT");
 204   string s1("AAAAAAAT");
 205   Sequence seq1(s1);
 206
 207   Mussa analysis;
 208   analysis.add_a_seq(s0);
 209   analysis.add_a_seq(s1);
 210   analysis.analyze(4,3);
 211   NwayPaths npath = analysis.paths();
 212   list<ConservedPath::path_type> result;
 213   list<vector<bool> > reversed;
 214   list<ConservedPath>::iterator pathz_i = npath.pathz.begin();
 215
 216   list<ConservedPath> selected_paths;
 217   selected_paths.push_back(*pathz_i);
 218   analysis.createLocalAlignment(selected_paths.begin(),
 219                                 selected_paths.end(),
 220                                 result,
 221                                 reversed);
 222
 223   local_align_test(analysis.sequences(), result, reversed);
 224
 225   ++pathz_i;
 226   result.clear();
 227   reversed.clear();
 228   selected_paths.clear();
 229   selected_paths.push_back(*pathz_i);
 230   analysis.createLocalAlignment(selected_paths.begin(),
 231                                 selected_paths.end(),
 232                                 result,
 233                                 reversed);
 234   local_align_test(analysis.sequences(), result, reversed);
 235
 236
 237 }
 238
 239