cryoSPARC调用motioncorr2程序运行报错

cryoSPARC调用motioncorr2程序运行报错

报错信息:

Traceback (most recent call last):
  File "cryosparc_master/cryosparc_compute/run.py", line 116, in cryosparc_master.cryosparc_compute.run.main
  File "/root/cryosparcuser/cryosparc/cryosparc_worker/cryosparc_compute/jobs/motioncorrection/run_motioncor2.py", line 398, in run_motioncor2_wrapper
    global_shifts = parse_global_alignment(full_frame_log_abs)
  File "/root/cryosparcuser/cryosparc/cryosparc_worker/cryosparc_compute/jobs/motioncorrection/run_motioncor2.py", line 22, in parse_global_alignment
    with open(full_frame_alignment_log) as logfile:
FileNotFoundError: [Errno 2] No such file or directory: '/data6/x/CS-x/J72/motioncor2_logs/00-Patch-Full.log'

原因

worker安装路径下/cryosparc_compute/jobs/motioncorrection/run_motioncor2.py 代码与motioncorr2_1.6.3版本调用命令不同

解决方案:

替换run_motioncor2.py代码内容

## ---------------------------------------------------------------------------
##    Copyright (c) 2024 Structura Biotechnology Inc. All rights reserved.
##         Do not reproduce or redistribute, in whole or in part.
##      Use of this code is permitted only under licence from Structura.
##                   Contact us at info@structura.bio.
## ---------------------------------------------------------------------------junquan


from builtins import str
from builtins import range
import os
import time
import numpy as n
from .. import runcommon as rc
import matplotlib.pyplot as plt
import matplotlib
from ...blobio import mrc
from ...blobio import tiff

def parse_global_alignment(full_frame_alignment_log):
    shifts = []
    with open(full_frame_alignment_log) as logfile:
        for line in logfile:
            stripped_line = line.strip()
            if '#' not in stripped_line and len(stripped_line) > 0:
                values = stripped_line.split()
                shifts.append([float(values[1]), float(values[2])])

    global_shifts = n.array([shifts])
    return global_shifts

def parse_patch_alignment(patch_patch_alignment_log):
    shifts = []
    patches = []
    with open(patch_patch_alignment_log) as logfile:
        for line in logfile:
            stripped_line = line.strip()
            if '#' in stripped_line or len(stripped_line) == 0: continue

            values = stripped_line.split()
            if int(values[0]) == 1:
                # e.g. 1   740.91   766.51     0.00     0.00     0.00     0.00
                # this means we've reached a new patch
                shifts.append([])
                # record the coordinates of the patch
                patches.append([float(values[1]), float(values[2])])

            shifts[-1].append([float(values[-2]), float(values[-1])])

    patch_shifts = n.array(shifts)
    patch_coords = n.array(patches)

    return patch_shifts, patch_coords

def find_extension(abs_path):
    ext = str(os.path.splitext(abs_path)[1])
    if ext in ['.mrc', '.mrcs', '.stk']:
        return 'mrc'
    elif ext in ['.tiff', '.tif']:
        return 'tif'
    else:
        assert False, "Cant read extension %s for input file path" % (ext)

