OneScience-Group/flex_ddG_tutorial

1#!/usr/bin/python2 3from __future__ import print_function4 5import socket6import sys7import os8import subprocess9 10script_dir = os.path.dirname(os.path.abspath(__file__))11project_root = os.path.dirname(script_dir)12conf_dir = os.path.join(project_root, 'conf')13inputs_dir = os.path.join(script_dir, 'inputs')14 15use_multiprocessing = True16if use_multiprocessing:17    import multiprocessing18    max_cpus = 2 # We might want to not run on the full number of cores, as Rosetta take about 2 Gb of memory per instance19 20###################################################################################################################################################################21# Important: The variables below are set to values that will make the run complete faster (as a tutorial example), but will not give scientifically valid results.22#            Please change them to the "normal" default values before a real run.23###################################################################################################################################################################24 25#rosetta_scripts_path = os.path.expanduser("~/rosetta/source/bin/rosetta_scripts")26rosetta_scripts_path = os.path.expanduser(27"/public/home/scnb9biwet/jiangqq/flex_ddG_tutorial-master/software/rosetta3.9/main/source/bin/rosetta_scripts.default.linuxgccrelease"28)29nstruct = 3 # Normally 3530max_minimization_iter = 5 # Normally 500031abs_score_convergence_thresh = 200.0 # Normally 1.032number_backrub_trials = 10 # Normally 3500033backrub_trajectory_stride = 5 # Can be whatever you want, if you would like to see results from earlier time points in the backrub trajectory. 7000 is a reasonable number, to give you three checkpoints for a 35000 step run, but you could also set it to 35000 for quickest run time (as the final minimization and packing steps will only need to be run one time).34path_to_script = os.path.join(conf_dir, 'ddG-backrub.xml')35residue_to_mutate = ('B', 49, '') # Residue position to perfrom saturation mutatagenesis. Format: (Chain, PDB residue number, insertion code).36 37if not os.path.isfile(rosetta_scripts_path):38    print('ERROR: "rosetta_scripts_path" variable must be set to the location of the "rosetta_scripts" binary executable')39    print('This file might look something like: "rosetta_scripts.linuxgccrelease"')40    raise Exception('Rosetta scripts missing')41 42def run_flex_ddg_saturation( name, input_path, input_pdb_path, chains_to_move, mut_aa, nstruct_i ):43    output_directory = os.path.join( project_root, 'output_saturation', os.path.join( '%s_%s' % (name, mut_aa), '%02d' % nstruct_i ) )44    if not os.path.isdir(output_directory):45        os.makedirs(output_directory)46 47    mutation_chain, mutation_resi, mutation_icode = residue_to_mutate48    resfile_path = os.path.join( output_directory, 'mutate_%s%d%s_to_%s.resfile' % (mutation_chain, mutation_resi, mutation_icode, mut_aa) )49    with open( resfile_path, 'w') as f:50        # Header must be NATAA, not NATRO. This resfile is used both to pick the mutated51        # position and as the task operation for the mutant PackRotamersMover. With NATRO,52        # the mutant branch freezes the 8 A neighbor shell while the wild type branch53        # repacks it, which biases every ddG upwards by several REU (a self-mutation to the54        # native amino acid comes out at +6 instead of 0).55        f.write( 'NATAA\nstart\n%d%s %s PIKAA %s\n' % (mutation_resi, mutation_icode, mutation_chain, mut_aa) )56 57    flex_ddg_args = [58        os.path.abspath(rosetta_scripts_path),59        "-s %s" % os.path.abspath(input_pdb_path),60        '-parser:protocol', os.path.abspath(path_to_script),61        '-parser:script_vars',62        'chainstomove=' + chains_to_move,63        'mutate_resfile_relpath=' + os.path.abspath( resfile_path ),64        'number_backrub_trials=%d' % number_backrub_trials,65        'max_minimization_iter=%d' % max_minimization_iter,66        'abs_score_convergence_thresh=%.1f' % abs_score_convergence_thresh,67        'backrub_trajectory_stride=%d' % backrub_trajectory_stride ,68        '-restore_talaris_behavior',69        '-in:file:fullatom',70        '-ignore_unrecognized_res',71        '-ignore_zero_occupancy false',72        '-ex1',73        '-ex2',74    ]75 76    log_path = os.path.join(output_directory, 'rosetta.out')77 78    print( 'Running Rosetta with args:' )79    print( ' '.join(flex_ddg_args) )80    print( 'Output logged to:', os.path.abspath(log_path) )81    print()82 83    outfile = open(log_path, 'w')84    process = subprocess.Popen(flex_ddg_args, stdout=outfile, stderr=subprocess.STDOUT, close_fds = True, cwd = output_directory)85    returncode = process.wait()86    outfile.close()87 88if __name__ == '__main__':89    mutation_chain, mutation_resi, mutation_icode = residue_to_mutate90    cases = []91    for nstruct_i in range(1, nstruct + 1 ):92        for case_name in os.listdir(inputs_dir):93            case_path = os.path.join( inputs_dir, case_name )94            for f in os.listdir(case_path):95                if f.endswith('.pdb'):96                    input_pdb_path = os.path.join( case_path, f )97                    break98 99            with open( os.path.join( case_path, 'chains_to_move.txt' ), 'r' ) as f:100                chains_to_move = f.readlines()[0].strip()101 102            for mut_aa in 'ACDEFGHIKLMNPQRSTVWY':103                cases.append( ('%s_%s%d%s' % (case_name, mutation_chain, mutation_resi, mutation_icode), case_path, input_pdb_path, chains_to_move, mut_aa, nstruct_i) )104 105    if use_multiprocessing:106        pool = multiprocessing.Pool( processes = min(max_cpus, multiprocessing.cpu_count()) )107 108    for args in cases:109        if use_multiprocessing:110            pool.apply_async( run_flex_ddg_saturation, args = args )111        else:112            run_flex_ddg_saturation( *args )113 114    if use_multiprocessing:115        pool.close()116        pool.join()117