from __future__ import with_statement from phaser.pipeline import mr_object from phaser.pipeline import engine from phaser.pipeline import expert from phaser.pipeline import phaser_ai from phaser.pipeline import style from phaser import tbx_utils from phaser import output from phaser import cli from phaser import sculptor import iotbx.phil from libtbx import runtime_utils from libtbx import easy_pickle from libtbx.utils import Sorry, Usage, null_out import subprocess import os.path import sys PROGRAM = "MRage" VERSION = "0.1.0" def read_homology_search_file(file_object): if file_object.name.endswith( ".hhr" ): from iotbx import bioinformatics data = file_object.read() try: result = bioinformatics.hhsearch_parser( output = data ) except ValueError, e: raise Sorry, e return result elif file_object.name.endswith( ".xml" ): from iotbx.bioinformatics import ncbi_blast_xml from iotbx.bioinformatics import ebi_wu_blast_xml parse_modules = [ ncbi_blast_xml, ebi_wu_blast_xml ] for pm in parse_modules: try: result = pm.parse( source = file_object ) except SyntaxError, e: file_object.seek( 0 ) continue return result else: raise Sorry, "No XML parser was able to process this file" elif any( file_object.name.endswith( ext ) for ext in [ ".pir", ".aln", ".clustal", ".fasta", ".fa", ".ali" ] ): from iotbx.bioinformatics import alignment_as_hsearch try: result = alignment_as_hsearch.parse( source = file_object ) except RuntimeError: raise Sorry, "No alignment parser was able to process this file" return result raise Sorry, "Unknown homology search file" def convert_to_file_root(file_name): return os.path.splitext( os.path.basename( file_name ) )[0] def perform_ncbi_homology_search(sequence, name, logger): from iotbx.pdb import fetch file_name = "ncbi_blast_%s.xml" % name logger.info( msg = "Starting BLAST search at NCBI" ) data = fetch.get_ncbi_pdb_blast( sequence = sequence.sequence, file_name = file_name ) logger.info( msg = "Search results saved in %s" % file_name ) from phaser.pipeline import homology parser = homology.get_ncbi_blast_parser() try: result = parser( data = data ) except SyntaxError: raise Sorry, "Invalid BLAST XML format received" return ( result, file_name ) perform_ncbi_homology_search.CAPTION = "NCBI_BLAST" def perform_local_blast_homology_search(sequence, name, logger): cmdline = ( "blastall", "-p", "blastp", "-d", "pdb", "-m", "7", ) logger.info( msg = "Starting local BLAST search" ) try: process = subprocess.Popen( cmdline, stdin = subprocess.PIPE, stdout = subprocess.PIPE, stderr = subprocess.PIPE ) except OSError, e: raise Sorry, "Cannot execute 'blastall': %s" % e ( out, err ) = process.communicate( input = str( sequence ) ) if err: raise Sorry, "Error exit from BLAST: %s" % err file_name = "local_blast_%s.xml" % name with open( file_name, "w" ) as fout: fout.write( out ) logger.info( msg = "Search results saved in %s" % file_name ) from phaser.pipeline import homology parser = homology.get_ncbi_blast_parser() try: result = parser( data = out ) except SyntaxError: raise Sorry, "Invalid BLAST XML format received" return ( result, file_name ) perform_local_blast_homology_search.CAPTION = "Local_BLAST" HOMOLOGY_SEARCH_ON = { "ncbi": perform_ncbi_homology_search, "local": perform_local_blast_homology_search, } def get_proxy_for_multiprocessing(params, logger): return engine.MPJob get_proxy_for_multiprocessing.CAPTION = "multiprocessing" def get_proxy_for_threading(params, logger): return engine.TJob get_proxy_for_threading.CAPTION = "threading" def get_proxy_for_queue(name, command, factory, qinterface_type, qslot_cpus, logger): if command is not None: logger.info( msg = "Using %s as submission command to %s" % ( command, name ) ) import shlex command = tuple( shlex.split( command ) ) if 1 < qslot_cpus: logger.info( msg = "Using %d cpus for single jobs" % qslot_cpus ) from functools import partial return partial( engine.QueueJob.Generic, qinterface = qinterface_type( command = command ), factory = factory, qslot_cpus = qslot_cpus, ) def get_proxy_for_lsf(params, logger): from libtbx.queuing_system_utils import generic return get_proxy_for_queue( name = "LSF", command = params.submission_command, factory = generic.Job, qinterface_type = generic.lsf_interface, qslot_cpus = params.qslot_cpus, logger = logger ) get_proxy_for_lsf.CAPTION = "Load_Sharing_Facility" def get_proxy_for_sge(params, logger): from libtbx.queuing_system_utils import generic return get_proxy_for_queue( name = "SGE", command = params.submission_command, factory = generic.Job, qinterface_type = generic.sge_interface, qslot_cpus = params.qslot_cpus, logger = logger ) get_proxy_for_sge.CAPTION = "Sun_Grid_Engine" def get_proxy_for_pbs(params, logger): from libtbx.queuing_system_utils import generic return get_proxy_for_queue( name = "PBS", command = params.submission_command, factory = generic.PBSJob, qinterface_type = generic.pbs_interface, qslot_cpus = params.qslot_cpus, logger = logger ) get_proxy_for_pbs.CAPTION = "Portable_Batch_System" PROXY_FOR = { "multiprocessing": get_proxy_for_multiprocessing, "threading": get_proxy_for_threading, "lsf": get_proxy_for_lsf, "sge": get_proxy_for_sge, "pbs": get_proxy_for_pbs, } def get_assembly_acceptance_always_policy(): return expert.AssemblyCreate.Always() get_assembly_acceptance_always_policy.CAPTION = "always" def get_assembly_acceptance_observed_policy(): return expert.AssemblyCreate.Observed() get_assembly_acceptance_observed_policy.CAPTION = "observed" def get_assembly_acceptance_never_policy(): return expert.AssemblyCreate.Never() get_assembly_acceptance_never_policy.CAPTION = "never" ASSEMBLY_ACCEPTANCE_POLICY_FOR = { "always": get_assembly_acceptance_always_policy, "observed": get_assembly_acceptance_observed_policy, "never": get_assembly_acceptance_never_policy, } def get_assembly_ensemble_creation_always_policy(): return expert.SearchModelMaintenance.Always() get_assembly_ensemble_creation_always_policy.CAPTION = "always" def get_assembly_ensemble_creation_observed_policy(): return expert.SearchModelMaintenance.Observed() get_assembly_ensemble_creation_observed_policy.CAPTION = "observed" def get_assembly_ensemble_creation_never_policy(): return expert.SearchModelMaintenance.Never() get_assembly_ensemble_creation_never_policy.CAPTION = "never" ASSEMBLY_ENSEMBLE_CREATION_POLICY_FOR = { "always": get_assembly_ensemble_creation_always_policy, "observed": get_assembly_ensemble_creation_observed_policy, "never": get_assembly_ensemble_creation_never_policy, } def get_ks_for_quick_mode(params): rot_expert = expert.RotationSearch( cutoff = params.rotation_peaks_cutoff ) rps_expert = expert.RotationPeakSalvage() level1 = [ expert.FindEnsemble(), expert.TranslationSearch(), expert.PeakEvaluate(), rps_expert, rot_expert, expert.Ensembler(), expert.Sculptor( protocols = params.sculptor_protocols ), expert.DataFetch( site = params.pdb_mirror ), ] level2 = [ expert.Coordinator(), expert.SuperposeSolve(), expert.PooledPackingCheck( pool = params.packing_pool ), expert.ApplyRotationSearchResults( rotation_search_expert = rot_expert ) ] postprocessors = [ expert.Amalgamation(), expert.LocalSymmetryExpert( remove_unsupported_symmetry = False ), expert.AssemblyRemove(), expert.AssemblyUpdate(), ASSEMBLY_ACCEPTANCE_POLICY_FOR[ params.assembly_acceptance_policy ](), ASSEMBLY_ENSEMBLE_CREATION_POLICY_FOR[ params.assembly_ensemble_creation_policy ](), expert.ObjectCompletion(), rps_expert, ] return ( level1, level2, postprocessors ) def get_ks_for_full_mode(params): rot_expert = expert.RotationSearch( cutoff = params.rotation_peaks_cutoff ) level2 = [ expert.FindEnsemble(), expert.PooledPackingCheck( pool = params.packing_pool ), expert.TranslationSearch(), expert.ApplyRotationSearchResults( rotation_search_expert = rot_expert ), rot_expert, expert.Sculptor( protocols = params.sculptor_protocols ), expert.DataFetch( site = params.pdb_mirror ), ] return ( [], level2, [] ) MODE_SETUP_FOR = { "quick": get_ks_for_quick_mode, "full": get_ks_for_full_mode, } SYMMETRY_OPTIONS_FOR = { "dataset": phaser_ai.SolutionProgress.Dataset, "enantiomorph": phaser_ai.SolutionProgress.Enantiomorph, "pointgroup": phaser_ai.SolutionProgress.Pointgroup, } def calculate_rmsd_with_phaser(identity, num_residues): import phaser return phaser.rms_estimate.rms( identity = identity, num_residues = num_residues, ) def calculate_rmsd_with_floored_chothia_lesk(identity, num_residues): import phaser return max( phaser.rms_estimate.chothia_and_lesk( identity = identity ), phaser.rms_estimate.chothia_and_lesk_floor, ) def calculate_rmsd_from_rmsd(identity, num_residues): # rms passed in as "identity" return identity CALCULATE_RMSD_WITH = { "default": calculate_rmsd_with_phaser, "rmsd": calculate_rmsd_from_rmsd, "floored_chothia_and_lesk": calculate_rmsd_with_floored_chothia_lesk, } PHIL_MASTER = iotbx.phil.parse( """ hklin = None .type = path .short_caption = Reflections .style = bold file_type:hkl OnChange:update_phaser_hkl_info input_file labin = None .type = str .short_caption = Data labels .input_size = 160 .style = bold renderer:draw_phaser_labin_widget resolution_cutoff = 2.5 .type = float .input_size = 80 .style = bold no_tncs = False .type = bool .optional = False .short_caption = Do not check for translational NCS crystal_symmetry { unit_cell = None .type=unit_cell .style = bold space_group = None .type = space_group .style = bold } mode = %(mode)s .type = choice .optional = False .short_caption = Search mode symmetry_exploration = %(symmetry)s .type = choice .optional = False .input_size = 160 output { root = "mrage" .type = str .optional = False .style = hidden max_solutions_to_write = 0 .type = int .optional = False .short_caption = Max. solutions to write .input_size = 50 .style = spinner gui_base_dir = None .type = path .short_caption = Output folder .help = Base directory in which the PHENIX GUI will create a \ folder containing the results. Ignored when running the \ program from the command line. .style = bold output_dir save_aniso_data = False .type = bool .expert_level = 3 include scope libtbx.phil.interface.tracking_params } composition { count = None .type = int .optional = True .help = Overall number of copies of the complete molecule (which \ may be multiple components). If left blank, this will be \ guessed based on the expected solvent content. .short_caption = Overall count .input_size = 50 .style = spinner component .multiple = True { sequence = None .type = path .optional = True .short_caption = Component sequence .style = bold file_type:seq input_file stoichiometry = 1 .type=int .optional = False .short_caption = Component stoichiometry .input_size = 50 .help = Number of copies of this component relative to the \ overal macromolecule. This is separate from the overall \ count parameter in the Basic Options tab, which defines the \ number of copies of the complete unit. For instance, if \ you expect the basic molecule to be two copies of this \ component and one copy of a second component, the \ stoichiometry is 2. If you are searching for an unknown \ number of copies of a single component, you should leave \ the stoichiometry set to 1. .style = bold spinner ensemble .multiple = True .optional=True .short_caption = Ensemble { coordinates .multiple = True .optional = True { pdb = None .optional = False .type = path .style = input_file identity = None .optional = False .type = float error_translation = %(error.translation)s .optional = False .type = choice .help = "Identity to RMS conversion" } } model_collection .multiple = True .optional = True { trim = False .type = bool .optional = False .short_caption = Trim superposed models coordinates .multiple = True .optional = True { pdb = None .optional = False .type = path .style = input_file identity = None .optional = False .type = float error_translation = %(error.translation)s .optional = False .type = choice .help = "Identity to RMS conversion" } } template .multiple = True .optional = True { pdb = None .optional = False .type = path .style = input_file alignment = None .optional = True .type = path } homology .multiple = True .optional = True { file_name = None .optional = False .type = path .style = input_file max_hits = 3 .optional = False .type = int } search .optional = False { services = %(search.service)s .type = choice(multi=True) .optional = True .caption = %(search.service.caption)s .short_caption = Automatically search PDB with .help = Performs a BLAST search of the PDB using either \ a locally installed version of BLAST and the PDB \ sequence database (requires separate installation), or \ the NCBI's BLAST server (requires Internet access). max_hits = 3 .optional = False .type = int .short_caption = Max. number of hits to use .input_size = 50 .style = spinner } } } assembly .multiple = True .optional = True .help = "Define known assembly" { use_assembled = False .optional = False .type = bool .help = "Use assembly as model in search or only for completion" local_search .optional = False .help = "Define local search extent" { sweep = 5 .optional = False .type = float .help = "Rotational search radius" extent = 5 .optional = False .type = float .help = "Translational search extent" } component .multiple = True .optional = True .help = "Define a component of the assembly" { file_name = None .optional = False .type = path .help = "Assembly component structure" operation .optional = False .help = "Model transformation to participate in the assembly" { euler = 0 0 0 .optional = False .type = floats( size = 3 ) .help = "Rotation" displacement = 0 0 0 .optional = False .type = floats( size = 3 ) .help = "Translation" } } } queue { technology = %(queue.technology)s .type = choice .optional = False .short_caption = Parallelization method .caption = %(queue.technology.caption)s .input_size = 160 .style = bold renderer:draw_mrage_technology_choice cpus = 1 .type = int .short_caption = Number of CPUs .style = renderer:draw_nproc_widget submission_command = None .type = str .short_caption = Queue submit command qslot_cpus = 1 .type = int .optional = False } packing_pool = 20 .type = int .optional = False rotation_peaks_cutoff = 0.75 .type = float .optional = False post_refinement_cutoff = 0.75 .type = float .optional = False final_selection_cutoff = 0.75 .type = float .optional = False b_factor_refinement = True .type = bool .optional = False .short_caption = Refine B-factors of placed model significant_peak_threshold = 7.0 .type = float .optional = False sculptor_protocols = %(sculptor.protocols)s .type = choice( multi = True ) .optional = False template_equivalence = True .type = bool .optional = False assembly_acceptance_policy = %(assembly.acceptance)s .type = choice .optional = False assembly_ensemble_creation_policy = %(assembly.ensemble.creation)s .type = choice .optional = False pdb_mirror = *rcsb pdbe pdbj .type = choice .optional = False exclude_pdb_ids = None .type = strings .short_caption = Exclude PDB IDs from search .help = Developer option to prevent specific structures from being \ used as search models, in order to test performance with lower \ sequence identity structures. .expert_level = 3 simple_run .style = noauto .help = Parameters for simplified GUI - not intended for command-line use. { enable = False .type = bool .help = GUI parameter - not intended for command-line use. .style = hidden .expert_level = 3 component_sequence = None .type = path .style = input_file file_type:seq .short_caption = Component sequence template_model = None .type = path .multiple = True .optional = True .style = input_file file_type:pdb .short_caption = Template model alignment = None .type = path .multiple = True .optional = True .style = input_file file_type:aln,hhr .short_caption = Sequence alignment blast_services = %(search.service)s .type = choice(multi=True) .optional = True .caption = %(search.service.caption)s .short_caption = Automatically search PDB with .help = Performs a BLAST search of the PDB using either \ a locally installed version of BLAST and the PDB \ sequence database (requires separate installation), or \ the NCBI's BLAST server (requires Internet access). max_blast_hits = 3 .optional = False .type = int .short_caption = Max. number of hits to use .input_size = 50 .style = spinner } """ % { "queue.technology": tbx_utils.choice_string( possible = PROXY_FOR, default = "multiprocessing" ), "queue.technology.caption": " ".join( f.CAPTION for f in PROXY_FOR.values() ), "symmetry": tbx_utils.choice_string( possible = SYMMETRY_OPTIONS_FOR, default = "dataset" ), "mode": tbx_utils.choice_string( possible = MODE_SETUP_FOR, default = "quick" ), "search.service": tbx_utils.choice_string( possible = HOMOLOGY_SEARCH_ON, default = None ), "search.service.caption": " ".join( f.CAPTION for f in HOMOLOGY_SEARCH_ON.values() ), "sculptor.protocols": tbx_utils.choice_string( possible = sculptor.get_known_protocols(), default = "all" ), "assembly.acceptance": tbx_utils.choice_string( possible = ASSEMBLY_ACCEPTANCE_POLICY_FOR, default = "always", ), "assembly.ensemble.creation": tbx_utils.choice_string( possible = ASSEMBLY_ENSEMBLE_CREATION_POLICY_FOR, default = "observed", ), "error.translation": tbx_utils.choice_string( possible = CALCULATE_RMSD_WITH, default = "default", ) }, process_includes=True ) def read_any_format_data(hklin, labin, crystal_symmetry, resolution, no_tncs, logger): from iotbx import file_reader hkl_file = file_reader.any_file( hklin ) hkl_file.assert_file_type( "hkl" ) input_array = None known_arrays = hkl_file.file_object.as_miller_arrays() if labin is not None: for miller_array in known_arrays: array_labels = miller_array.info().label_string() if array_labels == labin: input_array = miller_array break else: logger.info( msg = "Arrays found in %s:" % hklin ) logger.indent() for array in known_arrays: logger.info( msg = array.info().label_string() ) logger.dedent() raise Sorry, "Couldn't find array %s in file %s." % ( labin, hklin ) else: suitable_arrays = [ array for array in known_arrays if ( array.is_xray_amplitude_array() or array.is_xray_intensity_array() or array.is_xray_reconstructed_amplitude_array() ) ] if len( suitable_arrays ) != 1: logger.info( msg = "Suitable arrays found in %s" % hklin ) logger.indent() for array in suitable_arrays: logger.info( msg = array.info().label_string() ) logger.dedent() raise Sorry, "Could not select reflection array: multiple candidates" else: input_array = suitable_arrays[0] logger.info( msg = "Reflection array selected: %s" % input_array.info().label_string() ) assert input_array is not None if crystal_symmetry.space_group is not None: logger.info( msg = "Changing space group to %s" % crystal_symmetry.space_group ) input_array = input_array.customized_copy( space_group_info = crystal_symmetry.space_group ) if crystal_symmetry.unit_cell is not None: logger.info( msg = "Changing unit cell to %s" % crystal_symmetry.unit_cell ) input_array = input_array.customized_copy( unit_cell = crystal_symmetry.unit_cell ) if input_array.space_group_info() is None or input_array.unit_cell() is None: logger.info( msg = "No crystal symmetry information for the data in %s" % hklin ) raise Sorry, "Please define crystal_symmetry parameters first" if input_array.is_xray_intensity_array(): logger.info( msg = "Running French-Wilson amplitude conversion" ) from cctbx import french_wilson input_array = french_wilson.french_wilson_scale( miller_array = input_array, log = logger.info ) logger.info( msg = output.section_close( text = "Amplitude conversion finished" ) ) if input_array.anomalous_flag(): logger.info( msg = "Averaging Bijvoet pairs" ) input_array = input_array.average_bijvoet_mates() if resolution is not None: input_array = input_array.resolution_filter( d_min = resolution ) logger.info( msg = "Data high resolution limit: %.2f A" % input_array.d_min() ) if not input_array.is_xray_amplitude_array(): raise Sorry, "%s is not X-ray amplitude data" % labin import phaser cell = input_array.unit_cell().parameters() hall = input_array.space_group_info().type().hall_symbol() inp = phaser.InputANO() inp.setREFL( input_array.indices(), input_array.data(), input_array.sigmas() ) inp.setCELL6( cell ) inp.setSPAC_HALL( hall ) inp.setMUTE( True ) inp.setHKLO( False ) result = phaser.runANO( inp ) if result.Failed(): logger.info( msg = result.logfile() ) raise Sorry, "Error during anisotropy correction" # Check for presence of TNCS and bail out if detected if not no_tncs: tncsinp = phaser.InputNCS() tncsinp.setREFL( result.getMiller(), result.getCorrectedF(), result.getCorrectedSIGF(), ) tncsinp.setCELL6( cell ) tncsinp.setSPAC_HALL( hall ) tncsinp.setMUTE( True ) tncsinp.setHKLO( False ) tncsres = phaser.runNCS( tncsinp ) if tncsres.hasTNCS(): raise Sorry, "TNCS detected in the data, and MRage does not handle TNCS yet" else: logger.info( "No TNCS detected" ) xray_data = mr_object.XrayData( miller = result.getMiller(), fp = result.getCorrectedF(), sigfp = result.getCorrectedSIGF(), cell = cell, resolution = input_array.d_min() ) return ( xray_data, hall ) def get_number_of_residues(file_name): if not os.path.isfile( file_name ): raise Sorry, "Cannot find PDB file: %s" % file_name import iotbx.pdb root = iotbx.pdb.input( file_name ).construct_hierarchy() if not root.models(): raise Sorry, "PDB file empty: %s" % file_name return len( list( root.models()[0].residue_groups() ) ) def run(args, logger): logger.info( msg = output.banner( text = "%s version %s" % ( PROGRAM, VERSION ) ) ) accumulator = tbx_utils.PhilAccumulator( master_phil = PHIL_MASTER ) accumulator.register_file_handler( handler = tbx_utils.ExtensionFileHandler( extensions = [ ".mtz" ], template = "hklin=%s" ) ) for argument in args: argument.process( accumulator = accumulator ) try: merged_phil = accumulator.merge() except RuntimeError, e: raise Sorry, e logger.info( msg = output.underlined( text = "All configuration options:" ) ) logger.info( msg = merged_phil.as_str() ) try: params = merged_phil.extract() except RuntimeError, e: raise Sorry, e # generate canonical parameters for simple GUI runs if (params.simple_run.enable) : params = combine_simple_params(params) params.simple_run.enable = False # Validate parameters validate_params_basic( params = params ) # Create current state problem = mr_object.Problem( hklin = params.hklin, labin = params.labin, crystal_symmetry = params.crystal_symmetry, resolution = params.resolution_cutoff, no_tncs = params.no_tncs, ) logger.info( msg = "Reading data..." ) problem.load( logger = logger ) if params.output.save_aniso_data: logger.info( msg = "Saving corrected data..." ) problem.enpickle( root = params.output.root ) components = [] models = [] if params.exclude_pdb_ids is not None: excludes = set( c.strip().upper() for c in params.exclude_pdb_ids ) logger.info( "Excluding PDB ids: %s" % ", ".join( excludes ) ) else: excludes = set() for comp in params.composition.component: if comp.sequence is None: logger.warning( msg = "No sequences specified" ) logger.warning( msg = "Molecular weight will be calculated from specified models" ) logger.warning( msg = ( "This is suboptimal for distant models " + "and better results could be obtained by " + "specifying the sequence" ) ) import iotbx.pdb from phaser import rsam from phaser import mmt assert comp.ensemble # check in validate_params_basic weights = [ [ rsam.structure_weight( root = iotbx.pdb.input( coor.pdb ).construct_hierarchy(), mmt = mmt.PROTEIN ) for coor in md.coordinates ] for md in comp.ensemble ] averages = [ sum( w ) / len( w ) for w in weights if w ] assert averages mw = sum( averages ) / len( averages ) c = mr_object.Component( mw = mw ) else: seqfile = tbx_utils.SequenceObject.from_file( file_name = comp.sequence ) if len( seqfile.object ) != 1: raise Sorry, "Sequence file '%s' contain multiple sequences" % comp.sequence sequence = mr_object.SequenceData( sequence = seqfile.object[0] ) from phaser.pipeline import domain_analysis c = mr_object.SequenceComponent( sequence = sequence, selection = domain_analysis.ChainSequenceSelection.from_selection( selection = [ True ] * len( sequence ) ) ) logger.info( msg = "Molecular weight for this component: %.1f" % c.mw ) assert 0 < comp.stoichiometry # check in validate_params_basic components.append( ( c, comp.stoichiometry ) ) for md in comp.ensemble: assert md.coordinates # check in validate_params_basic pdbs = [] rmsds = [] for coor in md.coordinates: # Also does sanity checks on file rms = CALCULATE_RMSD_WITH[ coor.error_translation ]( identity = coor.identity, num_residues = get_number_of_residues( file_name = coor.pdb ), ) logger.info( msg = "%s: quality = %.2f, method = %s => rmsd = %.3f" % ( coor.pdb, coor.identity, coor.error_translation, rms, ) ) rmsds.append( rms ) pdbs.append( coor.pdb ) m = mr_object.MultifileEnsemble( pdbs = pdbs, rmsds = rmsds, composition = mr_object.Composition( components = [ c ] ) ) models.append( m ) for mc in comp.model_collection: assert mc.coordinates # check in validate_params_basic pdbs = [] rmsds = [] for coor in mc.coordinates: # Also does sanity checks on file rms = CALCULATE_RMSD_WITH[ coor.error_translation ]( identity = coor.identity, num_residues = get_number_of_residues( file_name = coor.pdb ), ) logger.info( msg = "%s: quality = %.2f, method = %s => rmsd = %.3f" % ( coor.pdb, coor.identity, coor.error_translation, rms, ) ) rmsds.append( rms ) pdbs.append( coor.pdb ) m = mr_object.ModelCollection( pdbs = pdbs, rmsds = rmsds, trim = mc.trim, composition = mr_object.Composition( components = [ c ] ) ) models.append( m ) for tm in comp.template: if not os.path.isfile( tm.pdb ): raise Sorry, "Cannot find PDB file: %s" % tm.pdb if tm.alignment is None: assert comp.sequence # check in validate_params_basic assert isinstance( c, mr_object.SequenceComponent ) logger.info( msg = "Generating alignment for template" ) from mmtbx.msa import get_muscle_alignment_ordered from iotbx import bioinformatics import iotbx.pdb from phaser import rsam from phaser import mmt root = iotbx.pdb.input( tm.pdb ).construct_hierarchy() ali = get_muscle_alignment_ordered( sequences = [ c.sequence.sequence, bioinformatics.sequence( name = tm.pdb, sequence = rsam.one_letter_sequence( rgs = root.models()[0].chains()[0].residue_groups(), one_letter_for = dict( zip( mmt.PROTEIN.three_letter_codes, mmt.PROTEIN.one_letter_codes ) ), unknown = mmt.PROTEIN.unknown_one_letter, ) ), ] ) logger.info( msg = ali ) file_name = "muscle_%s_%s.aln" % ( convert_to_file_root( file_name = comp.sequence ), convert_to_file_root( file_name = tm.pdb ), ) with open( file_name, "w" ) as outfile: outfile.write( str( ali ) ) logger.info( msg = "Saved alignment as %s" % file_name ) alignment = tbx_utils.AlignmentObject( alignment = ali, name = file_name ) else: alignment = tbx_utils.AlignmentObject.from_file( file_name = tm.alignment ) if comp.sequence: assert isinstance( c, mr_object.SequenceComponent ) t_comp = sequence_component_selection( sequence_data = c.sequence, alignment = alignment.object.copy() ) else: t_comp = c t = mr_object.Template( pdb_file = tm.pdb, alignment = alignment, composition = mr_object.Composition( components = [ t_comp ] ) ) models.append( t ) for hs in comp.homology: try: file_object = open( hs.file_name ) except IOError, e: raise Sorry, "Error while opening: %s" % e search = read_homology_search_file( file_object = file_object ) hitgen = ( ( i, h ) for (i, h ) in enumerate( search.hits(), start = 1 ) if h.identifier.upper() not in excludes ) for (ctr, ( index, hit ) ) in zip( range( hs.max_hits ), hitgen ): if comp.sequence: assert isinstance( c, mr_object.SequenceComponent ) h_comp = sequence_component_selection( sequence_data = c.sequence, alignment = hit.alignment.copy() ) else: h_comp = c h = mr_object.HomologySearchHit( index = index, search = hs.file_name, hit = hit, composition = mr_object.Composition( components = [ h_comp ] ) ) models.append( h ) if comp.search.services: assert comp.sequence # check in validate_params_basic assert isinstance( c, mr_object.SequenceComponent ) for service in comp.search.services: ( search, file_name ) = HOMOLOGY_SEARCH_ON[ service ]( sequence = c.sequence.sequence, name = convert_to_file_root( file_name = comp.sequence ), logger = logger ) hitgen = ( ( i, h ) for (i, h ) in enumerate( search.hits(), start = 1 ) if h.identifier.upper() not in excludes ) for ( ctr, ( index, hit ) ) in zip( range( comp.search.max_hits ), hitgen ): s_comp = sequence_component_selection( sequence_data = c.sequence, alignment = hit.alignment.copy() ) h = mr_object.HomologySearchHit( index = index, search = file_name, hit = hit, composition = mr_object.Composition( components = [ s_comp ] ) ) models.append( h ) if not models: raise Sorry, "No models have been defined" for e in models: problem.data.write( obj = e ) if params.composition.count is None: logger.info( msg = "Guessing number of copies based on Matthews coefficient" ) from phaser import mmt from mmtbx.scaling import matthews count_probabilities = matthews.number_table( components = [ matthews.component( mw = comp.mw * stoich, rho_spec = mmt.PROTEIN.specific_volume ) for ( comp, stoich ) in components ], density_calculator = matthews.density_calculator( crystal = problem.crystal() ) ) if not count_probabilities: logger.warning( msg = "Components do not fit into asymmetric unit" ) raise Sorry, "Asymmetric unit too full" params.composition.count = max( count_probabilities, key = lambda p: p[1] )[0] logger.info( msg = "Most likely count: %d" % params.composition.count ) multifile_ensembles = [ m for m in models if isinstance( m, mr_object.MultifileEnsemble ) ] import scitbx.math import scitbx.matrix for ( i_ass, assembly ) in enumerate( params.assembly, start = 1 ): logger.info( msg = "%s. assembly:" % i_ass ) logger.indent() definitions = [] for ( i_comp, comp ) in enumerate( assembly.component, start = 1 ): logger.info( msg = "%s. component: %s" % ( i_comp, comp.file_name ) ) ense_find = ( mfe for mfe in multifile_ensembles if comp.file_name in mfe.pdbs ) try: ense = ense_find.next() except StopIteration: raise Sorry, "Cannot find component: %s is not defined as an ensemble" % comp.file_name logger.info( msg = "Corresponding ensemble: %s" % ense ) op = scitbx.matrix.rt( ( scitbx.math.euler_angles_zyz_matrix( *comp.operation.euler ), comp.operation.displacement, ) ) definitions.append( ( ense, op ) ) logger.dedent() assert definitions # check in validate_params_basic a = mr_object.Assembly( definitions = definitions, symmetry = mr_object.PointGroup(), strategy = mr_object.SearchStrategy( sweep = assembly.local_search.sweep, extent = assembly.local_search.extent, ), ) problem.data.write( obj = a ) if assembly.use_assembled: regense = mr_object.RegularEnsemble( assembly = a ) problem.data.write( obj = regense ) if params.queue.cpus == 1: logger.info( msg = "Running all calculations on the main thread" ) queue = engine.Queue( nproc = params.queue.cpus, proxy = engine.Job ) else: logger.info( msg = "Calculations will be run in parallel using %d CPUs" % params.queue.cpus ) queue = engine.Queue( nproc = params.queue.cpus, proxy = PROXY_FOR[ params.queue.technology ]( params = params.queue, logger = logger ) ) import operator from phaser.pipeline import superposition state = SYMMETRY_OPTIONS_FOR[ params.symmetry_exploration ]( composition = mr_object.Composition( components = reduce( operator.add, [ [ comp ] * ( stoich * params.composition.count ) for ( comp, stoich ) in components ] ) ), problem = problem, queue = queue, output = logger, superposition = superposition.Proxy( method = superposition.simple_ssm_superposition ), template_equivalence = params.template_equivalence, identity_to_rms = calculate_rmsd_with_phaser, ) # Create solution process control ( level1, level2, postprocessors ) = MODE_SETUP_FOR[ params.mode ]( params = params ) plan = phaser_ai.Plan( level1 = level1, level2 = level2, postprocessors = postprocessors, criterion = lambda peak: params.significant_peak_threshold <= peak.tfz, b_factor_refinement = params.b_factor_refinement, threshold_postrefinement = params.post_refinement_cutoff, ) # Solve state.solve( plan = plan ) # Output logger.info( msg = style.title( text = "Solutions found" ) ) assert state.history assert len( state.cases ) == len( state.history[-1] ) import pickle result_files = [] for ( case, results ) in zip( state.cases, state.history[-1] ): if results is None: continue if all( st == case.problem.root for st in case.score_for ): continue selection = phaser_ai.SolutionSelection( case = case, threshold = params.final_selection_cutoff, template_equivalence = params.template_equivalence, ) print_solutions_as_text( selection = selection, logger = logger ) logger.info( "Top solution:" ) logger.info( selection.solutions[0] ) selection.solutions[0].probability = calculate_success_probability( structure = selection.solutions[0], case = case, stream = logger.info, ) prefix = "%s_%s" % ( params.output.root, selection.case.space_group_suffix_symbol(), ) with open( "%s_results.pkl" % prefix, "w" ) as fout: logger.info( msg = "Writing results to %s" % fout.name ) pickle.dump( ( prefix, selection ), fout ) result_files.append( os.path.abspath( fout.name ) ) for ( index, structure ) in enumerate( selection.solutions, start = 1 ): if params.output.max_solutions_to_write < index: break logger.info( msg = output.underlined( text = "Solution #%s" % index ) ) ( pdb_file, mtz_file ) = run_llg_job( data = selection.case.calculation_data(), structure = structure, root = "%s_%d" % ( prefix, index ), stream = logger.info, ) logger.indent() logger.info( msg = "PDB file: %s" % pdb_file ) logger.info( msg = "Map file: %s" % mtz_file ) logger.dedent() logger.info( msg = output.blank() ) logger.info( msg = "Finished" ) return result_files def calculate_success_probability(structure, case, stream): if getattr( structure, "probability", None ) is not None: probability = structure.probability else: from phaser.pipeline import evaluation stream.write( "Evaluating probability of success\n" ) peaks = [] if 1 < len( structure.peaks ): stream.write( "Multi-molecule structure: recalculate TFZs for molecules (except last)\n" ) from phaser.pipeline import calculation data = case.calculation_data() for ( index, peak ) in enumerate( structure.peaks[:-1], start = 1 ): calc = calculation.TFZCalculation( data = data, partial = mr_object.Structure( peaks = [ p for p in structure.peaks if p != peak ] ), peak = peak, ) raw = calc.run() result = calc.process( results = raw ) tfz = result.data.tfz stream.write( "%d. %s: TFZ=%.2f\n" % ( index, peak, tfz ) ) peaks.append( evaluation.TFZBasedProbabilityPeakDescription( tfz = tfz ) ) # Last peak tfz = structure.tfz stream.write( "%d. %s: TFZ=%.2f\n" % ( len( structure.peaks ), structure.peaks[-1], tfz ) ) peaks.append( evaluation.TFZBasedProbabilityPeakDescription( tfz = tfz ) ) prob = evaluation.TFZBasedProbability( collection = evaluation.TFZ_POLAR_NONPOLAR, description = evaluation.TFZBasedProbabilityStructureDescription( space_group_info = case.space_group_info(), peaks = peaks, ), ) stream.write( "Evaluation for probability of solution being correct:\n" ) stream.write( str( prob ) + "\n" ) probability = prob.total return probability def sequence_component_selection(sequence_data, alignment): from phaser import rsam rsam.alignment_stich( sequence = sequence_data.sequence, alignment = alignment, index = 0 ) from phaser.pipeline import domain_analysis selection = domain_analysis.ChainSequenceSelection.from_selection( selection = [ m != alignment.gap for ( t, m ) in zip( alignment.alignments[0], alignment.alignments[-1] ) if t != alignment.gap ] ) return mr_object.SequenceComponent( sequence = sequence_data, selection = selection ) def print_solutions_as_text(selection, logger): logger.info( msg = output.subtitle( text = "Solution list for space group %s" % selection.case.space_group_symbol() ) ) if selection.solutions: logger.info( msg = "%d solution(s) found" % len( selection.solutions ) ) logger.info( msg = output.blank() ) for ( index, structure ) in enumerate( selection.solutions, start = 1 ): print_solution_info( case = selection.case, index = index, structure = structure, logger = logger ) else: logger.info( msg = "No solutions found" ) def print_solutions_as_xml(selection, logger): from xml.etree.ElementTree import Element, ElementTree root = Element( tag = "solutions", attrib = { "space_group": selection.case.space_group_symbol(), "count": str( len( selection.solutions ) ), } ) for ( index, structure ) in enumerate( selection.solutions, start = 1 ): solution = Element( tag = "solution", attrib = { "index": str( index ), "LLG": "%.2f" % selection.case.score_for[ structure ], "significant": str( structure.significant ) } ) path = search_history_for( structure = structure ) previous = path.next() # there is at least one item annotation = Element( tag = "history" ) solution.append( annotation ) for current in path: extra = len( current.peaks ) - len( previous.peaks ) assert 0 < extra extension = Element( tag = "extension", attrib = { "RFZ": ( "%.1f" % current.rfz ) if current.rfz is not None else "n/a", "TFZ": "%.2f" % current.tfz, "LLG": "%.2f" % selection.case.score_for[ current ], "DLLG": "%.2f" % ( selection.case.score_for[ current ] - selection.case.score_for[ previous ] ), } ) for comp in [ p.ensemble for p in current.peaks[ -extra: ] ]: found = Element( tag = "ensemble", text = str( comp ) ) extension.append( found ) previous = current annotation.append( extension ) root.append( solution ) ElementTree( element = root ).write( logger.info ) logger.info( msg = output.blank() ) def print_solutions_as_script(selection, logger): logger.info( msg = "SPACEGROUP %s" % selection.case.space_group_symbol() ) logger.info( msg = output.blank() ) import operator ensembles = set( reduce( operator.add, [ [ peak.ensemble for peak in structure.peaks ] for structure in selection.solutions ] ) ) logger.info( msg = "# Ensemble definitions" ) for ense in ensembles: logger.info( msg = ense.phaser_keyword_format() ) logger.info( msg = output.blank() ) logger.info( msg = "# Solutions" ) for structure in selection.solutions: mr_set = structure.mr_set() path = search_history_for( structure = structure ) previous = path.next() # there is at least one item annotation = [ "RFZ=%s TFZ=%.2f LLG=%.2f" % ( ( "%.1f" % current.rfz ) if current.rfz is not None else "n/a", current.tfz, selection.case.score_for[ current ], ) for current in path ] mr_set.ANNOTATION =" %s" % " ".join( annotation ) logger.info( msg = mr_set.unparse() ) logger.info( msg = output.blank() ) def print_solution_info(case, index, structure, logger): logger.info( msg = "Solution %d, score=%.3f" % ( index, case.score_for[ structure ] ) ) logger.indent() logger.info( msg = str( structure ) ) logger.info( msg = style.underlined( text = "Search history" ) ) history = search_history_for( structure = structure ) previous = history.next() # there is at least one item if previous.peaks: logger.info( msg = "Starting point, score=%.3f" % ( index, case.score_for[ previous ] ) ) logger.indent() logger.info( msg = str( previous ) ) logger.dedent() logger.info( msg = "" ) lines = [] for current in history: extra = len( current.peaks ) - len( previous.peaks ) assert 0 < extra comps = [ p.ensemble for p in current.peaks[ -extra: ] ] mentry = ( "[ %s ]" % " + ".join( str( e ) for e in comps ) if len( comps ) != 1 else str( comps[0] ) ) lines.append( [ mentry, ( "%.1f" % current.rfz ) if current.rfz is not None else "n/a", current.tfz, case.score_for[ current ], case.score_for[ current ] - case.score_for[ previous ], ] ) previous = current formatted = tbx_utils.format_table_entries( captions = [ "Model", "RFZ", "TFZ", "LLG", "DLLG" ], formats = [ "%s", "%s", "%.1f", "%.1f", "%.1f" ], lines = lines ) logger.info( msg = tbx_utils.assemble_table( lines = formatted ) ) logger.dedent() logger.info( msg = "" ) # XXX for PHENIX GUI def format_solution_history (case, structure) : history = search_history_for( structure = structure ) previous = history.next() # there is at least one item lines = [] for current in history: extra = len( current.peaks ) - len( previous.peaks ) assert 0 < extra comps = [ p.ensemble for p in current.peaks[ -extra: ] ] mentry = ( "[ %s ]" % " + ".join( str( e ) for e in comps ) if len( comps ) != 1 else str( comps[0] ) ) lines.append( [ mentry, current.tfz, case.score_for[ current ], case.score_for[ current ] - case.score_for[ previous ], ] ) previous = current formatted = tbx_utils.format_table_entries( captions = [ "Model", "TFZ", "LLG", "DLLG" ], formats = [ "%s", "%.1f", "%.1f", "%.1f" ], lines = lines ) return formatted[1:] def search_history_for(structure): history = [ structure ] while history[-1].predecessor: history.append( history[-1].predecessor ) return reversed( history ) def format_solution_list (results) : table = [] results = [ r for p,r in results ] # get rid of prefix results_ = sorted(results, lambda x,y: cmp(y.best_score, x.best_score)) for sg_result in results_ : case = sg_result.case space_group = case.space_group_symbol() for index, solution in enumerate(sg_result.solutions, start=1) : row = [ space_group, str(index), "%.2f" % case.score_for[solution] ] table.append(row) return table def get_solution_history (results, space_group, solution_id) : assert (solution_id > 0) assert isinstance(results, list) and isinstance(space_group, str) for prefix, sg_result in results : sg = sg_result.case.space_group_symbol() if (sg == space_group) : for index, solution in enumerate(sg_result.solutions, start=1) : if (index == solution_id) : return format_solution_history(sg_result.case, solution) return None # brunett.solutions functions def structure_output_single(structure, data, mtz_file, stream): import iotbx.pdb from phaser import chisel model = iotbx.pdb.hierarchy.model() chain_ids = chisel.chain_ids() cell = data.uctbx_unit_cell() for p in structure.peaks: r = p.root( cell = cell ) assert r.models() for c in r.models()[0].chains(): copy = c.detached_copy() copy.id = chain_ids.next() model.append_chain( copy ) root = iotbx.pdb.hierarchy.root() root.append_model( model ) root.atoms_reset_serial() return root def structure_output_ensemble(structure, data, mtz_file, stream): import iotbx.pdb from phaser import chisel import itertools import string model = iotbx.pdb.hierarchy.model() chain_ids = chisel.chain_ids() cell = data.uctbx_unit_cell() for p in structure.peaks: altlocs = itertools.cycle( string.ascii_uppercase ) for r in p.roots( cell = cell ): for m in r.models(): for c in m.chains(): copy = c.detached_copy() copy.id = chain_ids.next() altloc = altlocs.next() for rg in copy.residue_groups(): ags = rg.atom_groups() assert ags ags[0].altloc = altloc for ag in ags[1:]: rg.remove_atom_group( ag ) model.append_chain( copy ) root = iotbx.pdb.hierarchy.root() root.append_model( model ) root.atoms_reset_serial() return root def structure_output_combined(structure, data, mtz_file, stream): from solve_resolve.resolve_python.resolve_utils import get_array_dict array_dict,high_resolution,crystal_symmetry=get_array_dict( map_coeffs_array = None, mtz_in = mtz_file, map_coeffs = "FWT,PHWT" ) map_coeffs_array = array_dict[ "map_coeffs" ] from phenix.command_line import combine_models combined_chains = [] cell = data.uctbx_unit_cell() for ( index, cmg ) in enumerate( structure.chain_model_groups( cell = cell ), start = 1 ): stream.write( "### Combining chain group %d ###" % index ) lines = [] for comp in cmg.composition.components: try: seq = comp.residue_sequence() except AttributeError: stream.write( "WARNING: no sequence available for %s\n" % comp ) continue lines.append( "".join( seq ) ) sequence = ">SEQUENCE\n%s\n" % "\n".join( lines ) assert cmg.models current = to_pdb_input_object( model = cmg.models[0].model ) for other in cmg.models[1:]: combination = combine_models.combine_models( args = [ "output_files.pdb_out=None", "output_files.log=None", "output_files.params_out= None", ], quiet = False, out = stream, pdb_input = current, second_pdb_input = to_pdb_input_object( model = other.model ), map_coeffs_array = map_coeffs_array, seq_file_as_string = sequence, ) current = combination.pdb_output root = current.construct_hierarchy() assert len( root.models() ) == 1 combined_chains.extend( root.models()[0].chains() ) from phaser import chisel model = iotbx.pdb.hierarchy.model() chain_ids = chisel.chain_ids() for chain in combined_chains: copy = chain.detached_copy() copy.id = chain_ids.next() model.append_chain( copy ) root = iotbx.pdb.hierarchy.root() root.append_model( model ) root.atoms_reset_serial() return root def structure_output_best_single_model(structure, data, mtz_file, stream): import iotbx.pdb from phaser import chisel model = iotbx.pdb.hierarchy.model() chain_ids = chisel.chain_ids() cell = data.uctbx_unit_cell() constituents = list( structure.atomic_constituents( cell = cell ) ) stream.write( "There are %s constituents in the structure\n" % len( constituents ) ) for ( index, current ) in enumerate( constituents, start = 1 ): inactives = [ c for c in constituents if c != current ] cmgs = list( current.ensemble.chain_model_groups() ) assert len( cmgs ) == 1 scoreds = [] stream.write( "%d. constituent: %s models\n" % ( index,len( cmgs[0].models ) ) ) for ( mindex, cmodel ) in enumerate( cmgs[0].models, start = 1 ): root = iotbx.pdb.hierarchy.root() root.append_model( cmodel.model ) ensemble = mr_object.MultimodelEnsemble( roots = [ root ], rmsds = [ cmodel.rmsd ], composition = cmgs[0].composition, ) peak = mr_object.Peak( ensemble = ensemble, rotation = current.rotation, translation = current.translation, bfactor = current.bfactor, ) result = calculate_llg_map( data = data, structure = mr_object.Structure( peaks = inactives + [ peak ] ), hklo = False, ) assert len( result.getValues() ) == 1 score = result.getValues()[0] scoreds.append( ( peak, score ) ) stream.write( " %d. model: score = %.2f\n" % ( mindex, score ) ) best_pair = max( scoreds, key = lambda p: p[1] ) stream.write( " Best: index = %d, score = %.2f\n" % ( scoreds.index( best_pair ) + 1, best_pair[1], ) ) r = best_pair[0].root( cell = cell ) assert len( r.models() ) == 1 for c in r.models()[0].chains(): copy = c.detached_copy() copy.id = chain_ids.next() model.append_chain( copy ) root = iotbx.pdb.hierarchy.root() root.append_model( model ) root.atoms_reset_serial() return root def to_pdb_input_object(model): import iotbx.pdb root = iotbx.pdb.hierarchy.root() root.append_model( model ) return iotbx.pdb.input( source_info = "string", lines = root.as_pdb_string(), ) def calculate_llg_map(data, structure, root = "PHASER", hklo = True): import phaser input = phaser.InputMR_LLG() data.apply( input = input ) input.setXYZO( False ) input.setHKLO( hklo ) input.setMUTE( True ) input.setROOT( root ) input.setMACA_PROT( "OFF" ) input.setJOBS( 1 ) ensembles = structure.ensembles() for ense in set( ensembles ): ense.apply_ensemble_definition( input ) input.setSOLU( [ structure.mr_set() ] ) result = phaser.runMR_LLG( input ) if not result.Success(): raise Sorry, result.logfile() return result def run_llg_job(data, structure, root, stream, method = structure_output_single): calculate_llg_map( data = data, structure = structure, root = root, hklo = True ) mtz_file = "%s.1.mtz" % root assert os.path.isfile( mtz_file ) pdb_file = "%s.1.pdb" % root pdb_data = method( structure = structure, data = data, mtz_file = mtz_file, stream = stream, ) pdb_data.write_pdb_file( file_name = pdb_file, crystal_symmetry = data.crystal_symmetry() ) assert os.path.isfile( pdb_file ) return ( pdb_file, mtz_file ) # Validation functions def validate_params (params) : if (params.simple_run.enable) : params = combine_simple_params(params) validate_params_basic( params = params ) if (params.labin is None) : raise Sorry("Please select labels for the data to use. Both amplitudes "+ "or intensities and sigmas are required.") if (params.crystal_symmetry.space_group is None) : raise Sorry("Space group undefined.") if (params.crystal_symmetry.unit_cell is None) : raise Sorry("Unit cell undefined.") if (params.output.gui_base_dir is None) : raise Sorry("Please specify an output directory.") return True def validate_params_basic (params) : if (params.hklin is None) : raise Sorry("No reflections file defined!") if (params.output.root is None) : raise Sorry("Please specify a base output file name.") if (len(params.composition.component) == 0) : raise Sorry("At least one component (including search model(s), "+ "template model(s), or homology result(s)) must be defined.") for i, component in enumerate(params.composition.component, start=1) : if (component.sequence is None) and (not component.ensemble): raise Sorry("Molecular weight needed for component %d." % i) if any( len( md.coordinates ) == 0 for md in component.ensemble ): raise Sorry("Empty ensemble definition for component %d" % i) if any( len( md.coordinates ) == 0 for md in component.model_collection ): raise Sorry("Empty model_collection definition for component %d" % i) if (component.stoichiometry is None) : raise Sorry("Must set stoichiometry for component %d." % i) if (component.stoichiometry < 0) : raise Sorry("Stoichiometry out-of-range: 0 <= stoichiometry") if (component.sequence is None) and any( tm.alignment is None for tm in component.template ) : raise Sorry("You must define the target sequence for component %d (required for alignment)." % i) if (component.sequence is None and component.search.services): raise Sorry("You must define the target sequence for component %d (required for homology_search)" % i) if not any( component_has_models( component = c ) for c in params.composition.component ): raise Sorry( "You must either at least one automatic homology search " + "to perform, or supply one or more ensembles, processed model collections, " + "template models, or homology search results for at least one component" ) if params.packing_pool <= 0: raise Sorry("Specify positive number for packing_pool") if params.rotation_peaks_cutoff < 0 or 1 < params.rotation_peaks_cutoff: raise Sorry("rotation_peaks_cutoff out-of-range: 0 <= RPC <=1") if any( len( a.component ) == 0 for a in params.assembly ): raise Sorry("Empty assembly definition") if (sys.platform == "win32") : if (params.queue.technology not in ["multiprocessing", "threading"]) : raise Sorry("Only multiprocessing and threading are allowed for "+ "parallel processing when running on Windows.") return True def component_has_models(component): return ( len( component.search.services ) != 0 or len( component.ensemble ) != 0 or len( component.model_collection ) != 0 or len( component.template ) != 0 or len( component.homology ) != 0 ) # XXX messy hack to allow automatic generation of simplified GUI with a single # list control for input files. since PHIL will generate a separate component # scope for each sub-parameter, this method merges their contents into a # single scope. the parameters are then converted back to PHIL objects, which # are then saved to the config file used to run the job. # # (this sounds uglier than it really is - since the GUI limits inputs to # sequence, template models, and BLAST search parameters, it is very simple to # assemble a single scope.) def combine_simple_params (params) : param_phil = iotbx.phil.parse(""" composition.component.sequence = \"%s\" """ % params.simple_run.component_sequence) new_params = PHIL_MASTER.fetch(source=param_phil).extract() params.composition.__phil_join__( new_params.composition ) assert (len(params.composition.component) == 1) params.composition.component[0].search.services = \ params.simple_run.blast_services params.composition.component[0].search.max_hits = \ params.simple_run.max_blast_hits for template in params.simple_run.template_model : param_phil = iotbx.phil.parse(""" composition.component.template.pdb = \"%s\" """ % template) new_params = PHIL_MASTER.fetch(source=param_phil).extract() params.composition.component[0].template.append( new_params.composition.component[0].template[0]) for alignment in params.simple_run.alignment : param_phil = iotbx.phil.parse(""" composition.component.homology.file_name = \"%s\" """ % alignment) new_params = PHIL_MASTER.fetch(source=param_phil).extract() params.composition.component[0].homology.append( new_params.composition.component[0].homology[0]) return params class launcher (runtime_utils.target_with_save_result) : def run (self) : if (not os.path.isdir(self.output_dir)) : os.makedirs(self.output_dir) os.chdir(self.output_dir) parser = cli.PCLIParser(master_phil=PHIL_MASTER) params = parser.parse_args(args=self.args) logger = output.SingleStream(stream=sys.stdout, level=1, gui=True) result = run(args=params.phils, logger=logger) if hasattr(sys.stdout, "_flush") : sys.stdout._flush() # XXX: gross # the result is just the name of the final pickle file return result class write_files (object) : def __init__ (self, results, space_group, solution_id, output_dir) : self.results = results self.space_group = space_group self.solution_id = solution_id self.output_dir = output_dir def __call__ (self, *args, **kwds) : os.chdir(self.output_dir) for prefix, sg_result in self.results : sg = sg_result.case.space_group_symbol() if (sg == self.space_group) : for index, structure in enumerate(sg_result.solutions, start=1) : if (index == self.solution_id) : logger = output.SingleStream(stream=sys.stdout) sg_result.case.problem.load(logger=logger) (pdb_file, mtz_file) = run_llg_job( data = sg_result.case.calculation_data(), structure = structure, root = "%s_%d" % (prefix, self.solution_id), stream = logger.info) return (os.path.abspath(pdb_file), os.path.abspath(mtz_file)) return (None, None)