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