opticall avatar opticall committed 397426f

changed outlier threshold on blocking to 5 sds, added option to change hwep threshold for rescue

Comments (0)

Files changed (1)

opticall/opticall.cpp

 #include <ctime>
 #include <string>
 #include <limits>
+#include <sstream>
 
 #include <time.h>
 
 	
 	for (int i = 0; i < distances.rows(); i++)
 	{
-		if ( distances(i,0) > mean_dist + 4.0*sd_dist  )
+		if ( distances(i,0) > mean_dist + 5.0*sd_dist  )
 		{
 			exc_indices.push_back(i);
 			cout << "excluding " << inputdata.row(i) << endl;
 	//char * output_prob_ext = ".probs";
 	//char * output_call_ext = ".calls";
 
-	string version = "0.1.5";
+	string version = "0.2.0";
 	cout << "opticall version " << version << endl;
 	cout << "thank you for choosing opticall for your genotyping needs" << endl;
 
 	string output_call_ext = ".calls";
 
 	int sample_block_size = 50000;
+	
+	double resc_threshhold = 5*1e-15;
 
 
 	if (argc < 5) { // Check the value of argc. If not enough parameters have been passed, inform user and exit.
 		} else if (string(argv[i]) == "-inblock") {
 			fname_block = string(argv[i + 1]);
 			i++;
+		} else if (string(argv[i]) == "-hwep") {
+			string hwep = string(argv[i + 1]);
+			istringstream os(hwep);
+			os >> resc_threshhold;
+			cout << "rescue threshold is: " << resc_threshhold << endl;
+			i++;
 		} else {
 			cout << "Not enough or invalid arguments, please try again.\n";
 			cout << i << endl;
 		/****
 		 * STEP 3 - remove sampled contour if HWE-p test failed
 		 */
-
-		if (  isnan(hwep) || hwep < 5*1e-15) {
+		
+		//5*1e-15
+
+		if (  isnan(hwep) || hwep < resc_threshhold) {
 			cout << "Out of HWE, running Rescue" << endl;
 
 			//ok so how will I do this?
 			hwep = hwepvalue(calls_from_posterior_probs(post_probs));
 			cout << "HWEP after Rescue" << snpinfo2[l][0] << " " << hwep << endl;
 
-			if (isnan(hwep) || hwep < 5*1e-15) {
+			if (isnan(hwep) || hwep < resc_threshhold) {
 				cout << " rescue can't fix it! Calling everything unknown" << endl;
 				post_probs = postprobstounknown(post_probs); 
 			}
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.