Skip to content
Snippets Groups Projects
Snakefile 10.1 KiB
Newer Older
#where are we (can be configured through -d/--directory flag)
JOBDIR=os.path.abspath(os.path.curdir)
if JOBDIR[-1] != "/": # this checks if jobdir ends with slash if not it adds a slash
#data specific config file, expected to be inside JOBDIR
configfile: "tomancak_czi.json"
datasets = glob.glob('*00.h5') # searches JOBDIR for files that match this wildcard expression
#TODO: this should go into a python module in this path
fre = re.compile(r'(?P<xml_base>\w+)-(?P<file_id>\d+)-00.h5')
xml_merge_in = []
for ds in datasets:
   bn = os.path.basename(ds)
   bn_res = fre.search(bn)
   xml_base,file_id = bn_res.group('xml_base'),bn_res.group('file_id')
   xml_merge_in.append("{xbase}.job_{fid}.xml".format(xbase=xml_base, fid=int(file_id)))
   
#TODO: this should go into a python module in this path
def produce_string(_fstring, *args, **kwargs):
   contents = dict()
   for item in args:
      if type(item) == type(kwargs):
         contents.update(item)
         
   contents.update(kwargs)
   return _fstring.format(**contents)

rule done:
    input: [ ds+"_deconvolved" for ds in datasets ]

# defining xml for czi dataset
rule define_xml_czi:
    input: glob.glob('*.czi')
    output: config["common"]["first_xml_filename"] + ".xml"
Christopher Schmied's avatar
Christopher Schmied committed
    message: "Execute define_xml_czi on the following files {input}"
    log: "define_xml_czi.log"
Christopher Schmied's avatar
Christopher Schmied committed
        cmd_string = produce_string(
        	"""{fiji-prefix} {fiji-app} \
        -Dimage_file_directory={jdir} \
	-Dfirst_czi={first_czi} \
	-Dangle_1={angle_1} \
	-Dangle_2={angle_2} \
	-Dangle_3={angle_3} \
	-Dangle_4={angle_4} \
	-Dangle_5={angle_5} \
	-Dchannel_1={channel_1} \
	-Dchannel_2={channel_2} \
	-Dillumination_1={illumination_1} \
	-Drotation_around={rotation_around} \
	-Dpixel_distance_x={pixel_distance_x} \
	-Dpixel_distance_y={pixel_distance_y} \
	-Dpixel_distance_z={pixel_distance_z} \
	-Dpixel_unit={pixel_unit} \
	-Dfirst_xml_filename={first_xml_filename} \
	-- --no-splash {path_bsh}""",
	config["common"],
	config["define_xml_czi"],
	jdir=JOBDIR,
	path_bsh=config["common"]["bsh_directory"] + config["define_xml_czi"]["bsh_file"])
	cmd_string += "> {log} 2>&1"
# defining xml for tif dataset
rule define_xml_tif:
    input: glob.glob('*.tif')
    output: config["common"]["first_xml_filename"] + ".xml"
    message: "Execute define_xml_czi on the following files {input}"
    log: "define_xml_tif.log"
    run:
        cmd_string = produce_string(
        		"""{fiji-prefix} {fiji-app} \
        	-Dimage_file_directory={jdir} \
        	-Dtimepoints={timepoints} \
        	-Dacquisition_angles={acquisition_angles} \
        	-Dchannels={channels} \
        	-Dimage_file_pattern={image_file_pattern} \
        	-Dpixel_distance_x={pixel_distance_x} \
        	-Dpixel_distance_y={pixel_distance_y} \
        	-Dpixel_distance_z={pixel_distance_z} \
        	-Dpixel_unit={pixel_unit} \
        	-Dxml_filename={first_xml_filename} \
        	-Dtype_of_dataset={type_of_dataset} \
        	-Dmultiple_timepoints={multiple_timepoints} \
        	-Dmultiple_channels={multiple_channels} \
        	-Dmultiple_illumination_directions={multiple_illumination_directions} \
        	-Dmultiple_angles={multiple_angles} \
        	-Dimglib_container={imglib_container} \
		-- --no-splash {path_bsh}""",
	config["common"],
	config["define_xml_tif"],
	jdir=JOBDIR,
	path_bsh=config["common"]["bsh_directory"] + config["define_xml_tif"]["bsh_file"])
	cmd_string +="> {log} 2>&1 && touch {output}"
	shell(cmd_string)
	
