Getting closer to a subanalysis mode
[mussa.git] / alg / mussa.cpp
index 4677a6c84e871c0db9cac12dfa0ae20eedf76f18..79078b2cb1447b3ba4c2ca46d0a3500116fc2337 100644 (file)
@@ -27,6 +27,7 @@ using namespace std;
 
 
 Mussa::Mussa()
+  : color_mapper(new AnnotationColors)
 {
   clear();
   connect(&the_paths, SIGNAL(progress(const std::string&, int, int)), 
@@ -60,7 +61,7 @@ Mussa::clear()
   win_append = false;
   thres_append = false;
   motif_sequences.clear();
-  color_mapper.clear();
+  if(color_mapper) color_mapper->clear();
   the_paths.clear();
 }
 
@@ -173,7 +174,7 @@ void Mussa::createLocalAlignment(std::list<ConservedPath>::iterator begin,
                                  std::list<ConservedPath::path_type>& result,
                                  std::list<std::vector<bool> >& reversed)
 {
-  const vector<Sequence>& raw_seq = the_seqs;
+  const vector_sequence_type& raw_seq = the_seqs;
   ConservedPath::path_type aligned_path;
   size_t i2, i3;
   int x_start, x_end;
@@ -213,15 +214,17 @@ void Mussa::createLocalAlignment(std::list<ConservedPath>::iterator begin,
         x_start = (abs(a_path[i2]-rc_1+win_i));
         x_end =   (abs(a_path[i2+1]-rc_2+win_i));
         
+        boost::shared_ptr<Sequence> cur(raw_seq[i2]) ;
+        boost::shared_ptr<Sequence> next(raw_seq[i2+1]);
         // RC case handling
         // ugh, and xor...only want rc coloring if just one of the nodes is rc
         // if both nodes are rc, then they are 'normal' relative to each other
         if((rc_list[i2] || rc_list[i2+1] )&&!(rc_list[i2] && rc_list[i2+1]))
         { //the hideous rc matching logic - not complex, but annoying
-          if(!(( (raw_seq[i2][x_start]=='A')&&(raw_seq[i2+1][x_end]=='T')) ||
-                ((raw_seq[i2][x_start]=='T')&&(raw_seq[i2+1][x_end]=='A')) ||
-                ((raw_seq[i2][x_start]=='G')&&(raw_seq[i2+1][x_end]=='C')) ||
-                ((raw_seq[i2][x_start]=='C')&&(raw_seq[i2+1][x_end]=='G'))) ) 
+          if(!(( ((*cur)[x_start]=='A')&&((*next)[x_end]=='T')) ||
+                (((*cur)[x_start]=='T')&&((*next)[x_end]=='A')) ||
+                (((*cur)[x_start]=='G')&&((*next)[x_end]=='C')) ||
+                (((*cur)[x_start]=='C')&&((*next)[x_end]=='G'))) ) 
           {
             full_match = false;
           } else {
@@ -230,9 +233,8 @@ void Mussa::createLocalAlignment(std::list<ConservedPath>::iterator begin,
         }
         else
         { // forward match
-          if (!( (raw_seq[i2][x_start] == raw_seq[i2+1][x_end]) &&
-                (raw_seq[i2][x_start] != 'N') &&
-                (raw_seq[i2+1][x_end] != 'N') ) ) {
+          if (!( ((*cur)[x_start] == (*next)[x_end]) &&
+                ((*cur)[x_start] != 'N') && ((*next)[x_end] != 'N') ) ) {
             full_match = false;
           } else {
             aligned_path.push_back(x_start);
@@ -253,12 +255,19 @@ void Mussa::createLocalAlignment(std::list<ConservedPath>::iterator begin,
 }
 
 
-void Mussa::append_sequence(Sequence a_seq)
+void Mussa::append_sequence(const Sequence& a_seq)
+{
+  boost::shared_ptr<Sequence> seq_copy(new Sequence(a_seq));
+  the_seqs.push_back(seq_copy);
+}
+
+void Mussa::append_sequence(boost::shared_ptr<Sequence> a_seq)
 {
   the_seqs.push_back(a_seq);
 }
 
-const vector<Sequence>& 
+
+const vector<boost::shared_ptr<Sequence> >& 
 Mussa::sequences() const
 {
   return the_seqs;
@@ -267,10 +276,10 @@ Mussa::sequences() const
 void Mussa::load_sequence(fs::path seq_file, fs::path annot_file, 
                           int fasta_index, int sub_seq_start, int sub_seq_end)
 {
-  Sequence aseq;
-  aseq.load_fasta(seq_file, fasta_index, sub_seq_start, sub_seq_end);
+  boost::shared_ptr<Sequence> aseq(new Sequence);
+  aseq->load_fasta(seq_file, fasta_index, sub_seq_start, sub_seq_end);
   if ( not annot_file.empty() ) {
-    aseq.load_annot(annot_file, sub_seq_start, sub_seq_end);
+    aseq->load_annot(annot_file, sub_seq_start, sub_seq_end);
   }
   the_seqs.push_back(aseq);
 }
@@ -451,7 +460,7 @@ Mussa::seqcomp()
       all_comps[i].push_back(dummy_comp);
   }
   for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++) {
-    seq_lens.push_back(the_seqs[i].size());
+    seq_lens.push_back(the_seqs[i]->size());
   }
   int seqcomps_done = 0;
   int seqcomps_todo = (the_seqs.size() * (the_seqs.size()-1)) / 2;
@@ -462,8 +471,8 @@ Mussa::seqcomp()
     {
       //cout << "seqcomping: " << i << " v. " << i2 << endl;
       all_comps[i][i2].setup(window, threshold);
-      all_comps[i][i2].seqcomp(the_seqs[i], the_seqs[i2], false);
-      all_comps[i][i2].seqcomp(the_seqs[i], the_seqs[i2].rev_comp(),true);
+      all_comps[i][i2].seqcomp(*the_seqs[i], *the_seqs[i2], false);
+      all_comps[i][i2].seqcomp(*the_seqs[i], the_seqs[i2]->rev_comp(),true);
       ++seqcomps_done;
       emit progress("seqcomp", seqcomps_done, seqcomps_todo);
     }
@@ -472,7 +481,6 @@ Mussa::seqcomp()
 void
 Mussa::nway()
 {
-  vector<string> some_Seqs;
 
   the_paths.set_soft_threshold(soft_thres);
 
@@ -484,11 +492,12 @@ Mussa::nway()
   }
   else if (ana_mode == EntropyNway)
   {
+    vector<string> some_Seqs;
     //unlike other methods, entropy needs to look at the sequence at this stage
     some_Seqs.clear();
     for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
     {
-      some_Seqs.push_back(the_seqs[i]);
+      some_Seqs.push_back(*the_seqs[i]);
     }
 
     the_paths.setup_ent(ent_thres, some_Seqs); // ent analysis extra setup
@@ -541,7 +550,7 @@ Mussa::save()
 
     for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
     {
-      the_seqs[i].save(save_file);
+      the_seqs[i]->save(save_file);
     }
 
     save_file << "</Mussa_Sequence>" << endl;
@@ -574,7 +583,6 @@ Mussa::load(fs::path ana_file)
   fs::path a_file_path; 
   fs::path ana_path(ana_file);
   bool parsing_path;
-  Sequence tmp_seq;
   string err_msg;
   ostringstream append_info;
   vector<FLPs> empty_FLP_vector;
@@ -600,9 +608,9 @@ Mussa::load(fs::path ana_file)
   // this is a bit of a hack due to C++ not acting like it should with files
   for (i = 1; i <= seq_num; i++)
   {
-    tmp_seq.clear();
+    boost::shared_ptr<Sequence> tmp_seq(new Sequence);
     //cout << "mussa_class: loading museq frag... " << a_file_path.string() << endl;
-    tmp_seq.load_museq(a_file_path, i);
+    tmp_seq->load_museq(a_file_path, i);
     the_seqs.push_back(tmp_seq);
   }
   
@@ -638,7 +646,7 @@ Mussa::save_old()
   save_file.open(analysis_name, ios::out);
 
   for(vector<Sequence>::size_type i = 0; i < the_seqs.size(); i++)
-    save_file << the_seqs[i] << endl;
+    save_file << *(the_seqs[i]) << endl;
 
   save_file << window << endl;
   save_file.close();
@@ -668,7 +676,7 @@ Mussa::load_old(char * load_file_path, int s_num)
   for(i = 0; i < seq_num; i++)
   {
     getline(save_file, file_data_line);
-    a_seq = file_data_line;
+    boost::shared_ptr<Sequence> a_seq(new Sequence(file_data_line));
     the_seqs.push_back(a_seq);
   }
 
@@ -712,7 +720,7 @@ void Mussa::add_motifs(const vector<string>& motifs,
   for(size_t i = 0; i != motifs.size(); ++i)
   {
     motif_sequences.insert(motifs[i]);
-    color_mapper.appendInstanceColor("motif", motifs[i], colors[i]);
+    color_mapper->appendInstanceColor("motif", motifs[i], colors[i]);
   }
   update_sequences_motifs();
 }
@@ -758,7 +766,7 @@ void Mussa::load_motifs(std::istream &in)
       // sequence wasn't found
       motif_sequences.insert(seq);
       Color c(red, green, blue);
-      color_mapper.appendInstanceColor("motif", seq, c);
+      color_mapper->appendInstanceColor("motif", seq, c);
     } else {
       clog << "sequence " << seq << " was already defined skipping" 
            << endl;
@@ -779,18 +787,18 @@ void Mussa::update_sequences_motifs()
 {
   // once we've loaded all the motifs from the file, 
   // lets attach them to the sequences
-  for(vector<Sequence>::iterator seq_i = the_seqs.begin();
+  for(vector<boost::shared_ptr<Sequence> >::iterator seq_i = the_seqs.begin();
       seq_i != the_seqs.end();
       ++seq_i)
   {
     // clear out old motifs
-    seq_i->clear_motifs();
+    (*seq_i)->clear_motifs();
     // for all the motifs in our set, attach them to the current sequence
     for(set<string>::iterator motif_i = motif_sequences.begin();
         motif_i != motif_sequences.end();
         ++motif_i)
     {
-      seq_i->add_motif(*motif_i);
+      (*seq_i)->add_motif(*motif_i);
     }
   }
 }
@@ -800,7 +808,7 @@ const set<string>& Mussa::motifs() const
   return motif_sequences;
 }
 
-AnnotationColors& Mussa::colorMapper()
+boost::shared_ptr<AnnotationColors> Mussa::colorMapper()
 {
   return color_mapper;
 }