bugfix for chrs with no non-zero obs?

ivazquez · May 27, 2014 · 533772e · 533772e
1 parent d37754f
commit 533772e
Show file tree

Hide file tree

Showing 3 changed files with 12 additions and 10 deletions.
diff --git a/src/clone-prior.cpp b/src/clone-prior.cpp
@@ -48,6 +48,8 @@ void Clone::set_cna_prior( gsl_vector * prior, int sample){
 //used in SNV-only mode, w/o correlation and w/o cn-info...
 void Clone::initialize_snv_prior_param(){// SNV prior, conditional on max-tcn
   if (initial_snv_prior_param != NULL) gsl_matrix_free(initial_snv_prior_param);
+  initial_snv_prior_param = NULL;
+  if (nClones == 0) return;
   initial_snv_prior_param = gsl_matrix_calloc( maxtcn+1, maxtcn+1);
   gsl_matrix_set( initial_snv_prior_param, 0, 0, snv_fpr);
   double p = snv_pen_high;// penalty for higher genotypes

diff --git a/src/cloneHD-inference.cpp b/src/cloneHD-inference.cpp
@@ -254,7 +254,7 @@ int infer_clones( gsl_matrix * Clones, gsl_vector * Mass, Clone * myClone, cmdl_
   if (best_mass[bestn] != NULL){
     myClone->set_mass(best_mass[bestn]);
   }
-  if (!cnaEmit->is_set && snvEmit->is_set && !snvEmit->connect){
+  if (!cnaEmit->is_set && snvEmit->is_set && !snvEmit->connect && bestn > 0){
     myClone->initialize_snv_prior_param();
   }
   if (best_priors[bestn] != NULL && bestn > 0){

diff --git a/src/common-functions.cpp b/src/common-functions.cpp
@@ -81,7 +81,7 @@ void get_data( const char * data_fn, Emission * myEmit){
   }
   int ct=0,l;
   int chr=0,old=-1, sample=0;
-  int d,r, keep=0;
+  int d,r, keep=0, wait=0;
   //now collect all data...
   while( data_ifs.good()){
     line.clear();
@@ -93,17 +93,17 @@ void get_data( const char * data_fn, Emission * myEmit){
     line_ss >> chr >> l;//chromosome and locus
     if (chr != old){
       if (myEmit->chrs.count(chr) == 0){
-	cout<<"ERROR in get_data(): chr "<<chr<<" was not expected"<<endl<<line<<endl;
-	for (int s=0; s<myEmit->nSamples; s++){
-	  printf("sample %i = chr %i, idx = %i\n", 
-		 s+1, myEmit->chr[s], myEmit->idx_of[myEmit->chr[s]]);
-	}
-	exit(1);
+	printf("WARNING: chr %2i in file %s will be ignored.\n", chr, data_fn);
+	wait = 1;
+      }
+      else{
+	sample = myEmit->idx_of[chr];
+	ct  = 0;
+	wait = 0;
       }
-      sample = myEmit->idx_of[chr];
-      ct  = 0;
       old = chr;
     }
+    if (wait) continue;
     if (ct >= myEmit->nSites[sample]) continue;
     keep = 0;
     for (int t=0; t<myEmit->nTimes; t++){