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