Merge branch 'bugfix-pick-match'

andrej-fischer · Apr 25, 2014 · e751fd5 · e751fd5
2 parents aabcd3a + c25725f
commit e751fd5
Showing 1 changed file with 28 additions and 20 deletions.
diff --git a/src/pre-filter.cpp b/src/pre-filter.cpp
@@ -65,7 +65,7 @@ void get_opts( int argc, const char ** argv, cmdl_opts& opts);
 void test_opts(cmdl_opts& opts);
 void default_opts(cmdl_opts& opts);
 void pre_filter( Emission& dataEmit, cmdl_opts& opts);  
-void pick_match( Emission& pickEmit, Emission& matchEmit, cmdl_opts& opts);
+void pick_match( Emission * pickEmit, Emission * matchEmit, cmdl_opts& opts);
 
 
 // *** MAIN START***
@@ -97,7 +97,7 @@ int main (int argc, const char * argv[]){
     get_data( opts.match_fn, &matchEmit);
   }
   if (dataEmit.is_set) pre_filter( dataEmit, opts); 
-  if (pickEmit.is_set) pick_match( pickEmit, matchEmit, opts);
+  if (pickEmit.is_set) pick_match( &pickEmit, &matchEmit, opts);
   //done
   return (0);
 }
@@ -325,43 +325,51 @@ void pre_filter( Emission& dataEmit, cmdl_opts& opts){
 }
 
 
-void pick_match( Emission& pickEmit, Emission& matchEmit, cmdl_opts& opts){
+void pick_match( Emission * pickEmit, Emission * matchEmit, cmdl_opts& opts){
   char buff[1024];
   sprintf( buff, "%s.pref.txt", opts.pre);
   FILE * pref_fp  = fopen(buff,"w"); 
-  for (int s=0; s < pickEmit.nSamples; s++){
-    int chr=pickEmit.chr[s];
-    if ( matchEmit.chrs.count(chr) == 0){
+  for (int s=0; s < pickEmit->nSamples; s++){
+    int chr=pickEmit->chr[s];
+    if ( matchEmit->chrs.count(chr) == 0){
       printf( "ERROR: chr %i in %s could not be found in %s\n",
 	      chr, opts.pick_fn, opts.match_fn
 	      );
       exit(1);
     }
-    int matchSample = matchEmit.idx_of[chr];
+    int matchSample = matchEmit->idx_of[chr];
     // get all bins and determine bin width...
     // !!!NOTE: assumes uniform bin width!!!
     std::map<int,int> diffs;
-    unsigned int * mloci = matchEmit.loci[matchSample];
-    for (int l=1; l<matchEmit.nSites[matchSample]; l++){
-      diffs[int(mloci[l] - mloci[l-1])]++;
+    unsigned int * mloci = matchEmit->loci[matchSample];
+    int binw=0;
+    for (int l=1; l<matchEmit->nSites[matchSample]; l++){
+      binw = int(mloci[l]) - int(mloci[l-1]);
+      if (diffs.count(binw) == 0) diffs.insert(std::pair<int,int>(binw,0));
+      diffs[binw] =+ 1;
     }
-    int binw=0, ct=0;
+    int ct = 0;
     std::map<int,int>::iterator it;
-    for (it = diffs.begin(); it != diffs.end(); ++it){
+    for (it = diffs.begin(); it != diffs.end(); it++){
       if (it->second > ct){
 	binw = it->first;
 	ct   = it->second;
       }
     }
     //now pick...
-    int midx=0;
-    unsigned int * ploci = pickEmit.loci[s];
-    for (int idx=0; idx<pickEmit.nSites[s]; idx++){
-      while (mloci[midx] < ploci[idx]) midx++; 
-      if (int(ploci[idx]) <= int(mloci[midx]) - binw) continue;
-      fprintf( pref_fp, "%-2i %12i", pickEmit.chr[s], pickEmit.loci[s][idx]);
-      for (int t=0;t<pickEmit.nTimes; t++){
-	fprintf( pref_fp, " %-3i %-3i", pickEmit.reads[t][s][idx], pickEmit.depths[t][s][idx]);
+    int midx = 0;
+    int mlocus = (int) mloci[0];
+    for (int idx=0; idx<pickEmit->nSites[s]; idx++){
+      int plocus = (int) pickEmit->loci[s][idx];
+      while (mlocus < plocus){
+	midx++;
+	if (midx == matchEmit->nSites[matchSample]) break;
+	mlocus = (int) mloci[midx];	
+      }
+      if ( mlocus < plocus || plocus <= mlocus - binw ) continue;
+      fprintf( pref_fp, "%-2i %12i", chr, plocus);
+      for ( int t=0; t<pickEmit->nTimes; t++){
+	fprintf( pref_fp, " %-3i %-3i", pickEmit->reads[t][s][idx], pickEmit->depths[t][s][idx]);
       }
       fprintf( pref_fp, "\n");
     }