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