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 //! Do simple operations work correctly?
  25 BOOST_AUTO_TEST_CASE( sequence_filter )
  26 {
  27   Sequence s1("AATTGGCC");
  28   BOOST_CHECK_EQUAL(s1.get_seq(), "AATTGGCC");
  29
  30   Sequence s2("aattggcc");
  31   BOOST_CHECK_EQUAL(s2.get_seq(), "AATTGGCC");
  32   BOOST_CHECK_EQUAL(s2.rev_comp(), "GGCCAATT");
  33   BOOST_CHECK_EQUAL(s2.size(), s2.get_seq().size());
  34   BOOST_CHECK_EQUAL(s2.c_seq(), s2.get_seq().c_str());
  35
  36   Sequence s3("asdfg");
  37   BOOST_CHECK_EQUAL(s3.get_seq(), "ANNNG");
  38   BOOST_CHECK_EQUAL(s3.subseq(0,2), "AN");
  39
  40   s3.set_filtered_sequence("AAGGCCTT", 0, 2);
  41   BOOST_CHECK_EQUAL(s3.get_seq(), "AA");
  42   s3.set_filtered_sequence("AAGGCCTT", 2, 2);
  43   BOOST_CHECK_EQUAL( s3.get_seq(), "GG");
  44   s3.set_filtered_sequence("AAGGCCTT", 4);
  45   BOOST_CHECK_EQUAL( s3.get_seq(), "CCTT");
  46
  47   s3.clear();
  48   BOOST_CHECK_EQUAL(s3.get_seq(), "");
  49
  50   s3.set_seq("AAGGFF");
  51   BOOST_CHECK_EQUAL(s3.get_seq(), "AAGGNN");
  52 }
  53
  54 //! Can we load data from a file
  55 BOOST_AUTO_TEST_CASE( sequence_load )
  56 {
  57   fs::path seq_path(fs::path(EXAMPLE_DIR)/ "seq" );
  58   seq_path /=  "human_mck_pro.fa";
  59   Sequence s;
  60   s.load_fasta(seq_path);
  61   BOOST_CHECK_EQUAL(s.subseq(0, 5), "GGATC"); // first few chars of fasta file
  62   BOOST_CHECK_EQUAL(s.get_header(), "gi|180579|gb|M21487.1|HUMCKMM1 Human "
  63                                     "muscle creatine kinase gene (CKMM), "
  64                                     "5' flank");
  65 }
  66
  67 BOOST_AUTO_TEST_CASE( annotation_load )
  68 {
  69   string annot_data = "human\n"
  70                       "0 10 name type\n"
  71                       "10 20 name\n"
  72                       "20 30\n"
  73                       "15 20 backward\n";
  74   string s('A',100);
  75   s += "GCTGCT";
  76   Sequence seq(s);
  77
  78   //istringstream annot_stream(annot_data);
  79   seq.parse_annot(annot_data, 0, 0);
  80   typedef std::list<annot> annot_list_t;
  81   annot_list_t annots = seq.annotations();
  82   for(annot_list_t::iterator annot_i = annots.begin();
  83       annot_i != annots.end();
  84       ++annot_i)
  85   {
  86     std::cout << "start " << annot_i->start << endl;
  87   }
  88 }
  89
  90 // ticket:83 when you try to load a sequence from a file that doesn't
  91 // have fasta headers it crashes.
  92 BOOST_AUTO_TEST_CASE( sequence_past_end )
  93 {
  94   fs::path seq_path(fs::path(EXAMPLE_DIR)/ "seq" );
  95   seq_path /=  "misformated_seq.fa";
  96   Sequence s;
  97   BOOST_CHECK_THROW( s.load_fasta(seq_path), mussa_load_error );
  98 }
  99
 100 BOOST_AUTO_TEST_CASE ( sequence_empty )
 101 {
 102   Sequence s;
 103   BOOST_CHECK_EQUAL( s.empty(), true );
 104   s = "AAAGGG";
 105   BOOST_CHECK_EQUAL( s.empty(), false );
 106 }
 107
 108 BOOST_AUTO_TEST_CASE ( sequence_iterators )
 109 {
 110   std::string seq_string = "AAGGCCTTNNTATA";
 111   Sequence s(seq_string);
 112   const Sequence cs(s);
 113   std::string::size_type count = 0;
 114
 115   std::string::iterator str_itor;
 116   Sequence::iterator s_itor;
 117   Sequence::const_iterator cs_itor;
 118
 119   for( str_itor = seq_string.begin(),
 120        s_itor   = s.begin(),
 121        cs_itor  = cs.begin();
 122        str_itor != seq_string.end() and
 123        s_itor   != s.end() and
 124        cs_itor  != cs.end();
 125        ++str_itor, ++s_itor, ++cs_itor, ++count)
 126   {
 127     BOOST_CHECK_EQUAL ( *str_itor, *s_itor );
 128     BOOST_CHECK_EQUAL ( *s_itor, *cs_itor );
 129     BOOST_CHECK_EQUAL ( *cs_itor, *str_itor );
 130   }
 131   BOOST_CHECK_EQUAL( seq_string.size(), count );
 132   BOOST_CHECK_EQUAL( s.size(), count );
 133   BOOST_CHECK_EQUAL( cs.size(), count );
 134 }
 135
 136 BOOST_AUTO_TEST_CASE( sequence_motifs )
 137 {
 138   string m("AAAA");
 139   string bogus("AATTGGAA");
 140   Sequence s1("AAAAGGGGCCCCTTTT");
 141
 142   list<motif>::const_iterator motif_i = s1.motifs().begin();
 143   list<motif>::const_iterator motif_end = s1.motifs().end();
 144
 145   // do our iterators work?
 146   BOOST_CHECK( motif_i == s1.motifs().begin() );
 147   BOOST_CHECK( motif_end == s1.motifs().end() );
 148   BOOST_CHECK( motif_i == motif_end );
 149
 150
 151   s1.add_motif(bogus);
 152   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 153   s1.add_motif(m);
 154   BOOST_CHECK( s1.motifs().begin() != s1.motifs().end() );
 155   BOOST_CHECK_EQUAL( s1.motifs().size(), 2 );
 156
 157   for(motif_i = s1.motifs().begin();
 158       motif_i != s1.motifs().end();
 159       ++motif_i)
 160   {
 161     BOOST_CHECK_EQUAL( motif_i->type, "motif" );
 162     BOOST_CHECK_EQUAL( motif_i->name, m);
 163     BOOST_CHECK_EQUAL( motif_i->sequence, m);
 164   }
 165
 166   s1.clear_motifs();
 167   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 168 }
 169
 170 BOOST_AUTO_TEST_CASE( annot_test )
 171 {
 172   annot a(0, 10, "test", "thing");
 173
 174   BOOST_CHECK_EQUAL( a.start, 0 );
 175   BOOST_CHECK_EQUAL( a.end,   10 );
 176   BOOST_CHECK_EQUAL( a.type,  "test" );
 177   BOOST_CHECK_EQUAL( a.name,  "thing" );
 178
 179   motif m(10, "AAGGCC");
 180   BOOST_CHECK_EQUAL( m.start, 10 );
 181   BOOST_CHECK_EQUAL( m.type, "motif" );
 182   BOOST_CHECK_EQUAL( m.name, "AAGGCC" );
 183   BOOST_CHECK_EQUAL( m.end,  10+6 );
 184 }
 185
 186 BOOST_AUTO_TEST_CASE( annotate_from_sequence )
 187 {
 188   string s("CCGCCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 189            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 190            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 191            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 192            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 193            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 194            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 195            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 196   string gc("GCCCCC");
 197   string gga("GGACACCTC");
 198   Sequence seq(s);
 199
 200   std::list<Sequence> query_list;
 201   std::list<string> string_list;
 202   query_list.push_back(Sequence(gc));
 203   string_list.push_back(gc);
 204   query_list.push_back(Sequence(gga));
 205   string_list.push_back(gga);
 206
 207   BOOST_CHECK_EQUAL( seq.annotations().size(), 0 );
 208   seq.find_sequences(query_list.begin(), query_list.end());
 209
 210   int count = 0;
 211   for(list<string>::iterator string_i = string_list.begin();
 212       string_i != string_list.end();
 213       ++string_i)
 214   {
 215     string::size_type pos=0;
 216     while(pos != string::npos) {
 217       pos = s.find(*string_i, pos);
 218       if (pos != string::npos) {
 219         ++count;
 220         ++pos;
 221       }
 222     }
 223   }
 224   BOOST_CHECK_EQUAL(seq.annotations().size(), count);
 225 }
 226