# resave .czi dataset as hdf5	
rule resave_czi_hdf5:
    input: rules.define_xml_czi.output, glob.glob('*.czi')
    output: "hdf5_Stock68.h5", "hdf5_Stock68.xml" 
    message: "Execute resave_hdf5 on {input}"
#    threads: int(config["resave_hdf5"]["parallel_timepoints"]) # parallel timepoints should tell me how many timepoints to expect 
    log: "export_czi.log"
    run:
        part_string = produce_string("""{fiji-prefix} {fiji-app} \
        -Dimage_file_directory={jdir} \
        -Dfirst_xml_filename={first_xml_filename} \
        -Dhdf5_xml_filename={hdf5_xml_filename} \
        -Dresave_angle={resave_angle} \
        -Dresave_channel={resave_channel} \
        -Dresave_illumination={resave_illumination} \
        -Dresave_timepoint={resave_timepoint} \
        -Dsubsampling_factors={subsampling_factors} \
        -Dhdf5_chunk_sizes={hdf5_chunk_sizes} \
        -Dtimepoints_per_partition={timepoints_per_partition} \
        -Dsetups_per_partition={setups_per_partition} \
        -Drun_only_job_number={job_number} \
        -- --no-splash {path_bsh} > {log}-{job_number}.log 2>&1 &\n""", # the & submits everyting at once
           config["common"],
           config["define_xml_czi"],
           config["resave_hdf5"],
           jdir=JOBDIR,
           path_bsh=config["common"]["bsh_directory"] + config["resave_hdf5"]["bsh_file"])
        cmd_string = ""
        for i in range(int(config["resave_hdf5"]["parallel_timepoints"])):
        	cmd_string += part_string.format(job_number=i) 
        shell(cmd_string) 
       	       	
    input: "{xml_base}-{file_id}-00.h5"
    #input: rules.resave_hdf5.output, "{xml_base}-{file_id}-00.h5"
    output: "{xml_base}-{file_id,\d+}-00.h5_registered", #"{xml_base}.job_{file_id,\d+}.xml"
    log: "{xml_base}-{file_id}-registration.log"
    run:
        cmd_string = produce_string(
           """{fiji-prefix} {fiji-app} \
-Dparallel_timepoints={file_id_w} \
-Dimage_file_directory={jdir} \
-Dxml_filename={input_xml} \
-Dprocess_timepoint={timepoint} \
-Dprocess_channel={channel} \
-Dprocess_illumination={illuminations} \
-Dprocess_angle={angle} \
-Dprocessing_channel={proc-ch} \
-Dlabel_interest_points={label_interest_points} \
-Dtype_of_registration={type_of_registration} \
-Dtype_of_detection={type_of_detection} \
-Dsubpixel_localization={subpixel_localization} \
-Dimglib_container={imglib_container} \
-Dradius_1={radius_1} \
-Dradius_2={radius_2} \
-Dthreshold={threshold} \
-Dinitial_sigma={initial_sigma} \
-Dthreshold_gaussian={threshold_gaussian} \
-Dregistration_algorithm={algorithm} \
-Dinterest_points_channel_0={interest_points_channel_0} \
-Dinterest_points_channel_1={interest_points_channel_1} \
-Dfix_tiles={fix_tiles} \
-Dmap_back_tiles={map_back_tiles} \
-Dtransformation_model={transformation_model} \
-Dmodel_to_regularize_with={model_to_regularize_with} \
-Dlambda={lambda} \
-Dallowed_error_for_ransac={allowed_error_for_ransac} \
-Ddetection_min_max={detection_min_max} \
-- --no-splash {bsh_file}""",
           config["common"], 
           config["registration"],
           file_id_w="{wildcards.file_id}",
           jdir=JOBDIR,
           input_xml="{wildcards.xml_base}")
        cmd_string += "> {log} 2>&1 && touch {output}"
        shell(cmd_string)
        #shell("touch {output}")

