Added newline char at end of file to remove gcc warning

[mussa.git] / alg / sequence.cpp

diff --git a/alg/sequence.cpp b/alg/sequence.cpp
index e1923c68d49e484ac769800c66fccff7952ef6a5..147fdce3b7e2c015eb20741964327954f70a924a 100644 (file)
--- a/alg/sequence.cpp
+++ b/alg/sequence.cpp
@@ -193,23 +193,29 @@ void Sequence::load_fasta(fs::path file_path, alphabet_ref a,
       errormsg << file_path.native_file_string()
                << " did not have any fasta sequences" << std::endl;
       throw sequence_empty_file_error(errormsg.str());
+    } catch(sequence_invalid_load_error e) {
+      std::ostringstream msg;
+      msg << file_path.native_file_string();
+      msg << " " << e.what();
+      throw sequence_invalid_load_error(msg.str());
     }
   }
 }
 
-void Sequence::load_fasta(std::iostream& file, 
+void Sequence::load_fasta(std::istream& file, 
                           int seq_num, int start_index, int end_index)
 {
   load_fasta(file, alphabet, seq_num, start_index, end_index);
 }
 
 void
-Sequence::load_fasta(std::iostream& data_file, alphabet_ref a, 
+Sequence::load_fasta(std::istream& data_file, alphabet_ref a, 
                      int seq_num, 
                      int start_index, int end_index)
 {
   std::string file_data_line;
   int header_counter = 0;
+  size_t line_counter = 0;
   bool read_seq = true;
   std::string rev_comp;
   std::string sequence_raw;
@@ -224,6 +230,7 @@ Sequence::load_fasta(std::iostream& data_file, alphabet_ref a,
   while ( (!data_file.eof()) && (header_counter < seq_num) )
   {
     multiplatform_getline(data_file, file_data_line);
+    ++line_counter;
     if (file_data_line.substr(0,1) == ">")
       header_counter++;
   }
@@ -235,6 +242,7 @@ Sequence::load_fasta(std::iostream& data_file, alphabet_ref a,
 
     while ( !data_file.eof() && read_seq ) {
       multiplatform_getline(data_file,file_data_line);
+      ++line_counter;
       if (file_data_line.substr(0,1) == ">")
         read_seq = false;
       else {
@@ -245,7 +253,10 @@ Sequence::load_fasta(std::iostream& data_file, alphabet_ref a,
            if(alpha.exists(*line_i)) {
              sequence_raw += *line_i;
            } else {
-            throw sequence_invalid_load_error("Unrecognized characters in fasta sequence");
+            std::ostringstream msg;
+            msg << "Unrecognized characters in fasta sequence at line ";
+            msg << line_counter;
+            throw sequence_invalid_load_error(msg.str());
            }
          }
       }
@@ -287,7 +298,7 @@ void Sequence::set_filtered_sequence(const std::string &in_seq,
   for(size_type i = 0; i != count; ++i, ++seq_i)
   {
     if (alpha_impl.exists(*seq_i)) {
-      new_seq->append(1, *seq_i);
+      new_seq->append(1, toupper(*seq_i));
     } else {
       new_seq->append(1, 'N');
     }
@@ -318,8 +329,16 @@ Sequence::load_annot(fs::path file_path, int start_index, int end_index)
     data.push_back(c);
   }
   data_stream.close();
-        
-  parse_annot(data, start_index, end_index);
+
+  try {  
+    parse_annot(data, start_index, end_index);
+  } catch(annotation_load_error e) {
+    std::ostringstream msg;
+    msg << file_path.native_file_string()
+        << " "
+        << e.what();
+    throw annotation_load_error(msg.str());
+  }
 }
 
 /* If this works, yikes, this is some brain hurting code.
@@ -418,7 +437,7 @@ Sequence::parse_annot(std::string data, int start_index, int end_index)
   std::string seq;
   std::list<annot> parsed_annots;
   std::list<Sequence> query_seqs;
-  int parsed=1;
+  int parsed=0;
 
   bool ok = spirit::parse(data.begin(), data.end(),
               (
@@ -463,7 +482,7 @@ Sequence::parse_annot(std::string data, int start_index, int end_index)
                     ((spirit::ch_p('>')|spirit::str_p("&gt;")) >> 
                        (*(spirit::print_p))[spirit::assign_a(name)] >>
                        spirit::eol_p >> 
-                       (+(spirit::chset<>(Alphabet::nucleic_alphabet.c_str())))[spirit::assign_a(seq)]
+                       (+(spirit::chset<>(Alphabet::nucleic_cstr)))[spirit::assign_a(seq)]
                      )[push_back_seq(query_seqs, name, seq, parsed)]
                     ) >>
                     *spirit::space_p
@@ -580,9 +599,9 @@ Sequence Sequence::rev_comp() const
   std::string rc_map = create_reverse_map();
 
   // reverse and convert
-  seq_string::const_reverse_iterator seq_i;
-  seq_string::const_reverse_iterator seq_end = seq->rend();  
-  for(seq_i = seq->rbegin(); 
+  Sequence::const_reverse_iterator seq_i;
+  Sequence::const_reverse_iterator seq_end = rend();  
+  for(seq_i = rbegin(); 
       seq_i != seq_end;
       ++seq_i)
   {
@@ -633,15 +652,15 @@ const Alphabet& Sequence::get_alphabet(alphabet_ref alpha) const
 {
   switch (alpha) {
     case reduced_dna_alphabet:
-      return Alphabet::reduced_dna_alphabet;
+      return Alphabet::reduced_dna_alphabet();
     case reduced_rna_alphabet:
-      return Alphabet::reduced_rna_alphabet;
+      return Alphabet::reduced_rna_alphabet();
     case reduced_nucleic_alphabet:
-      return Alphabet::reduced_nucleic_alphabet;
+      return Alphabet::reduced_nucleic_alphabet();
     case nucleic_alphabet:
-      return Alphabet::nucleic_alphabet;
+      return Alphabet::nucleic_alphabet();
     case protein_alphabet:
-      return Alphabet::protein_alphabet;    
+      return Alphabet::protein_alphabet();    
     default:
       throw std::runtime_error("unrecognized alphabet type");
       break;
@@ -711,11 +730,52 @@ Sequence::const_iterator Sequence::end() const
   }
 }
 
+Sequence::const_reverse_iterator Sequence::rbegin() const
+{
+  if (seq and seq_count != 0)
+    return seq->rbegin()+(seq->size()-(seq_start+seq_count));
+  else 
+    return Sequence::const_reverse_iterator();
+}
+
+Sequence::const_reverse_iterator Sequence::rend() const
+{
+  if (seq and seq_count != 0) {
+    return rbegin() + seq_count;
+  } else {
+    return Sequence::const_reverse_iterator();
+  }
+}
+
 bool Sequence::empty() const
 {
   return (seq_count == 0) ? true : false;
 }
 
+Sequence::size_type Sequence::find_first_not_of(
+  const std::string& query, 
+  Sequence::size_type index)
+{
+  typedef std::set<std::string::value_type> sequence_set;
+  sequence_set match_set;
+  
+  for(const_iterator query_item = query.begin();
+      query_item != query.end();
+      ++query_item)
+  {
+    match_set.insert(*query_item);
+  }  
+  for(const_iterator base = begin();
+      base != end();
+      ++base)
+  {
+    if(match_set.find(*base) == match_set.end()) {
+      return base-begin();
+    } 
+  }
+  return Sequence::npos;
+}
+ 
 Sequence::size_type Sequence::start() const
 {
   if (parent)
@@ -904,12 +964,12 @@ Sequence::motif_scan(const Sequence& a_motif, std::vector<int> * motif_match_sta
   Sequence::value_type motif_char;
   Sequence::value_type seq_char;
 
-  while (seq_i < seq->length())
+  while (seq_i < size())
   {
     // this is pretty much a straight translation of Nora's python code
     // to match iupac letter codes
     motif_char = toupper(a_motif[motif_i]);
-    seq_char = toupper(seq->at(seq_i));
+    seq_char = toupper(seq->at(seq_start+seq_i));
     if (motif_char =='N')
       motif_i++;
     else if (motif_char == seq_char)