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   Sequence s1("AATTGGCC");
  56   BOOST_CHECK_EQUAL(s1.get_seq(), "AATTGGCC");
  57
  58   Sequence s2("aattggcc");
  59   BOOST_CHECK_EQUAL(s2.get_seq(), "AATTGGCC");
  60   BOOST_CHECK_EQUAL(s2.rev_comp(), "GGCCAATT");
  61   BOOST_CHECK_EQUAL(s2.size(), s2.get_seq().size());
  62   BOOST_CHECK_EQUAL(s2.c_seq(), s2.get_seq().c_str());
  63
  64   Sequence s3("asdfg");
  65   BOOST_CHECK_EQUAL(s3.get_seq(), "ANNNG");
  66   BOOST_CHECK_EQUAL(s3.subseq(0,2), "AN");
  67
  68   s3.set_filtered_sequence("AAGGCCTT", 0, 2);
  69   BOOST_CHECK_EQUAL(s3.get_seq(), "AA");
  70   s3.set_filtered_sequence("AAGGCCTT", 2, 2);
  71   BOOST_CHECK_EQUAL( s3.get_seq(), "GG");
  72   s3.set_filtered_sequence("AAGGCCTT", 4);
  73   BOOST_CHECK_EQUAL( s3.get_seq(), "CCTT");
  74
  75   s3.clear();
  76   BOOST_CHECK_EQUAL(s3.get_seq(), "");
  77
  78   s3.set_seq("AAGGFF");
  79   BOOST_CHECK_EQUAL(s3.get_seq(), "AAGGNN");
  80 }
  81
  82 //! Can we load data from a file
  83 BOOST_AUTO_TEST_CASE( sequence_load )
  84 {
  85   fs::path seq_path(fs::path(EXAMPLE_DIR)/ "seq" );
  86   seq_path /=  "human_mck_pro.fa";
  87   Sequence s;
  88   s.load_fasta(seq_path);
  89   BOOST_CHECK_EQUAL(s.subseq(0, 5), "GGATC"); // first few chars of fasta file
  90   BOOST_CHECK_EQUAL(s.get_header(), "gi|180579|gb|M21487.1|HUMCKMM1 Human "
  91                                     "muscle creatine kinase gene (CKMM), "
  92                                     "5' flank");
  93 }
  94
  95 BOOST_AUTO_TEST_CASE( annotation_load )
  96 {
  97   string annot_data = "human\n"
  98                       "0 10 name   type\n"
  99                       "10 20 myf7\n"
 100                       "20 30 myod\n"
 101                       "50\t55 anothername\n"
 102                       "60 50 backward\n"
 103                       ">ident3 asdf\n"
 104                       "GCT\n"
 105                       "gCTn\n"
 106                       "75\t90\tname2\ttype2\n"
 107                       ;
 108   string s(100, 'A');
 109   s += "GCTGCTAATT";
 110   Sequence seq(s);
 111
 112   //istringstream annot_stream(annot_data);
 113   seq.parse_annot(annot_data, 0, 0);
 114   std::list<annot> annots_list = seq.annotations();
 115   std::vector<annot> annots(annots_list.begin(), annots_list.end());
 116   BOOST_REQUIRE_EQUAL( annots.size(), 7);
 117   BOOST_CHECK_EQUAL( annots[0].start, 0 );
 118   BOOST_CHECK_EQUAL( annots[0].end, 10 );
 119   BOOST_CHECK_EQUAL( annots[0].type, "type");
 120   BOOST_CHECK_EQUAL( annots[0].name, "name");
 121   BOOST_CHECK_EQUAL( annots[1].name, "myf7");
 122   BOOST_CHECK_EQUAL( annots[2].name, "myod");
 123   BOOST_CHECK_EQUAL( annots[3].name, "anothername");
 124   BOOST_CHECK_EQUAL( annots[4].name, "backward");
 125   BOOST_CHECK_EQUAL( annots[5].name, "name2");
 126   BOOST_CHECK_EQUAL( annots[5].end, 90);
 127   // sequence defined annotations will always be after the
 128   // absolute positions
 129   BOOST_CHECK_EQUAL( annots[6].name, "ident3 asdf");
 130   BOOST_CHECK_EQUAL( annots[6].start, 100);
 131
 132   //BOOST_CHECK_EQUAL( annots
 133 }
 134
 135 // ticket:83 when you try to load a sequence from a file that doesn't
 136 // have fasta headers it crashes.
 137 BOOST_AUTO_TEST_CASE( sequence_past_end )
 138 {
 139   fs::path seq_path(fs::path(EXAMPLE_DIR)/ "seq" );
 140   seq_path /=  "misformated_seq.fa";
 141   Sequence s;
 142   BOOST_CHECK_THROW( s.load_fasta(seq_path), mussa_load_error );
 143 }
 144
 145 BOOST_AUTO_TEST_CASE ( sequence_empty )
 146 {
 147   Sequence s;
 148   BOOST_CHECK_EQUAL( s.empty(), true );
 149   s = "AAAGGG";
 150   BOOST_CHECK_EQUAL( s.empty(), false );
 151 }
 152
 153 BOOST_AUTO_TEST_CASE ( sequence_iterators )
 154 {
 155   std::string seq_string = "AAGGCCTTNNTATA";
 156   Sequence s(seq_string);
 157   const Sequence cs(s);
 158   std::string::size_type count = 0;
 159
 160   std::string::iterator str_itor;
 161   Sequence::iterator s_itor;
 162   Sequence::const_iterator cs_itor;
 163
 164   for( str_itor = seq_string.begin(),
 165        s_itor   = s.begin(),
 166        cs_itor  = cs.begin();
 167        str_itor != seq_string.end() and
 168        s_itor   != s.end() and
 169        cs_itor  != cs.end();
 170        ++str_itor, ++s_itor, ++cs_itor, ++count)
 171   {
 172     BOOST_CHECK_EQUAL ( *str_itor, *s_itor );
 173     BOOST_CHECK_EQUAL ( *s_itor, *cs_itor );
 174     BOOST_CHECK_EQUAL ( *cs_itor, *str_itor );
 175   }
 176   BOOST_CHECK_EQUAL( seq_string.size(), count );
 177   BOOST_CHECK_EQUAL( s.size(), count );
 178   BOOST_CHECK_EQUAL( cs.size(), count );
 179 }
 180
 181 BOOST_AUTO_TEST_CASE( sequence_motifs )
 182 {
 183   string m("AAAA");
 184   string bogus("AATTGGAA");
 185   Sequence s1("AAAAGGGGCCCCTTTT");
 186
 187   list<motif>::const_iterator motif_i = s1.motifs().begin();
 188   list<motif>::const_iterator motif_end = s1.motifs().end();
 189
 190   // do our iterators work?
 191   BOOST_CHECK( motif_i == s1.motifs().begin() );
 192   BOOST_CHECK( motif_end == s1.motifs().end() );
 193   BOOST_CHECK( motif_i == motif_end );
 194
 195
 196   s1.add_motif(bogus);
 197   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 198   s1.add_motif(m);
 199   BOOST_CHECK( s1.motifs().begin() != s1.motifs().end() );
 200   BOOST_CHECK_EQUAL( s1.motifs().size(), 2 );
 201
 202   for(motif_i = s1.motifs().begin();
 203       motif_i != s1.motifs().end();
 204       ++motif_i)
 205   {
 206     BOOST_CHECK_EQUAL( motif_i->type, "motif" );
 207     BOOST_CHECK_EQUAL( motif_i->name, m);
 208     BOOST_CHECK_EQUAL( motif_i->sequence, m);
 209   }
 210
 211   s1.clear_motifs();
 212   BOOST_CHECK( s1.motifs().begin() == s1.motifs().end() );
 213 }
 214
 215 BOOST_AUTO_TEST_CASE( annot_test )
 216 {
 217   annot a(0, 10, "test", "thing");
 218
 219   BOOST_CHECK_EQUAL( a.start, 0 );
 220   BOOST_CHECK_EQUAL( a.end,   10 );
 221   BOOST_CHECK_EQUAL( a.type,  "test" );
 222   BOOST_CHECK_EQUAL( a.name,  "thing" );
 223
 224   motif m(10, "AAGGCC");
 225   BOOST_CHECK_EQUAL( m.start, 10 );
 226   BOOST_CHECK_EQUAL( m.type, "motif" );
 227   BOOST_CHECK_EQUAL( m.name, "AAGGCC" );
 228   BOOST_CHECK_EQUAL( m.end,  10+6 );
 229 }
 230
 231 BOOST_AUTO_TEST_CASE( annotate_from_sequence )
 232 {
 233   string s("CCGCCCCCCATCATCGCGGCTCTCCGAGAGTCCCGCGCCCCACTCCCGGC"
 234            "ACCCACCTGACCGCGGGCGGCTCCGGCCCCGCTTCGCCCCACTGCGATCA"
 235            "GTCGCGTCCCGCAGGCCAGGCACGCCCCGCCGCTCCCGCTGCGCCGGGCG"
 236            "TCTGGGACCTCGGGCGGCTCCTCCGAGGGGCGGGGCAGCCGGGAGCCACG"
 237            "CCCCCGCAGGTGAGCCGGCCACGCCCACCGCCCGTGGGAAGTTCAGCCTC"
 238            "GGGGCTCCAGCCCCGCGGGAATGGCAGAACTTCGCACGCGGAACTGGTAA"
 239            "CCTCCAGGACACCTCGAATCAGGGTGATTGTAGCGCAGGGGCCTTGGCCA"
 240            "AGCTAAAACTTTGGAAACTTTAGATCCCAGACAGGTGGCTTTCTTGCAGT");
 241   string gc("GCCCCC");
 242   string gga("GGACACCTC");
 243   Sequence seq(s);
 244
 245   std::list<Sequence> query_list;
 246   std::list<string> string_list;
 247   query_list.push_back(Sequence(gc));
 248   string_list.push_back(gc);
 249   query_list.push_back(Sequence(gga));
 250   string_list.push_back(gga);
 251
 252   BOOST_CHECK_EQUAL( seq.annotations().size(), 0 );
 253   seq.find_sequences(query_list.begin(), query_list.end());
 254
 255   int count = 0;
 256   for(list<string>::iterator string_i = string_list.begin();
 257       string_i != string_list.end();
 258       ++string_i)
 259   {
 260     string::size_type pos=0;
 261     while(pos != string::npos) {
 262       pos = s.find(*string_i, pos);
 263       if (pos != string::npos) {
 264         ++count;
 265         ++pos;
 266       }
 267     }
 268   }
 269   BOOST_CHECK_EQUAL(seq.annotations().size(), count);
 270 }
 271