alg/sequence.cpp

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