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