# LIBTBX_SET_DISPATCHER_NAME phenix.table_one from __future__ import division from mmtbx.validation import dummy_validation from mmtbx.validation import molprobity import mmtbx.command_line.molprobity from iotbx import file_reader import iotbx.table_one from libtbx import object_oriented_patterns as oop from libtbx.str_utils import make_header from libtbx.utils import Sorry, null_out from libtbx import runtime_utils import libtbx.phil.command_line from libtbx import easy_pickle from libtbx import easy_mp from libtbx import Auto import libtbx.phil import os.path import time import sys structure_params_str = """ structure .multiple = True .optional = True .short_caption = Structure files .style = auto_align box { name = None .type = str .input_size = 400 .short_caption = Structure name .help = Structure name (will become column label) .style = bold pdb_file = None .type = path .short_caption = PDB file .style = bold file_type:pdb .help = PDB file mtz_file = None .type = path .short_caption = MTZ file .style = bold file_type:mtz OnChange:extract_table_one_labels .help = MTZ file data_labels = None .type = str .input_size = 160 .style = renderer:draw_table_one_label_widget r_free_flags.label = None .type = str .short_caption = R-free flags label .input_size = 160 .style = renderer:draw_table_one_label_widget .help = R-free flags label wavelength = None .type = float cif_file = None .type = path .multiple = True .optional = True .short_caption = CIF file .help = Restraints file cif_directory = None .type = path .short_caption = CIF directory .style = directory .help = Directory containing restraints - all CIF files found will be \ used. data_type = *xray neutron .type = choice unmerged_data = None .type = path .help = Secondary reflections file with unmerged intensities, for \ calculation of merging statistics. .style = bold file_type:hkl OnChange:extract_unmerged_intensities \ help_page:unmerged_data.htm unmerged_labels = None .type = str .input_size = 160 .style = renderer:draw_unmerged_intensities_widget } """ master_phil_str = """ table_one { %s processing { re_compute_r_factors = True .type = bool .short_caption = Always re-compute R-factors .style = bold n_bins = 10 .type = int .short_caption = Number of resolution bins count_anomalous_pairs_separately = False .type = bool .short_caption = Count anomalous pairs separately .help = If true, the program will treat F+ and F- (if present) as \ independent reflections when calculating data statistics. (Not \ recommended.) ligand_selection = None .type = str .short_caption = Ligand atom selection .help = If specified, will determine which atoms or residues are \ counted as ligands (instead of automatic behavior). .input_size = 400 } multiprocessing { include scope libtbx.easy_mp.parallel_phil_str_no_threading } output { directory = None .type = path .help = This is only used by the PHENIX GUI. .short_caption = Output directory .style = bold output_dir job_title = None .type = str .input_size = 400 .style = bold show_missing_fields = True .type = bool format = txt csv *rtf .type = choice(multi=True) .caption = Text CSV RTF .short_caption = Output formats .style = bold base_name = Table1 .type = str .short_caption = Base file name .style = bold verbose = True .type = str text_field_separation = 2 .type = int } }""" % structure_params_str class _ (oop.injector, molprobity.molprobity) : """ Injector dummy class to add as_table1_column() method to the main molprobity object. This extracts statistics from the various validation objects and unmerged data, checks for consistency, and returns an iotbx.table_one.column object. """ def as_table1_column (self, label, wavelength, log, re_compute_r_factors=Auto) : """ Extract information for display in the traditional 'Table 1' of crystallographic statistics in structure articles. """ outer_shell = None data_stats = self.data_stats if (data_stats is None) : data_stats = dummy_validation() merging_stats = dummy_validation() merging_outer = dummy_validation() n_refl_uniq = data_stats.n_refl n_refl_refine = data_stats.n_refl_refine n_free = data_stats.n_free completeness = data_stats.completeness d_max_min = self.d_max_min() d_max, d_min = d_max_min if (self.merging is not None) : merging_stats = self.merging.overall merging_outer = self.merging.bins[-1] n_refl_uniq = merging_stats.n_uniq epsilon = 0.001 if ((merging_stats.d_min > d_min + 2*epsilon) or (merging_stats.d_max < d_max - 2*epsilon)) : raise Sorry(("Resolution limits for unmerged data in the structure "+ "'%s' do not cover the "+ "full range present in the merged data: %g - %g (merged) versus "+ "%g - %g (unmerged)") % (label, d_max, d_min, merging_stats.d_max, merging_stats.d_min)) r_work = self.r_work() r_free = self.r_free() n_tls_groups = None if (self.header_info is not None) : if (self.header_info.n_tls_groups > 0) : n_tls_groups = self.header_info.n_tls_groups use_header_values = (not re_compute_r_factors or (not self.header_info.is_phenix_refinement() and (re_compute_r_factors is Auto))) r_work, r_free, warned = rfactor_sanity_check( r_work_pdb=self.header_info.r_work, r_free_pdb=self.header_info.r_free, r_work_fmodel=r_work, r_free_fmodel=r_free, out=log, structure_name=label, re_compute_r_factors=not use_header_values) if (use_header_values) : n_refl_refine = data_stats.n_refl return iotbx.table_one.column( label=label, space_group=self.space_group_info(), unit_cell=self.unit_cell().parameters(), # data properties wavelength=wavelength, d_max_min=d_max_min, n_refl_all=merging_stats.n_obs, n_refl=n_refl_uniq, multiplicity=merging_stats.mean_redundancy, completeness=completeness, i_over_sigma=merging_stats.i_over_sigma_mean, wilson_b=data_stats.wilson_b, r_sym=merging_stats.r_merge, r_meas=merging_stats.r_meas, cc_one_half=merging_stats.cc_one_half, cc_star=merging_stats.cc_star, # refinement n_refl_refine=n_refl_refine, n_free=n_free, r_work=r_work, r_free=r_free, cc_work=merging_stats.cc_work, cc_free=merging_stats.cc_free, # model properties n_atoms=self.model_stats.all.n_non_hd, n_macro_atoms=self.model_stats.macromolecules.n_non_hd, n_ligand_atoms=self.model_stats.ligands.n_non_hd, n_waterss=self.model_stats.water.n_non_hd, n_residues=self.model_stats.n_protein, bond_rmsd=self.rms_bonds(), angle_rmsd=self.rms_angles(), rama_favored=self.rama_favored(), rama_allowed=self.rama_allowed(), rama_outliers=self.rama_outliers(), rota_outliers=self.rota_outliers(), clashscore=self.clashscore(), adp_mean=self.model_stats.all.b_mean, adp_mean_mm=self.model_stats.macromolecules.b_mean, adp_mean_lig=self.model_stats.ligands.b_mean, adp_mean_wat=self.model_stats.water.b_mean, n_tls_groups=n_tls_groups, anomalous_flag=data_stats.anomalous_flag, ).add_outer_shell( # XXX we need a consistency check here as well d_max_min=(data_stats.d_max_outer, data_stats.d_min_outer), n_refl=data_stats.n_refl_outer, n_refl_all=merging_outer.n_obs, n_refl_refine=data_stats.n_refl_refine_outer, n_free=data_stats.n_free_outer, cc_one_half=merging_outer.cc_one_half, cc_star=merging_outer.cc_star, r_sym=merging_outer.r_merge, r_meas=merging_outer.r_meas, i_over_sigma=merging_outer.i_over_sigma_mean, multiplicity=merging_outer.mean_redundancy, completeness=merging_outer.completeness, cc_work=merging_outer.cc_work, cc_free=merging_outer.cc_free, r_work=data_stats.r_work_outer, r_free=data_stats.r_free_outer) # XXX This is possibly problematic. We should encourage the deposition of as # much data as possible, regardless of what was actually used in the final # refinement. def resolution_sanity_check ( d_max_pdb, d_min_pdb, d_max_f_obs, d_min_f_obs, structure_name, out) : warned = False if (d_max_pdb is not None) : (d_max, d_min) = (d_max_pdb, d_min_pdb) print >> out, \ "Using resolution limits reported in PDB file for structure %s" % \ structure_name if ((d_max != d_max_pdb) or (d_min != d_min_pdb)) : warned = True print >> out, """\ *** WARNING: Resolution limits in the PDB file for structure '%s' are inconsistent with resolution limits in the MTZ file: %.2f - %.2f (in PDB file) %.2f - %.2f (in MTZ file) This is not a fatal error, since the data processing may output more data than are used in refinement, but you should check to make sure that you have not accidentally specified the wrong model file. """ % (structure_name, d_max_pdb, d_min_pdb, d_max_f_obs, d_min_f_obs) else : (d_max, d_min) = (d_max_f_obs, d_min_f_obs) return (d_max, d_min, warned) def rfactor_sanity_check ( r_work_pdb, r_free_pdb, r_work_fmodel, r_free_fmodel, out, structure_name, re_compute_r_factors, tolerance_ignore=0.001, tolerance_warn=0.004) : warned = False if (r_work_pdb is not None) and (r_free_pdb is not None) : if (not re_compute_r_factors) : print >> out, "Using R-factors in PDB header" (r_work, r_free) = (r_work_pdb, r_free_pdb) else : print >> out, "Using re-computed R-factors" (r_work, r_free) = (r_work_fmodel, r_free_fmodel) delta_r_work = abs(r_work_pdb - r_work_fmodel) delta_r_free = abs(r_free_pdb - r_free_fmodel) if ((delta_r_work > tolerance_warn) or (delta_r_free > tolerance_warn)) : warned = True print >> out, """\ *** WARNING: R-factors reported in the PDB file do not match those calculated by phenix.model_vs_data: r_work=%.4f r_free=%.4f (in PDB file) r_work=%.4f r_free=%.4f (from phenix.model_vs_data) This is not a fatal error, but if we can't recalculate the reported R-factors, others probably won't be able to either. This may be due to inconsistent use of a twin operator; phenix.model_vs_data may apply one automatically, but only if it significantly improves R-work. Alternately, it can indicate that the model was further modified after the last round of refinement, that the wrong MTZ file was used, or that a different resolution limit was used for refinement. """ % (r_work_pdb, r_free_pdb, r_work_fmodel, r_free_fmodel) elif ((delta_r_work > tolerance_ignore) or (delta_r_free > tolerance_ignore)) : print >> out, """\ NOTE: small discrepancy in R-factors reported in PDB file versus those calculated by phenix.model_vs_data: r_work=%.4f r_free=%.4f (in PDB file) r_work=%.4f r_free=%.4f (from phenix.model_vs_data) This is probably not a big deal, but please double-check that you used the correct data file and have not changed the resolution limit for refinement. """ % (r_work_pdb, r_free_pdb, r_work_fmodel, r_free_fmodel) else : print >> out, "Using re-computed R-factors for structure %s" % \ structure_name (r_work, r_free) = (r_work_fmodel, r_free_fmodel) return (r_work, r_free, warned) def run_single_structure (params, n_bins, count_anomalous_pairs_separately, ligand_selection=None, log=None) : if (log is None) : log = null_out() molprobity_args = [ "pdb.file_name=\"%s\"" % params.pdb_file, "xray_data.file_name=\"%s\"" % params.mtz_file, "xray_data.labels=\"%s\"" % params.data_labels, "xray_data.r_free_flags.file_name=\"%s\"" % params.mtz_file, "xray_data.r_free_flags.label=\"%s\"" % params.r_free_flags.label, "n_bins=%d" % n_bins, "count_anomalous_pairs_separately=%s" % count_anomalous_pairs_separately, "coot=False", "maps=False", "probe_dots=False", "use_pdb_header_resolution_cutoffs=True", "--quiet", ] + [ "monomers.file_name=\"%s\"" % cif for cif in params.cif_file ] if (ligand_selection is not None) : molprobity_args.append("ligand_selection=\"%s\"" % ligand_selection) if (params.unmerged_data is not None) : molprobity_args.extend([ "unmerged_data.file_name=\"%s\"" % params.unmerged_data, "unmerged_data.labels=\"%s\"" % params.unmerged_labels, ]) if (params.data_type == "neutron") : molprobity_args.append(["scattering_table=neutron", "keep_hydrogens=True"]) if (params.cif_directory is not None) : files = os.listdir(params.cif_directory) for file_name in files : full_path = os.path.join(params.cif_directory, file_name) if (os.path.isdir(full_path)) : continue f = file_reader.any_file(full_path) if (f.file_type == "cif") : molprobity_args.append("monomers.file_name=\"%s\"" % f.file_name) return mmtbx.command_line.molprobity.run(args=molprobity_args, out=log) class table_one (iotbx.table_one.table) : __slots__ = iotbx.table_one.table.__slots__ + [ "output_dir", "params", "output_files", ] #"n_warnings"] def __init__ (self, params, out=sys.stdout) : iotbx.table_one.table.__init__(self, text_field_separation=params.output.text_field_separation, count_anomalous_pairs_separately=\ params.processing.count_anomalous_pairs_separately) self.output_dir = os.getcwd() self.params = params self.output_files = [] make_header("Running data analysis and validation", out=out) results = easy_mp.parallel_map( iterable=range(len(self.params.structure)), func=self.run_single_structure, processes=params.multiprocessing.nproc, method=params.multiprocessing.technology) for structure, result in zip(params.structure, results) : print >> out, "" print >> out, "Collecting stats for structure %s" % structure.name column = result.validation.as_table1_column( label=structure.name, wavelength=structure.wavelength, re_compute_r_factors=params.processing.re_compute_r_factors, log=out) self.add_column(column) def run_single_structure (self, i_struct) : return run_single_structure( params=self.params.structure[i_struct], n_bins=self.params.processing.n_bins, count_anomalous_pairs_separately=self.count_anomalous_pairs_separately, ligand_selection=self.params.processing.ligand_selection) def save_multiple (self, file_base, formats) : for format in formats : file_name = "%s.%s" % (file_base, format) method = getattr(self, "save_%s" % format) method(file_name) self.output_files.append((file_name, "Table as '%s'" % format)) def finish_job (self, job=None) : return (self.output_files, []) def extract_labels (params, out, parameter_scope="structure") : """ Guess MTZ file column labels for experimental data and R-free flags. Only invoked when this program is run from the command line, but the Phenix GUI does something similar. """ for i, structure in enumerate(params.structure) : if (structure.mtz_file is None) : raise Sorry("Missing MTZ file for structure #%d." % (i+1)) if ([structure.data_labels, structure.r_free_flags.label].count(None)>0) : mtz_file = file_reader.any_file(structure.mtz_file, force_type="hkl") mtz_file.assert_file_type("hkl") server = mtz_file.file_server file_name = mtz_file.file_name if (structure.data_labels is None) : print >>out, "Attempting to guess labels for %s..." % file_name data = server.get_xray_data( file_name=file_name, labels=None, ignore_all_zeros=True, parameter_scope=parameter_scope, parameter_name="data_labels") structure.data_labels = data.info().label_string() if (structure.r_free_flags.label is None) : print >>out, "Attempting to guess R-free label for %s..." % file_name rfree = server.get_r_free_flags( file_name=file_name, label=None, test_flag_value=None, disable_suitability_test=False, parameter_scope=parameter_scope+".r_free_flags") structure.r_free_flags.label = rfree[0].info().label_string() def run (args, out=sys.stdout, auto_extract_labels=True, use_current_directory_if_not_specified=False, warn=True) : master_params = libtbx.phil.parse(master_phil_str, process_includes=True) if (len(args) == 0) : print >> out, """\ ************************************************************************ phenix.table_one - statistics harvesting for publication ************************************************************************ note: this is somewhat difficult to configure on the command line at present; you may find it more convenient to use the PHENIX GUI. """ print >> out, "# Parameter template for phenix.table_one:" master_params.show(out=out) print >> out, "# (the 'structure' scope may be copied as many times as " print >> out, "# necessary to handle multiple datasets.)" print >> out, "# Alternate usage:" print >> out, "# phenix.table_one model.pdb data.mtz [logfile]*" return None if (warn) : print >> out, """ note: this is somewhat difficult to configure on the command line at present; you may find it more convenient to use the PHENIX GUI. """ time.sleep(2) master_parmas = libtbx.phil.parse(master_phil_str) interpreter = libtbx.phil.command_line.argument_interpreter( master_phil=master_params, home_scope="table_one") file_phil = [] cmdline_phil = [] pdb_file = None mtz_file = None unmerged_data = None log_files = [] for arg in args : if os.path.isfile(arg) : f = file_reader.any_file(arg) if (f.file_type == "phil") : file_phil.append(f.file_object) elif (f.file_type == "pdb") : pdb_file = f.file_name elif (f.file_type == "hkl") : mtz_file = f.file_name elif (f.file_type == "txt") : log_files.append(f.file_name) else : if arg.startswith("unmerged_data=") : unmerged_data = os.path.abspath("=".join(arg.split("=")[1:])) continue if arg.startswith("--") : arg = arg[2:] + "=True" try : arg_phil = interpreter.process(arg=arg) except RuntimeError : print >> out, "Ignoring unknown argument %s" % arg else : cmdline_phil.append(arg_phil) working_phil = master_params.fetch(sources=file_phil+cmdline_phil) params = working_phil.extract() if (pdb_file is not None) : if (len(params.table_one.structure) > 0) : raise Sorry("You already have a structure defined in the parameter "+ "file; to add structures, you should edit the parameters instead of "+ "specifying additional PDB and data files on the command line.") if (mtz_file is None) : raise Sorry("You have supplied a PDB file, but no corresponding MTZ "+ "file.") log_file_str = "\n".join([ "log_file=%s" % f for f in log_files ]) structure_params = libtbx.phil.parse(structure_params_str) new_structure = structure_params.extract().structure[0] new_structure.pdb_file = pdb_file new_structure.mtz_file = mtz_file new_structure.unmerged_data = unmerged_data params.table_one.structure.append(new_structure) if auto_extract_labels : extract_labels(params.table_one, out=out) if use_current_directory_if_not_specified : if (params.table_one.output.directory is None) : params.table_one.output.directory = os.getcwd() validate_params(params) if (params.table_one.multiprocessing.nproc is None) : params.table_one.multiprocessing.nproc = 1 final_phil = master_params.format(python_object=params) if params.table_one.output.verbose : print >> out, "" print >> out, "#Final effective parameters:" final_phil.show(out=out) print >> out, "#---end" print >> out, "" final_phil.show(out=open("table_one.eff", "w")) table1 = table_one(params.table_one, out=out) easy_pickle.dump("%s.pkl" % params.table_one.output.base_name, table1) table1.save_multiple( file_base=params.table_one.output.base_name, formats=params.table_one.output.format) return table1 def validate_params (params) : if (len(params.table_one.structure) == 0) : raise Sorry("No structures defined.") if (not os.path.isdir(params.table_one.output.directory)) : raise Sorry("Please specify a valid output directory.") for i, struct in enumerate(params.table_one.structure) : if (None in [struct.pdb_file, struct.mtz_file]) : raise Sorry(("Structure #%d is missing either a PDB file or an MTZ "+ "file.") % (i+1)) # elif (None in [struct.data_labels, struct.r_free_flags.label]) : # raise Sorry(("Need both data labels and R-free flag label for MTZ file "+ # "%s (structure #%d).") % (struct.mtz_file, i+1)) if (params.table_one.output.text_field_separation < 1) : raise Sorry("Text field separation must be at least one character.") if (len(params.table_one.output.format) == 0) : raise Sorry("No formats selected for output!") return True class launcher (runtime_utils.target_with_save_result) : def run (self) : if not os.path.exists(self.output_dir) : os.makedirs(self.output_dir) os.chdir(self.output_dir) return run(args=list(self.args), out=sys.stdout, warn=False) if (__name__ == "__main__") : run(args=sys.argv[1:], auto_extract_labels=True, use_current_directory_if_not_specified=True)