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