move strand into seqspan

[mussa.git] / alg / sequence.cpp

diff --git a/alg/sequence.cpp b/alg/sequence.cpp
index 47cca2a7dfdfb46c831491ee276f20ec5dd6f4a0..68fa01460542956a6dabf5864786a50f47080737 100644 (file)
--- a/alg/sequence.cpp
+++ b/alg/sequence.cpp
@@ -80,8 +80,7 @@ motif::~motif()
 
 
 Sequence::Sequence(AlphabetRef alphabet)
-  : seq(new SeqSpan("", alphabet)),
-    strand(UnknownStrand)
+  : seq(new SeqSpan("", alphabet, SeqSpan::PlusStrand))
 {
 }
 
@@ -89,25 +88,24 @@ Sequence::~Sequence()
 {
 }
 
-Sequence::Sequence(const char *seq, AlphabetRef alphabet_)
-  : strand(UnknownStrand),
-    header(""),
+Sequence::Sequence(const char *seq, AlphabetRef alphabet_, SeqSpan::strand_type strand_)
+  : header(""),
     species("")
 {
-  set_filtered_sequence(seq, alphabet_);
+  set_filtered_sequence(seq, alphabet_, 0, npos, strand_);
 }
 
-Sequence::Sequence(const std::string& seq, AlphabetRef alphabet_) 
-  : strand(UnknownStrand),
-    header(""),
+Sequence::Sequence(const std::string& seq, 
+                   AlphabetRef alphabet_,
+                   SeqSpan::strand_type strand_) 
+  : header(""),
     species("")
 {
-  set_filtered_sequence(seq, alphabet_);
+  set_filtered_sequence(seq, alphabet_, 0, seq.size(), strand_);
 }
 
 Sequence::Sequence(const Sequence& o)
   : seq(o.seq),
-    strand(o.strand),
     header(o.header),
     species(o.species),
     annots(o.annots),
@@ -117,7 +115,6 @@ Sequence::Sequence(const Sequence& o)
 
 Sequence::Sequence(const Sequence* o)
   : seq(o->seq),
-    strand(o->strand),
     header(o->header),
     species(o->species),
     annots(o->annots),
@@ -127,7 +124,6 @@ Sequence::Sequence(const Sequence* o)
 
 Sequence::Sequence(const SeqSpanRef& seq_ref)
   : seq(seq_ref),
-    strand(UnknownStrand),
     header(""),
     species("")
 {
@@ -137,7 +133,6 @@ Sequence &Sequence::operator=(const Sequence& s)
 {
   if (this != &s) {
     seq = s.seq;
-    strand = s.strand;
     header = s.header;
     species = s.species;
     annots = s.annots;
@@ -220,7 +215,7 @@ Sequence::load_fasta(std::istream& data_file, AlphabetRef a,
   std::string rev_comp;
   std::string sequence_raw;
   std::string seq_tmp;      // holds sequence during basic filtering
-  const Alphabet &alpha = get_alphabet();
+  const Alphabet &alpha = Alphabet::get_alphabet(a);
 
   if (seq_num == 0) {
     throw mussa_load_error("fasta sequence number is 1 based (can't be 0)");
@@ -273,7 +268,7 @@ Sequence::load_fasta(std::istream& data_file, AlphabetRef a,
       std::string msg("The selected sequence appears to be empty"); 
       throw sequence_empty_error(msg);
     }
-    set_filtered_sequence(sequence_raw, a, start_index, end_index-start_index);
+    set_filtered_sequence(sequence_raw, a, start_index, end_index-start_index, SeqSpan::PlusStrand);
   } else {
     std::string errormsg("There were no fasta sequences");
     throw sequence_empty_file_error(errormsg);
@@ -284,7 +279,7 @@ void Sequence::set_filtered_sequence(const std::string &in_seq,
                                      AlphabetRef alphabet_,
                                      size_type start,
                                      size_type count,
-                                     strand_type strand_)
+                                     SeqSpan::strand_type strand_)
 {
   if ( count == npos)
     count = in_seq.size() - start;
@@ -302,9 +297,8 @@ void Sequence::set_filtered_sequence(const std::string &in_seq,
       new_seq.append(1, 'N');
     }
   }
-  SeqSpanRef new_seq_ref(new SeqSpan(new_seq, alphabet_)); 
+  SeqSpanRef new_seq_ref(new SeqSpan(new_seq, alphabet_, strand_)); 
   seq = new_seq_ref;
-  strand = strand_;
 }
 
 void
@@ -547,15 +541,16 @@ void Sequence::copy_children(Sequence &new_seq, size_type start, size_type count
   
 }
 Sequence
-Sequence::subseq(size_type start, size_type count) const
+Sequence::subseq(size_type start, size_type count, SeqSpan::strand_type strand) const
 {
+  // FIXME: should i really allow a subsequence of an empty sequence?
   if (!seq) {
     Sequence new_seq;
     return new_seq;
   }
 
   Sequence new_seq = *this;
-  new_seq.seq = seq->subseq(start, count);
+  new_seq.seq = seq->subseq(start, count, strand);
 
   copy_children(new_seq, start, count);
   
@@ -566,21 +561,8 @@ Sequence::subseq(size_type start, size_type count) const
 // FIXME: This needs to be moved into SeqSpan
 Sequence Sequence::rev_comp() const
 {
-  std::string rev_comp;
-  rev_comp.reserve(length());
-  
-  std::string rc_map = seq->seq->create_reverse_complement_map();
-
-  // reverse and convert
-  Sequence::const_reverse_iterator seq_i;
-  Sequence::const_reverse_iterator seq_end = rend();  
-  for(seq_i = rbegin(); 
-      seq_i != seq_end;
-      ++seq_i)
-  {
-    rev_comp.append(1, rc_map[*seq_i]);
-  }
-  return Sequence(rev_comp, seq->seq->get_alphabet_ref()); 
+  // a reverse complement is the whole opposite strand 
+  return subseq(0, npos, SeqSpan::OppositeStrand);
 }
 
 const Alphabet& Sequence::get_alphabet() const
@@ -623,7 +605,7 @@ Sequence::get_name() const
 
 void Sequence::set_sequence(const std::string& s, AlphabetRef a) 
 {
-  set_filtered_sequence(s, a);
+  set_filtered_sequence(s, a, 0, s.size(), SeqSpan::PlusStrand);
 }
 
 std::string Sequence::get_sequence() const
@@ -640,7 +622,6 @@ void
 Sequence::clear()
 {
   seq.reset();
-  strand = UnknownStrand;
   header.clear();
   species.clear();
   annots.clear();
@@ -696,7 +677,7 @@ Sequence::load_museq(fs::path load_file_path, int seq_num)
   }
   getline(load_file, file_data_line);
   // looks like the sequence is written as a single line
-  set_filtered_sequence(file_data_line, reduced_dna_alphabet);
+  set_filtered_sequence(file_data_line, reduced_dna_alphabet, 0, file_data_line.size(), SeqSpan::PlusStrand);
   getline(load_file, file_data_line);
   getline(load_file, file_data_line);
   if (file_data_line == "<Annotations>")