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_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                       ;
 110   string s(100, 'A');
 111   s += "GCTGCTAATT";
 112   Sequence seq(s);
 113
 114   //istringstream annot_stream(annot_data);
 115   seq.parse_annot(annot_data, 0, 0);
 116   std::list<annot> annots_list = seq.annotations();
 117   std::vector<annot> annots(annots_list.begin(), annots_list.end());
 118   BOOST_REQUIRE_EQUAL( annots.size(), 7);
 119   BOOST_CHECK_EQUAL( annots[0].start, 0 );
 120   BOOST_CHECK_EQUAL( annots[0].end, 10 );
 121   BOOST_CHECK_EQUAL( annots[0].type, "type");
 122   BOOST_CHECK_EQUAL( annots[0].name, "name");
 123   BOOST_CHECK_EQUAL( annots[1].name, "myf7");
 124   BOOST_CHECK_EQUAL( annots[2].name, "myod");
 125   BOOST_CHECK_EQUAL( annots[3].name, "anothername");
 126   BOOST_CHECK_EQUAL( annots[4].name, "backward");
 127   BOOST_CHECK_EQUAL( annots[5].name, "name2");
 128   BOOST_CHECK_EQUAL( annots[5].end, 90);
 129   // sequence defined annotations will always be after the
 130   // absolute positions
 131   BOOST_CHECK_EQUAL( annots[6].name, "ident3 asdf");
 132   BOOST_CHECK_EQUAL( annots[6].start, 100);
 133
 134   //BOOST_CHECK_EQUAL( annots
 135 }
 136
 137 // ticket:83 when you try to load a sequence from a file that doesn't
 138 // have fasta headers it crashes.
 139 BOOST_AUTO_TEST_CASE( sequence_past_end )
 140 {
 141   fs::path seq_path(fs::path(EXAMPLE_DIR)/ "seq" );
 142   seq_path /=  "misformated_seq.fa";
 143   Sequence s;
 144   BOOST_CHECK_THROW( s.load_fasta(seq_path), mussa_load_error );
 145 }
 146
 147 BOOST_AUTO_TEST_CASE ( sequence_empty )
 148 {
 149   Sequence s;
 150   BOOST_CHECK_EQUAL( s.empty(), true );
 151   s = "AAAGGG";
 152   BOOST_CHECK_EQUAL( s.empty(), false );
 153 }
 154
 155 BOOST_AUTO_TEST_CASE ( sequence_iterators )
 156 {
 157   std::string seq_string = "AAGGCCTTNNTATA";
 158   Sequence s(seq_string);
 159   const Sequence cs(s);
 160   std::string::size_type count = 0;
 161
 162   std::string::iterator str_itor;
 163   Sequence::iterator s_itor;
 164   Sequence::const_iterator cs_itor;
 165
 166   for( str_itor = seq_string.begin(),
 167        s_itor   = s.begin(),
 168        cs_itor  = cs.begin();
 169        str_itor != seq_string.end() and
 170        s_itor   != s.end() and
 171        cs_itor  != cs.end();
 172        ++str_itor, ++s_itor, ++cs_itor, ++count)
 173   {
 174     BOOST_CHECK_EQUAL ( *str_itor, *s_itor );
 175     BOOST_CHECK_EQUAL ( *s_itor, *cs_itor );
 176     BOOST_CHECK_EQUAL ( *cs_itor, *str_itor );
 177   }
 178   BOOST_CHECK_EQUAL( seq_string.size(), count );
 179   BOOST_CHECK_EQUAL( s.size(), count );
 180   BOOST_CHECK_EQUAL( cs.size(), count );
 181 }
 182
 183 BOOST_AUTO_TEST_CASE( sequence_motifs )
 184 {
 185   string m("AAAA");
 186   string bogus("AATTGGAA");
 187   Sequence s1("AAAAGGGGCCCCTTTT");
 188
 189   list<motif>::const_iterator motif_i = s1.motifs().begin();
 190   list<motif>::const_iterator motif_end = s1.motifs().end();
 191
 192   // do our iterators work?
 193   BOOST_CHECK( motif_i == s1.motifs().begin() );
 194   BOOST_CHECK( motif_end == s1.motifs().end() );
 195   BOOST_CHECK( motif_i == motif_end );
 196
 197
 198   s1.add_motif(bogus);
 199   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 200   s1.add_motif(m);
 201   BOOST_CHECK( s1.motifs().begin() != s1.motifs().end() );
 202   BOOST_CHECK_EQUAL( s1.motifs().size(), 2 );
 203
 204   for(motif_i = s1.motifs().begin();
 205       motif_i != s1.motifs().end();
 206       ++motif_i)
 207   {
 208     BOOST_CHECK_EQUAL( motif_i->type, "motif" );
 209     BOOST_CHECK_EQUAL( motif_i->name, m);
 210     BOOST_CHECK_EQUAL( motif_i->sequence, m);
 211   }
 212
 213   s1.clear_motifs();
 214   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 215 }
 216
 217 BOOST_AUTO_TEST_CASE( annot_test )
 218 {
 219   annot a(0, 10, "test", "thing");
 220
 221   BOOST_CHECK_EQUAL( a.start, 0 );
 222   BOOST_CHECK_EQUAL( a.end,   10 );
 223   BOOST_CHECK_EQUAL( a.type,  "test" );
 224   BOOST_CHECK_EQUAL( a.name,  "thing" );
 225
 226   motif m(10, "AAGGCC");
 227   BOOST_CHECK_EQUAL( m.start, 10 );
 228   BOOST_CHECK_EQUAL( m.type, "motif" );
 229   BOOST_CHECK_EQUAL( m.name, "AAGGCC" );
 230   BOOST_CHECK_EQUAL( m.end,  10+6 );
 231 }
 232
 233 BOOST_AUTO_TEST_CASE( annotate_from_sequence )
 234 {
 235   string s("CCGCCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 236            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 237            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 238            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 239            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 240            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 241            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 242            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 243   string gc("GCCCCC");
 244   string gga("GGACACCTC");
 245   Sequence seq(s);
 246
 247   std::list<Sequence> query_list;
 248   std::list<string> string_list;
 249   query_list.push_back(Sequence(gc));
 250   string_list.push_back(gc);
 251   query_list.push_back(Sequence(gga));
 252   string_list.push_back(gga);
 253
 254   BOOST_CHECK_EQUAL( seq.annotations().size(), 0 );
 255   seq.find_sequences(query_list.begin(), query_list.end());
 256
 257   int count = 0;
 258   for(list<string>::iterator string_i = string_list.begin();
 259       string_i != string_list.end();
 260       ++string_i)
 261   {
 262     string::size_type pos=0;
 263     while(pos != string::npos) {
 264       pos = s.find(*string_i, pos);
 265       if (pos != string::npos) {
 266         ++count;
 267         ++pos;
 268       }
 269     }
 270   }
 271   BOOST_CHECK_EQUAL(seq.annotations().size(), count);
 272 }
 273
 274 BOOST_AUTO_TEST_CASE( subseq_annotation_test )
 275 {
 276   string s("CCGCCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 277            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 278            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 279            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 280            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 281            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 282            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 283            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 284   Sequence seq(s);
 285
 286
 287   seq.add_annotation(annot(0, 10, "0-10", "0-10"));
 288   seq.add_annotation(annot(10, 20, "10-20", "10-20"));
 289   seq.add_annotation(annot(0, 20, "0-20", "0-20"));
 290   seq.add_annotation(annot(8, 12, "8-12", "8-12"));
 291   seq.add_annotation(annot(100, 5000, "100-5000", "100-5000"));
 292
 293   Sequence subseq = seq.subseq(5, 10);
 294   const list<annot> annots = subseq.annotations();
 295   // generate some ground truth
 296   list<annot> correct;
 297   correct.push_back(annot(0, 5, "0-10",  "0-10"));
 298   correct.push_back(annot(5,10, "10-20", "10-20"));
 299   correct.push_back(annot(0,10, "0-20",  "0-20"));
 300   correct.push_back(annot(3, 7, "8-12",  "8-12"));
 301   BOOST_REQUIRE_EQUAL( annots.size(), correct.size() );
 302
 303   list<annot>::iterator correct_i = correct.begin();
 304   list<annot>::const_iterator annot_i = annots.begin();
 305   for(; annot_i != annots.end(); ++annot_i, ++correct_i)
 306   {
 307     BOOST_CHECK( *annot_i == *correct_i );
 308   }
 309 }
 310
 311 BOOST_AUTO_TEST_CASE( out_operator )
 312 {
 313   string s("AAGGCCTT");
 314   Sequence seq(s);
 315
 316   ostringstream buf;
 317   buf << s;
 318   BOOST_CHECK_EQUAL( s, buf.str() );
 319 }
 320