Commit b7e75e53 authored by Maximilian Beckers's avatar Maximilian Beckers
Browse files

restructuring of the programm to make it callable as a python package

parent 77a4ff64
No preview for this file type
<?xml version="1.0" encoding="UTF-8"?>
<module type="PYTHON_MODULE" version="4">
<component name="NewModuleRootManager">
<content url="file://$MODULE_DIR$" />
<orderEntry type="inheritedJdk" />
<orderEntry type="sourceFolder" forTests="false" />
<orderEntry type="library" name="R Skeletons" level="application" />
<orderEntry type="library" name="R User Library" level="project" />
</component>
<component name="TestRunnerService">
<option name="projectConfiguration" value="pytest" />
<option name="PROJECT_TEST_RUNNER" value="pytest" />
</component>
</module>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="Encoding" addBOMForNewFiles="with NO BOM" />
</project>
\ No newline at end of file
<component name="libraryTable">
<library name="R User Library">
<CLASSES />
<SOURCES />
</library>
</component>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectRootManager" version="2" project-jdk-name="Python 2.7 (2)" project-jdk-type="Python SDK" />
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ProjectModuleManager">
<modules>
<module fileurl="file://$PROJECT_DIR$/.idea/FDRthresholding.iml" filepath="$PROJECT_DIR$/.idea/FDRthresholding.iml" />
</modules>
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="VcsDirectoryMappings">
<mapping directory="$PROJECT_DIR$" vcs="Git" />
</component>
</project>
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<project version="4">
<component name="ChangeListManager">
<list default="true" id="77e0d2d5-4887-4fc9-9c52-109dc2916798" name="Default Changelist" comment="">
<change beforePath="$PROJECT_DIR$/FDRcontrol.py" beforeDir="false" afterPath="$PROJECT_DIR$/FDRcontrol.py" afterDir="false" />
<change beforePath="$PROJECT_DIR$/FDRutil.py" beforeDir="false" />
<change beforePath="$PROJECT_DIR$/locscaleUtil.py" beforeDir="false" />
<change beforePath="$PROJECT_DIR$/mapUtil.py" beforeDir="false" />
</list>
<option name="EXCLUDED_CONVERTED_TO_IGNORED" value="true" />
<option name="SHOW_DIALOG" value="false" />
<option name="HIGHLIGHT_CONFLICTS" value="true" />
<option name="HIGHLIGHT_NON_ACTIVE_CHANGELIST" value="false" />
<option name="LAST_RESOLUTION" value="IGNORE" />
</component>
<component name="FileEditorManager">
<leaf SIDE_TABS_SIZE_LIMIT_KEY="300">
<file pinned="false" current-in-tab="true">
<entry file="file://$PROJECT_DIR$/FDRcontrol.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="2320">
<caret line="168" column="45" selection-start-line="168" selection-start-column="45" selection-end-line="168" selection-end-column="45" />
<folding>
<element signature="e#86#108#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/confidenceMapUtil/FDRutil.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="1230">
<caret line="82" selection-start-line="82" selection-end-line="82" />
<folding>
<element signature="e#0#18#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</file>
<file pinned="false" current-in-tab="false">
<entry file="file://$PROJECT_DIR$/confidenceMapUtil/mapUtil.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="240">
<caret line="16" column="73" selection-start-line="16" selection-start-column="73" selection-end-line="16" selection-end-column="73" />
<folding>
<element signature="e#0#21#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</file>
</leaf>
</component>
<component name="Git.Settings">
<option name="RECENT_GIT_ROOT_PATH" value="$PROJECT_DIR$" />
</component>
<component name="IdeDocumentHistory">
<option name="CHANGED_PATHS">
<list>
<option value="$PROJECT_DIR$/mapUtil.py" />
<option value="$PROJECT_DIR$/FDRutil.py" />
<option value="$PROJECT_DIR$/FDRcontrol.py" />
<option value="$PROJECT_DIR$/confidenceMapUtil/FDRutil.py" />
</list>
</option>
</component>
<component name="ProjectFrameBounds">
<option name="x" value="1536" />
<option name="y" value="72" />
<option name="width" value="1547" />
<option name="height" value="924" />
</component>
<component name="ProjectLevelVcsManager">
<ConfirmationsSetting value="1" id="Add" />
</component>
<component name="ProjectView">
<navigator proportions="" version="1">
<foldersAlwaysOnTop value="true" />
</navigator>
<panes>
<pane id="Scope" />
<pane id="ProjectPane">
<subPane>
<expand>
<path>
<item name="FDRthresholding" type="b2602c69:ProjectViewProjectNode" />
<item name="FDRthresholding" type="462c0819:PsiDirectoryNode" />
</path>
</expand>
<select />
</subPane>
</pane>
</panes>
</component>
<component name="PropertiesComponent">
<property name="last_opened_file_path" value="$PROJECT_DIR$" />
</component>
<component name="RunDashboard">
<option name="ruleStates">
<list>
<RuleState>
<option name="name" value="ConfigurationTypeDashboardGroupingRule" />
</RuleState>
<RuleState>
<option name="name" value="StatusDashboardGroupingRule" />
</RuleState>
</list>
</option>
</component>
<component name="RunManager">
<configuration name="Unnamed" type="PythonConfigurationType" factoryName="Python" nameIsGenerated="true">
<module name="FDRthresholding" />
<option name="INTERPRETER_OPTIONS" value="" />
<option name="PARENT_ENVS" value="true" />
<envs>
<env name="PYTHONUNBUFFERED" value="1" />
</envs>
<option name="SDK_HOME" value="$USER_HOME$/PycharmProjects/GUI/venv/bin/python" />
<option name="WORKING_DIRECTORY" value="" />
<option name="IS_MODULE_SDK" value="false" />
<option name="ADD_CONTENT_ROOTS" value="true" />
<option name="ADD_SOURCE_ROOTS" value="true" />
<option name="SCRIPT_NAME" value="" />
<option name="PARAMETERS" value="" />
<option name="SHOW_COMMAND_LINE" value="false" />
<option name="EMULATE_TERMINAL" value="false" />
<option name="MODULE_MODE" value="false" />
<option name="REDIRECT_INPUT" value="false" />
<option name="INPUT_FILE" value="" />
<method v="2" />
</configuration>
</component>
<component name="SvnConfiguration">
<configuration />
</component>
<component name="TaskManager">
<task active="true" id="Default" summary="Default task">
<changelist id="77e0d2d5-4887-4fc9-9c52-109dc2916798" name="Default Changelist" comment="" />
<created>1547640142735</created>
<option name="number" value="Default" />
<option name="presentableId" value="Default" />
<updated>1547640142735</updated>
</task>
<servers />
</component>
<component name="ToolWindowManager">
<frame x="1536" y="72" width="1547" height="924" extended-state="0" />
<editor active="true" />
<layout>
<window_info active="true" content_ui="combo" id="Project" order="0" visible="true" weight="0.2564784" />
<window_info id="Structure" order="1" side_tool="true" weight="0.25" />
<window_info id="Favorites" order="2" side_tool="true" />
<window_info anchor="bottom" id="Message" order="0" />
<window_info anchor="bottom" id="Find" order="1" />
<window_info anchor="bottom" id="Run" order="2" weight="0.32932693" />
<window_info anchor="bottom" id="Debug" order="3" weight="0.4" />
<window_info anchor="bottom" id="Cvs" order="4" weight="0.25" />
<window_info anchor="bottom" id="Inspection" order="5" weight="0.4" />
<window_info anchor="bottom" id="TODO" order="6" />
<window_info anchor="bottom" id="Version Control" order="7" />
<window_info anchor="bottom" id="Terminal" order="8" />
<window_info anchor="bottom" id="Event Log" order="9" side_tool="true" />
<window_info anchor="bottom" id="Python Console" order="10" />
<window_info anchor="right" id="Commander" internal_type="SLIDING" order="0" type="SLIDING" weight="0.4" />
<window_info anchor="right" id="Ant Build" order="1" weight="0.25" />
<window_info anchor="right" content_ui="combo" id="Hierarchy" order="2" weight="0.25" />
<window_info anchor="right" id="R Graphics" order="3" />
<window_info anchor="right" id="R Packages" order="4" />
</layout>
</component>
<component name="editorHistoryManager">
<entry file="file://$PROJECT_DIR$/locscaleUtil.py" />
<entry file="file://$PROJECT_DIR$/mapUtil.py" />
<entry file="file://$PROJECT_DIR$/FDRutil.py" />
<entry file="file://$PROJECT_DIR$/confidenceMapUtil/mapUtil.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="240">
<caret line="16" column="73" selection-start-line="16" selection-start-column="73" selection-end-line="16" selection-end-column="73" />
<folding>
<element signature="e#0#21#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/confidenceMapUtil/FDRutil.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="1230">
<caret line="82" selection-start-line="82" selection-end-line="82" />
<folding>
<element signature="e#0#18#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
<entry file="file://$PROJECT_DIR$/FDRcontrol.py">
<provider selected="true" editor-type-id="text-editor">
<state relative-caret-position="2320">
<caret line="168" column="45" selection-start-line="168" selection-start-column="45" selection-end-line="168" selection-end-column="45" />
<folding>
<element signature="e#86#108#0" expanded="true" />
</folding>
</state>
</provider>
</entry>
</component>
</project>
\ No newline at end of file
#Author: Maximilian Beckers, EMBL Heidelberg, Sachse Group (2017)
#import some stuff
from FDRutil import *
from mapUtil import *
from locscaleUtil import *
from confidenceMapUtil import mapUtil, locscaleUtil, FDRutil
import numpy as np
import argparse, os, sys
import subprocess
......@@ -66,6 +64,7 @@ cmdl_parser.add_argument('-stepSize', '--stepSize', metavar="stepSize_locScale",
#********************** main function ***********************
#************************************************************
def main():
start = time.time();
......@@ -80,7 +79,7 @@ def main():
print("Step Size cannot be bigger than the window_size. Job is killed ...")
return;
launch_amplitude_scaling(args);
locscaleUtil.launch_amplitude_scaling(args);
else:
#no ampltidue scaling will be done
print('************************************************');
......@@ -111,25 +110,6 @@ def main():
map = mrcfile.open(filename, mode='r+');
mapData = np.copy(map.data);
apix = args.apix;
#get boxCoordinates
if args.noiseBox is None:
boxCoord = 0;
else:
boxCoord = args.noiseBox;
#set test procdure
if args.testProc is not None:
testProc = args.testProc;
else:
testProc = 'rightSided';
#set ECDF
if args.ecdf:
ECDF = 1;
else:
ECDF = 0;
#set output filename
if args.outputFilename is not None:
......@@ -137,139 +117,71 @@ def main():
else:
splitFilename = os.path.splitext(os.path.basename(filename));
sizeMap = mapData.shape;
if args.lowPassFilter is not None:
frequencyMap = calculate_frequency_map(mapData);
providedRes = apix/float(args.lowPassFilter);
mapData = lowPassFilter(np.fft.rfftn(mapData), frequencyMap, providedRes, mapData.shape);
#handle FDR correction procedure
if args.method is not None:
method = args.method;
else:
#default is Benjamini-Yekutieli
method = 'BY';
if args.window_size is not None:
wn = args.window_size;
wn = int(wn);
if wn < 20:
print("Provided window size is quite small. Please think about potential inaccuracies of your noise estimates!");
else:
wn = max(int(0.05*sizeMap[0]),10);
#generate a circular Mask
sphere_radius = (np.min(sizeMap) // 2);
circularMaskData = makeCircularMask( np.copy(mapData), sphere_radius);
#if mask is provided, take it
if args.mask is not None:
mask = mrcfile.open(args.mask, mode='r+');
maskData = np.copy(mask.data);
else:
maskData = circularMaskData;
#do some diagnostics and check for normality of map
makeDiagnosticPlot(mapData, wn, 0, False, boxCoord);
checkNormality(mapData, wn, boxCoord);
maskData = None;
#estimate noise statistics
if args.locResMap is None: #if no local Resolution map is given,don't do any filtration
mean, var, _ = estimateNoiseFromMap(mapData, wn, boxCoord);
#if varianceMap is given, use it
if args.varianceMap is not None:
varMap = mrcfile.open(args.varianceMap, mode='r+');
var = np.copy(varMap.data);
#if varianceMap is given, use it
if args.varianceMap is not None:
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+');
mean = np.copy(meanMap.data);
if np.isscalar(mean) and np.isscalar(var):
output = "Estimated noise statistics: mean: " + repr(mean) + " and variance: " + repr(var);
else:
output = "Using user provided noise statistics";
print(output);
#if meanMap is given, use it
if args.meanMap is not None:
meanMap = mrcfile.open(args.meanMap, mode='r+');
meanMapData = np.copy(meanMap.data);
else:
meanMapData = None;
else: #do localFiltration and estimate statistics from this map
#if local resolutions are given, use them
if args.locResMap is not None:
locResMap = mrcfile.open(args.locResMap, mode='r+');
locResMapData = np.copy(locResMap.data);
mapData, mean, var, ECDF = localFiltration(mapData, locResMapData, apix, True, wn, boxCoord, args.mask, maskData, ECDF);
locFiltMap = mrcfile.new(splitFilename[0] + '_locFilt.mrc', overwrite=True)
mapData = np.float32(mapData);
locFiltMap.set_data(mapData);
locFiltMap.voxel_size = apix;
locFiltMap.close();
#calculate the qMap
qMap = calcQMap(mapData, mean, var, ECDF, wn, boxCoord, circularMaskData, method, testProc);
else:
locResMapData = None;
#if a explicit thresholding is wished, do so
if args.fdr is not None:
#run the actual analysis
confidenceMap, locFiltMap, binMap, maskedMap = FDRutil.calculateConfidenceMap(mapData, args.apix, args.noiseBox, args.testProc, args.ecdf, args.lowPassFilter, args.method, args.window_size, locResMapData, meanMapData, varMapData, args.fdr);
fdr = args.fdr;
#threshold the qMap
binMap = binarizeMap(qMap, fdr);
#apply the thresholded qMap to data
maskedMap = np.multiply(binMap, mapData);
minMapValue = np.min(maskedMap[np.nonzero(maskedMap)]);
maskedMap = np.multiply(maskedMap, maskData);
if args.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);
print(output);
binMapMRC = mrcfile.new(splitFilename[0] + '_FDR' + str(fdr) + '_binMap.mrc', overwrite=True);
if locFiltMap is not None:
locFiltMapMRC = mrcfile.new(splitFilename[0] + '_locFilt.mrc', overwrite=True);
locFiltMap = np.float32(locFiltMap);
locFiltMapMRC.set_data(locFiltMap);
locFiltMapMRC.voxel_size = args.apix;
locFiltMapMRC.close();
if binMap is not None:
binMapMRC = mrcfile.new(splitFilename[0] + '_FDR' + str(args.fdr) + '_binMap.mrc', overwrite=True);
binMap = np.float32(binMap);
binMapMRC.set_data(binMap);
binMapMRC.voxel_size = apix;
binMapMRC.voxel_size = args.apix;
binMapMRC.close();
maskedMapMRC = mrcfile.new(splitFilename[0] + '_FDR'+ str(fdr) + '_maskedMap.mrc', overwrite=True);
if maskedMap is not None:
maskedMapMRC = mrcfile.new(splitFilename[0] + '_FDR'+ str(args.fdr) + '_maskedMap.mrc', overwrite=True);
maskedMap = np.float32(maskedMap);
maskedMapMRC.set_data(maskedMap);
maskedMapMRC.voxel_size = apix;
maskedMapMRC.voxel_size = args.apix;
maskedMapMRC.close();
else:
#threshold the qMap
fdr = 0.01;
binMap = binarizeMap(qMap, fdr);
#apply the thresholded qMap to data
maskedMap = np.multiply(binMap, np.copy(mapData));
minMapValue = np.min(maskedMap[np.nonzero(maskedMap)]);
if args.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);
print(output);
#invert qMap for visualization tools
confidenceMap = np.subtract(1.0, qMap);
#apply lowpass-filtered mask to maps
confidenceMap = np.multiply(confidenceMap, circularMaskData);
#write the confidence Maps
confidenceMapMRC = mrcfile.new(splitFilename[0] + '_confidenceMap.mrc', overwrite=True);
confidenceMap = np.float32(confidenceMap);
confidenceMapMRC.set_data(confidenceMap);
confidenceMapMRC.voxel_size = apix;
confidenceMapMRC.voxel_size = args.apix;
confidenceMapMRC.close();
end = time.time();
totalRuntime = end -start;
printSummary(args, totalRuntime);
mapUtil.printSummary(args, totalRuntime);
if (__name__ == "__main__"):
main()
......
......@@ -4,9 +4,132 @@ import math
import gc
import os
import sys
import mapUtil
#Author: Maximilian Beckers, EMBL Heidelberg, Sachse Group
def calculateConfidenceMap(em_map, apix, noiseBox, testProc, ecdf, lowPassFilter_resolution, method, window_size, locResMap,
meanMap, varMap, fdr):
#*********************************************
#******* this function calc. confMaps ********
#*********************************************
# get boxCoordinates
if noiseBox is None:
boxCoord = 0;
else:
boxCoord = noiseBox;
# set test procdure
if testProc is not None:
testProc = testProc;
else:
testProc = 'rightSided';
# set ECDF
if ecdf:
ECDF = 1;
else:
ECDF = 0;
sizeMap = em_map.shape;
if lowPassFilter is not None:
frequencyMap = calculate_frequency_map(em_map);
providedRes = apix/float(lowPassFilter_resolution);
em_map = lowPassFilter(np.fft.rfftn(em_map), frequencyMap, providedRes, em_map.shape);
# handle FDR correction procedure
if method is not None:
method = method;
else:
# default is Benjamini-Yekutieli
method = 'BY';
if window_size is not None:
wn = window_size;
wn = int(wn);
if wn < 20:
print("Provided window size is quite small. Please think about potential inaccuracies of your noise estimates!");
else:
wn = max(int(0.05 * sizeMap[0]), 10);
# generate a circular Mask
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);
checkNormality(em_map, wn, boxCoord);
# estimate noise statistics
if locResMap is None: # if no local Resolution map is given,don't do any filtration
mean, var, _ = estimateNoiseFromMap(em_map, wn, boxCoord);
if varMap is not None:
var = varMap;
if meanMap is not None:
mean = meanMap;
if np.isscalar(mean) and np.isscalar(var):
output = "Estimated noise statistics: mean: " + repr(mean) + " and variance: " + repr(var);
else:
output = "Using user provided noise statistics";
print(output);
locFiltMap = None;
else: # do localFiltration and estimate statistics from this map
em_map, mean, var, ECDF = mapUtil.localFiltration(em_map, locResMap, apix, True, wn, boxCoord, ECDF);
locFiltMap = em_map;
# calculate the qMap
qMap = calcQMap(em_map, mean, var, ECDF, wn, boxCoord, circularMaskData, method, testProc);
# if a explicit thresholding is wished, do so
if fdr is not None:
fdr = fdr;
# threshold the qMap
binMap = binarizeMap(qMap, fdr);
# apply the thresholded qMap to data
maskedMap = np.multiply(binMap, em_map);
minMapValue = np.min(maskedMap[np.nonzero(maskedMap)]);
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);
print(output);
else:
# threshold the qMap
fdr = 0.01;
binMap = binarizeMap(qMap, fdr);
# apply the thresholded qMap to data
maskedMap = np.multiply(binMap, np.copy(em_map));
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);
print(output);
binMap = None;
maskedMap = None;