X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=mussa.git;a=blobdiff_plain;f=alg%2Fsequence.cpp;h=55ff5bf91bb575f9322102fafc78d4259ac1bd62;hp=6cb4e2cf001b443be1bb5f10d0ceda96be9eb7aa;hb=67888dae3b16b9d69aa846e393f11e7ff3633f16;hpb=5c3dc8c42679629c19ece07c0e63a53b1c69663f diff --git a/alg/sequence.cpp b/alg/sequence.cpp index 6cb4e2c..55ff5bf 100644 --- a/alg/sequence.cpp +++ b/alg/sequence.cpp @@ -40,27 +40,7 @@ namespace spirit = boost::spirit; #include #include -annot::annot() - : begin(0), - end(0), - type(""), - name("") -{ -} - -annot::annot(int begin, int end, std::string type, std::string name) - : begin(begin), - end(end), - type(type), - name(name) -{ -} - -annot::~annot() -{ -} - -bool operator==(const annot& left, const annot& right) +bool operator==(const motif& left, const motif& right) { return ((left.begin== right.begin) and (left.end == right.end) and @@ -68,8 +48,20 @@ bool operator==(const annot& left, const annot& right) (left.name == right.name)); } -motif::motif(int begin, std::string motif) - : annot(begin, begin+motif.size(), "motif", motif), +motif::motif() + : begin(0), + end(0), + type("motif"), + name(""), + sequence("") +{ +} + +motif::motif(int begin_, std::string motif) + : begin(begin_), + end(begin_+motif.size()), + type("motif"), + name(motif), sequence(motif) { } @@ -78,9 +70,9 @@ motif::~motif() { } - Sequence::Sequence(AlphabetRef alphabet) - : seq(new SeqSpan("", alphabet, SeqSpan::PlusStrand)), + : seq(new SeqSpan("", alphabet, SeqSpan::PlusStrand)), + annotation_list(new SeqSpanRefList), motif_list(new MotifList) { } @@ -92,6 +84,7 @@ Sequence::~Sequence() Sequence::Sequence(const char *seq, AlphabetRef alphabet_, SeqSpan::strand_type strand_) : header(""), species(""), + annotation_list(new SeqSpanRefList), motif_list(new MotifList) { set_filtered_sequence(seq, alphabet_, 0, npos, strand_); @@ -102,6 +95,7 @@ Sequence::Sequence(const std::string& seq, SeqSpan::strand_type strand_) : header(""), species(""), + annotation_list(new SeqSpanRefList), motif_list(new MotifList) { set_filtered_sequence(seq, alphabet_, 0, seq.size(), strand_); @@ -111,7 +105,7 @@ Sequence::Sequence(const Sequence& o) : seq(o.seq), header(o.header), species(o.species), - annots(o.annots), + annotation_list(o.annotation_list), motif_list(o.motif_list) { } @@ -120,7 +114,7 @@ Sequence::Sequence(const Sequence* o) : seq(o->seq), header(o->header), species(o->species), - annots(o->annots), + annotation_list(o->annotation_list), motif_list(o->motif_list) { } @@ -129,7 +123,7 @@ Sequence::Sequence(const SequenceRef o) : seq(new SeqSpan(o->seq)), header(o->header), species(o->species), - annots(o->annots), + annotation_list(o->annotation_list), motif_list(o->motif_list) { } @@ -138,6 +132,7 @@ Sequence::Sequence(const SeqSpanRef& seq_ref) : seq(seq_ref), header(""), species(""), + annotation_list(new SeqSpanRefList), motif_list(new MotifList) { } @@ -148,7 +143,7 @@ Sequence &Sequence::operator=(const Sequence& s) seq = s.seq; header = s.header; species = s.species; - annots = s.annots; + annotation_list = s.annotation_list; motif_list = s.motif_list; } return *this; @@ -371,20 +366,23 @@ Sequence::load_annot(fs::path file_path, int start_index, int end_index) */ struct push_back_annot { - std::list& annot_list; + Sequence* parent; + SeqSpanRefListRef children; int& begin; int& end; std::string& name; std::string& type; int &parsed; - push_back_annot(std::list& annot_list_, + push_back_annot(Sequence* parent_seq, + SeqSpanRefListRef children_list, int& begin_, int& end_, std::string& name_, std::string& type_, int &parsed_) - : annot_list(annot_list_), + : parent(parent_seq), + children(children_list), begin(begin_), end(end_), name(name_), @@ -396,8 +394,7 @@ struct push_back_annot { void operator()(std::string::const_iterator, std::string::const_iterator) const { - //std::cout << "adding annot: " << begin << "|" << end << "|" << name << "|" << type << std::endl; - annot_list.push_back(annot(begin, end, name, type)); + children->push_back(parent->make_annotation(name, type, begin, end)); ++parsed; }; }; @@ -446,8 +443,8 @@ Sequence::parse_annot(std::string data, int start_index, int end_index) int end=0; std::string name; std::string type; - std::string seq; - std::list parsed_annots; + std::string seqstr; + SeqSpanRefListRef parsed_annots(new SeqSpanRefList); std::list query_seqs; int parsed=0; @@ -489,13 +486,13 @@ Sequence::parse_annot(std::string data, int start_index, int end_index) ) // to understand how this group gets set // read the comment above struct push_back_annot - )[push_back_annot(parsed_annots, start, end, type, name, parsed)] + )[push_back_annot(this, parsed_annots, start, end, name, type, parsed)] | ((spirit::ch_p('>')|spirit::str_p(">")) >> (*(spirit::print_p))[spirit::assign_a(name)] >> spirit::eol_p >> - (+(spirit::chset<>(Alphabet::nucleic_cstr)))[spirit::assign_a(seq)] - )[push_back_seq(query_seqs, name, seq, parsed)] + (+(spirit::chset<>(Alphabet::nucleic_cstr)))[spirit::assign_a(seqstr)] + )[push_back_seq(query_seqs, name, seqstr, parsed)] ) >> *spirit::space_p ) @@ -506,33 +503,56 @@ Sequence::parse_annot(std::string data, int start_index, int end_index) msg << "Error parsing annotation #" << parsed; throw annotation_load_error(msg.str()); } + // If everything loaded correctly add the sequences to our annotation list // add newly parsed annotations to our sequence - std::copy(parsed_annots.begin(), parsed_annots.end(), std::back_inserter(annots)); - // go seearch for query sequences + std::copy(parsed_annots->begin(), parsed_annots->end(), std::back_inserter(*annotation_list)); + // go search for query sequences find_sequences(query_seqs.begin(), query_seqs.end()); } -void Sequence::add_annotation(const annot& a) +void Sequence::add_annotation(const SeqSpanRef a) { - annots.push_back(a); + annotation_list->push_back(a); } -const std::list& Sequence::annotations() const +void Sequence::add_annotation(std::string name, std::string type, size_type start, size_type stop) { - return annots; + add_annotation(make_annotation(name, type, start, stop)); +} + +SeqSpanRef +Sequence::make_annotation(std::string name, std::string type, size_type start, size_type stop) const +{ + // we want things to be in the positive direction + if (stop < start) { + size_type tmp = start; + start = stop; + stop = tmp; + } + size_type count = stop - start; + SeqSpanRef new_annot(seq->subseq(start, count, SeqSpan::UnknownStrand)); + AnnotationsRef metadata(new Annotations(name)); + metadata->set("type", type); + new_annot->setAnnotations(metadata); + return new_annot; +} + +const SeqSpanRefList& Sequence::annotations() const +{ + return *annotation_list; } void Sequence::copy_children(Sequence &new_seq, size_type start, size_type count) const { new_seq.motif_list = motif_list; - new_seq.annots.clear(); + new_seq.annotation_list.reset(new SeqSpanRefList); - for(std::list::const_iterator annot_i = annots.begin(); - annot_i != annots.end(); + for(SeqSpanRefList::const_iterator annot_i = annotation_list->begin(); + annot_i != annotation_list->end(); ++annot_i) { - size_type annot_begin= annot_i->begin; - size_type annot_end = annot_i->end; + size_type annot_begin= (*annot_i)->start(); + size_type annot_end = (*annot_i)->stop(); if (annot_begin < start+count) { if (annot_begin >= start) { @@ -547,8 +567,9 @@ void Sequence::copy_children(Sequence &new_seq, size_type start, size_type count annot_end = count; } - annot new_annot(annot_begin, annot_end, annot_i->type, annot_i->name); - new_seq.annots.push_back(new_annot); + SeqSpanRef new_annot(seq->subseq(annot_begin, annot_end)); + new_annot->setAnnotations((*annot_i)->annotations()); + new_seq.annotation_list->push_back(new_annot); } } } @@ -562,7 +583,7 @@ Sequence::subseq(size_type start, size_type count, SeqSpan::strand_type strand) return new_seq; } - Sequence new_seq = *this; + Sequence new_seq(*this); new_seq.seq = seq->subseq(start, count, strand); if (seq->annotations()) { AnnotationsRef a(new Annotations(*(seq->annotations()))); @@ -640,15 +661,18 @@ Sequence::clear() seq.reset(); header.clear(); species.clear(); - annots.clear(); + annotation_list.reset(new SeqSpanRefList); motif_list.reset(new MotifList); } void Sequence::save(fs::fstream &save_file) { + std::string type("type"); + std::string empty_str(""); //fstream save_file; - std::list::iterator annots_i; + SeqSpanRefList::iterator annots_i; + AnnotationsRef metadata; // not sure why, or if i'm doing something wrong, but can't seem to pass // file pointers down to this method from the mussa control class @@ -661,10 +685,14 @@ Sequence::save(fs::fstream &save_file) save_file << "" << std::endl; save_file << species << std::endl; - for (annots_i = annots.begin(); annots_i != annots.end(); ++annots_i) + for (annots_i = annotation_list->begin(); + annots_i != annotation_list->end(); + ++annots_i) { - save_file << annots_i->begin << " " << annots_i->end << " " ; - save_file << annots_i->name << " " << annots_i->type << std::endl; + metadata = (*annots_i)->annotations(); + save_file << (*annots_i)->parentStart() << " " << (*annots_i)->parentStop() << " " ; + save_file << metadata->name() << " " + << metadata->getdefault(type, empty_str) << std::endl; } save_file << "" << std::endl; //save_file.close(); @@ -676,11 +704,17 @@ Sequence::load_museq(fs::path load_file_path, int seq_num) fs::fstream load_file; std::string file_data_line; int seq_counter; - annot an_annot; + //annot an_annot; + int annot_begin; + int annot_end; + std::string annot_name; + std::string annot_type; + std::string::size_type space_split_i; std::string annot_value; - annots.clear(); + annotation_list.reset(new SeqSpanRefList); + load_file.open(load_file_path, std::ios::in); seq_counter = 0; @@ -709,39 +743,41 @@ Sequence::load_museq(fs::path load_file_path, int seq_num) // get annot start index space_split_i = file_data_line.find(" "); annot_value = file_data_line.substr(0,space_split_i); - an_annot.begin = atoi (annot_value.c_str()); + annot_begin = atoi (annot_value.c_str()); file_data_line = file_data_line.substr(space_split_i+1); // get annot end index space_split_i = file_data_line.find(" "); annot_value = file_data_line.substr(0,space_split_i); - an_annot.end = atoi (annot_value.c_str()); + annot_end = atoi (annot_value.c_str()); if (space_split_i == std::string::npos) // no entry for type or name { std::cout << "seq, annots - no type or name\n"; - an_annot.type = ""; - an_annot.name = ""; + annot_name = ""; + annot_type = ""; } else // else get annot type { file_data_line = file_data_line.substr(space_split_i+1); space_split_i = file_data_line.find(" "); annot_value = file_data_line.substr(0,space_split_i); - an_annot.type = annot_value; + //an_annot.type = annot_value; + annot_type = annot_value; if (space_split_i == std::string::npos) // no entry for name { std::cout << "seq, annots - no name\n"; - an_annot.name = ""; + annot_name = ""; } else // get annot name { file_data_line = file_data_line.substr(space_split_i+1); space_split_i = file_data_line.find(" "); annot_value = file_data_line.substr(0,space_split_i); - an_annot.type = annot_value; + // this seems like its wrong? + annot_type = annot_value; } } - annots.push_back(an_annot); // don't forget to actually add the annot + add_annotation(annot_name, annot_type, annot_begin, annot_end); } //std::cout << "seq, annots: " << an_annot.start << ", " << an_annot.end // << "-->" << an_annot.type << "::" << an_annot.name << std::endl; @@ -856,7 +892,6 @@ Sequence::motif_scan(const Sequence& a_motif, std::vector * motif_match_sta // end Nora stuff, now we see if a match is found this pass if (motif_i == motif_len) { - annot new_motif; motif_match_starts->push_back(seq_i - motif_len + 1); motif_i = 0; } @@ -871,16 +906,11 @@ void Sequence::add_string_annotation(std::string a_seq, { std::vector seq_starts = find_motif(a_seq); - //std::cout << "searching for " << a_seq << " found " << seq_starts.size() << std::endl; - for(std::vector::iterator seq_start_i = seq_starts.begin(); seq_start_i != seq_starts.end(); ++seq_start_i) { - annots.push_back(annot(*seq_start_i, - *seq_start_i+a_seq.size(), - "", - name)); + add_annotation(name, "", *seq_start_i, *seq_start_i+a_seq.size()); } }