Commits

opticall committed 9415f36

added option to alter meanintensity standard deviation limit

  • Participants
  • Parent commits 9a9fce7

Comments (0)

Files changed (1)

File opticall/opticall.cpp

 	bool noblank = false;
 	bool meanintfilter = false;
 	bool outliers_only = false;
+	bool X_chrom = false;
 
 	if (argc < 5) { // Check the value of argc. If not enough parameters have been passed, inform user and exit.
 		cout << "Usage is -in <infile> -out <outdir>\n"; // Inform the user of how to use the program
 			os >> G_CALL_THRESHOLD;
 			cout << "probability call threshold is: " << G_CALL_THRESHOLD << endl;
 			i++;
+		} else if (string(argv[i]) == "-meanintsds") {
+			string meanintsds = string(argv[i + 1]);
+			istringstream os(meanintsds);
+			os >> initial_outliers_numsds;
+			cout << "meanintensity filtering limit is: " << initial_outliers_numsds << " standard deviations" << endl;
+			i++;
 		} else if (string(argv[i]) == "-nointcutoff") {
 			cout << "turning off internal max-intensity cutoff" << endl;
 			nonorm = true;
 		} else if (string(argv[i]) == "-noblank") {
 			cout << "turning off HWE based SNP blanking" << endl;
 			noblank = true;
+		} else if (string(argv[i]) == "-X") {
+			cout << "X chromosome calling mode" << endl;
+			X_chrom = true;
 		} else {
 			cout << "Not enough or invalid arguments, please try again.\n";
 			cout << i << endl;
 		//analyse these probs for HWE-p
 		//go from probs to calls
 
-		double hwep = hwepvalue(calls_from_posterior_probs(post_probs));
-		cout << "HWEP " << snpinfo2[l][0] << " " << hwep << endl;
-
-
-		/****
-		 * STEP 3 - remove sampled contour if HWE-p test failed
-		 */
+		if (!X_chrom)
+		{
 		
-		//5*1e-15
-
-		if (  isnan(hwep) || hwep < resc_threshhold) {
-			cout << "Out of HWE, running Rescue" << endl;
-
-			//ok so how will I do this?
-			//=============================
-			//BEGINNING OF HACKY HACK HACK - to fix shifted intensities
-
-
-			/*
-			int data_count_ill_fix = lineblock.rows()*lineblock.cols()/2;
-			//fix the data from a line
-			MatrixXd line_data_T(2,data_count_ill_fix);
-			MatrixXd line_data;
-			line_data_T = Map<MatrixXd>(lineblock.data(),2,data_count_ill_fix);
-			line_data = line_data_T.transpose();
-			*/
-			MatrixXd line_data = transform_data(int_data);
-
-			MatrixXd blockmus_ill_fix;
-			vector<MatrixXd> blockcovars_ill_fix;
-		    MatrixXd blockjs_ill_fix;
-		    MatrixXd blockvs_ill_fix;
-
-			MatrixXd vs_ill_fix(4,1);
-			vector<MatrixXd> custommus_ill_fix;
-			vector<MatrixXd> customcovars_ill_fix;
-			vs_ill_fix << 1,1,1.3,1;
-			custommus_ill_fix.push_back(Matrix<double, 1, 2>::Zero());
-			customcovars_ill_fix.push_back(Matrix<double, 2, 2>::Identity() * 100);
-
-
-			//run a mixture model to get back the distributions
-
-
-			//genotype_mm_block_caller(line_data, vs_ill_fix, custommus_ill_fix, customcovars_ill_fix, blockmus_ill_fix, blockcovars_ill_fix, blockjs_ill_fix, blockvs_ill_fix);
-
-
-			clock_t tStart = clock();
-			//cout << "rescue calling step" << endl;
-			genotype_mm_rescue_caller(line_data, vs_ill_fix, custommus_ill_fix, customcovars_ill_fix,
+			double hwep = hwepvalue(calls_from_posterior_probs(post_probs));
+			cout << "HWEP " << snpinfo2[l][0] << " " << hwep << endl;
+
+
+			/****
+			 * STEP 3 - remove sampled contour if HWE-p test failed
+			 */
+		
+			//5*1e-15
+
+			if (  isnan(hwep) || hwep < resc_threshhold) {
+				cout << "Out of HWE, running Rescue" << endl;
+
+				//ok so how will I do this?
+				//=============================
+				//BEGINNING OF HACKY HACK HACK - to fix shifted intensities
+
+
+				/*
+				 int data_count_ill_fix = lineblock.rows()*lineblock.cols()/2;
+				 //fix the data from a line
+				 MatrixXd line_data_T(2,data_count_ill_fix);
+				 MatrixXd line_data;
+				 line_data_T = Map<MatrixXd>(lineblock.data(),2,data_count_ill_fix);
+				 line_data = line_data_T.transpose();
+				 */
+				MatrixXd line_data = transform_data(int_data);
+
+				MatrixXd blockmus_ill_fix;
+				vector<MatrixXd> blockcovars_ill_fix;
+				MatrixXd blockjs_ill_fix;
+				MatrixXd blockvs_ill_fix;
+
+				MatrixXd vs_ill_fix(4,1);
+				vector<MatrixXd> custommus_ill_fix;
+				vector<MatrixXd> customcovars_ill_fix;
+				vs_ill_fix << 1,1,1.3,1;
+				custommus_ill_fix.push_back(Matrix<double, 1, 2>::Zero());
+				customcovars_ill_fix.push_back(Matrix<double, 2, 2>::Identity() * 100);
+
+
+				//run a mixture model to get back the distributions
+
+
+				//genotype_mm_block_caller(line_data, vs_ill_fix, custommus_ill_fix, customcovars_ill_fix, blockmus_ill_fix, blockcovars_ill_fix, blockjs_ill_fix, blockvs_ill_fix);
+
+
+				clock_t tStart = clock();
+				//cout << "rescue calling step" << endl;
+				genotype_mm_rescue_caller(line_data, vs_ill_fix, custommus_ill_fix, customcovars_ill_fix,
 					blockmus_ill_fix,  blockcovars_ill_fix,  blockjs_ill_fix,  blockvs_ill_fix,
 					blockmus, blockcovars, blockjs, blockvs, sd_sample);
-			cout << "rescue step took " << (double)(clock() - tStart)/CLOCKS_PER_SEC  << endl;
-
-			//cout << "===============" << endl;
-			//cout << blockmus_ill_fix << endl;
-			//cout << blockjs_ill_fix << endl;
-			//cout << "===============" << endl;
-
-
-			arrange_mixture_components(line_data, blockmus_ill_fix,  blockcovars_ill_fix,  blockjs_ill_fix,  blockvs_ill_fix);
-
-			//cout << "mixture components arranged" << endl;
-
-
-
-			post_probs = probs_j_rescue(line_data,blockmus_ill_fix, blockcovars_ill_fix, blockjs_ill_fix, blockvs_ill_fix);
-
-			//cout << "mixture components arranged and post_probs calculated" << endl;
-
-			hwep = hwepvalue(calls_from_posterior_probs(post_probs));
-			cout << "HWEP after Rescue " << snpinfo2[l][0] << " " << hwep << endl;
-
-			if ((isnan(hwep) || hwep < resc_threshhold)) {
+				cout << "rescue step took " << (double)(clock() - tStart)/CLOCKS_PER_SEC  << endl;
+
+				//cout << "===============" << endl;
+				//cout << blockmus_ill_fix << endl;
+				//cout << blockjs_ill_fix << endl;
+				//cout << "===============" << endl;
+
+
+				arrange_mixture_components(line_data, blockmus_ill_fix,  blockcovars_ill_fix,  blockjs_ill_fix,  blockvs_ill_fix);
+
+				//cout << "mixture components arranged" << endl;
+
+
+
+				post_probs = probs_j_rescue(line_data,blockmus_ill_fix, blockcovars_ill_fix, blockjs_ill_fix, blockvs_ill_fix);
+
+				//cout << "mixture components arranged and post_probs calculated" << endl;
+
+				hwep = hwepvalue(calls_from_posterior_probs(post_probs));
+				cout << "HWEP after Rescue " << snpinfo2[l][0] << " " << hwep << endl;
+
+				if ((isnan(hwep) || hwep < resc_threshhold)) {
 				
-				if ( noblank == false) {
-					cout << " rescue can't fix it! Calling everything unknown" << endl;
-					post_probs = postprobstounknown(post_probs);
+					if ( noblank == false) {
+						cout << " rescue can't fix it! Calling everything unknown" << endl;
+						post_probs = postprobstounknown(post_probs);
+					}
+					else {
+						cout << "SNP would be set to unknown if -noblank was not applied" << endl;
+					}
 				}
-				else {
-					cout << "SNP would be set to unknown if -noblank was not applied" << endl;
-				}
+
 			}
 
 		}
 
 
-
-
 		//reinflate data to set the nan rows to unknown
 		//cout << "reinflating outputting post probs " << endl;
 		post_probs = probs_with_nans(post_probs,nan_samples, excluded_samples, int_data_w_nans.rows() );