def make_command(
        exec_path,
        abs_input_path,
        abs_output_path,
        abs_logfile_path,
        gpu_resources,
        gpu_mem_usage,
        patch_x,
        patch_y,
        patch_overlap,
        abs_gain_ref_path,
        gain_rotate_num,
        gain_flip_num,
        abs_dark_ref_path,
        abs_defect_correct_path,
        alignment_config,
        iter,
        tol,
        throw,
        trunc,
        align_to_frame_num,
        dose_weight,
        psize_A,
        accel_kv,
        frame_dose,
        tilt_start,
        tilt_step,
        selected_sum,
        sum_range_a,
        sum_range_b,
        major_scale, 
        minor_scale,
        distortion_angle,
        fourier_cropping,
        low_snr,
        b_factor_global,
        b_factor_local,
        group_num,
        frame_motion_blurring,
        phase_only,
        extra
        ):

        path_to_motioncor2_executable = exec_path

        command = [
            path_to_motioncor2_executable,
        ]

        #find extension of input movie
        input_type = find_extension(abs_input_path)

        if input_type == "tif": command.extend(['-InTiff', abs_input_path])
        elif input_type == "mrc": command.extend(['-InMrc', abs_input_path])

        if abs_output_path is not None: command.extend(['-OutMrc', abs_output_path])

        if abs_gain_ref_path is not None: command.extend(['-Gain', abs_gain_ref_path])

        if gain_rotate_num != 0: command.extend(['-RotGain', gain_rotate_num])
        
        if gain_flip_num != 0: command.extend(['-FlipGain', gain_flip_num])

        if abs_dark_ref_path is not None: 
            dark_ext = str(os.path.splitext(abs_dark_ref_path)[1])
            assert dark_ext in ['.mrc', '.mrcs', '.stk'], "Dark Reference file must be a .mrc file"
            command.extend(['-Dark', abs_dark_ref_path])
        
        if abs_defect_correct_path is not None:
            defect_ext = str(os.path.splitext(abs_defect_correct_path)[1])
            assert defect_ext in ['.txt'], "Defect file must be a .txt file"
            command.extend(['-DefectFile', abs_defect_correct_path])

        if patch_x is not None and patch_y is not None: 
            if patch_overlap is None:
                command.extend(['-Patch', patch_x, patch_y])
            else:
                command.extend(['-Patch', patch_x, patch_y, patch_overlap])
        
        if alignment_config: command.extend(['-Iter', iter, '-Tol', tol])
        
        if throw is not None: command.extend(['-Throw', throw])
        
        if trunc is not None: command.extend(['-Trunc', trunc])
        
        if align_to_frame_num is not None: command.extend(['-FmRef', align_to_frame_num])

        if dose_weight: command.extend(['-Kv', accel_kv, '-PixSize', psize_A, '-FmDose', frame_dose])
        
        if (tilt_start is not None and tilt_step is not None): command.extend(['-Tilt', tilt_start, tilt_step])
        
        if selected_sum: command.extend(['-SumRange', sum_range_a, sum_range_b])

        if major_scale is not None and minor_scale is not None and distortion_angle is not None: command.extend(['-Mag', major_scale, minor_scale, distortion_angle])

        if fourier_cropping is not None: command.extend(['-FtBin', fourier_cropping])

        if low_snr: 
            if (b_factor_global is not None) & (b_factor_local is not None): command.extend(['-Bft', b_factor_global, b_factor_local])
            if group_num is not None: command.extend (['-Group', group_num])
            
        if frame_motion_blurring: command.extend(['-InFmMotion 1'])
        
        if phase_only: command.extend(['-PhaseOnly'])

        if gpu_resources != 0:
            command.extend(['-Gpu', gpu_resources])

        command.extend(['-GpuMemUsage', gpu_mem_usage])

        command.extend(['-LogDir', abs_logfile_path])  #---------------------------------junquan-----------------------

        if extra is not None: command.extend([extra])

        command = [str(c) for c in command] # convert everything to string

        return(command)

