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   // so i should probably be passing the parse function some iterators
 328   // but the annotations files are (currently) small, so i think i can
 329   // get away with loading the whole file into memory
 330   std::string data;
 331   char c;
 332   while(data_stream.good()) {
 333     data_stream.get(c);
 334     data.push_back(c);
 335   }
 336   data_stream.close();
 337
 338   try {
 339     parse_annot(data, start_index, end_index);
 340   } catch(annotation_load_error e) {
 341     std::ostringstream msg;
 342     msg << file_path.native_file_string()
 343         << " "
 344         << e.what();
 345     throw annotation_load_error(msg.str());
 346   }
 347 }
 348
 349 /* If this works, yikes, this is some brain hurting code.
 350  *
 351  * what's going on is that when pb_annot is instantiated it stores references
 352  * to begin, end, name, type, declared in the parse function, then
 353  * when operator() is called it grabs values from those references
 354  * and uses that to instantiate an annot object and append that to our
 355  * annotation list.
 356  *
 357  * This weirdness is because the spirit library requires that actions
 358  * conform to a specific prototype operator()(IteratorT, IteratorT)
 359  * which doesn't provide any useful opportunity for me to actually
 360  * grab the results of our parsing.
 361  *
 362  * so I instantiate this structure in order to have a place to grab
 363  * my data from.
 364  */
 365
 366 struct push_back_annot {
 367   Sequence* parent;
 368   SeqSpanRefListRef children;
 369   int& begin;
 370   int& end;
 371   std::string& name;
 372   std::string& type;
 373   int &parsed;
 374
 375   push_back_annot(Sequence* parent_seq,
 376                   SeqSpanRefListRef children_list,
 377                   int& begin_,
 378                   int& end_,
 379                   std::string& name_,
 380                   std::string& type_,
 381                   int &parsed_)
 382   : parent(parent_seq),
 383     children(children_list),
 384     begin(begin_),
 385     end(end_),
 386     name(name_),
 387     type(type_),
 388     parsed(parsed_)
 389   {
 390   }
 391
 392   void operator()(std::string::const_iterator,
 393                   std::string::const_iterator) const
 394   {
 395     children->push_back(parent->make_annotation(name, type, begin, end));
 396     ++parsed;
 397   };
 398 };
 399
 400 struct push_back_seq {
 401   std::list<Sequence>& seq_list;
 402   std::string& name;
 403   std::string& seq;
 404   int &parsed;
 405
 406   push_back_seq(std::list<Sequence>& seq_list_,
 407                 std::string& name_,
 408                 std::string& seq_,
 409                 int &parsed_)
 410   : seq_list(seq_list_),
 411     name(name_),
 412     seq(seq_),
 413     parsed(parsed_)
 414   {
 415   }
 416
 417   void operator()(std::string::const_iterator,
 418                   std::string::const_iterator) const
 419   {
 420     // filter out newlines from our sequence
 421     std::string new_seq;
 422     for(std::string::const_iterator seq_i = seq.begin();
 423         seq_i != seq.end();
 424         ++seq_i)
 425     {
 426       if (*seq_i != '\015' && *seq_i != '\012') new_seq += *seq_i;
 427     }
 428     //std::cout << "adding seq: " << name << " " << new_seq << std::endl;
 429
 430     Sequence s(new_seq);
 431     s.set_fasta_header(name);
 432     seq_list.push_back(s);
 433     ++parsed;
 434   };
 435 };
 436
 437 void
 438 Sequence::parse_annot(std::string data, int start_index, int end_index)
 439 {
 440   int start=0;
 441   int end=0;
 442   std::string name;
 443   std::string type;
 444   std::string seqstr;
 445   SeqSpanRefListRef parsed_annots(new SeqSpanRefList);
 446   std::list<Sequence> query_seqs;
 447   int parsed=0;
 448
 449   bool ok = spirit::parse(data.begin(), data.end(),
 450               (
 451                //begin grammar
 452                  !(
 453                     (
 454                       spirit::alpha_p >>
 455                       +(spirit::graph_p)
 456                     )[spirit::assign_a(species)] >>
 457                     +(spirit::space_p)
 458                   ) >>
 459                   *(
 460                      ( // ignore html tags
 461                        *(spirit::space_p) >>
 462                        spirit::ch_p('<') >>
 463                        +(~spirit::ch_p('>')) >>
 464                        spirit::ch_p('>') >>
 465                        *(spirit::space_p)
 466                      )
 467                    |
 468                     ( // parse an absolute location name
 469                      (spirit::uint_p[spirit::assign_a(start)] >>
 470                       +spirit::space_p >>
 471                       spirit::uint_p[spirit::assign_a(end)] >>
 472                       +spirit::space_p >>
 473                       (
 474                          spirit::alpha_p >>
 475                          *spirit::graph_p
 476                       )[spirit::assign_a(name)] >>
 477                       // optional type
 478                       !(
 479                           +spirit::space_p >>
 480                           (
 481                             spirit::alpha_p >>
 482                             *spirit::graph_p
 483                           )[spirit::assign_a(type)]
 484                       )
 485                       // to understand how this group gets set
 486                       // read the comment above struct push_back_annot
 487                      )[push_back_annot(this, parsed_annots, start, end, name, type, parsed)]
 488                    |
 489                     ((spirit::ch_p('>')|spirit::str_p("&gt;")) >>
 490                        (*(spirit::print_p))[spirit::assign_a(name)] >>
 491                        spirit::eol_p >>
 492                        (+(spirit::chset<>(Alphabet::nucleic_cstr)))[spirit::assign_a(seqstr)]
 493                      )[push_back_seq(query_seqs, name, seqstr, parsed)]
 494                     ) >>
 495                     *spirit::space_p
 496                    )
 497               //end grammar
 498               )).full;
 499   if (not ok) {
 500     std::stringstream msg;
 501     msg << "Error parsing annotation #" << parsed;
 502     throw annotation_load_error(msg.str());
 503   }
 504   // If everything loaded correctly add the sequences to our annotation list
 505   // add newly parsed annotations to our sequence
 506   std::copy(parsed_annots->begin(), parsed_annots->end(), std::back_inserter(*annotation_list));
 507   // go search for query sequences
 508   find_sequences(query_seqs.begin(), query_seqs.end());
 509 }
 510
 511 void Sequence::add_annotation(const SeqSpanRef a)
 512 {
 513   annotation_list->push_back(a);
 514 }
 515
 516 void Sequence::add_annotation(std::string name, std::string type, size_type start, size_type stop)
 517 {
 518   add_annotation(make_annotation(name, type, start, stop));
 519 }
 520
 521 SeqSpanRef
 522 Sequence::make_annotation(std::string name, std::string type, size_type start, size_type stop) const
 523 {
 524   // we want things to be in the positive direction
 525   if (stop < start) {
 526     size_type tmp = start;
 527     start = stop;
 528     stop = tmp;
 529   }
 530   size_type count = stop - start;
 531   SeqSpanRef new_annot(seq->subseq(start, count,  SeqSpan::UnknownStrand));
 532   AnnotationsRef metadata(new Annotations(name));
 533   metadata->set("type", type);
 534   new_annot->setAnnotations(metadata);
 535   return new_annot;
 536 }
 537
 538 const SeqSpanRefList& Sequence::annotations() const
 539 {
 540   return *annotation_list;
 541 }
 542
 543 void Sequence::copy_children(Sequence &new_seq, size_type start, size_type count) const
 544 {
 545   new_seq.motif_list = motif_list;
 546   new_seq.annotation_list.reset(new SeqSpanRefList);
 547
 548   for(SeqSpanRefList::const_iterator annot_i = annotation_list->begin();
 549       annot_i != annotation_list->end();
 550       ++annot_i)
 551   {
 552     size_type annot_begin= (*annot_i)->start();
 553     size_type annot_end = (*annot_i)->stop();
 554
 555     if (annot_begin < start+count) {
 556       if (annot_begin >= start) {
 557         annot_begin -= start;
 558       } else {
 559         annot_begin = 0;
 560       }
 561
 562       if (annot_end < start+count) {
 563         annot_end -= start;
 564       } else {
 565         annot_end = count;
 566       }
 567       SeqSpanRef new_annot(new_seq.seq->subseq(annot_begin, annot_end));
 568       new_annot->setAnnotations((*annot_i)->annotations());
 569       new_seq.annotation_list->push_back(new_annot);
 570     }
 571   }
 572 }
 573
 574 Sequence
 575 Sequence::subseq(size_type start, size_type count, SeqSpan::strand_type strand) const
 576 {
 577   // FIXME: should i really allow a subsequence of an empty sequence?
 578   if (!seq) {
 579     Sequence new_seq;
 580     return new_seq;
 581   }
 582
 583   Sequence new_seq(*this);
 584   new_seq.seq = seq->subseq(start, count, strand);
 585   if (seq->annotations()) {
 586     AnnotationsRef a(new Annotations(*(seq->annotations())));
 587     new_seq.seq->setAnnotations(a);
 588   }
 589   copy_children(new_seq, start, count);
 590
 591   return new_seq;
 592 }
 593
 594
 595 // FIXME: This needs to be moved into SeqSpan
 596 Sequence Sequence::rev_comp() const
 597 {
 598   // a reverse complement is the whole opposite strand
 599   return subseq(0, npos, SeqSpan::OppositeStrand);
 600 }
 601
 602 const Alphabet& Sequence::get_alphabet() const
 603 {
 604   return (seq) ? seq->get_alphabet() : Alphabet::empty_alphabet();
 605 }
 606
 607 void Sequence::set_fasta_header(std::string header_)
 608 {
 609   header = header_;
 610 }
 611
 612 void Sequence::set_species(const std::string& name)
 613 {
 614   species = name;
 615 }
 616
 617 std::string Sequence::get_species() const
 618 {
 619   return species;
 620 }
 621
 622
 623 std::string
 624 Sequence::get_fasta_header() const
 625 {
 626   return header;
 627 }
 628
 629 std::string
 630 Sequence::get_name() const
 631 {
 632   if (header.size() > 0)
 633     return header;
 634   else if (species.size() > 0)
 635     return species;
 636   else
 637     return "";
 638 }
 639
 640 void Sequence::set_sequence(const std::string& s, AlphabetRef a)
 641 {
 642   set_filtered_sequence(s, a, 0, s.size(), SeqSpan::PlusStrand);
 643 }
 644
 645 std::string Sequence::get_sequence() const
 646 {
 647   return seq->sequence();
 648 }
 649
 650 Sequence::const_reference Sequence::operator[](Sequence::size_type i) const
 651 {
 652   return at(i);
 653 }
 654
 655 void
 656 Sequence::clear()
 657 {
 658   seq.reset();
 659   header.clear();
 660   species.clear();
 661   annotation_list.reset(new SeqSpanRefList);
 662   motif_list.reset(new MotifList);
 663 }
 664
 665 void
 666 Sequence::save(fs::fstream &save_file)
 667 {
 668   std::string type("type");
 669   std::string empty_str("");
 670   //fstream save_file;
 671   SeqSpanRefList::iterator annots_i;
 672   AnnotationsRef metadata;
 673
 674   // not sure why, or if i'm doing something wrong, but can't seem to pass
 675   // file pointers down to this method from the mussa control class
 676   // so each call to save a sequence appends to the file started by mussa_class
 677   //save_file.open(save_file_path.c_str(), std::ios::app);
 678
 679   save_file << "<Sequence>" << std::endl;
 680   save_file << *this << std::endl;
 681   save_file << "</Sequence>" << std::endl;
 682
 683   save_file << "<Annotations>" << std::endl;
 684   save_file << species << std::endl;
 685   for (annots_i = annotation_list->begin();
 686        annots_i != annotation_list->end();
 687        ++annots_i)
 688   {
 689     metadata = (*annots_i)->annotations();
 690     save_file << (*annots_i)->parentStart() << " " << (*annots_i)->parentStop() << " " ;
 691     save_file << metadata->name() << " "
 692               << metadata->getdefault(type, empty_str) << std::endl;
 693   }
 694   save_file << "</Annotations>" << std::endl;
 695   //save_file.close();
 696 }
 697
 698 //void
 699 //Sequence::load_museq(fs::path load_file_path, int seq_num)
 700 //{
 701 //  fs::fstream load_file;
 702 //  std::string file_data_line;
 703 //  int seq_counter;
 704 //  //annot an_annot;
 705 //  int annot_begin;
 706 //  int annot_end;
 707 //  std::string annot_name;
 708 //  std::string annot_type;
 709 //
 710 //  std::string::size_type space_split_i;
 711 //  std::string annot_value;
 712 //
 713 //  annotation_list.reset(new SeqSpanRefList);
 714 //
 715 //  load_file.open(load_file_path, std::ios::in);
 716 //
 717 //  seq_counter = 0;
 718 //  // search for the seq_num-th sequence
 719 //  while ( (!load_file.eof()) && (seq_counter < seq_num) )
 720 //  {
 721 //    getline(load_file,file_data_line);
 722 //    if (file_data_line == "<Sequence>")
 723 //      seq_counter++;
 724 //  }
 725 //  getline(load_file, file_data_line);
 726 //  // looks like the sequence is written as a single line
 727 //  set_filtered_sequence(file_data_line, reduced_dna_alphabet, 0, file_data_line.size(), SeqSpan::PlusStrand);
 728 //  getline(load_file, file_data_line);
 729 //  getline(load_file, file_data_line);
 730 //  if (file_data_line == "<Annotations>")
 731 //  {
 732 //    getline(load_file, file_data_line);
 733 //    species = file_data_line;
 734 //    while ( (!load_file.eof())  && (file_data_line != "</Annotations>") )
 735 //    {
 736 //      getline(load_file,file_data_line);
 737 //      if ((file_data_line != "") && (file_data_line != "</Annotations>"))
 738 //      {
 739 //        // need to get 4 values...almost same code 4 times...
 740 //        // get annot start index
 741 //        space_split_i = file_data_line.find(" ");
 742 //        annot_value = file_data_line.substr(0,space_split_i);
 743 //        annot_begin = atoi (annot_value.c_str());
 744 //        file_data_line = file_data_line.substr(space_split_i+1);
 745 //        // get annot end index
 746 //        space_split_i = file_data_line.find(" ");
 747 //        annot_value = file_data_line.substr(0,space_split_i);
 748 //        annot_end = atoi (annot_value.c_str());
 749 //
 750 //        if (space_split_i == std::string::npos)  // no entry for type or name
 751 //        {
 752 //          std::cout << "seq, annots - no type or name\n";
 753 //          annot_name = "";
 754 //          annot_type = "";
 755 //        }
 756 //        else   // else get annot type
 757 //        {
 758 //          file_data_line = file_data_line.substr(space_split_i+1);
 759 //          space_split_i = file_data_line.find(" ");
 760 //          annot_value = file_data_line.substr(0,space_split_i);
 761 //          //an_annot.type = annot_value;
 762 //          annot_type = annot_value;
 763 //          if (space_split_i == std::string::npos)  // no entry for name
 764 //          {
 765 //            std::cout << "seq, annots - no name\n";
 766 //            annot_name = "";
 767 //          }
 768 //          else          // get annot name
 769 //          {
 770 //            file_data_line = file_data_line.substr(space_split_i+1);
 771 //            space_split_i = file_data_line.find(" ");
 772 //            annot_value = file_data_line.substr(0,space_split_i);
 773 //            // this seems like its wrong?
 774 //            annot_type = annot_value;
 775 //          }
 776 //        }
 777 //        add_annotation(annot_name, annot_type, annot_begin, annot_end);
 778 //      }
 779 //      //std::cout << "seq, annots: " << an_annot.start << ", " << an_annot.end
 780 //      //     << "-->" << an_annot.type << "::" << an_annot.name << std::endl;
 781 //    }
 782 //  }
 783 //  load_file.close();
 784 //}
 785
 786 SequenceRef Sequence::load_museq(boost::filesystem::fstream& load_file)
 787 {
 788   boost::shared_ptr<Sequence> seq(new Sequence);
 789   std::string file_data_line;
 790   int seq_counter;
 791   //annot an_annot;
 792   int annot_begin;
 793   int annot_end;
 794   std::string annot_name;
 795   std::string annot_type;
 796
 797   std::string::size_type space_split_i;
 798   std::string annot_value;
 799
 800   //seq->annotation_list.reset(new SeqSpanRefList);
 801
 802   seq_counter = 0;
 803   // search for the next sequence
 804   int seq_num = 1;
 805   while ( (!load_file.eof()) && (seq_counter < seq_num) )
 806   {
 807     getline(load_file,file_data_line);
 808     if (file_data_line == "<Sequence>")
 809       seq_counter++;
 810   }
 811
 812   // Could not find next sequence
 813   if (load_file.eof())
 814   {
 815     seq.reset();
 816     return seq;
 817   }
 818
 819   getline(load_file, file_data_line);
 820   // looks like the sequence is written as a single line
 821   seq->set_filtered_sequence(file_data_line, reduced_dna_alphabet, 0, file_data_line.size(), SeqSpan::PlusStrand);
 822   getline(load_file, file_data_line);
 823   getline(load_file, file_data_line);
 824   if (file_data_line == "<Annotations>")
 825   {
 826     getline(load_file, file_data_line);
 827     seq->set_species(file_data_line);
 828     while ( (!load_file.eof())  && (file_data_line != "</Annotations>") )
 829     {
 830       getline(load_file,file_data_line);
 831       if ((file_data_line != "") && (file_data_line != "</Annotations>"))
 832       {
 833         // need to get 4 values...almost same code 4 times...
 834         // get annot start index
 835         space_split_i = file_data_line.find(" ");
 836         annot_value = file_data_line.substr(0,space_split_i);
 837         annot_begin = atoi (annot_value.c_str());
 838         file_data_line = file_data_line.substr(space_split_i+1);
 839         // get annot end index
 840         space_split_i = file_data_line.find(" ");
 841         annot_value = file_data_line.substr(0,space_split_i);
 842         annot_end = atoi (annot_value.c_str());
 843
 844         if (space_split_i == std::string::npos)  // no entry for type or name
 845         {
 846           std::cout << "seq, annots - no type or name\n";
 847           annot_name = "";
 848           annot_type = "";
 849         }
 850         else   // else get annot type
 851         {
 852           file_data_line = file_data_line.substr(space_split_i+1);
 853           space_split_i = file_data_line.find(" ");
 854           annot_value = file_data_line.substr(0,space_split_i);
 855           //an_annot.type = annot_value;
 856           annot_type = annot_value;
 857           if (space_split_i == std::string::npos)  // no entry for name
 858           {
 859             std::cout << "seq, annots - no name\n";
 860             annot_name = "";
 861           }
 862           else          // get annot name
 863           {
 864             file_data_line = file_data_line.substr(space_split_i+1);
 865             space_split_i = file_data_line.find(" ");
 866             annot_value = file_data_line.substr(0,space_split_i);
 867             // this seems like its wrong?
 868             annot_type = annot_value;
 869           }
 870         }
 871         seq->add_annotation(annot_name, annot_type, annot_begin, annot_end);
 872       }
 873       //std::cout << "seq, annots: " << an_annot.start << ", " << an_annot.end
 874       //     << "-->" << an_annot.type << "::" << an_annot.name << std::endl;
 875     }
 876   }
 877   //load_file.close();
 878   return seq;
 879 }
 880
 881
 882 void Sequence::add_motif(const Sequence& a_motif)
 883 {
 884   std::vector<int> motif_starts = find_motif(a_motif);
 885
 886   for(std::vector<int>::iterator motif_start_i = motif_starts.begin();
 887       motif_start_i != motif_starts.end();
 888       ++motif_start_i)
 889   {
 890     motif_list->push_back(motif(*motif_start_i, a_motif.get_sequence()));
 891   }
 892 }
 893
 894 void Sequence::clear_motifs()
 895 {
 896   if (motif_list)
 897     motif_list->clear();
 898   else
 899     motif_list.reset(new MotifList);
 900 }
 901
 902 const Sequence::MotifList& Sequence::motifs() const
 903 {
 904   return *motif_list;
 905 }
 906
 907 std::vector<int>
 908 Sequence::find_motif(const Sequence& a_motif) const
 909 {
 910   std::vector<int> motif_match_starts;
 911   Sequence norm_motif_rc;
 912
 913   motif_match_starts.clear();
 914   // std::cout << "motif is: " << norm_motif << std::endl;
 915
 916   if (a_motif.size() > 0)
 917   {
 918     //std::cout << "Sequence: none blank motif\n";
 919     motif_scan(a_motif, &motif_match_starts);
 920
 921     norm_motif_rc = a_motif.rev_comp();;
 922     // make sure not to do search again if it is a palindrome
 923     if (norm_motif_rc != a_motif) {
 924       motif_scan(norm_motif_rc, &motif_match_starts);
 925     }
 926   }
 927   return motif_match_starts;
 928 }
 929
 930 void
 931 Sequence::motif_scan(const Sequence& a_motif, std::vector<int> * motif_match_starts) const
 932 {
 933   // if there's no sequence we can't scan for it?
 934   // should this throw an exception?
 935   if (!seq) return;
 936
 937   std::string::size_type seq_i = 0;
 938   Sequence::size_type motif_i = 0;
 939   Sequence::size_type motif_len = a_motif.length();
 940   Sequence::value_type motif_char;
 941   Sequence::value_type seq_char;
 942
 943   while (seq_i < size())
 944   {
 945     // this is pretty much a straight translation of Nora's python code
 946     // to match iupac letter codes
 947     motif_char = toupper(a_motif[motif_i]);
 948     seq_char = toupper(seq->at(seq_i));
 949     if (motif_char =='N')
 950       motif_i++;
 951     else if (motif_char == seq_char)
 952       motif_i++;
 953     else if ((motif_char =='M') && ((seq_char=='A') || (seq_char=='C')))
 954       motif_i++;
 955     else if ((motif_char =='R') && ((seq_char=='A') || (seq_char=='G')))
 956       motif_i++;
 957     else if ((motif_char =='W') && ((seq_char=='A') || (seq_char=='T')))
 958       motif_i++;
 959     else if ((motif_char =='S') && ((seq_char=='C') || (seq_char=='G')))
 960       motif_i++;
 961     else if ((motif_char =='Y') && ((seq_char=='C') || (seq_char=='T')))
 962       motif_i++;
 963     else if ((motif_char =='K') && ((seq_char=='G') || (seq_char=='T')))
 964       motif_i++;
 965     else if ((motif_char =='V') &&
 966              ((seq_char=='A') || (seq_char=='C') || (seq_char=='G')))
 967       motif_i++;
 968     else if ((motif_char =='H') &&
 969              ((seq_char=='A') || (seq_char=='C') || (seq_char=='T')))
 970       motif_i++;
 971     else if ((motif_char =='D') &&
 972              ((seq_char=='A') || (seq_char=='G') || (seq_char=='T')))
 973       motif_i++;
 974     else if ((motif_char =='B') &&
 975              ((seq_char=='C') || (seq_char=='G') || (seq_char=='T')))
 976       motif_i++;
 977     else
 978     {
 979       // if a motif doesn't match, erase our current trial and try again
 980       seq_i -= motif_i;
 981       motif_i = 0;
 982     }
 983
 984     // end Nora stuff, now we see if a match is found this pass
 985     if (motif_i == motif_len)
 986     {
 987       motif_match_starts->push_back(seq_i - motif_len + 1);
 988       motif_i = 0;
 989     }
 990
 991     seq_i++;
 992   }
 993   //std::cout << std::endl;
 994 }
 995
 996 void Sequence::add_string_annotation(std::string a_seq,
 997                                      std::string name)
 998 {
 999   std::vector<int> seq_starts = find_motif(a_seq);
1000
1001   for(std::vector<int>::iterator seq_start_i = seq_starts.begin();
1002       seq_start_i != seq_starts.end();
1003       ++seq_start_i)
1004   {
1005     add_annotation(name, "", *seq_start_i, *seq_start_i+a_seq.size());
1006   }
1007 }
1008
1009 void Sequence::find_sequences(std::list<Sequence>::iterator start,
1010                               std::list<Sequence>::iterator end)
1011 {
1012   while (start != end) {
1013     add_string_annotation(start->get_sequence(), start->get_fasta_header());
1014     ++start;
1015   }
1016 }
1017
1018
1019 std::ostream& operator<<(std::ostream& out, const Sequence& s)
1020 {
1021   if (s.seq) {
1022     for(Sequence::const_iterator s_i = s.begin(); s_i != s.end(); ++s_i) {
1023       out << *s_i;
1024     }
1025   }
1026   return out;
1027 }
1028
1029 bool operator<(const Sequence& x, const Sequence& y)
1030 {
1031   Sequence::const_iterator x_i = x.begin();
1032   Sequence::const_iterator y_i = y.begin();
1033   // for sequences there's some computation associated with computing .end
1034   // so lets cache it.
1035   Sequence::const_iterator xend = x.end();
1036   Sequence::const_iterator yend = y.end();
1037   while(1) {
1038     if( x_i == xend and y_i == yend ) {
1039       return false;
1040     } else if ( x_i == xend ) {
1041       return true;
1042     } else if ( y_i == yend ) {
1043       return false;
1044     } else if ( (*x_i) < (*y_i)) {
1045       return true;
1046     } else if ( (*x_i) > (*y_i) ) {
1047       return false;
1048     } else {
1049       ++x_i;
1050       ++y_i;
1051     }
1052   }
1053 }
1054
1055 template <typename Iter1, typename Iter2>
1056 static
1057 bool sequence_insensitive_equality(Iter1 abegin, Iter1 aend, Iter2 bbegin, Iter2 bend)
1058 {
1059   Iter1 aseq_i = abegin;
1060   Iter2 bseq_i = bbegin;
1061   if (aend-abegin == bend-bbegin) {
1062     // since the length of the two sequences is equal, we only need to
1063     // test one.
1064     for(; aseq_i != aend; ++aseq_i, ++bseq_i) {
1065       if (toupper(*aseq_i) != toupper(*bseq_i)) {
1066         return false;
1067       }
1068     }
1069     return true;
1070   } else {
1071     return false;
1072   }
1073 }
1074
1075 bool operator==(const Sequence& x, const Sequence& y)
1076 {
1077   if (x.seq and y.seq) {
1078     // both x and y are defined
1079     if (SeqSpan::isFamily(x.seq, y.seq)) {
1080       // both are part of the same SeqSpan tree
1081       return *(x.seq) == *(y.seq);
1082     } else {
1083       // we'll have to do a real comparison
1084       return sequence_insensitive_equality<SeqSpan::const_iterator, SeqSpan::const_iterator>(
1085                x.begin(), x.end(),
1086                y.begin(), y.end()
1087              );
1088     }
1089   } else {
1090     // true if they're both empty (with either a null SeqSpanRef or
1091     // a zero length string
1092     return (x.size() == y.size());
1093   }
1094 }
1095
1096 bool operator!=(const Sequence& x, const Sequence& y)
1097 {
1098   return not operator==(x, y);
1099 }
1100