alg/test/test_sequence.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
   6 #include <list>
   7 #include <iostream>
   8 #include <sstream>
   9
  10 #include "alg/sequence.hpp"
  11 #include "mussa_exceptions.hpp"
  12
  13 using namespace std;
  14
  15 //! when we try to load a missing file, do we get an error?
  16 BOOST_AUTO_TEST_CASE( sequence_load_exception )
  17 {
  18   Sequence s;
  19   // there should be errors when we try to load something that doesn't exist
  20   BOOST_CHECK_THROW( s.load_fasta("alkejralk", 1, 0, 0), mussa_load_error);
  21   BOOST_CHECK_THROW( s.load_annot("alkejralk", 0, 0), mussa_load_error);
  22 }
  23
  24 BOOST_AUTO_TEST_CASE( sequence_eol_conventions )
  25 {
  26   string header(">Header");
  27   string line1("AAAAGGGGCCCCTTTTT");
  28   string line2("AAAAGGGGCCCCTTTTT");
  29   int seq_len = line1.size() + line2.size();
  30
  31   stringstream cr;
  32   cr << header << "\015" << line1 << "\015" << line2 << "\015";
  33   Sequence seq_cr;
  34   seq_cr.load_fasta(cr);
  35
  36   stringstream crlf;
  37   crlf << header << "\015\012" << line1 << "\015\012" << line2 << "\015\012";
  38   Sequence seq_crlf;
  39   seq_crlf.load_fasta(crlf);
  40
  41   stringstream lf;
  42   lf << header << "\012" << line1 << "\012" << line2 << "\012";
  43   Sequence seq_lf;
  44   seq_lf.load_fasta(lf);
  45
  46   BOOST_CHECK_EQUAL(seq_cr.size(),   seq_len);
  47   BOOST_CHECK_EQUAL(seq_crlf.size(), seq_len);
  48   BOOST_CHECK_EQUAL(seq_cr.size(),   seq_len);
  49 }
  50
  51
  52 //! Do simple operations work correctly?
  53 BOOST_AUTO_TEST_CASE( sequence_filter )
  54 {
  55   const char *core_seq = "AATTGGCC";
  56   Sequence s1(core_seq);
  57   BOOST_CHECK_EQUAL(s1, core_seq);
  58
  59   Sequence s2("aattggcc");
  60   BOOST_CHECK_EQUAL(s2, "AATTGGCC");
  61   BOOST_CHECK_EQUAL(s2.rev_comp(), "GGCCAATT");
  62   BOOST_CHECK_EQUAL(s2.size(), s2.size());
  63   BOOST_CHECK_EQUAL(s2.c_str(), core_seq);
  64
  65   Sequence s3("asdfg");
  66   BOOST_CHECK_EQUAL(s3, "ANNNG");
  67   BOOST_CHECK_EQUAL(s3.subseq(0,2), "AN");
  68
  69   s3.set_filtered_sequence("AAGGCCTT", 0, 2);
  70   BOOST_CHECK_EQUAL(s3, "AA");
  71   s3.set_filtered_sequence("AAGGCCTT", 2, 2);
  72   BOOST_CHECK_EQUAL( s3, "GG");
  73   s3.set_filtered_sequence("AAGGCCTT", 4);
  74   BOOST_CHECK_EQUAL( s3, "CCTT");
  75
  76   s3.clear();
  77   BOOST_CHECK_EQUAL(s3, "");
  78
  79   s3 = "AAGGFF";
  80   BOOST_CHECK_EQUAL(s3, "AAGGNN");
  81 }
  82
  83 //! Can we load data from a file
  84 BOOST_AUTO_TEST_CASE( sequence_load )
  85 {
  86   fs::path seq_path(fs::path(EXAMPLE_DIR)/ "seq" );
  87   seq_path /=  "human_mck_pro.fa";
  88   Sequence s;
  89   s.load_fasta(seq_path);
  90   BOOST_CHECK_EQUAL(s.subseq(0, 5), "GGATC"); // first few chars of fasta file
  91   BOOST_CHECK_EQUAL(s.subseq(2, 3), "ATC");
  92   BOOST_CHECK_EQUAL(s.get_fasta_header(), "gi|180579|gb|M21487.1|HUMCKMM1 Human "
  93                                     "muscle creatine kinase gene (CKMM), "
  94                                     "5' flank");
  95 }
  96
  97 BOOST_AUTO_TEST_CASE( annotation_load )
  98 {
  99   string annot_data = "human\n"
 100                       "0 10 name   type\n"
 101                       "10 20 myf7\n"
 102                       "20 30 myod\n"
 103                       "50\t55 anothername\n"
 104                       "60 50 backward\n"
 105                       ">ident3 asdf\n"
 106                       "GCT\n"
 107                       "gCTn\n"
 108                       "75\t90\tname2\ttype2\n"
 109                       "100 120 name-asdf type!@#$%\n"
 110                       ;
 111   string s(100, 'A');
 112   s += "GCTGCTAATT";
 113   Sequence seq(s);
 114
 115   //istringstream annot_stream(annot_data);
 116   seq.parse_annot(annot_data, 0, 0);
 117   std::list<annot> annots_list = seq.annotations();
 118   std::vector<annot> annots(annots_list.begin(), annots_list.end());
 119   BOOST_REQUIRE_EQUAL( annots.size(), 8);
 120   BOOST_CHECK_EQUAL( annots[0].start, 0 );
 121   BOOST_CHECK_EQUAL( annots[0].end, 10 );
 122   BOOST_CHECK_EQUAL( annots[0].type, "type");
 123   BOOST_CHECK_EQUAL( annots[0].name, "name");
 124   BOOST_CHECK_EQUAL( annots[1].name, "myf7");
 125   BOOST_CHECK_EQUAL( annots[2].name, "myod");
 126   BOOST_CHECK_EQUAL( annots[3].name, "anothername");
 127   BOOST_CHECK_EQUAL( annots[4].name, "backward");
 128   BOOST_CHECK_EQUAL( annots[5].name, "name2");
 129   BOOST_CHECK_EQUAL( annots[5].end, 90);
 130   BOOST_CHECK_EQUAL( annots[6].start, 100);
 131   BOOST_CHECK_EQUAL( annots[6].end, 120);
 132   BOOST_CHECK_EQUAL( annots[6].name, "name-asdf");
 133   BOOST_CHECK_EQUAL( annots[6].type, "type!@#$%");
 134   // sequence defined annotations will always be after the
 135   // absolute positions
 136   BOOST_CHECK_EQUAL( annots[7].name, "ident3 asdf");
 137   BOOST_CHECK_EQUAL( annots[7].start, 100);
 138
 139   //BOOST_CHECK_EQUAL( annots
 140 }
 141
 142 BOOST_AUTO_TEST_CASE( annotation_load_no_species_name )
 143 {
 144   string annot_data = "0 10 name   type\n"
 145                       "10 20 myf7\n"
 146                       "20 30 myod\n"
 147                       "50\t55 anothername\n"
 148                       "60 50 backward\n"
 149                       ">ident3 asdf\n"
 150                       "GCT\n"
 151                       "gCTn\n"
 152                       "75\t90\tname2\ttype2\n"
 153                       "100 120 name-asdf type!@#$%\n"
 154                       ;
 155   string s(100, 'A');
 156   s += "GCTGCTAATT";
 157   Sequence seq(s);
 158
 159   //istringstream annot_stream(annot_data);
 160   seq.parse_annot(annot_data, 0, 0);
 161   std::list<annot> annots_list = seq.annotations();
 162   std::vector<annot> annots(annots_list.begin(), annots_list.end());
 163   BOOST_REQUIRE_EQUAL( annots.size(), 8);
 164   BOOST_CHECK_EQUAL( annots[0].start, 0 );
 165   BOOST_CHECK_EQUAL( annots[0].end, 10 );
 166   BOOST_CHECK_EQUAL( annots[0].type, "type");
 167 }
 168
 169 // ticket:83 when you try to load a sequence from a file that doesn't
 170 // have fasta headers it crashes.
 171 BOOST_AUTO_TEST_CASE( sequence_past_end )
 172 {
 173   fs::path seq_path(fs::path(EXAMPLE_DIR)/ "seq" );
 174   seq_path /=  "misformated_seq.fa";
 175   Sequence s;
 176   BOOST_CHECK_THROW( s.load_fasta(seq_path), mussa_load_error );
 177 }
 178
 179 BOOST_AUTO_TEST_CASE ( sequence_empty )
 180 {
 181   Sequence s;
 182   BOOST_CHECK_EQUAL( s.empty(), true );
 183   s = "AAAGGG";
 184   BOOST_CHECK_EQUAL( s.empty(), false );
 185 }
 186
 187 BOOST_AUTO_TEST_CASE ( sequence_iterators )
 188 {
 189   std::string seq_string = "AAGGCCTTNNTATA";
 190   Sequence s(seq_string);
 191   const Sequence cs(s);
 192   std::string::size_type count = 0;
 193
 194   std::string::iterator str_itor;
 195   Sequence::iterator s_itor;
 196   Sequence::const_iterator cs_itor;
 197
 198   for( str_itor = seq_string.begin(),
 199        s_itor   = s.begin(),
 200        cs_itor  = cs.begin();
 201        str_itor != seq_string.end() and
 202        s_itor   != s.end() and
 203        cs_itor  != cs.end();
 204        ++str_itor, ++s_itor, ++cs_itor, ++count)
 205   {
 206     BOOST_CHECK_EQUAL ( *str_itor, *s_itor );
 207     BOOST_CHECK_EQUAL ( *s_itor, *cs_itor );
 208     BOOST_CHECK_EQUAL ( *cs_itor, *str_itor );
 209   }
 210   BOOST_CHECK_EQUAL( seq_string.size(), count );
 211   BOOST_CHECK_EQUAL( s.size(), count );
 212   BOOST_CHECK_EQUAL( cs.size(), count );
 213 }
 214
 215 BOOST_AUTO_TEST_CASE( sequence_motifs )
 216 {
 217   string m("AAAA");
 218   string bogus("AATTGGAA");
 219   Sequence s1("AAAAGGGGCCCCTTTT");
 220
 221   list<motif>::const_iterator motif_i = s1.motifs().begin();
 222   list<motif>::const_iterator motif_end = s1.motifs().end();
 223
 224   // do our iterators work?
 225   BOOST_CHECK( motif_i == s1.motifs().begin() );
 226   BOOST_CHECK( motif_end == s1.motifs().end() );
 227   BOOST_CHECK( motif_i == motif_end );
 228
 229   // this shouldn't show up
 230   s1.add_motif(bogus);
 231   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 232   BOOST_CHECK_EQUAL( s1.motifs().size(), 0 );
 233
 234   s1.add_motif(m);
 235   BOOST_CHECK( s1.motifs().begin() != s1.motifs().end() );
 236   BOOST_CHECK_EQUAL( s1.motifs().size(), 2 );
 237
 238   for(motif_i = s1.motifs().begin();
 239       motif_i != s1.motifs().end();
 240       ++motif_i)
 241   {
 242     BOOST_CHECK_EQUAL( motif_i->type, "motif" );
 243     BOOST_CHECK_EQUAL( motif_i->name, m);
 244     BOOST_CHECK_EQUAL( motif_i->sequence, m);
 245   }
 246
 247   s1.clear_motifs();
 248   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 249
 250   /* FIXME: enable this when i find a way of passing storing the motif name
 251   // does our annotation travel?
 252   Sequence motif_seq(m);
 253   motif_seq.set_fasta_header("hi");
 254   s1.add_motif(motif_seq);
 255
 256   BOOST_CHECK_EQUAL(s1.motifs().size(), 2);
 257   for(motif_i = s1.motifs().begin();
 258       motif_i != s1.motifs().end();
 259       ++motif_i)
 260   {
 261     BOOST_CHECK_EQUAL( motif_i->type, "motif" );
 262     BOOST_CHECK_EQUAL( motif_i->name, "hi");
 263     BOOST_CHECK_EQUAL( motif_i->sequence, m);
 264   }
 265   */
 266 }
 267
 268 BOOST_AUTO_TEST_CASE( annot_test )
 269 {
 270   annot a(0, 10, "test", "thing");
 271
 272   BOOST_CHECK_EQUAL( a.start, 0 );
 273   BOOST_CHECK_EQUAL( a.end,   10 );
 274   BOOST_CHECK_EQUAL( a.type,  "test" );
 275   BOOST_CHECK_EQUAL( a.name,  "thing" );
 276
 277   motif m(10, "AAGGCC");
 278   BOOST_CHECK_EQUAL( m.start, 10 );
 279   BOOST_CHECK_EQUAL( m.type, "motif" );
 280   BOOST_CHECK_EQUAL( m.name, "AAGGCC" );
 281   BOOST_CHECK_EQUAL( m.end,  10+6 );
 282 }
 283
 284 BOOST_AUTO_TEST_CASE( annotate_from_sequence )
 285 {
 286   string s("CCGCCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 287            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 288            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 289            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 290            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 291            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 292            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 293            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 294   string gc("GCCCCC");
 295   string gga("GGACACCTC");
 296   Sequence seq(s);
 297
 298   std::list<Sequence> query_list;
 299   std::list<string> string_list;
 300   query_list.push_back(Sequence(gc));
 301   string_list.push_back(gc);
 302   query_list.push_back(Sequence(gga));
 303   string_list.push_back(gga);
 304
 305   BOOST_CHECK_EQUAL( seq.annotations().size(), 0 );
 306   seq.find_sequences(query_list.begin(), query_list.end());
 307
 308   int count = 0;
 309   for(list<string>::iterator string_i = string_list.begin();
 310       string_i != string_list.end();
 311       ++string_i)
 312   {
 313     string::size_type pos=0;
 314     while(pos != string::npos) {
 315       pos = s.find(*string_i, pos);
 316       if (pos != string::npos) {
 317         ++count;
 318         ++pos;
 319       }
 320     }
 321   }
 322   BOOST_CHECK_EQUAL(seq.annotations().size(), count);
 323 }
 324
 325 BOOST_AUTO_TEST_CASE( subseq_annotation_test )
 326 {
 327   string s("CCGCCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 328            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 329            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 330            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 331            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 332            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 333            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 334            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 335   Sequence seq(s);
 336
 337
 338   seq.add_annotation(annot(0, 10, "0-10", "0-10"));
 339   seq.add_annotation(annot(10, 20, "10-20", "10-20"));
 340   seq.add_annotation(annot(0, 20, "0-20", "0-20"));
 341   seq.add_annotation(annot(8, 12, "8-12", "8-12"));
 342   seq.add_annotation(annot(100, 5000, "100-5000", "100-5000"));
 343
 344   Sequence subseq = seq.subseq(5, 10);
 345   const list<annot> annots = subseq.annotations();
 346   // generate some ground truth
 347   list<annot> correct;
 348   correct.push_back(annot(0, 5, "0-10",  "0-10"));
 349   correct.push_back(annot(5,10, "10-20", "10-20"));
 350   correct.push_back(annot(0,10, "0-20",  "0-20"));
 351   correct.push_back(annot(3, 7, "8-12",  "8-12"));
 352   BOOST_REQUIRE_EQUAL( annots.size(), correct.size() );
 353
 354   list<annot>::iterator correct_i = correct.begin();
 355   list<annot>::const_iterator annot_i = annots.begin();
 356   for(; annot_i != annots.end(); ++annot_i, ++correct_i)
 357   {
 358     BOOST_CHECK( *annot_i == *correct_i );
 359   }
 360 }
 361
 362 BOOST_AUTO_TEST_CASE( motif_annotation_update )
 363 {
 364   string s("CCGTCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 365            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 366            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 367            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 368            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 369            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 370            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 371            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 372   Sequence seq(s);
 373
 374   // starting conditions
 375   BOOST_CHECK_EQUAL(seq.annotations().size(), 0);
 376   BOOST_CHECK_EQUAL(seq.motifs().size(), 0);
 377   seq.add_annotation(annot(0, 10, "0-10", "0-10"));
 378   seq.add_annotation(annot(10, 20, "10-20", "10-20"));
 379   seq.add_annotation(annot(0, 20, "0-20", "0-20"));
 380   BOOST_CHECK_EQUAL(seq.annotations().size(), 3);
 381   BOOST_CHECK_EQUAL(seq.motifs().size(), 0);
 382   seq.add_motif("CCGTCCC");
 383   BOOST_CHECK_EQUAL(seq.annotations().size(), 3);
 384   BOOST_CHECK_EQUAL(seq.motifs().size(), 1);
 385   seq.clear_motifs();
 386   BOOST_CHECK_EQUAL(seq.annotations().size(), 3);
 387   BOOST_CHECK_EQUAL(seq.motifs().size(), 0);
 388 }
 389
 390 BOOST_AUTO_TEST_CASE( out_operator )
 391 {
 392   string s("AAGGCCTT");
 393   Sequence seq(s);
 394
 395   ostringstream buf;
 396   buf << s;
 397   BOOST_CHECK_EQUAL( s, buf.str() );
 398 }
 399
 400 BOOST_AUTO_TEST_CASE( get_name )
 401 {
 402   Sequence seq("AAGGCCTT");
 403
 404   BOOST_CHECK_EQUAL( seq.get_name(), "" );
 405   seq.set_species("hooman"); // anyone remember tradewars?
 406   BOOST_CHECK_EQUAL( seq.get_name(), "hooman");
 407   seq.set_fasta_header("fasta human");
 408   BOOST_CHECK_EQUAL( seq.get_name(), "fasta human");
 409 }