rule xml_merge:
    input: [ item+"_registered" for item in datasets ] #xml_merge_in 
    output: "{xml_base}_merge.xml"
    log: "{xml_base}_merge.log"
    run:
        cmd_string = produce_string("{fiji-prefix} {fiji-app} -Dimage_file_directory={jdir} -Dmerged_xml={output} -- --no-splash {bsh_file} > {log} 2>&1",
                                    config["common"],
                                    config["xml_merge"],
                                    log="{log}",
        shell(cmd_string)

rule external_trafo:
    input: rules.xml_merge.output
    output: rules.xml_merge.output[0] + "_external_trafo"
    log: "external_trafo.log"
    run:
        cmd_string = produce_string("""{fiji-prefix} {fiji-app} \
-Dimage_file_directory={jdir} \
Christopher Schmied's avatar
Christopher Schmied committed
	-Dmerged_xml={input} \
-Dtransform_angle={angle} \
-Dtransform_channel={channel} \
-Dtransform_illumination={illumination} \
-Dtransform_timepoint={timepoint} \
-Dtransformation={transformation} \
-Dapply_transformation={apply_transformation} \
-Ddefine_mode_transform={define_mode_transform} \
-Dmatrix_transform={matrix_transform} \
-- --no-splash {bsh_file}""",
                                    config["common"], 
                                    config["external_transform"],
        cmd_string += "> {log} 2>&1 && touch {output}"
        shell(cmd_string)

rule deconvolution:
    input: rules.external_trafo.output, "{xml_base}-{file_id,\d+}-00.h5", merged_xml="{xml_base}_merge.xml"
    output: "{xml_base}-{file_id,\d+}-00.h5_deconvolved"
    log: "{xml_base}-{file_id,\d+}-00-deconvolution.log"
    run:
        cmd_string = produce_string("""{fiji-prefix} {fiji-app} \
-Dimage_file_directory={jdir} \
-Ddeco_output_file_directory={jdir} \
-Dparallel_timepoints={file_id_w} \
-Dprocess_timepoint={process_timepoint} \
-Dprocess_channel={process_channel} \
-Dprocess_illumination={process_illumination} \
-Dprocess_angle={process_angle} \
-Dminimal_x_deco={minimal_x} \
-Dminimal_y_deco={minimal_y} \
-Dminimal_z_deco={minimal_z} \
-Dmaximal_x_deco={maximal_x} \
-Dmaximal_y_deco={maximal_y} \
-Dmaximal_z_deco={maximal_z} \
-Dimglib2_container_deco={imglib2_container} \
-Dtype_of_iteration={type_of_iteration} \
-Dosem_acceleration={osem_acceleration} \
-DTikhonov_parameter={Tikhonov_parameter} \
-Dcompute={compute} \
-Dpsf_estimation={psf_estimation} \
-Ddirectory_cuda={directory_cuda} \
-Ddetections_to_extract_psf_for_channel_0={detections_to_extract_psf_for_channel_0} \
-Ddetections_to_extract_psf_for_channel_1={detections_to_extract_psf_for_channel_1} \
-Dpsf_size_x={psf_size_x} \
-Dpsf_size_y={psf_size_y} \
-Dpsf_size_z={psf_size_z} \
-Diterations={iterations} \
                                    config["common"], 
                                    config["deconvolution"],
                                    file_id_w="{wildcards.file_id}",
                                    jdir=JOBDIR,
                                    merged_xml_file="{input.merged_xml}"
        )
        cmd_string += "> {log} 2>&1 && touch {output}"
        shell(cmd_string)
  
    shell : "rm -rf *registered *tif *log *_deconvolved *.xml~* interestpoints "+str(" ".join(xml_merge_in))