def run_motioncor2_wrapper(job):
    """ function that runs the actual job. """

    params = rc.com.get_merged_params(job)

    proj_dir_abs = rc.get_project_dir_abs()
    job_dir_rel = job['job_dir']
    job_dir_abs = os.path.join(proj_dir_abs, job_dir_rel)
    #puid = job['project_uid']
    #juid = job['uid']

    res_alloc = job['resources_allocated']
    hostname = res_alloc['hostname']
    print("Running job on hostname %s", hostname)
    print("Allocated Resources : ", res_alloc)
    cuda_devs = res_alloc['slots']['GPU']
    cuda_dev = ("".join(str(dev) + " " for dev in cuda_devs)).strip() #converts [0,1,2] into '0 1 2'

    #terms
    rc.log("CryoSPARC does not distribute MotionCor2 binaries. Please ensure you have your own copy of MotionCor2 installed under the terms of the MotionCor2 Non-Commercial Software License Agreement available at: https://docs.google.com/forms/d/e/1FAIpQLSfAQm5MA81qTx90W9JL6ClzSrM77tytsvyyHh1ZZWrFByhmfQ/formResponse. For-profit users must contact David Agard for licensing information prior to download. Structura Biotechnology Inc. makes no warranty regarding MotionCor2.")
    
    movies_dset = rc.load_input_group(input_group_name='movies')
    movies_dset.add_fields(rc.com.get_dataset_columns('exposure.micrograph_blob', 'micrograph_blob_non_dw')) 
    movies_dset.add_fields(rc.com.get_dataset_columns('exposure.motion', 'rigid_motion'))
    movies_dset.add_fields(rc.com.get_dataset_columns('exposure.motion', 'local_motion'))
    #if the user selected the dose weighting option (which applies to all micrographs outputted)
    do_dose_weighting = params.get('dose_weighting', False)
    if do_dose_weighting:
        movies_dset.add_fields(rc.com.get_dataset_columns('exposure.micrograph_blob', 'micrograph_blob')) 

    os.makedirs(os.path.join(proj_dir_abs, job_dir_rel, 'motioncorrected'))
    rc.log("Created dirctory %s" %(os.path.join(proj_dir_abs, job_dir_rel, 'motioncorrected')))
    os.makedirs(os.path.join(proj_dir_abs, job_dir_rel, 'motioncor2_logs'))
    rc.log("Created directory %s"% (os.path.join(proj_dir_abs, job_dir_rel, 'motioncor2_logs')))

    #TODO: Find a better way to store the path to motioncor2 executable (maybe in config.sh?)
    motioncor2_exec_path = params.get('exec_path', None)
    assert motioncor2_exec_path is not None, "Must provide path to motioncor2 executable."

    movies = movies_dset.rows()
    keep_idxs = []

    #TODO: BATCH PROCESSING MODE
    # SINGLE PROCESSING MODE
    num_mov_processed = 0
    checkpoint_freq = 5
    num_mov_total = len(movies[:params.get('random_num')])
    fulltic = time.time()

    rc.log("Starting to process exposures...")

    for movie in movies[:params.get('random_num')]:
        
        if num_mov_processed % checkpoint_freq == 0: 
            rc.log_checkpoint()
            rc.log("--------------------------------------------------------------")
            rc.log("Processed %d of %d movies in %.2fs " % (num_mov_processed, num_mov_total, time.time()-fulltic))
        keep_idxs.append(movie.idx)
        
        movietic = time.time()
        tic = time.time()

        #find the full path to the input movie
        mov_path_rel = movie['movie_blob/path']
        mov_path_abs = os.path.join(proj_dir_abs, mov_path_rel)
        rc.log('Raw movie filepath located at: %s - creating MotionCor2 command string...' % (mov_path_rel))

        input_type = find_extension(mov_path_abs)
        if input_type == "mrc":
            mov_raw = mrc.read_mrc(mov_path_abs)
        elif input_type == "tif":
            mov_raw = tiff.read_tiff(mov_path_abs)

        #full path to output micrograph
        basename = rc.com.get_movie_basename(movie['movie_blob/path'])
        output_path_rel = os.path.join(job_dir_rel, 'motioncorrected', basename + '_motioncor2_aligned.mrc' )
        output_path_abs = os.path.join(proj_dir_abs, output_path_rel)
        #full path to dose weighted micrograph
        if do_dose_weighting:
            dw_output_path_rel = os.path.join(job_dir_rel, 'motioncorrected', basename + '_motioncor2_aligned_DW.mrc' )
            dw_output_path_abs = os.path.join(proj_dir_abs, dw_output_path_rel)

        #motioncor2 log prefix file
        motioncor2_log_file_prefix = os.path.join(proj_dir_abs, job_dir_rel, 'motioncor2_logs') #, str(movie.idx)

        psize_raw = movie['movie_blob/psize_A']
        psize_final = psize_raw
        fourier_cropping = params.get('fourier_cropping', None)
        if fourier_cropping is not None:
            psize_final = psize_raw * fourier_cropping

        #obtain movie's microscope params (its a requirement that these are provided)
        if not params['override_dose_weighting']:
            accel_kv = movie['mscope_params/accel_kv']
            movie_shape = tuple(movie['movie_blob/shape'])
            frame_dose = movie['mscope_params/total_dose_e_per_A2'] / movie_shape[0] #divide frame dose by total number of frames
        else:
            accel_kv = params['accel_kv']
            frame_dose = params['fm_dose_e_per_A2']

        #get gain reference details
        gain_path_abs = None; flip_x = False; flip_y = False; rotate_num = 0; flip_num = 0
        if not movie['movie_blob/is_gain_corrected']:
            gain_path_rel = movie['gain_ref_blob/path']
            gain_path_abs = os.path.join(proj_dir_abs, gain_path_rel)
            flip_x = movie['gain_ref_blob/flip_x']
            flip_y = movie['gain_ref_blob/flip_y']

            if flip_x:
                flip_num = 2
            elif flip_y:
                flip_num = 1

            gain_rotate_num = movie['gain_ref_blob/rotate_num']

            if (gain_rotate_num >= 1) & (gain_rotate_num <= 91):
                rotate_num = 1
            elif (gain_rotate_num > 91) & (gain_rotate_num <= 181):
                rotate_num = 2
            elif (gain_rotate_num > 181) & (gain_rotate_num <= 271):
                rotate_num = 3
        
        command = make_command( exec_path = motioncor2_exec_path,
                                abs_input_path = mov_path_abs,
                                abs_output_path = output_path_abs,
                                abs_logfile_path = motioncor2_log_file_prefix, 
                                gpu_resources = cuda_dev,
                                gpu_mem_usage = params.get('gpu_mem_usage', 0.5),
                                patch_x = params.get('patch_x', 5.0),
                                patch_y = params.get('patch_y', 5.0),
                                patch_overlap = params.get('patch_overlap', None),
                                abs_gain_ref_path = gain_path_abs,
                                gain_rotate_num = rotate_num,
                                gain_flip_num = flip_num,
                                abs_dark_ref_path = params.get('dark_ref_path', None),
                                abs_defect_correct_path = params.get('defect_path', None),
                                alignment_config = params.get('alignment_config', False),
                                iter = params.get('max_iterations', None),
                                tol = params.get('alignment_tolerance', None),
                                throw = params.get('throw', None),
                                trunc = params.get('trunc', None),
                                align_to_frame_num = params.get('frame_index', None),
                                dose_weight = params.get('dose_weighting', False),
                                psize_A = psize_raw,
                                accel_kv = accel_kv,
                                frame_dose = frame_dose,
                                tilt_start = params.get('tilt_start', None),
                                tilt_step = params.get('tilt_step', None),
                                selected_sum = params.get('selected_sum', False),
                                sum_range_a = params.get('sum_range_a', None),
                                sum_range_b = params.get('sum_range_b', None),
                                major_scale = params.get('major_scale', None),
                                minor_scale = params.get('minor_scale', None),
                                distortion_angle = params.get('distortion_angle', None),
                                fourier_cropping = fourier_cropping,
                                low_snr = params.get('low_snr', False),
                                b_factor_global = params.get('b_factor_global', None),
                                b_factor_local = params.get('b_factor_local', None),
                                group_num = params.get('group_num', 1),
                                frame_motion_blurring = params.get('frame_motion_blurring', False),
                                phase_only = params.get('phase_only', False),
                                extra = params.get('extra', None),
        )

        rc.log('Finished creating MotionCor2 command string in %.2fs' % (time.time() - tic))
        
        rc.log('Starting MotionCor2 process...')
        tic = time.time()
        try:
            rc.log('Running MotionCor2 command:', (' ').join(command))
            joblog = job_dir_abs + '/job.log'
            process = rc.start_external_process(command, log_file_path=joblog)
            rc.log('Running process %i' % process.pid)
            process.communicate()
            rc.remove_subprocess(process.pid)
            assert os.path.exists(output_path_abs), 'motioncor2 failed to produce output file %s' % output_path_abs
        except Exception as e:
            rc.log("ERROR", e)

        rc.log('Finished MotionCor2 process in %.2fs' % (time.time() - tic))

        tic = time.time()

        with open(output_path_abs) as mrc_file:
            output_mic_header = mrc.read_mrc_header(mrc_file, debug_fname=output_path_abs)

        movie['micrograph_blob_non_dw/path'] = output_path_rel
        movie['micrograph_blob_non_dw/idx'] = 0
        movie['micrograph_blob_non_dw/shape'] = [output_mic_header['ny'], output_mic_header['nx']]
        movie['micrograph_blob_non_dw/psize_A'] = psize_final
        movie['micrograph_blob_non_dw/format'] = 'MRC/2'
        movie['micrograph_blob_non_dw/is_background_subtracted'] = False

        if do_dose_weighting:
            with open(dw_output_path_abs) as mrc_file:
                dw_output_mic_header = mrc.read_mrc_header(mrc_file, debug_fname=dw_output_path_abs)

            movie['micrograph_blob/path'] = dw_output_path_rel
            movie['micrograph_blob/idx'] = 0
            movie['micrograph_blob/shape'] = [dw_output_mic_header['ny'], dw_output_mic_header['nx']]
            movie['micrograph_blob/psize_A'] = psize_final
            movie['micrograph_blob/format'] = 'MRC/2'
            movie['micrograph_blob/is_background_subtracted'] = False

        frame_start = params.get('throw', 0)
        frame_end = params.get('trunc', None)
        zero_shift_frame = params.get('align_to_frame_num', 0)

        movie['rigid_motion/type'] = 'rigid'
        movie['rigid_motion/idx'] = 0
        movie['rigid_motion/frame_start'] = 0 if frame_start is None else frame_start
        movie['rigid_motion/frame_end'] = len(mov_raw) if frame_end is None else frame_end
        movie['rigid_motion/zero_shift_frame'] = 0 if zero_shift_frame is None else zero_shift_frame  # BEFORE frame start is applied
        rigid_motion_traj_rel = os.path.join(job_dir_rel, 'motioncorrected', basename + '_rigid_motion_traj.npy' )
        rigid_motion_traj_abs = os.path.join(proj_dir_abs, rigid_motion_traj_rel)
        full_frame_log_abs = motioncor2_log_file_prefix+'/'+basename+'-Patch-Full.log'
        global_shifts = parse_global_alignment(full_frame_log_abs)
        rc.com.write_array_to_numpy_file(global_shifts, rigid_motion_traj_abs)
        movie['rigid_motion/path'] = rigid_motion_traj_rel

        movie['local_motion/type'] = 'patch'
        movie['local_motion/idx'] = 0
        movie['local_motion/frame_start'] = 0 if frame_start is None else frame_start
        movie['local_motion/frame_end'] = len(mov_raw) if frame_end is None else frame_end
        movie['local_motion/zero_shift_frame'] = 0 if zero_shift_frame is None else zero_shift_frame  # BEFORE frame start is applied
        local_patch_traj_rel = os.path.join(job_dir_rel, 'motioncorrected', basename + '_local_patch_traj.npy' )
        local_patch_traj_abs = os.path.join(proj_dir_abs, local_patch_traj_rel)
        patch_patch_log_abs = motioncor2_log_file_prefix+'/'+basename+'-Patch-Patch.log'
        patch_shifts, patch_coords = parse_patch_alignment(patch_patch_log_abs)
        rc.com.write_array_to_numpy_file(patch_shifts, local_patch_traj_abs)
        movie['local_motion/path'] = local_patch_traj_rel

        if num_mov_processed < params['num_plots']:

            fit_cumulative_patch_shifts = []
            for patch in range(0,len(patch_coords)):
                fit_cumulative_patch_shifts.append([])
                for shift in patch_shifts[patch]:
                    fit_cumulative_patch_shifts[-1].append([patch_coords[patch][0]+shift[0], patch_coords[patch][1]+shift[1]])
            fit_cumulative_patch_shifts = n.array(fit_cumulative_patch_shifts)
            
            matplotlib.rcParams.update({'font.size': 8})
            motioncorrected_figures = plt.figure(figsize=(8,4))
            plt.subplot(121)
            plt.title('Local Motion Correction (Patch 0)')
            plt.plot(fit_cumulative_patch_shifts[0][:,0],fit_cumulative_patch_shifts[0][:,1],'-k.', label='Patch Motion Alignment')
            plt.axis('equal')
            plt.grid(True)
            plt.subplot(122)
            plt.title('Global Motion Correction')
            plt.plot(global_shifts[0][:,0],global_shifts[0][:,1],'-k.', label='Global Motion Alignment')
            plt.axis('equal')
            plt.grid(True)
            plt.tight_layout()
            rc.log_plot(motioncorrected_figures, 'Local and Global Motion for %s' % basename)

            if num_mov_processed < 1:
                rc.set_tile_image('img%d'%num_mov_processed, motioncorrected_figures, 1, 2)
            if num_mov_processed < 1:
                rc.set_output_group_image('micrographs', motioncorrected_figures)

            rc.log('Finished plotting in %.2fs' % (time.time() - tic))

        rc.log('================ Completed motion correction of exposure in %.2fs' % (time.time() - movietic))
        num_mov_processed += 1

    outpath_rel = os.path.join(job_dir_rel, 'movies_motioncor2_aligned.cs')
    movies_dset = movies_dset.take(keep_idxs)
    movies_dset.save(os.path.join(proj_dir_abs, outpath_rel))
    
    rc.output('micrographs', 'micrograph_blob_non_dw', outpath_rel, 0, len(movies_dset))
    if do_dose_weighting:
        rc.output('micrographs', 'micrograph_blob', outpath_rel, 0, len(movies_dset))
    rc.output('micrographs', 'rigid_motion', outpath_rel, 0, len(movies_dset))
    rc.output('micrographs', 'local_motion', outpath_rel, 0, len(movies_dset))
    
    rc.log("Finished processing all exposures!")


评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值