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