alg/test/test_sequence.cpp

   1 #define BOOST_AUTO_TEST_MAIN
   2 #include <boost/test/auto_unit_test.hpp>
   3 #include <boost/filesystem/path.hpp>
   4 #include <boost/filesystem/operations.hpp>
   5 namespace fs=boost::filesystem;
   6
   7 #include <boost/algorithm/string/case_conv.hpp>
   8
   9 #include <list>
  10 #include <iostream>
  11 #include <sstream>
  12
  13 #include <boost/archive/text_oarchive.hpp>
  14 #include <boost/archive/text_iarchive.hpp>
  15 #include <boost/archive/xml_oarchive.hpp>
  16 #include <boost/archive/xml_iarchive.hpp>
  17
  18 #include "alg/sequence.hpp"
  19 #include "mussa_exceptions.hpp"
  20
  21 using namespace std;
  22
  23 BOOST_AUTO_TEST_CASE( sequence_get_sequence )
  24 {
  25         Sequence s;
  26         // make sure that retrieving the sequence doesn't throw an error
  27         BOOST_CHECK_EQUAL(s.get_sequence(), std::string() );
  28 }
  29
  30 //! when we try to load a missing file, do we get an error?
  31 BOOST_AUTO_TEST_CASE( sequence_load_exception )
  32 {
  33   Sequence s;
  34   // there should be errors when we try to load something that doesn't exist
  35   BOOST_CHECK_THROW( s.load_fasta("alkejralk", 1, 0, 0), mussa_load_error);
  36   BOOST_CHECK_THROW( s.load_annot("alkejralk", 0, 0), mussa_load_error);
  37 }
  38
  39 BOOST_AUTO_TEST_CASE( sequence_eol_conventions )
  40 {
  41   string header(">Header");
  42   string line1("AAAAGGGGCCCCTTTTT");
  43   string line2("AAAAGGGGCCCCTTTTT");
  44   int seq_len = line1.size() + line2.size();
  45
  46   stringstream cr;
  47   cr << header << "\015" << line1 << "\015" << line2 << "\015";
  48   Sequence seq_cr;
  49   seq_cr.load_fasta(cr);
  50
  51   stringstream crlf;
  52   crlf << header << "\015\012" << line1 << "\015\012" << line2 << "\015\012";
  53   Sequence seq_crlf;
  54   seq_crlf.load_fasta(crlf);
  55
  56   stringstream lf;
  57   lf << header << "\012" << line1 << "\012" << line2 << "\012";
  58   Sequence seq_lf;
  59   seq_lf.load_fasta(lf);
  60
  61   BOOST_CHECK_EQUAL(seq_cr.size(),   seq_len);
  62   BOOST_CHECK_EQUAL(seq_crlf.size(), seq_len);
  63   BOOST_CHECK_EQUAL(seq_cr.size(),   seq_len);
  64 }
  65
  66
  67 //! Do simple operations work correctly?
  68 BOOST_AUTO_TEST_CASE( sequence_filter )
  69 {
  70   const char *core_seq = "AATTGGCC";
  71   Sequence s1(core_seq, Sequence::reduced_dna_alphabet);
  72   BOOST_CHECK_EQUAL(s1, core_seq);
  73
  74   Sequence s2("aattggcc", Sequence::reduced_dna_alphabet);
  75   BOOST_CHECK_EQUAL(s2, "AATTGGCC");
  76   BOOST_CHECK_EQUAL(s2.rev_comp(), "GGCCAATT");
  77   BOOST_CHECK_EQUAL(s2.rev_comp(), "ggccaatt"); // we should be case insensitive now
  78   BOOST_CHECK_EQUAL(s2.size(), s2.size());
  79   //We're currently forcing sequences to uppercase
  80   BOOST_CHECK_EQUAL(s2.get_sequence(), "AATTGGCC");
  81
  82   Sequence s3("asdfg", Sequence::reduced_dna_alphabet);
  83   BOOST_CHECK_EQUAL(s3, "ANNNG");
  84   BOOST_CHECK_EQUAL(s3.subseq(0,2), "AN");
  85
  86   s3.set_filtered_sequence("AAGGCCTT", Sequence::reduced_dna_alphabet, 0, 2);
  87   BOOST_CHECK_EQUAL(s3, "AA");
  88   s3.set_filtered_sequence("AAGGCCTT", Sequence::reduced_dna_alphabet, 2, 2);
  89   BOOST_CHECK_EQUAL( s3, "GG");
  90   s3.set_filtered_sequence("AAGGCCTT", Sequence::reduced_dna_alphabet, 4);
  91   BOOST_CHECK_EQUAL( s3, "CCTT");
  92
  93   s3 = "AAGGFF";
  94   BOOST_CHECK_EQUAL(s3, "AAGGNN");
  95 }
  96
  97 BOOST_AUTO_TEST_CASE( sequence_nucleic_alphabet )
  98 {
  99   std::string agct("AGCT");
 100   Sequence seq(agct, Sequence::nucleic_alphabet);
 101   BOOST_CHECK_EQUAL(seq.size(), agct.size());
 102   BOOST_CHECK_EQUAL(seq.get_sequence(), agct);
 103
 104   std::string bdv("BDv");
 105   Sequence seq_bdv(bdv, Sequence::nucleic_alphabet);
 106   BOOST_CHECK_EQUAL(seq_bdv.size(), bdv.size());
 107   // forcing sequence to upper case
 108   BOOST_CHECK_EQUAL(seq_bdv.get_sequence(),
 109                     boost::algorithm::to_upper_copy(bdv));
 110
 111 }
 112
 113 BOOST_AUTO_TEST_CASE( sequence_default_alphabet )
 114 {
 115   std::string agct("AGCT");
 116   Sequence seq(agct);
 117   BOOST_CHECK_EQUAL(seq.size(), agct.size());
 118   BOOST_CHECK_EQUAL(seq.get_sequence(), agct);
 119   BOOST_CHECK_EQUAL(seq[0], agct[0]);
 120   BOOST_CHECK_EQUAL(seq[1], agct[1]);
 121   BOOST_CHECK_EQUAL(seq[2], agct[2]);
 122   BOOST_CHECK_EQUAL(seq[3], agct[3]);
 123
 124   std::string bdv("BDv");
 125   Sequence seq_bdv(bdv);
 126   BOOST_CHECK_EQUAL(seq_bdv.size(), bdv.size());
 127   // default alphabet only allows AGCTUN
 128   BOOST_CHECK_EQUAL(seq_bdv.get_sequence(), "NNN");
 129 }
 130
 131 BOOST_AUTO_TEST_CASE( subseq_names )
 132 {
 133   Sequence s1("AAGGCCTT", Sequence::reduced_dna_alphabet);
 134   s1.set_species("species");
 135   s1.set_fasta_header("a fasta header");
 136   Sequence s2 = s1.subseq(2,2);
 137   BOOST_CHECK_EQUAL(s2, "GG");
 138   BOOST_CHECK_EQUAL(s2.get_species(), s1.get_species());
 139   BOOST_CHECK_EQUAL(s2.get_fasta_header(), s1.get_fasta_header());
 140 }
 141
 142 BOOST_AUTO_TEST_CASE( sequence_start_stop )
 143 {
 144   Sequence s1;
 145   BOOST_CHECK_EQUAL( s1.start(), 0 );
 146   BOOST_CHECK_EQUAL( s1.stop(), 0 );
 147
 148   std::string seq_string("AAGGCCTT");
 149   Sequence s2(seq_string, Sequence::reduced_dna_alphabet);
 150   BOOST_CHECK_EQUAL( s2.start(), 0 );
 151   BOOST_CHECK_EQUAL( s2.stop(), seq_string.size() );
 152
 153   std::string s3seq_string = seq_string.substr(2,3);
 154   Sequence s3 = s2.subseq(2,3);
 155   BOOST_CHECK_EQUAL( s3.start(), 2);
 156   BOOST_CHECK_EQUAL( s3.stop(), 2+3);
 157   BOOST_CHECK_EQUAL( s3.size(), 3);
 158   BOOST_CHECK_EQUAL( s3, s3seq_string);
 159
 160   std::string s4seq_string = s3seq_string.substr(1,1);
 161   Sequence s4 = s3.subseq(1,1);
 162   BOOST_CHECK_EQUAL( s4.start(), 1 );
 163   BOOST_CHECK_EQUAL( s4.stop(), 1+1);
 164   BOOST_CHECK_EQUAL( s4.size(), 1);
 165   BOOST_CHECK_EQUAL( s4, s4seq_string);
 166 }
 167
 168 //! Can we load data from a file
 169 BOOST_AUTO_TEST_CASE( sequence_load )
 170 {
 171   fs::path seq_path(fs::path(EXAMPLE_DIR, fs::native)/ "seq");
 172   seq_path /=  "human_mck_pro.fa";
 173   Sequence s;
 174   s.load_fasta(seq_path, Sequence::reduced_dna_alphabet);
 175   BOOST_CHECK_EQUAL(s.subseq(0, 5), "GGATC"); // first few chars of fasta file
 176   BOOST_CHECK_EQUAL(s.subseq(2, 3), "ATC");
 177   BOOST_CHECK_EQUAL(s.get_fasta_header(), "gi|180579|gb|M21487.1|HUMCKMM1 Human "
 178                                     "muscle creatine kinase gene (CKMM), "
 179                                     "5' flank");
 180 }
 181
 182 BOOST_AUTO_TEST_CASE( sequence_load_fasta_error )
 183 {
 184   fs::path seq_path(fs::path(EXAMPLE_DIR, fs::native)/"seq");
 185   seq_path /= "broken.fa";
 186   bool exception_thrown = false;
 187   std::string exception_filename;
 188   Sequence s;
 189   try {
 190     s.load_fasta(seq_path);
 191   } catch(sequence_invalid_load_error e) {
 192     exception_thrown = true;
 193     size_t native_string_size = seq_path.native_file_string().size();
 194     std:string estr(e.what());
 195     BOOST_REQUIRE(estr.size() > native_string_size);
 196     std::copy(estr.begin(), estr.begin()+native_string_size,
 197               std::back_inserter(exception_filename));
 198   }
 199   BOOST_CHECK_EQUAL(exception_thrown, true);
 200   BOOST_CHECK_EQUAL(seq_path.native_file_string(), exception_filename);
 201 }
 202
 203 BOOST_AUTO_TEST_CASE( sequence_load_annot_error )
 204 {
 205   fs::path seq_path(fs::path(EXAMPLE_DIR, fs::native)/"seq");
 206   seq_path /= "mouse_mck_pro.fa";
 207   fs::path annot_path(fs::path(EXAMPLE_DIR, fs::native));
 208   annot_path /= "broken.annot";
 209   bool exception_thrown = false;
 210   Sequence s;
 211   s.load_fasta(seq_path);
 212
 213   std::string exception_filename;
 214   try {
 215     s.load_annot(annot_path, 0, 0);
 216   } catch(annotation_load_error e) {
 217     exception_thrown = true;
 218     std:string estr(e.what());
 219     size_t native_string_size = annot_path.native_file_string().size();
 220     BOOST_REQUIRE(estr.size() > native_string_size);
 221     std::copy(estr.begin(), estr.begin()+native_string_size,
 222               std::back_inserter(exception_filename));
 223   }
 224   BOOST_CHECK_EQUAL(exception_thrown, true);
 225   BOOST_CHECK_EQUAL(annot_path.native_file_string(), exception_filename);
 226 }
 227
 228 BOOST_AUTO_TEST_CASE( sequence_load_dna_reduced )
 229 {
 230   std::string reduced_dna_fasta_string(">foo\nAAGGCCTTNN\n");
 231   std::stringstream reduced_dna_fasta(reduced_dna_fasta_string);
 232   std::string invalid_dna_fasta_string(">wrong\nAUSSI\n");
 233   std::stringstream invalid_dna_fasta(invalid_dna_fasta_string);
 234   std::string reduced_rna_fasta_string(">foo\nAAGGCCUUNN-\n");
 235   std::stringstream reduced_rna_fasta(reduced_rna_fasta_string);
 236   std::string garbage_string(">foo\na34ralk3547oilk,jual;(*&^#%-\n");
 237   std::stringstream garbage_fasta(garbage_string);
 238
 239   Sequence s;
 240   s.load_fasta(reduced_dna_fasta, Sequence::reduced_dna_alphabet);
 241   BOOST_CHECK_THROW(s.load_fasta(invalid_dna_fasta,
 242                                  Sequence::reduced_dna_alphabet),
 243                     sequence_invalid_load_error);
 244   BOOST_CHECK_THROW(s.load_fasta(reduced_rna_fasta,
 245                                  Sequence::reduced_dna_alphabet),
 246                     sequence_invalid_load_error);
 247   BOOST_CHECK_THROW(s.load_fasta(garbage_fasta,
 248                                  Sequence::reduced_dna_alphabet),
 249                     sequence_invalid_load_error);
 250
 251 }
 252
 253 BOOST_AUTO_TEST_CASE( sequence_load_rna_reduced )
 254 {
 255   std::string reduced_rna_fasta_string(">foo\nAAGGCCUUNN\n");
 256   std::stringstream reduced_rna_fasta(reduced_rna_fasta_string);
 257   std::string invalid_rna_fasta_string(">wrong\nATSSI\n");
 258   std::stringstream invalid_rna_fasta(invalid_rna_fasta_string);
 259   std::string reduced_dna_fasta_string(">foo\nAAGGCCTTNN\n");
 260   std::stringstream reduced_dna_fasta(reduced_dna_fasta_string);
 261   std::string garbage_string(">foo\na34ralk3547oilk,jual;(*&^#%-\n");
 262   std::stringstream garbage_fasta(garbage_string);
 263
 264   Sequence s;
 265   s.load_fasta(reduced_rna_fasta, Sequence::reduced_rna_alphabet);
 266   BOOST_CHECK_THROW(s.load_fasta(invalid_rna_fasta,
 267                                  Sequence::reduced_rna_alphabet),
 268                     sequence_invalid_load_error);
 269   BOOST_CHECK_THROW(s.load_fasta(reduced_dna_fasta,
 270                                  Sequence::reduced_rna_alphabet),
 271                     sequence_invalid_load_error);
 272   BOOST_CHECK_THROW(s.load_fasta(garbage_fasta,
 273                                  Sequence::reduced_rna_alphabet),
 274                     sequence_invalid_load_error);
 275 }
 276
 277 BOOST_AUTO_TEST_CASE( sequence_load_fasta_default )
 278 {
 279   std::string reduced_rna_fasta_string(">foo\nAAGGCCUUNN\n");
 280   std::stringstream reduced_rna_fasta1(reduced_rna_fasta_string);
 281   std::stringstream reduced_rna_fasta2(reduced_rna_fasta_string);
 282   std::string invalid_nucleotide_fasta_string(">wrong\nATSSI\n");
 283   std::stringstream invalid_nucleotide_fasta(invalid_nucleotide_fasta_string);
 284   std::string reduced_dna_fasta_string(">foo\nAAGGCCTTNN\n");
 285   std::stringstream reduced_dna_fasta1(reduced_dna_fasta_string);
 286   std::stringstream reduced_dna_fasta2(reduced_dna_fasta_string);
 287   std::string garbage_string(">foo\na34ralk3547oilk,jual;(*&^#%-\n");
 288   std::stringstream garbage_fasta(garbage_string);
 289
 290   Sequence s;
 291   Sequence specific;
 292   // there's two copies of reduced_rna_fasta because i didn't feel like
 293   // figuring out how to properly reset the read pointer in a stringstream
 294   s.load_fasta(reduced_rna_fasta1);
 295   specific.load_fasta(reduced_rna_fasta2, Sequence::reduced_nucleic_alphabet);
 296   BOOST_CHECK_EQUAL(s, specific);
 297
 298   s.load_fasta(reduced_dna_fasta1);
 299   specific.load_fasta(reduced_dna_fasta2, Sequence::reduced_nucleic_alphabet);
 300   BOOST_CHECK_EQUAL(s, specific);
 301
 302   BOOST_CHECK_THROW(s.load_fasta(invalid_nucleotide_fasta),
 303                     sequence_invalid_load_error);
 304   BOOST_CHECK_THROW(s.load_fasta(garbage_fasta),
 305                     sequence_invalid_load_error);
 306 }
 307
 308 BOOST_AUTO_TEST_CASE( sequence_load_multiple_sequences_one_fasta )
 309 {
 310   std::string fasta_file(
 311     ">gi|10129974|gb|AF188002.1|AF188002\n"
 312     "AAAAGGCTCCTGTCATATTGTGTCCTGCTCTGGTCTGC\n"
 313     ">gi|180579|gb|M21487.1|HUMCKMM1\n"
 314     "CGCCGAGAGCGCTTGCTCTGCCCAGATCTCGGCGAGTC\n"
 315     ">gi|1621|emb|X55146.1|OCMCK1\n"
 316     "CTCCCTGAGGGGAGTGCCCCGCTTAGCCC\n"
 317   );
 318   istringstream seq1_file(fasta_file);
 319   Sequence seq1;
 320   seq1.load_fasta(seq1_file, 1, 0, 0);
 321   BOOST_CHECK_EQUAL(seq1.get_sequence(), "AAAAGGCTCCTGTCATATTGTGTCCTGCTCTGGTCTGC");
 322
 323   istringstream seq2_file(fasta_file);
 324   Sequence seq2;
 325   seq2.load_fasta(seq2_file, 2, 0, 0);
 326   BOOST_CHECK_EQUAL(seq2.get_sequence(), "CGCCGAGAGCGCTTGCTCTGCCCAGATCTCGGCGAGTC");
 327
 328   istringstream seq3_file(fasta_file);
 329   Sequence seq3;
 330   seq3.load_fasta(seq3_file, 3, 0, 0);
 331   BOOST_CHECK_EQUAL(seq3.get_sequence(), "CTCCCTGAGGGGAGTGCCCCGCTTAGCCC");
 332 }
 333
 334 BOOST_AUTO_TEST_CASE( sequence_reverse_complement )
 335 {
 336   std::string iupac_symbols("AaCcGgTtRrYySsWwKkMmBbDdHhVvNn-~.?");
 337   Sequence seq(iupac_symbols, Sequence::nucleic_alphabet);
 338   Sequence seqr = seq.rev_comp();
 339
 340   BOOST_CHECK( seq != seqr );
 341   BOOST_CHECK_EQUAL( seq, seqr.rev_comp() );
 342   // forcing sequence to upper case
 343   BOOST_CHECK_EQUAL( seq.get_sequence(),
 344                      boost::algorithm::to_upper_copy(iupac_symbols) );
 345 }
 346
 347 BOOST_AUTO_TEST_CASE( sequence_reverse_complement_dna )
 348 {
 349   std::string dna_str("AGCTN");
 350   Sequence dna_seq(dna_str, Sequence::reduced_dna_alphabet);
 351   BOOST_CHECK_EQUAL(dna_seq.rev_comp(), "NAGCT");
 352   BOOST_CHECK_EQUAL(dna_seq, dna_seq.rev_comp().rev_comp());
 353 }
 354
 355 BOOST_AUTO_TEST_CASE( sequence_reverse_complement_rna )
 356 {
 357   std::string rna_str("AGCUN");
 358   Sequence rna_seq(rna_str, Sequence::reduced_rna_alphabet);
 359   BOOST_CHECK_EQUAL(rna_seq.rev_comp(), "NAGCU");
 360   BOOST_CHECK_EQUAL(rna_seq, rna_seq.rev_comp().rev_comp());
 361 }
 362
 363 BOOST_AUTO_TEST_CASE( sequence_reverse_complement_subseq )
 364 {
 365   std::string dna_str("AAAAAAAAAAGGGGGGGGGGG");
 366   Sequence seq(dna_str, Sequence::reduced_dna_alphabet);
 367   Sequence subseq = seq.subseq(8,4);
 368   BOOST_CHECK_EQUAL( subseq, "AAGG");
 369   Sequence rev_subseq = subseq.rev_comp();
 370   BOOST_CHECK_EQUAL( rev_subseq.size(), subseq.size() );
 371   BOOST_CHECK_EQUAL( rev_subseq.get_sequence(), "CCTT");
 372 }
 373
 374 BOOST_AUTO_TEST_CASE( sequence_reverse_iterator )
 375 {
 376   std::string dna_str("AAAAAAAAAAGGGGGGGGGGG");
 377   std::string dna_str_reversed(dna_str.rbegin(), dna_str.rend());
 378   Sequence seq(dna_str);
 379   std::string seq_reversed(seq.rbegin(), seq.rend());
 380   BOOST_CHECK_EQUAL(seq_reversed, dna_str_reversed);
 381
 382   std::string substr = dna_str.substr(8,4);
 383   Sequence subseq = seq.subseq(8,4);
 384   BOOST_CHECK_EQUAL(substr, subseq);
 385
 386   std::string substr_reversed(substr.rbegin(), substr.rend());
 387   std::string subseq_reversed(subseq.rbegin(), subseq.rend());
 388   BOOST_CHECK_EQUAL(substr_reversed, subseq_reversed);
 389 }
 390
 391 BOOST_AUTO_TEST_CASE( sequence_empty_reverse_iterator)
 392 {
 393   // so what happens with reverse interators when we have no sequence?
 394   Sequence seq1;
 395   Sequence seq2;
 396   Sequence seq3("AGCT");
 397
 398   // all the empty sequences should have equal iterators
 399   BOOST_CHECK(seq1.rbegin() == seq1.rend());
 400   BOOST_CHECK(seq1.rbegin() == seq2.rend());
 401
 402   // none of the seq1 iterators should equal any of the seq3 iterators
 403   BOOST_CHECK(seq1.rbegin() != seq3.rbegin());
 404   BOOST_CHECK(seq1.rbegin() != seq3.rend());
 405   BOOST_CHECK(seq1.rend() != seq3.rbegin());
 406   BOOST_CHECK(seq1.rend() != seq3.rend());
 407
 408   // seq3 iterators should work
 409   BOOST_CHECK(seq3.rbegin() != seq3.rend());
 410
 411 }
 412
 413 BOOST_AUTO_TEST_CASE( annotation_load )
 414 {
 415   string annot_data = "human\n"
 416                       "0 10 name   type\n"
 417                       "10 20 myf7\n"
 418                       "20 30 myod\n"
 419                       "50\t55 anothername\n"
 420                       "60 50 backward\n"
 421                       ">ident3 asdf\n"
 422                       "GCT\n"
 423                       "gCTn\n"
 424                       "75\t90\tname2\ttype2\n"
 425                       "100 120 name-asdf type!@#$%\n"
 426                       ;
 427   string s(100, 'A');
 428   s += "GCTGCTAATT";
 429   Sequence seq(s, Sequence::reduced_dna_alphabet);
 430
 431   //istringstream annot_stream(annot_data);
 432   seq.parse_annot(annot_data, 0, 0);
 433   std::list<annot> annots_list = seq.annotations();
 434   std::vector<annot> annots(annots_list.begin(), annots_list.end());
 435   BOOST_REQUIRE_EQUAL( annots.size(), 8);
 436   BOOST_CHECK_EQUAL( annots[0].begin, 0 );
 437   BOOST_CHECK_EQUAL( annots[0].end, 10 );
 438   BOOST_CHECK_EQUAL( annots[0].type, "type");
 439   BOOST_CHECK_EQUAL( annots[0].name, "name");
 440   BOOST_CHECK_EQUAL( annots[1].name, "myf7");
 441   BOOST_CHECK_EQUAL( annots[2].name, "myod");
 442   BOOST_CHECK_EQUAL( annots[3].name, "anothername");
 443   BOOST_CHECK_EQUAL( annots[4].name, "backward");
 444   BOOST_CHECK_EQUAL( annots[5].name, "name2");
 445   BOOST_CHECK_EQUAL( annots[5].end, 90);
 446   BOOST_CHECK_EQUAL( annots[6].begin, 100);
 447   BOOST_CHECK_EQUAL( annots[6].end, 120);
 448   BOOST_CHECK_EQUAL( annots[6].name, "name-asdf");
 449   BOOST_CHECK_EQUAL( annots[6].type, "type!@#$%");
 450   // sequence defined annotations will always be after the
 451   // absolute positions
 452   BOOST_CHECK_EQUAL( annots[7].name, "ident3 asdf");
 453   BOOST_CHECK_EQUAL( annots[7].begin, 100);
 454
 455   //BOOST_CHECK_EQUAL( annots
 456 }
 457
 458 BOOST_AUTO_TEST_CASE( annotation_broken_load )
 459 {
 460   string annot_data = "human\n"
 461                       "0 10 name   type\n"
 462                       "blah60 50 backward\n"
 463                       ">ident3 asdf\n"
 464                       "GCT\n"
 465                       "gCTn\n"
 466                       ;
 467   string s(100, 'A');
 468   s += "GCTGCTAATT";
 469   Sequence seq(s, Sequence::reduced_dna_alphabet);
 470
 471   BOOST_CHECK_THROW(seq.parse_annot(annot_data, 0, 0), annotation_load_error);
 472   BOOST_CHECK_EQUAL(seq.annotations().size(), 0);
 473   }
 474
 475 BOOST_AUTO_TEST_CASE(annotation_ucsc_html_load)
 476 {
 477   // this actually is basically what's returned by UCSC
 478   // (well actually with some of the sequence and copies of fasta blocks
 479   // removed to make the example shorter
 480   string annot_data = "\n"
 481     "<PRE>\n"
 482     ">hg17_knownGene_NM_001824_0 range=chr19:50517919-50517974 5'pad=0 3'pad=0 revComp=TRUE strand=- repeatMasking=none\n"
 483     "GGGTCAGTGTCACCTCCAGGATACAGACAG\n"
 484     "&gt;hg17_knownGene_NM_001824_3 range=chr19:50510563-50510695 5'pad=0 3'pad=0 revComp=TRUE strand=- repeatMasking=none\n"
 485     "GGTGGAGACGACCTGGACCCTAACTACGT\n"
 486     "</PRE>\n"
 487     "\n"
 488     "</BODY>\n"
 489     "</HTML>\n"
 490     ;
 491
 492   string s =
 493     "TGGGTCAGTGTCACCTCCAGGATACAGACAGCCCCCCTTCAGCCCAGCCCAGCCAG"
 494     "AAAAA"
 495     "GGTGGAGACGACCTGGACCCTAACTACGTGCTCAGCAGCCGCGTCCGCAC";
 496   Sequence seq(s, Sequence::reduced_dna_alphabet);
 497   seq.parse_annot(annot_data);
 498   std::list<annot> annots = seq.annotations();
 499   BOOST_CHECK_EQUAL( annots.size(), 2);
 500 }
 501
 502 BOOST_AUTO_TEST_CASE( annotation_load_no_species_name )
 503 {
 504   string annot_data = "0 10 name   type\n"
 505                       "10 20 myf7\n"
 506                       "20 30 myod\n"
 507                       "50\t55 anothername\n"
 508                       "60 50 backward\n"
 509                       ">ident3 asdf\n"
 510                       "GCT\n"
 511                       "gCTn\n"
 512                       "75\t90\tname2\ttype2\n"
 513                       "100 120 name-asdf type!@#$%\n"
 514                       ;
 515   string s(100, 'A');
 516   s += "GCTGCTAATT";
 517   Sequence seq(s, Sequence::reduced_dna_alphabet);
 518
 519   //istringstream annot_stream(annot_data);
 520   seq.parse_annot(annot_data, 0, 0);
 521   std::list<annot> annots_list = seq.annotations();
 522   std::vector<annot> annots(annots_list.begin(), annots_list.end());
 523   BOOST_REQUIRE_EQUAL( annots.size(), 8);
 524   BOOST_CHECK_EQUAL( annots[0].begin, 0 );
 525   BOOST_CHECK_EQUAL( annots[0].end, 10 );
 526   BOOST_CHECK_EQUAL( annots[0].type, "type");
 527 }
 528
 529 // ticket:83 when you try to load a sequence from a file that doesn't
 530 // have fasta headers it crashes.
 531 BOOST_AUTO_TEST_CASE( sequence_past_end )
 532 {
 533   fs::path seq_path(fs::path(EXAMPLE_DIR, fs::native)/ "seq" );
 534   seq_path /=  "misformated_seq.fa";
 535   Sequence s;
 536   BOOST_CHECK_THROW( s.load_fasta(seq_path), mussa_load_error );
 537 }
 538
 539 BOOST_AUTO_TEST_CASE ( sequence_empty )
 540 {
 541
 542   Sequence s;
 543   BOOST_CHECK_EQUAL( s.empty(), true );
 544   s = "AAAGGG";
 545   BOOST_CHECK_EQUAL( s.empty(), false );
 546   s.clear();
 547   BOOST_CHECK_EQUAL( s.empty(), true);
 548   s = "";
 549   BOOST_CHECK_EQUAL( s.empty(), true);
 550 }
 551
 552 BOOST_AUTO_TEST_CASE ( sequence_size )
 553 {
 554
 555   Sequence s;
 556   BOOST_CHECK_EQUAL( s.size(), 0);
 557   std::string seq_string("AAAGGG");
 558   s = seq_string;
 559   BOOST_CHECK_EQUAL( s.size(), seq_string.size() );
 560   s.clear();
 561   BOOST_CHECK_EQUAL( s.size(), 0);
 562   s = "";
 563   BOOST_CHECK_EQUAL( s.size(), 0);
 564 }
 565
 566 BOOST_AUTO_TEST_CASE( sequence_empty_equality )
 567 {
 568   Sequence szero("", Sequence::reduced_dna_alphabet);
 569   BOOST_CHECK_EQUAL(szero.empty(), true);
 570   BOOST_CHECK_EQUAL(szero, szero);
 571   BOOST_CHECK_EQUAL(szero, "");
 572
 573   Sequence sclear("AGCT", Sequence::reduced_dna_alphabet);
 574   sclear.clear();
 575   BOOST_CHECK_EQUAL(sclear.empty(), true);
 576   BOOST_CHECK_EQUAL(sclear, sclear);
 577   BOOST_CHECK_EQUAL(sclear, szero);
 578   BOOST_CHECK_EQUAL(sclear, "");
 579
 580 }
 581 BOOST_AUTO_TEST_CASE ( sequence_iterators )
 582 {
 583   std::string seq_string = "AAGGCCTTNNTATA";
 584   Sequence s(seq_string, Sequence::reduced_dna_alphabet);
 585   const Sequence cs(s);
 586   std::string::size_type count = 0;
 587
 588   std::string::iterator str_itor;
 589   Sequence::const_iterator s_itor;
 590   Sequence::const_iterator cs_itor;
 591
 592   for( str_itor = seq_string.begin(),
 593        s_itor   = s.begin(),
 594        cs_itor  = cs.begin();
 595        str_itor != seq_string.end() and
 596        s_itor   != s.end() and
 597        cs_itor  != cs.end();
 598        ++str_itor, ++s_itor, ++cs_itor, ++count)
 599   {
 600     BOOST_CHECK_EQUAL ( *str_itor, *s_itor );
 601     BOOST_CHECK_EQUAL ( *s_itor, *cs_itor );
 602     BOOST_CHECK_EQUAL ( *cs_itor, *str_itor );
 603   }
 604   BOOST_CHECK_EQUAL( seq_string.size(), count );
 605   BOOST_CHECK_EQUAL( s.size(), count );
 606   BOOST_CHECK_EQUAL( cs.size(), count );
 607 }
 608
 609 BOOST_AUTO_TEST_CASE( sequence_motifs )
 610 {
 611   string m("AAAA");
 612   string bogus("AATTGGAA");
 613   Sequence s1("AAAAGGGGCCCCTTTT", Sequence::reduced_dna_alphabet);
 614
 615   list<motif>::const_iterator motif_i = s1.motifs().begin();
 616   list<motif>::const_iterator motif_end = s1.motifs().end();
 617
 618   // do our iterators work?
 619   BOOST_CHECK( motif_i == s1.motifs().begin() );
 620   BOOST_CHECK( motif_end == s1.motifs().end() );
 621   BOOST_CHECK( motif_i == motif_end );
 622
 623   // this shouldn't show up
 624   s1.add_motif(bogus);
 625   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 626   BOOST_CHECK_EQUAL( s1.motifs().size(), 0 );
 627
 628   s1.add_motif(m);
 629   BOOST_CHECK( s1.motifs().begin() != s1.motifs().end() );
 630   BOOST_CHECK_EQUAL( s1.motifs().size(), 2 );
 631
 632   for(motif_i = s1.motifs().begin();
 633       motif_i != s1.motifs().end();
 634       ++motif_i)
 635   {
 636     BOOST_CHECK_EQUAL( motif_i->type, "motif" );
 637     BOOST_CHECK_EQUAL( motif_i->name, m);
 638     BOOST_CHECK_EQUAL( motif_i->sequence, m);
 639   }
 640
 641   s1.clear_motifs();
 642   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 643
 644   /* FIXME: enable this when i find a way of passing storing the motif name
 645   // does our annotation travel?
 646   Sequence motif_seq(m);
 647   motif_seq.set_fasta_header("hi");
 648   s1.add_motif(motif_seq);
 649
 650   BOOST_CHECK_EQUAL(s1.motifs().size(), 2);
 651   for(motif_i = s1.motifs().begin();
 652       motif_i != s1.motifs().end();
 653       ++motif_i)
 654   {
 655     BOOST_CHECK_EQUAL( motif_i->type, "motif" );
 656     BOOST_CHECK_EQUAL( motif_i->name, "hi");
 657     BOOST_CHECK_EQUAL( motif_i->sequence, m);
 658   }
 659   */
 660 }
 661
 662 BOOST_AUTO_TEST_CASE( sequence_motif_subseq)
 663 {
 664   // when searching for a motif on a subsequence we should
 665   // only search the subsequence ticket:199
 666   string aaaa("AAAA");
 667   string cccc("CCCC");
 668   Sequence s1("AAAANCCCC", Sequence::reduced_dna_alphabet);
 669
 670   // this shouldn't show up
 671   s1.add_motif(cccc);
 672   BOOST_CHECK_EQUAL( s1.motifs().size(), 1 );
 673
 674   s1.add_motif(aaaa);
 675   BOOST_CHECK_EQUAL( s1.motifs().size(), 2 );
 676
 677   Sequence subseq1 = s1.subseq(4,5);
 678   BOOST_CHECK_EQUAL(subseq1.motifs().size(), 2);
 679   subseq1.clear_motifs();
 680   BOOST_CHECK_EQUAL(subseq1.motifs().size(), 0);
 681   // this is outside of our subsequence, and so shouldn't be found
 682   subseq1.add_motif(aaaa);
 683   BOOST_CHECK_EQUAL( subseq1.motifs().size(), 0 );
 684
 685   subseq1.add_motif(cccc);
 686   BOOST_CHECK_EQUAL( subseq1.motifs().size(), 1);
 687   std::list<motif>::const_iterator motif_i = subseq1.motifs().begin();
 688   BOOST_REQUIRE(motif_i != subseq1.motifs().end());
 689   BOOST_CHECK_EQUAL(motif_i->begin, 1);
 690   BOOST_CHECK_EQUAL(motif_i->end, 5);
 691 }
 692
 693 BOOST_AUTO_TEST_CASE( annot_test )
 694 {
 695   annot a(0, 10, "test", "thing");
 696
 697   BOOST_CHECK_EQUAL( a.begin, 0 );
 698   BOOST_CHECK_EQUAL( a.end,   10 );
 699   BOOST_CHECK_EQUAL( a.type,  "test" );
 700   BOOST_CHECK_EQUAL( a.name,  "thing" );
 701
 702   motif m(10, "AAGGCC");
 703   BOOST_CHECK_EQUAL( m.begin, 10 );
 704   BOOST_CHECK_EQUAL( m.type, "motif" );
 705   BOOST_CHECK_EQUAL( m.name, "AAGGCC" );
 706   BOOST_CHECK_EQUAL( m.end,  10+6 );
 707 }
 708
 709 BOOST_AUTO_TEST_CASE( annotate_from_sequence )
 710 {
 711   string s("CCGCCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 712            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 713            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 714            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 715            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 716            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 717            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 718            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 719   string gc("GCCCCC");
 720   string gga("GGACACCTC");
 721   Sequence seq(s, Sequence::reduced_dna_alphabet);
 722
 723   std::list<Sequence> query_list;
 724   std::list<string> string_list;
 725   query_list.push_back(Sequence(gc));
 726   string_list.push_back(gc);
 727   query_list.push_back(Sequence(gga));
 728   string_list.push_back(gga);
 729
 730   BOOST_CHECK_EQUAL( seq.annotations().size(), 0 );
 731   seq.find_sequences(query_list.begin(), query_list.end());
 732
 733   int count = 0;
 734   for(list<string>::iterator string_i = string_list.begin();
 735       string_i != string_list.end();
 736       ++string_i)
 737   {
 738     string::size_type pos=0;
 739     while(pos != string::npos) {
 740       pos = s.find(*string_i, pos);
 741       if (pos != string::npos) {
 742         ++count;
 743         ++pos;
 744       }
 745     }
 746   }
 747   BOOST_CHECK_EQUAL(seq.annotations().size(), count);
 748   const std::list<annot> &a = seq.annotations();
 749   for (std::list<annot>::const_iterator annot_i = a.begin();
 750        annot_i != a.end();
 751        ++annot_i)
 752   {
 753     int count = annot_i->end - annot_i->begin ;
 754   }
 755 }
 756
 757 BOOST_AUTO_TEST_CASE( subseq_annotation_test )
 758 {
 759   string s("CCGCCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 760            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 761            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 762            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 763            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 764            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 765            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 766            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 767   Sequence seq(s, Sequence::reduced_dna_alphabet);
 768
 769
 770   seq.add_annotation(annot(0, 10, "0-10", "0-10"));
 771   seq.add_annotation(annot(10, 20, "10-20", "10-20"));
 772   seq.add_annotation(annot(0, 20, "0-20", "0-20"));
 773   seq.add_annotation(annot(8, 12, "8-12", "8-12"));
 774   seq.add_annotation(annot(100, 5000, "100-5000", "100-5000"));
 775
 776   Sequence subseq = seq.subseq(5, 10);
 777   const list<annot> annots = subseq.annotations();
 778   // generate some ground truth
 779   list<annot> correct;
 780   correct.push_back(annot(0, 5, "0-10",  "0-10"));
 781   correct.push_back(annot(5,10, "10-20", "10-20"));
 782   correct.push_back(annot(0,10, "0-20",  "0-20"));
 783   correct.push_back(annot(3, 7, "8-12",  "8-12"));
 784   BOOST_REQUIRE_EQUAL( annots.size(), correct.size() );
 785
 786   list<annot>::iterator correct_i = correct.begin();
 787   list<annot>::const_iterator annot_i = annots.begin();
 788   for(; annot_i != annots.end(); ++annot_i, ++correct_i)
 789   {
 790     BOOST_CHECK( *annot_i == *correct_i );
 791   }
 792 }
 793
 794 BOOST_AUTO_TEST_CASE( motif_annotation_update )
 795 {
 796   string s("CCGTCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 797            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 798            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 799            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 800            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 801            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 802            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 803            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 804   Sequence seq(s, Sequence::reduced_dna_alphabet);
 805
 806   // starting conditions
 807   BOOST_CHECK_EQUAL(seq.annotations().size(), 0);
 808   BOOST_CHECK_EQUAL(seq.motifs().size(), 0);
 809   seq.add_annotation(annot(0, 10, "0-10", "0-10"));
 810   seq.add_annotation(annot(10, 20, "10-20", "10-20"));
 811   seq.add_annotation(annot(0, 20, "0-20", "0-20"));
 812   BOOST_CHECK_EQUAL(seq.annotations().size(), 3);
 813   BOOST_CHECK_EQUAL(seq.motifs().size(), 0);
 814   seq.add_motif("CCGTCCC");
 815   BOOST_CHECK_EQUAL(seq.annotations().size(), 3);
 816   BOOST_CHECK_EQUAL(seq.motifs().size(), 1);
 817   seq.clear_motifs();
 818   BOOST_CHECK_EQUAL(seq.annotations().size(), 3);
 819   BOOST_CHECK_EQUAL(seq.motifs().size(), 0);
 820 }
 821
 822 BOOST_AUTO_TEST_CASE( out_operator )
 823 {
 824   string s("AAGGCCTT");
 825   Sequence seq(s, Sequence::reduced_dna_alphabet);
 826
 827   ostringstream buf;
 828   buf << s;
 829   BOOST_CHECK_EQUAL( s, buf.str() );
 830 }
 831
 832 BOOST_AUTO_TEST_CASE( get_name )
 833 {
 834   Sequence seq("AAGGCCTT", Sequence::reduced_dna_alphabet);
 835
 836   BOOST_CHECK_EQUAL( seq.get_name(), "" );
 837   seq.set_species("hooman"); // anyone remember tradewars?
 838   BOOST_CHECK_EQUAL( seq.get_name(), "hooman");
 839   seq.set_fasta_header("fasta human");
 840   BOOST_CHECK_EQUAL( seq.get_name(), "fasta human");
 841 }
 842
 843 BOOST_AUTO_TEST_CASE( serialize_simple )
 844 {
 845   std::string seq_string = "AAGGCCTT";
 846   Sequence seq(seq_string, Sequence::reduced_dna_alphabet);
 847   seq.set_species("ribbet");
 848   std::ostringstream oss;
 849   // allocate/deallocate serialization components
 850   {
 851     boost::archive::text_oarchive oarchive(oss);
 852     const Sequence& const_seq(seq);
 853     BOOST_CHECK_EQUAL(seq, const_seq);
 854     oarchive << const_seq;
 855   }
 856   Sequence seq_loaded;
 857   {
 858     std::istringstream iss(oss.str());
 859     boost::archive::text_iarchive iarchive(iss);
 860     iarchive >> seq_loaded;
 861   }
 862   BOOST_CHECK_EQUAL(seq_loaded, seq);
 863   BOOST_CHECK_EQUAL(seq.get_species(), "ribbet");
 864 }
 865
 866 BOOST_AUTO_TEST_CASE( serialize_tree )
 867 {
 868   std::string seq_string = "AAGGCCTT";
 869   Sequence seq(seq_string, Sequence::reduced_dna_alphabet);
 870   seq.set_species("ribbet");
 871   seq.add_motif("AA");
 872   seq.add_motif("GC");
 873   annot a1(6,7,"t","t");
 874   seq.add_annotation(a1);
 875
 876   std::ostringstream oss;
 877   // allocate/deallocate serialization components
 878   {
 879     boost::archive::text_oarchive oarchive(oss);
 880     const Sequence& const_seq(seq);
 881     BOOST_CHECK_EQUAL(seq, const_seq);
 882     oarchive << const_seq;
 883   }
 884
 885   Sequence seq_loaded;
 886   {
 887     std::istringstream iss(oss.str());
 888     boost::archive::text_iarchive iarchive(iss);
 889     iarchive >> seq_loaded;
 890   }
 891   BOOST_CHECK_EQUAL(seq_loaded, seq);
 892 }
 893
 894 // this writes out an "old" style annotated sequence
 895 // with annotations attached as "motifs" and "annots"
 896 BOOST_AUTO_TEST_CASE( serialize_xml_sequence )
 897 {
 898   std::string seq_string = "AAGGCCTT";
 899   Sequence seq(seq_string, Sequence::reduced_dna_alphabet);
 900   seq.set_species("ribbet");
 901   seq.add_motif("AA");
 902   seq.add_motif("GC");
 903   annot a1(6,7,"t","t");
 904   seq.add_annotation(a1);
 905
 906   std::ostringstream oss;
 907   // allocate/deallocate serialization components
 908   {
 909     boost::archive::xml_oarchive oarchive(oss);
 910     const Sequence& const_seq(seq);
 911     BOOST_CHECK_EQUAL(seq, const_seq);
 912     oarchive << boost::serialization::make_nvp("root", const_seq);
 913   }
 914   Sequence seq_loaded;
 915   {
 916     std::istringstream iss(oss.str());
 917     boost::archive::xml_iarchive iarchive(iss);
 918     iarchive >> boost::serialization::make_nvp("root", seq_loaded);
 919   }
 920   BOOST_CHECK_EQUAL(seq_loaded, seq);
 921 }
 922
 923 BOOST_AUTO_TEST_CASE( serialize_xml_two )
 924 {
 925   std::string seq_string = "AAGGCCTT";
 926   Sequence seq1(seq_string, Sequence::reduced_dna_alphabet);
 927   Sequence seq2(seq1);
 928
 929   std::ostringstream oss;
 930   // allocate/deallocate serialization components
 931   {
 932     boost::archive::xml_oarchive oarchive(oss);
 933     const Sequence& const_seq1(seq1);
 934     const Sequence& const_seq2(seq2);
 935     oarchive << boost::serialization::make_nvp("seq1", const_seq1);
 936     oarchive << boost::serialization::make_nvp("seq2", const_seq2);
 937   }
 938   //std::cout << "xml: " << oss.str() << std::endl;
 939   Sequence seq1_loaded;
 940   Sequence seq2_loaded;
 941   {
 942     std::istringstream iss(oss.str());
 943     boost::archive::xml_iarchive iarchive(iss);
 944     iarchive >> boost::serialization::make_nvp("seq1", seq1_loaded);
 945     iarchive >> boost::serialization::make_nvp("seq2", seq2_loaded);
 946   }
 947   BOOST_CHECK_EQUAL(seq1_loaded, seq1);
 948   BOOST_CHECK_EQUAL(seq2_loaded, seq2);
 949   // test if our pointers are the same
 950   BOOST_CHECK_EQUAL(seq1_loaded.c_str(), seq2_loaded.c_str());
 951 }