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\n" // type\n"
  71                       "10 20 myf\n"
  72                       "20 30 myod\n"
  73                       "50\t55 anothername\n"
  74                       "60 50 backward\n"
  75                       ">ident\n"
  76                       "GCT\n"
  77                       "GCT\n"
  78                       ;
  79   string s('A',100);
  80   s += "GCTGCT";
  81   Sequence seq(s);
  82
  83   //istringstream annot_stream(annot_data);
  84   seq.parse_annot(annot_data, 0, 0);
  85   std::list<annot> annots_list = seq.annotations();
  86   std::vector<annot> annots(annots_list.begin(), annots_list.end());
  87   BOOST_REQUIRE_EQUAL( annots.size(), 5);
  88   BOOST_CHECK_EQUAL( annots[0].start, 0 );
  89   BOOST_CHECK_EQUAL( annots[0].end, 10 );
  90   //BOOST_CHECK_EQUAL( annots[0].type, "type");
  91   BOOST_CHECK_EQUAL( annots[0].name, "name");
  92   //BOOST_CHECK_EQUAL( annots[1].name, "myf7");
  93
  94   //BOOST_CHECK_EQUAL( annots
  95 }
  96
  97 // ticket:83 when you try to load a sequence from a file that doesn't
  98 // have fasta headers it crashes.
  99 BOOST_AUTO_TEST_CASE( sequence_past_end )
 100 {
 101   fs::path seq_path(fs::path(EXAMPLE_DIR)/ "seq" );
 102   seq_path /=  "misformated_seq.fa";
 103   Sequence s;
 104   BOOST_CHECK_THROW( s.load_fasta(seq_path), mussa_load_error );
 105 }
 106
 107 BOOST_AUTO_TEST_CASE ( sequence_empty )
 108 {
 109   Sequence s;
 110   BOOST_CHECK_EQUAL( s.empty(), true );
 111   s = "AAAGGG";
 112   BOOST_CHECK_EQUAL( s.empty(), false );
 113 }
 114
 115 BOOST_AUTO_TEST_CASE ( sequence_iterators )
 116 {
 117   std::string seq_string = "AAGGCCTTNNTATA";
 118   Sequence s(seq_string);
 119   const Sequence cs(s);
 120   std::string::size_type count = 0;
 121
 122   std::string::iterator str_itor;
 123   Sequence::iterator s_itor;
 124   Sequence::const_iterator cs_itor;
 125
 126   for( str_itor = seq_string.begin(),
 127        s_itor   = s.begin(),
 128        cs_itor  = cs.begin();
 129        str_itor != seq_string.end() and
 130        s_itor   != s.end() and
 131        cs_itor  != cs.end();
 132        ++str_itor, ++s_itor, ++cs_itor, ++count)
 133   {
 134     BOOST_CHECK_EQUAL ( *str_itor, *s_itor );
 135     BOOST_CHECK_EQUAL ( *s_itor, *cs_itor );
 136     BOOST_CHECK_EQUAL ( *cs_itor, *str_itor );
 137   }
 138   BOOST_CHECK_EQUAL( seq_string.size(), count );
 139   BOOST_CHECK_EQUAL( s.size(), count );
 140   BOOST_CHECK_EQUAL( cs.size(), count );
 141 }
 142
 143 BOOST_AUTO_TEST_CASE( sequence_motifs )
 144 {
 145   string m("AAAA");
 146   string bogus("AATTGGAA");
 147   Sequence s1("AAAAGGGGCCCCTTTT");
 148
 149   list<motif>::const_iterator motif_i = s1.motifs().begin();
 150   list<motif>::const_iterator motif_end = s1.motifs().end();
 151
 152   // do our iterators work?
 153   BOOST_CHECK( motif_i == s1.motifs().begin() );
 154   BOOST_CHECK( motif_end == s1.motifs().end() );
 155   BOOST_CHECK( motif_i == motif_end );
 156
 157
 158   s1.add_motif(bogus);
 159   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 160   s1.add_motif(m);
 161   BOOST_CHECK( s1.motifs().begin() != s1.motifs().end() );
 162   BOOST_CHECK_EQUAL( s1.motifs().size(), 2 );
 163
 164   for(motif_i = s1.motifs().begin();
 165       motif_i != s1.motifs().end();
 166       ++motif_i)
 167   {
 168     BOOST_CHECK_EQUAL( motif_i->type, "motif" );
 169     BOOST_CHECK_EQUAL( motif_i->name, m);
 170     BOOST_CHECK_EQUAL( motif_i->sequence, m);
 171   }
 172
 173   s1.clear_motifs();
 174   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 175 }
 176
 177 BOOST_AUTO_TEST_CASE( annot_test )
 178 {
 179   annot a(0, 10, "test", "thing");
 180
 181   BOOST_CHECK_EQUAL( a.start, 0 );
 182   BOOST_CHECK_EQUAL( a.end,   10 );
 183   BOOST_CHECK_EQUAL( a.type,  "test" );
 184   BOOST_CHECK_EQUAL( a.name,  "thing" );
 185
 186   motif m(10, "AAGGCC");
 187   BOOST_CHECK_EQUAL( m.start, 10 );
 188   BOOST_CHECK_EQUAL( m.type, "motif" );
 189   BOOST_CHECK_EQUAL( m.name, "AAGGCC" );
 190   BOOST_CHECK_EQUAL( m.end,  10+6 );
 191 }
 192
 193 BOOST_AUTO_TEST_CASE( annotate_from_sequence )
 194 {
 195   string s("CCGCCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 196            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 197            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 198            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 199            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 200            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 201            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 202            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 203   string gc("GCCCCC");
 204   string gga("GGACACCTC");
 205   Sequence seq(s);
 206
 207   std::list<Sequence> query_list;
 208   std::list<string> string_list;
 209   query_list.push_back(Sequence(gc));
 210   string_list.push_back(gc);
 211   query_list.push_back(Sequence(gga));
 212   string_list.push_back(gga);
 213
 214   BOOST_CHECK_EQUAL( seq.annotations().size(), 0 );
 215   seq.find_sequences(query_list.begin(), query_list.end());
 216
 217   int count = 0;
 218   for(list<string>::iterator string_i = string_list.begin();
 219       string_i != string_list.end();
 220       ++string_i)
 221   {
 222     string::size_type pos=0;
 223     while(pos != string::npos) {
 224       pos = s.find(*string_i, pos);
 225       if (pos != string::npos) {
 226         ++count;
 227         ++pos;
 228       }
 229     }
 230   }
 231   BOOST_CHECK_EQUAL(seq.annotations().size(), count);
 232 }
 233