Commits

opticall committed e134fb7

update fixing bug when all internsities for a SNP are zero

Comments (0)

Files changed (1)

opticall/opticall.cpp

 	MatrixXd scaling_sd = (sd_int_data.array() / prior_sd.array());
 	MatrixXd scaling_matrix = scaling_sd.transpose() * scaling_sd;
 
-	/*
-	cout << "Scaling matrix " << endl;
-	cout << scaling_matrix << endl;
-	*/
+	
+	//cout << "Scaling matrix " << endl;
+	//cout << scaling_matrix << endl;
+	
 
 	for (int j = 0; j < numclasses-numcustomclasses; j++ ) {
 		priorcovars_scaled.push_back( priorcovars[j].array()*scaling_matrix.array() );
 
 		if (nan_samples.size() == data_count) {
 			//all the data for this SNP is NAN
-			cout << "all data for snp is NaN, calling everything unknown"  << endl;;
+			cout << "all intensity data for snp is NaN, calling everything unknown"  << endl;;
 			MatrixXd post_probs = probs_with_nans(4, int_data_w_nans.rows() );
 			//cout << post_probs << endl;
 			//then output the posterior probs to file
 		}
 
 
+		if (int_data.norm() == 0 || sd_of_data(int_data).norm() == 0) {
+			//all the data for this SNP is 0
+			cout << "all intensity data for snp is zero, calling everything unknown"  << endl;;
+			MatrixXd post_probs = probs_with_nans(4, int_data_w_nans.rows() );
+			//cout << post_probs << endl;
+			//then output the posterior probs to file
+			outputpostprobs(post_probs,out_probsfile,snpinfo2[l]);
+			outputcalls(post_probs,out_callsfile,snpinfo2[l]);
+			
+			cout << "line step took" << (double)(clock() - linetStart)/CLOCKS_PER_SEC  << endl;
+			
+			//go to the next SNP
+			continue;
+		}
+		
+		
+		
+		
+		
+		
 
 		MatrixXd post_probs = genotype_mm_line_caller_w_priors(int_data, line_vs, line_custommus, line_customcovars, blockmus, blockcovars, blockjs, blockvs, sd_sample);