diagnostic plot generation now callable as function

FDRcontrol.py

@@ -27,7 +27,7 @@ cmdl_parser.add_argument('-halfmap2', '--halfmap2', metavar="halfmap2.mrc",  typ
                          help='Input filename halfmap 2');
 cmdl_parser.add_argument('-mask', '--mask',metavar="mask",  type=str,
                          help='Input filename mask');
-cmdl_parser.add_argument('-p', '--apix', metavar="apix",  type=float, required=True,
+cmdl_parser.add_argument('-p', '--apix', metavar="apix",  type=float, required=False,
                          help='pixel Size of input map');
 cmdl_parser.add_argument('-fdr', '--fdr', metavar="fdr",  type=float, required=False,
                          help='False Discovery Rate');
@@ -94,10 +94,11 @@ def main():
 			else:
 				#load the maps
 				filename = args.em_map;
-				map1 = mrcfile.open(args.em_map, mode='r+');
+				map1 = mrcfile.open(args.em_map, mode='r');
+				args.apix = map1.voxel_size.x;
 				halfMapData1 = np.copy(map1.data);
-
-				map2 = mrcfile.open(args.halfmap2, mode='r+');
+				
+				map2 = mrcfile.open(args.halfmap2, mode='r');
 				halfMapData2 = np.copy(map2.data);
 
 				mapData = (halfMapData1 + halfMapData2)*0.5;
@@ -107,10 +108,11 @@ def main():
 		else:
 			#load single map
 			filename = args.em_map;
-			map = mrcfile.open(filename, mode='r+');
+			map = mrcfile.open(filename, mode='r');
+			args.apix = float(map.voxel_size.x);
 			mapData = np.copy(map.data);
-
-
+			
+			
 		#set output filename
 		if args.outputFilename is not None:
 			splitFilename = os.path.splitext(os.path.basename(args.outputFilename));
@@ -120,28 +122,28 @@ def main():
 
 		#if mask is provided, take it
 		if args.mask is not None:	
-			mask = mrcfile.open(args.mask, mode='r+');
+			mask = mrcfile.open(args.mask, mode='r');
 			maskData = np.copy(mask.data);
 		else:
 			maskData = None;
 
 		#if varianceMap is given, use it
 		if args.varianceMap is not None:
-			varMap = mrcfile.open(args.varianceMap, mode='r+');
+			varMap = mrcfile.open(args.varianceMap, mode='r');
 			varMapData = np.copy(varMap.data);
 		else:
 			varMapData = None;
 			
 		#if meanMap is given, use it
 		if args.meanMap is not None:
-			meanMap = mrcfile.open(args.meanMap, mode='r+');
+			meanMap = mrcfile.open(args.meanMap, mode='r');
 			meanMapData = np.copy(meanMap.data);
 		else:
 			meanMapData = None;
 
 		#if local resolutions are given, use them
 		if args.locResMap is not None:
-			locResMap = mrcfile.open(args.locResMap, mode='r+');
+			locResMap = mrcfile.open(args.locResMap, mode='r');
 			locResMapData = np.copy(locResMap.data);
 		else:
 			locResMapData = None;

confidenceMapUtil/FDRutil.py

@@ -16,6 +16,10 @@ def calculateConfidenceMap(em_map, apix, noiseBox, testProc, ecdf, lowPassFilter
 	#******* this function calc. confMaps ********
 	#*********************************************
 
+	#print pixelSize and give feedback to the user
+	output = "Pixel size was read as " + "%.2f" %apix + " Angstrom. If this is incorrect, please specify with -p pixelSize";
+	print(output);
+
 	# get boxCoordinates
 	if noiseBox is None:
 		boxCoord = 0;
@@ -60,8 +64,11 @@ def calculateConfidenceMap(em_map, apix, noiseBox, testProc, ecdf, lowPassFilter
 	sphere_radius = (np.min(sizeMap) // 2);
 	circularMaskData = mapUtil.makeCircularMask(np.copy(em_map), sphere_radius);
 
-	# do some diagnostics and check for normality of map
-	mapUtil.makeDiagnosticPlot(em_map, wn, 0, False, boxCoord);
+	# plot locations of noise estimation
+	pp = mapUtil.makeDiagnosticPlot(em_map, wn, 0, False, boxCoord);
+	pp.savefig("diag_image.pdf");
+	pp.close();
+
 	checkNormality(em_map, wn, boxCoord);
 
 	# estimate noise statistics
@@ -105,7 +112,7 @@ def calculateConfidenceMap(em_map, apix, noiseBox, testProc, ecdf, lowPassFilter
 		maskedMap = np.multiply(maskedMap, circularMaskData);
 
 		if locResMap is None:  # if no local Resolution map is give, then give the correspoding threshold, not usefule with local filtration
-			output = "Calculated map threshold: " + repr(minMapValue) + " at a FDR of " + repr(fdr);
+			output = "Calculated map threshold: " + repr(minMapValue) + " at a FDR of " + repr(fdr*100) + "%.";
 			print(output);
 	else:
 		# threshold the qMap
@@ -117,7 +124,7 @@ def calculateConfidenceMap(em_map, apix, noiseBox, testProc, ecdf, lowPassFilter
 		minMapValue = np.min(maskedMap[np.nonzero(maskedMap)]);
 
 		if locResMap is None:  # if no local Resolution map is give, then give the correspoding threshold, not usefule with local filtration
-			output = "Calculated map threshold: " + repr(minMapValue) + " at a FDR of " + repr(fdr);
+			output = "Calculated map threshold: " + repr(minMapValue) + " at a FDR of " + repr(fdr*100) + "%.";
 			print(output);
 
 		binMap = None;
@@ -721,7 +728,7 @@ def printSummary(args, time):
 			print(output);
 
 	#print pixel size
-	output = "Pixel size: " + repr(args.apix);
+	output = "Pixel size: " + "%.2f" %args.apix;
 	print(output);
 
 	#print method used for FDR-control

confidenceMapUtil/locscaleUtil.py

@@ -85,10 +85,11 @@ def prepare_mask_and_maps_for_scaling(args):
         else:
             #load the maps
 	    filename = args.em_map;
-	    map1 = mrcfile.open(args.em_map, mode='r+');
+	    map1 = mrcfile.open(args.em_map, mode='r');
+            apix = float(map1.voxel_size.x);
             halfMapData1 = np.copy(map1.data);
 
-            map2 = mrcfile.open(args.halfmap2, mode='r+');
+            map2 = mrcfile.open(args.halfmap2, mode='r');
             halfMapData2 = np.copy(map2.data);
 
             emmap = (halfMapData1 + halfMapData2)*0.5;
@@ -98,9 +99,15 @@ def prepare_mask_and_maps_for_scaling(args):
     else:
             #load single map
             filename = args.em_map;
-            map = mrcfile.open(filename, mode='r+');
+            map = mrcfile.open(filename, mode='r');
+            apix = float(map.voxel_size.x);
             emmap = np.copy(map.data);
 
+    if args.apix is None:
+        output = "Pixel size was read from file as " + "%.2f" %apix + " Angstrom. If this is incorrect, run the program with the flag -p pixelSize";
+        print(output);
+        args.apix = apix;       
+
     modmap = np.copy(mrcfile.open(args.model_map).data);
     
     if args.mask is None:
@@ -539,7 +546,7 @@ def write_out_final_volume_window_back_if_required(args, wn, window_bleed_and_pa
 def launch_amplitude_scaling(args):
 
     startTime = time.time();
-    emmap, modmap, mask, wn, wn_locscale, window_bleed_and_pad, method, locFilt, locResMap, boxCoord = prepare_mask_and_maps_for_scaling(args); 
+    emmap, modmap, mask, wn, wn_locscale, window_bleed_and_pad, method, locFilt, locResMap, boxCoord= prepare_mask_and_maps_for_scaling(args); 
     meanNoise, varNoise, sample = estimateNoiseFromMap(emmap, wn, boxCoord);
 
     #set output filenames
@@ -585,7 +592,10 @@ def launch_amplitude_scaling(args):
 
         endTime = time.time()
         runTime = endTime - startTime
-        makeDiagnosticPlot(emmap, wn, wn_locscale, True, boxCoord);
+       
+        pp = makeDiagnosticPlot(emmap, wn, wn_locscale, True, boxCoord);
+        pp.savefig("diag_image.pdf");
+        pp.close();
         printSummary(args, runTime)
 
     elif args.mpi:
@@ -605,6 +615,8 @@ def launch_amplitude_scaling(args):
 
             endTime = time.time()
             runTime = endTime - startTime;
-            makeDiagnosticPlot(emmap, wn, wn_locscale, True, boxCoord);
+            pp = makeDiagnosticPlot(emmap, wn, wn_locscale, True, boxCoord);
+            pp.savefig("diag_image.pdf");
+            pp.close();
             printSummary(args, runTime);
 

confidenceMapUtil/mapUtil.py

@@ -70,8 +70,13 @@ def localFiltration(map, locResMap, apix, localVariance, windowSize, boxCoord, E
 	#printProgressBar(counter, numRes, prefix = 'Progress:', suffix = 'Complete', bar_length = 50)
 	print("Start local filtering. This might take a few minutes ...");
 
+	counterRes = 0;
 	for tmpRes in locResArray:   
-		#print(tmpRes);
+		counterRes = counterRes + 1;
+		progress = counterRes/float(numRes);
+		if counterRes%(int(numRes/20.0)) == 0:
+			output = "%.1f" %(progress*100) + "% finished ..." ;
+			print(output);
 		
 		#get indices of voxels with the current resolution	
 		indices = np.where(locResMapData == tmpRes);
@@ -179,7 +184,6 @@ def makeDiagnosticPlot(map, windowSize, padded, singleBox, boxCoord):
 	plt.gray(); #make grayscale images
 	plt.rc('xtick', labelsize=8);    # fontsize of the tick labels
 	plt.rc('ytick', labelsize=8);    # fontsize of the tick labels
-	pp = PdfPages('diag_image.pdf');
 	gs = gridspec.GridSpec(1, 3);
 
 	#add image of y-z slice
@@ -202,11 +206,8 @@ def makeDiagnosticPlot(map, windowSize, padded, singleBox, boxCoord):
 	ax3.set_xlabel('Y');
 	ax3.set_ylabel('X');
 	ax3.imshow(sliceMapXY);
-
-	pp.savefig();
-	pp.close();
-	plt.close();
-
+		
+	return plt;
 #------------------------------------------------------------------------------------------------
 def printProgressBar (iteration, total, prefix = '', suffix = '', decimals = 1, bar_length = 100):