alg/sequence.cpp

   1 //  This file is part of the Mussa source distribution.
   2 //  http://mussa.caltech.edu/
   3 //  Contact author: Tristan  De Buysscher, tristan@caltech.edu
   4
   5 // This program and all associated source code files are Copyright (C) 2005
   6 // the California Institute of Technology, Pasadena, CA, 91125 USA.  It is
   7 // under the GNU Public License; please see the included LICENSE.txt
   8 // file for more information, or contact Tristan directly.
   9
  10
  11 //  This file is part of the Mussa source distribution.
  12 //  http://mussa.caltech.edu/
  13 //  Contact author: Tristan  De Buysscher, tristan@caltech.edu
  14
  15 // This program and all associated source code files are Copyright (C) 2005
  16 // the California Institute of Technology, Pasadena, CA, 91125 USA.  It is
  17 // under the GNU Public License; please see the included LICENSE.txt
  18 // file for more information, or contact Tristan directly.
  19
  20
  21 //                        ----------------------------------------
  22 //                           ---------- sequence.cc -----------
  23 //                        ----------------------------------------
  24 #include <boost/filesystem/fstream.hpp>
  25 #include <boost/filesystem/operations.hpp>
  26 namespace fs = boost::filesystem;
  27
  28 #include <boost/spirit/core.hpp>
  29 #include <boost/spirit/actor/push_back_actor.hpp>
  30 #include <boost/spirit/iterator/file_iterator.hpp>
  31 #include <boost/spirit/utility/chset.hpp>
  32 namespace spirit = boost::spirit;
  33
  34 #include "alg/sequence.hpp"
  35 #include "mussa_exceptions.hpp"
  36
  37 #include <string>
  38 #include <stdexcept>
  39 #include <iostream>
  40 #include <sstream>
  41 #include <set>
  42
  43 annot::annot()
  44  : begin(0),
  45    end(0),
  46    type(""),
  47    name("")
  48 {
  49 }
  50
  51 annot::annot(int begin, int end, std::string type, std::string name)
  52  : begin(begin),
  53    end(end),
  54    type(type),
  55    name(name)
  56 {
  57 }
  58
  59 annot::~annot()
  60 {
  61 }
  62
  63 bool operator==(const annot& left, const annot& right)
  64 {
  65   return ((left.begin== right.begin) and
  66           (left.end == right.end) and
  67           (left.type == right.type) and
  68           (left.name == right.name));
  69 }
  70
  71 motif::motif(int begin, std::string motif)
  72  : annot(begin, begin+motif.size(), "motif", motif),
  73    sequence(motif)
  74 {
  75 }
  76
  77 motif::~motif()
  78 {
  79 }
  80
  81
  82 Sequence::Sequence(alphabet_ref alphabet_)
  83   : seq(new SeqSpan("")),
  84     alphabet(alphabet_),
  85     strand(UnknownStrand)
  86 {
  87 }
  88
  89 Sequence::~Sequence()
  90 {
  91 }
  92
  93 Sequence::Sequence(const char *seq, alphabet_ref alphabet_)
  94   : alphabet(alphabet_),
  95     strand(UnknownStrand),
  96     header(""),
  97     species("")
  98 {
  99   set_filtered_sequence(seq, alphabet);
 100 }
 101
 102 Sequence::Sequence(const std::string& seq, alphabet_ref alphabet_)
 103   : alphabet(alphabet_),
 104     strand(UnknownStrand),
 105     header(""),
 106     species("")
 107 {
 108   set_filtered_sequence(seq, alphabet);
 109 }
 110
 111 Sequence::Sequence(const Sequence& o)
 112   : seq(o.seq),
 113     alphabet(o.alphabet),
 114     strand(o.strand),
 115     header(o.header),
 116     species(o.species),
 117     annots(o.annots),
 118     motif_list(o.motif_list)
 119 {
 120 }
 121
 122 Sequence::Sequence(const Sequence* o)
 123   : seq(o->seq),
 124     alphabet(o->alphabet),
 125     strand(o->strand),
 126     header(o->header),
 127     species(o->species),
 128     annots(o->annots),
 129     motif_list(o->motif_list)
 130 {
 131 }
 132
 133 Sequence::Sequence(const SeqSpanRef& seq_ref, alphabet_ref alphabet_)
 134   : seq(seq_ref),
 135     alphabet(alphabet),
 136     strand(UnknownStrand),
 137     header(""),
 138     species("")
 139 {
 140 }
 141
 142 Sequence &Sequence::operator=(const Sequence& s)
 143 {
 144   if (this != &s) {
 145     seq = s.seq;
 146     alphabet = s.alphabet;
 147     strand = s.strand;
 148     header = s.header;
 149     species = s.species;
 150     annots = s.annots;
 151     motif_list = s.motif_list;
 152   }
 153   return *this;
 154 }
 155
 156 static void multiplatform_getline(std::istream& in, std::string& line)
 157 {
 158   line.clear();
 159   char c;
 160   in.get(c);
 161   while(in.good() and !(c == '\012' or c == '\015') ) {
 162     line.push_back(c);
 163     in.get(c);
 164   }
 165   // if we have cr-lf eat it
 166   c = in.peek();
 167   if (c=='\012' or c == '\015') {
 168     in.get();
 169   }
 170 }
 171
 172 void Sequence::load_fasta(fs::path file_path, int seq_num, int start_index, int end_index)
 173 {
 174   load_fasta(file_path, alphabet, seq_num, start_index, end_index);
 175 }
 176
 177 //! load a fasta file into a sequence
 178 void Sequence::load_fasta(fs::path file_path, alphabet_ref a,
 179                           int seq_num, int start_index, int end_index)
 180 {
 181   fs::fstream data_file;
 182   data_file.open(file_path, std::ios::in);
 183
 184   if (!data_file.good())
 185   {
 186     throw mussa_load_error("Sequence File: "+file_path.string()+" not found");
 187   } else {
 188     try {
 189       load_fasta(data_file, a, seq_num, start_index, end_index);
 190     } catch(sequence_empty_error e) {
 191       // there doesn't appear to be any sequence
 192       // catch and rethrow to include the filename
 193       std::stringstream msg;
 194       msg << "The selected sequence in "
 195           << file_path.native_file_string()
 196           << " appears to be empty";
 197       throw sequence_empty_error(msg.str());
 198     } catch(sequence_empty_file_error e) {
 199       std::stringstream errormsg;
 200       errormsg << file_path.native_file_string()
 201                << " did not have any fasta sequences" << std::endl;
 202       throw sequence_empty_file_error(errormsg.str());
 203     } catch(sequence_invalid_load_error e) {
 204       std::ostringstream msg;
 205       msg << file_path.native_file_string();
 206       msg << " " << e.what();
 207       throw sequence_invalid_load_error(msg.str());
 208     }
 209   }
 210 }
 211
 212 void Sequence::load_fasta(std::istream& file,
 213                           int seq_num, int start_index, int end_index)
 214 {
 215   load_fasta(file, alphabet, seq_num, start_index, end_index);
 216 }
 217
 218 void
 219 Sequence::load_fasta(std::istream& data_file, alphabet_ref a,
 220                      int seq_num,
 221                      int start_index, int end_index)
 222 {
 223   std::string file_data_line;
 224   int header_counter = 0;
 225   size_t line_counter = 0;
 226   bool read_seq = true;
 227   std::string rev_comp;
 228   std::string sequence_raw;
 229   std::string seq_tmp;      // holds sequence during basic filtering
 230   const Alphabet &alpha = get_alphabet(a);
 231
 232   if (seq_num == 0) {
 233     throw mussa_load_error("fasta sequence number is 1 based (can't be 0)");
 234   }
 235
 236   // search for the header of the fasta sequence we want
 237   while ( (!data_file.eof()) && (header_counter < seq_num) )
 238   {
 239     multiplatform_getline(data_file, file_data_line);
 240     ++line_counter;
 241     if (file_data_line.substr(0,1) == ">")
 242       header_counter++;
 243   }
 244
 245   if (header_counter > 0) {
 246     header = file_data_line.substr(1);
 247
 248     sequence_raw = "";
 249
 250     while ( !data_file.eof() && read_seq ) {
 251       multiplatform_getline(data_file,file_data_line);
 252       ++line_counter;
 253       if (file_data_line.substr(0,1) == ">")
 254         read_seq = false;
 255       else {
 256         for (std::string::const_iterator line_i = file_data_line.begin();
 257              line_i != file_data_line.end();
 258              ++line_i)
 259          {
 260            if(alpha.exists(*line_i)) {
 261              sequence_raw += *line_i;
 262            } else {
 263             std::ostringstream msg;
 264             msg << "Unrecognized characters in fasta sequence at line ";
 265             msg << line_counter;
 266             throw sequence_invalid_load_error(msg.str());
 267            }
 268          }
 269       }
 270     }
 271
 272     // Lastly, if subselection of the sequence was specified we keep cut out
 273     // and only keep that part
 274     // end_index = 0 means no end was specified, so cut to the end
 275     if (end_index == 0)
 276       end_index = sequence_raw.size();
 277
 278     // sequence filtering for upcasing agctn and convert non AGCTN to N
 279     if (end_index-start_index <= 0) {
 280       std::string msg("The selected sequence appears to be empty");
 281       throw sequence_empty_error(msg);
 282     }
 283     set_filtered_sequence(sequence_raw, a, start_index, end_index-start_index);
 284   } else {
 285     std::string errormsg("There were no fasta sequences");
 286     throw sequence_empty_file_error(errormsg);
 287   }
 288 }
 289
 290 void Sequence::set_filtered_sequence(const std::string &in_seq,
 291                                      alphabet_ref alphabet_,
 292                                      size_type start,
 293                                      size_type count,
 294                                      strand_type strand_)
 295 {
 296   alphabet = alphabet_;
 297   if ( count == npos)
 298     count = in_seq.size() - start;
 299   std::string new_seq;
 300   new_seq.reserve(count);
 301
 302   // finally, the actual conversion loop
 303   const Alphabet& alpha_impl = get_alphabet(); // go get one of our actual alphabets
 304   std::string::const_iterator seq_i = in_seq.begin()+start;
 305   for(size_type i = 0; i != count; ++i, ++seq_i)
 306   {
 307     if (alpha_impl.exists(*seq_i)) {
 308       new_seq.append(1, toupper(*seq_i));
 309     } else {
 310       new_seq.append(1, 'N');
 311     }
 312   }
 313   SeqSpanRef new_seq_ref(new SeqSpan(new_seq));
 314   seq = new_seq_ref;
 315   strand = strand_;
 316 }
 317
 318 void
 319 Sequence::load_annot(fs::path file_path, int start_index, int end_index)
 320 {
 321   if (not fs::exists(file_path)) {
 322     throw mussa_load_error("Annotation File " + file_path.string() + " was not found");
 323   }
 324   if (fs::is_directory(file_path)) {
 325     throw mussa_load_error(file_path.string() +
 326             " is a directory, please provide a file for annotations."
 327           );
 328   }
 329   fs::fstream data_stream(file_path, std::ios::in);
 330   if (!data_stream)
 331   {
 332     throw mussa_load_error("Error loading annotation file " + file_path.string());
 333   }
 334
 335   // so i should probably be passing the parse function some iterators
 336   // but the annotations files are (currently) small, so i think i can
 337   // get away with loading the whole file into memory
 338   std::string data;
 339   char c;
 340   while(data_stream.good()) {
 341     data_stream.get(c);
 342     data.push_back(c);
 343   }
 344   data_stream.close();
 345
 346   try {
 347     parse_annot(data, start_index, end_index);
 348   } catch(annotation_load_error e) {
 349     std::ostringstream msg;
 350     msg << file_path.native_file_string()
 351         << " "
 352         << e.what();
 353     throw annotation_load_error(msg.str());
 354   }
 355 }
 356
 357 /* If this works, yikes, this is some brain hurting code.
 358  *
 359  * what's going on is that when pb_annot is instantiated it stores references
 360  * to begin, end, name, type, declared in the parse function, then
 361  * when operator() is called it grabs values from those references
 362  * and uses that to instantiate an annot object and append that to our
 363  * annotation list.
 364  *
 365  * This weirdness is because the spirit library requires that actions
 366  * conform to a specific prototype operator()(IteratorT, IteratorT)
 367  * which doesn't provide any useful opportunity for me to actually
 368  * grab the results of our parsing.
 369  *
 370  * so I instantiate this structure in order to have a place to grab
 371  * my data from.
 372  */
 373
 374 struct push_back_annot {
 375   std::list<annot>& annot_list;
 376   int& begin;
 377   int& end;
 378   std::string& name;
 379   std::string& type;
 380   int &parsed;
 381
 382   push_back_annot(std::list<annot>& annot_list_,
 383                   int& begin_,
 384                   int& end_,
 385                   std::string& name_,
 386                   std::string& type_,
 387                   int &parsed_)
 388   : annot_list(annot_list_),
 389     begin(begin_),
 390     end(end_),
 391     name(name_),
 392     type(type_),
 393     parsed(parsed_)
 394   {
 395   }
 396
 397   void operator()(std::string::const_iterator,
 398                   std::string::const_iterator) const
 399   {
 400     //std::cout << "adding annot: " << begin << "|" << end << "|" << name << "|" << type << std::endl;
 401     annot_list.push_back(annot(begin, end, name, type));
 402     ++parsed;
 403   };
 404 };
 405
 406 struct push_back_seq {
 407   std::list<Sequence>& seq_list;
 408   std::string& name;
 409   std::string& seq;
 410   int &parsed;
 411
 412   push_back_seq(std::list<Sequence>& seq_list_,
 413                 std::string& name_,
 414                 std::string& seq_,
 415                 int &parsed_)
 416   : seq_list(seq_list_),
 417     name(name_),
 418     seq(seq_),
 419     parsed(parsed_)
 420   {
 421   }
 422
 423   void operator()(std::string::const_iterator,
 424                   std::string::const_iterator) const
 425   {
 426     // filter out newlines from our sequence
 427     std::string new_seq;
 428     for(std::string::const_iterator seq_i = seq.begin();
 429         seq_i != seq.end();
 430         ++seq_i)
 431     {
 432       if (*seq_i != '\015' && *seq_i != '\012') new_seq += *seq_i;
 433     }
 434     //std::cout << "adding seq: " << name << " " << new_seq << std::endl;
 435
 436     Sequence s(new_seq);
 437     s.set_fasta_header(name);
 438     seq_list.push_back(s);
 439     ++parsed;
 440   };
 441 };
 442
 443 void
 444 Sequence::parse_annot(std::string data, int start_index, int end_index)
 445 {
 446   int start=0;
 447   int end=0;
 448   std::string name;
 449   std::string type;
 450   std::string seq;
 451   std::list<annot> parsed_annots;
 452   std::list<Sequence> query_seqs;
 453   int parsed=0;
 454
 455   bool ok = spirit::parse(data.begin(), data.end(),
 456               (
 457                //begin grammar
 458                  !(
 459                     (
 460                       spirit::alpha_p >>
 461                       +(spirit::graph_p)
 462                     )[spirit::assign_a(species)] >>
 463                     +(spirit::space_p)
 464                   ) >>
 465                   *(
 466                      ( // ignore html tags
 467                        *(spirit::space_p) >>
 468                        spirit::ch_p('<') >>
 469                        +(~spirit::ch_p('>')) >>
 470                        spirit::ch_p('>') >>
 471                        *(spirit::space_p)
 472                      )
 473                    |
 474                     ( // parse an absolute location name
 475                      (spirit::uint_p[spirit::assign_a(start)] >>
 476                       +spirit::space_p >>
 477                       spirit::uint_p[spirit::assign_a(end)] >>
 478                       +spirit::space_p >>
 479                       (
 480                          spirit::alpha_p >>
 481                          *spirit::graph_p
 482                       )[spirit::assign_a(name)] >>
 483                       // optional type
 484                       !(
 485                           +spirit::space_p >>
 486                           (
 487                             spirit::alpha_p >>
 488                             *spirit::graph_p
 489                           )[spirit::assign_a(type)]
 490                       )
 491                       // to understand how this group gets set
 492                       // read the comment above struct push_back_annot
 493                      )[push_back_annot(parsed_annots, start, end, type, name, parsed)]
 494                    |
 495                     ((spirit::ch_p('>')|spirit::str_p("&gt;")) >>
 496                        (*(spirit::print_p))[spirit::assign_a(name)] >>
 497                        spirit::eol_p >>
 498                        (+(spirit::chset<>(Alphabet::nucleic_cstr)))[spirit::assign_a(seq)]
 499                      )[push_back_seq(query_seqs, name, seq, parsed)]
 500                     ) >>
 501                     *spirit::space_p
 502                    )
 503               //end grammar
 504               )).full;
 505   if (not ok) {
 506     std::stringstream msg;
 507     msg << "Error parsing annotation #" << parsed;
 508     throw annotation_load_error(msg.str());
 509   }
 510   // add newly parsed annotations to our sequence
 511   std::copy(parsed_annots.begin(), parsed_annots.end(), std::back_inserter(annots));
 512   // go seearch for query sequences
 513   find_sequences(query_seqs.begin(), query_seqs.end());
 514 }
 515
 516 void Sequence::add_annotation(const annot& a)
 517 {
 518   annots.push_back(a);
 519 }
 520
 521 const std::list<annot>& Sequence::annotations() const
 522 {
 523   return annots;
 524 }
 525
 526 void Sequence::copy_children(Sequence &new_seq, size_type start, size_type count) const
 527 {
 528   new_seq.motif_list = motif_list;
 529   new_seq.annots.clear();
 530
 531   for(std::list<annot>::const_iterator annot_i = annots.begin();
 532       annot_i != annots.end();
 533       ++annot_i)
 534   {
 535     size_type annot_begin= annot_i->begin;
 536     size_type annot_end = annot_i->end;
 537
 538     if (annot_begin < start+count) {
 539       if (annot_begin >= start) {
 540         annot_begin -= start;
 541       } else {
 542         annot_begin = 0;
 543       }
 544
 545       if (annot_end < start+count) {
 546         annot_end -= start;
 547       } else {
 548         annot_end = count;
 549       }
 550
 551       annot new_annot(annot_begin, annot_end, annot_i->type, annot_i->name);
 552       new_seq.annots.push_back(new_annot);
 553     }
 554   }
 555
 556 }
 557 Sequence
 558 Sequence::subseq(size_type start, size_type count) const
 559 {
 560   if (!seq) {
 561     Sequence new_seq;
 562     return new_seq;
 563   }
 564
 565   Sequence new_seq = *this;
 566   new_seq.seq = seq->subseq(start, count);
 567
 568   copy_children(new_seq, start, count);
 569
 570   return new_seq;
 571 }
 572
 573 std::string Sequence::create_reverse_map() const
 574 {
 575   std::string rc_map(256, '~');
 576   // if we're rna, use U instead of T
 577   // we might want to add an "is_rna" to sequence at somepoint
 578   char TU = (alphabet == reduced_rna_alphabet) ? 'U' : 'T';
 579   char tu = (alphabet == reduced_rna_alphabet) ? 'u' : 't';
 580   rc_map['A'] = TU ; rc_map['a'] = tu ;
 581   rc_map['T'] = 'A'; rc_map['t'] = 'a';
 582   rc_map['U'] = 'A'; rc_map['u'] = 'a';
 583   rc_map['G'] = 'C'; rc_map['g'] = 'c';
 584   rc_map['C'] = 'G'; rc_map['c'] = 'g';
 585   rc_map['M'] = 'K'; rc_map['m'] = 'k';
 586   rc_map['R'] = 'Y'; rc_map['r'] = 'y';
 587   rc_map['W'] = 'W'; rc_map['w'] = 'w';
 588   rc_map['S'] = 'S'; rc_map['s'] = 's';
 589   rc_map['Y'] = 'R'; rc_map['y'] = 'r';
 590   rc_map['K'] = 'M'; rc_map['k'] = 'm';
 591   rc_map['V'] = 'B'; rc_map['v'] = 'b';
 592   rc_map['H'] = 'D'; rc_map['h'] = 'd';
 593   rc_map['D'] = 'H'; rc_map['d'] = 'h';
 594   rc_map['B'] = 'V'; rc_map['b'] = 'v';
 595   rc_map['N'] = 'N'; rc_map['n'] = 'n';
 596   rc_map['X'] = 'X'; rc_map['x'] = 'x';
 597   rc_map['?'] = '?';
 598   rc_map['.'] = '.';
 599   rc_map['-'] = '-';
 600   rc_map['~'] = '~'; // not really needed, but perhaps it's clearer.
 601   return rc_map;
 602 }
 603
 604 Sequence Sequence::rev_comp() const
 605 {
 606   std::string rev_comp;
 607   rev_comp.reserve(length());
 608
 609   std::string rc_map = create_reverse_map();
 610
 611   // reverse and convert
 612   Sequence::const_reverse_iterator seq_i;
 613   Sequence::const_reverse_iterator seq_end = rend();
 614   for(seq_i = rbegin();
 615       seq_i != seq_end;
 616       ++seq_i)
 617   {
 618     rev_comp.append(1, rc_map[*seq_i]);
 619   }
 620   return Sequence(rev_comp, alphabet);
 621 }
 622
 623 void Sequence::set_fasta_header(std::string header_)
 624 {
 625   header = header_;
 626 }
 627
 628 void Sequence::set_species(const std::string& name)
 629 {
 630   species = name;
 631 }
 632
 633 std::string Sequence::get_species() const
 634 {
 635   return species;
 636 }
 637
 638
 639 std::string
 640 Sequence::get_fasta_header() const
 641 {
 642   return header;
 643 }
 644
 645 std::string
 646 Sequence::get_name() const
 647 {
 648   if (header.size() > 0)
 649     return header;
 650   else if (species.size() > 0)
 651     return species;
 652   else
 653     return "";
 654 }
 655
 656 const Alphabet& Sequence::get_alphabet() const
 657 {
 658   return get_alphabet(alphabet);
 659 }
 660
 661 const Alphabet& Sequence::get_alphabet(alphabet_ref alpha) const
 662 {
 663   switch (alpha) {
 664     case reduced_dna_alphabet:
 665       return Alphabet::reduced_dna_alphabet();
 666     case reduced_rna_alphabet:
 667       return Alphabet::reduced_rna_alphabet();
 668     case reduced_nucleic_alphabet:
 669       return Alphabet::reduced_nucleic_alphabet();
 670     case nucleic_alphabet:
 671       return Alphabet::nucleic_alphabet();
 672     case protein_alphabet:
 673       return Alphabet::protein_alphabet();
 674     default:
 675       throw std::runtime_error("unrecognized alphabet type");
 676       break;
 677   }
 678 }
 679
 680 void Sequence::set_sequence(const std::string& s, alphabet_ref a)
 681 {
 682   set_filtered_sequence(s, a);
 683 }
 684
 685 std::string Sequence::get_sequence() const
 686 {
 687   return seq->sequence();
 688 }
 689
 690 Sequence::const_reference Sequence::operator[](Sequence::size_type i) const
 691 {
 692   return at(i);
 693 }
 694
 695 void
 696 Sequence::clear()
 697 {
 698   seq.reset();
 699   strand = UnknownStrand;
 700   header.clear();
 701   species.clear();
 702   annots.clear();
 703   motif_list.clear();
 704 }
 705
 706 void
 707 Sequence::save(fs::fstream &save_file)
 708 {
 709   //fstream save_file;
 710   std::list<annot>::iterator annots_i;
 711
 712   // not sure why, or if i'm doing something wrong, but can't seem to pass
 713   // file pointers down to this method from the mussa control class
 714   // so each call to save a sequence appends to the file started by mussa_class
 715   //save_file.open(save_file_path.c_str(), std::ios::app);
 716
 717   save_file << "<Sequence>" << std::endl;
 718   save_file << *this << std::endl;
 719   save_file << "</Sequence>" << std::endl;
 720
 721   save_file << "<Annotations>" << std::endl;
 722   save_file << species << std::endl;
 723   for (annots_i = annots.begin(); annots_i != annots.end(); ++annots_i)
 724   {
 725     save_file << annots_i->begin << " " << annots_i->end << " " ;
 726     save_file << annots_i->name << " " << annots_i->type << std::endl;
 727   }
 728   save_file << "</Annotations>" << std::endl;
 729   //save_file.close();
 730 }
 731
 732 void
 733 Sequence::load_museq(fs::path load_file_path, int seq_num)
 734 {
 735   fs::fstream load_file;
 736   std::string file_data_line;
 737   int seq_counter;
 738   annot an_annot;
 739   std::string::size_type space_split_i;
 740   std::string annot_value;
 741
 742   annots.clear();
 743   load_file.open(load_file_path, std::ios::in);
 744
 745   seq_counter = 0;
 746   // search for the seq_num-th sequence
 747   while ( (!load_file.eof()) && (seq_counter < seq_num) )
 748   {
 749     getline(load_file,file_data_line);
 750     if (file_data_line == "<Sequence>")
 751       seq_counter++;
 752   }
 753   getline(load_file, file_data_line);
 754   // looks like the sequence is written as a single line
 755   set_filtered_sequence(file_data_line, reduced_dna_alphabet);
 756   getline(load_file, file_data_line);
 757   getline(load_file, file_data_line);
 758   if (file_data_line == "<Annotations>")
 759   {
 760     getline(load_file, file_data_line);
 761     species = file_data_line;
 762     while ( (!load_file.eof())  && (file_data_line != "</Annotations>") )
 763     {
 764       getline(load_file,file_data_line);
 765       if ((file_data_line != "") && (file_data_line != "</Annotations>"))
 766       {
 767         // need to get 4 values...almost same code 4 times...
 768         // get annot start index
 769         space_split_i = file_data_line.find(" ");
 770         annot_value = file_data_line.substr(0,space_split_i);
 771         an_annot.begin = atoi (annot_value.c_str());
 772         file_data_line = file_data_line.substr(space_split_i+1);
 773         // get annot end index
 774         space_split_i = file_data_line.find(" ");
 775         annot_value = file_data_line.substr(0,space_split_i);
 776         an_annot.end = atoi (annot_value.c_str());
 777
 778         if (space_split_i == std::string::npos)  // no entry for type or name
 779         {
 780           std::cout << "seq, annots - no type or name\n";
 781           an_annot.type = "";
 782           an_annot.name = "";
 783         }
 784         else   // else get annot type
 785         {
 786           file_data_line = file_data_line.substr(space_split_i+1);
 787           space_split_i = file_data_line.find(" ");
 788           annot_value = file_data_line.substr(0,space_split_i);
 789           an_annot.type = annot_value;
 790           if (space_split_i == std::string::npos)  // no entry for name
 791           {
 792             std::cout << "seq, annots - no name\n";
 793             an_annot.name = "";
 794           }
 795           else          // get annot name
 796           {
 797             file_data_line = file_data_line.substr(space_split_i+1);
 798             space_split_i = file_data_line.find(" ");
 799             annot_value = file_data_line.substr(0,space_split_i);
 800             an_annot.type = annot_value;
 801           }
 802         }
 803         annots.push_back(an_annot);  // don't forget to actually add the annot
 804       }
 805       //std::cout << "seq, annots: " << an_annot.start << ", " << an_annot.end
 806       //     << "-->" << an_annot.type << "::" << an_annot.name << std::endl;
 807     }
 808   }
 809   load_file.close();
 810 }
 811
 812
 813 void Sequence::add_motif(const Sequence& a_motif)
 814 {
 815   std::vector<int> motif_starts = find_motif(a_motif);
 816
 817   for(std::vector<int>::iterator motif_start_i = motif_starts.begin();
 818       motif_start_i != motif_starts.end();
 819       ++motif_start_i)
 820   {
 821     motif_list.push_back(motif(*motif_start_i, a_motif.get_sequence()));
 822   }
 823 }
 824
 825 void Sequence::clear_motifs()
 826 {
 827   motif_list.clear();
 828 }
 829
 830 const std::list<motif>& Sequence::motifs() const
 831 {
 832   return motif_list;
 833 }
 834
 835 std::vector<int>
 836 Sequence::find_motif(const Sequence& a_motif) const
 837 {
 838   std::vector<int> motif_match_starts;
 839   Sequence norm_motif_rc;
 840
 841   motif_match_starts.clear();
 842   // std::cout << "motif is: " << norm_motif << std::endl;
 843
 844   if (a_motif.size() > 0)
 845   {
 846     //std::cout << "Sequence: none blank motif\n";
 847     motif_scan(a_motif, &motif_match_starts);
 848
 849     norm_motif_rc = a_motif.rev_comp();;
 850     // make sure not to do search again if it is a palindrome
 851     if (norm_motif_rc != a_motif) {
 852       motif_scan(norm_motif_rc, &motif_match_starts);
 853     }
 854   }
 855   return motif_match_starts;
 856 }
 857
 858 void
 859 Sequence::motif_scan(const Sequence& a_motif, std::vector<int> * motif_match_starts) const
 860 {
 861   // if there's no sequence we can't scan for it?
 862   // should this throw an exception?
 863   if (!seq) return;
 864
 865   std::string::size_type seq_i = 0;
 866   Sequence::size_type motif_i = 0;
 867   Sequence::size_type motif_len = a_motif.length();
 868   Sequence::value_type motif_char;
 869   Sequence::value_type seq_char;
 870
 871   while (seq_i < size())
 872   {
 873     // this is pretty much a straight translation of Nora's python code
 874     // to match iupac letter codes
 875     motif_char = toupper(a_motif[motif_i]);
 876     seq_char = toupper(seq->at(seq_i));
 877     if (motif_char =='N')
 878       motif_i++;
 879     else if (motif_char == seq_char)
 880       motif_i++;
 881     else if ((motif_char =='M') && ((seq_char=='A') || (seq_char=='C')))
 882       motif_i++;
 883     else if ((motif_char =='R') && ((seq_char=='A') || (seq_char=='G')))
 884       motif_i++;
 885     else if ((motif_char =='W') && ((seq_char=='A') || (seq_char=='T')))
 886       motif_i++;
 887     else if ((motif_char =='S') && ((seq_char=='C') || (seq_char=='G')))
 888       motif_i++;
 889     else if ((motif_char =='Y') && ((seq_char=='C') || (seq_char=='T')))
 890       motif_i++;
 891     else if ((motif_char =='K') && ((seq_char=='G') || (seq_char=='T')))
 892       motif_i++;
 893     else if ((motif_char =='V') &&
 894              ((seq_char=='A') || (seq_char=='C') || (seq_char=='G')))
 895       motif_i++;
 896     else if ((motif_char =='H') &&
 897              ((seq_char=='A') || (seq_char=='C') || (seq_char=='T')))
 898       motif_i++;
 899     else if ((motif_char =='D') &&
 900              ((seq_char=='A') || (seq_char=='G') || (seq_char=='T')))
 901       motif_i++;
 902     else if ((motif_char =='B') &&
 903              ((seq_char=='C') || (seq_char=='G') || (seq_char=='T')))
 904       motif_i++;
 905     else
 906     {
 907       // if a motif doesn't match, erase our current trial and try again
 908       seq_i -= motif_i;
 909       motif_i = 0;
 910     }
 911
 912     // end Nora stuff, now we see if a match is found this pass
 913     if (motif_i == motif_len)
 914     {
 915       annot new_motif;
 916       motif_match_starts->push_back(seq_i - motif_len + 1);
 917       motif_i = 0;
 918     }
 919
 920     seq_i++;
 921   }
 922   //std::cout << std::endl;
 923 }
 924
 925 void Sequence::add_string_annotation(std::string a_seq,
 926                                      std::string name)
 927 {
 928   std::vector<int> seq_starts = find_motif(a_seq);
 929
 930   //std::cout << "searching for " << a_seq << " found " << seq_starts.size() << std::endl;
 931
 932   for(std::vector<int>::iterator seq_start_i = seq_starts.begin();
 933       seq_start_i != seq_starts.end();
 934       ++seq_start_i)
 935   {
 936     annots.push_back(annot(*seq_start_i,
 937                            *seq_start_i+a_seq.size(),
 938                            "",
 939                            name));
 940   }
 941 }
 942
 943 void Sequence::find_sequences(std::list<Sequence>::iterator start,
 944                               std::list<Sequence>::iterator end)
 945 {
 946   while (start != end) {
 947     add_string_annotation(start->get_sequence(), start->get_fasta_header());
 948     ++start;
 949   }
 950 }
 951
 952
 953 std::ostream& operator<<(std::ostream& out, const Sequence& s)
 954 {
 955   if (s.seq) {
 956     for(Sequence::const_iterator s_i = s.begin(); s_i != s.end(); ++s_i) {
 957       out << *s_i;
 958     }
 959   }
 960   return out;
 961 }
 962
 963 bool operator<(const Sequence& x, const Sequence& y)
 964 {
 965   Sequence::const_iterator x_i = x.begin();
 966   Sequence::const_iterator y_i = y.begin();
 967   // for sequences there's some computation associated with computing .end
 968   // so lets cache it.
 969   Sequence::const_iterator xend = x.end();
 970   Sequence::const_iterator yend = y.end();
 971   while(1) {
 972     if( x_i == xend and y_i == yend ) {
 973       return false;
 974     } else if ( x_i == xend ) {
 975       return true;
 976     } else if ( y_i == yend ) {
 977       return false;
 978     } else if ( (*x_i) < (*y_i)) {
 979       return true;
 980     } else if ( (*x_i) > (*y_i) ) {
 981       return false;
 982     } else {
 983       ++x_i;
 984       ++y_i;
 985     }
 986   }
 987 }
 988
 989 template <typename Iter1, typename Iter2>
 990 static
 991 bool sequence_insensitive_equality(Iter1 abegin, Iter1 aend, Iter2 bbegin, Iter2 bend)
 992 {
 993   Iter1 aseq_i = abegin;
 994   Iter2 bseq_i = bbegin;
 995   if (aend-abegin == bend-bbegin) {
 996     // since the length of the two sequences is equal, we only need to
 997     // test one.
 998     for(; aseq_i != aend; ++aseq_i, ++bseq_i) {
 999       if (toupper(*aseq_i) != toupper(*bseq_i)) {
1000         return false;
1001       }
1002     }
1003     return true;
1004   } else {
1005     return false;
1006   }
1007 }
1008
1009 bool operator==(const Sequence& x, const Sequence& y)
1010 {
1011   if (x.seq and y.seq) {
1012     // both x and y are defined
1013     if (SeqSpan::isFamily(x.seq, y.seq)) {
1014       // both are part of the same SeqSpan tree
1015       return *(x.seq) == *(y.seq);
1016     } else {
1017       // we'll have to do a real comparison
1018       return sequence_insensitive_equality<SeqSpan::const_iterator, SeqSpan::const_iterator>(
1019                x.begin(), x.end(),
1020                y.begin(), y.end()
1021              );
1022     }
1023   } else {
1024     // true if they're both empty (with either a null SeqSpanRef or
1025     // a zero length string
1026     return (x.size() == y.size());
1027   }
1028 }
1029
1030 bool operator!=(const Sequence& x, const Sequence& y)
1031 {
1032   return not operator==(x, y);
1033 }
1034