diff --git a/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.py b/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.py index 7adb2dd..d23ea85 100644 --- a/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.py +++ b/scripts/multi_component_hydrogen_bond_propensity/multi_component_hydrogen_bond_propensity_report.py @@ -15,13 +15,13 @@ - Performs a multi-component HBP calculation for a given library of co-formers """ -import sys import os +import sys +import csv import glob +import json import argparse import tempfile -import subprocess -import json import matplotlib @@ -47,12 +47,6 @@ """.format(sys.executable) raise ImportError(error_message) -SCRIPT_DIR = os.path.dirname(__file__) -TEMPLATE_FILENAME = 'multi_component_hydrogen_bond_propensity_report.docx' -TEMPLATE_FILE = os.path.join(SCRIPT_DIR, TEMPLATE_FILENAME) -PAIR_TEMPLATE_FILENAME = 'multi_component_pair_hbp_report.docx' -PAIR_TEMPLATE_FILE = os.path.join(SCRIPT_DIR, PAIR_TEMPLATE_FILENAME) - ############################################################################### class PropensityCalc: @@ -81,7 +75,7 @@ def calculate(self): print(self.hbp.functional_groups) # Generate Training Dataset - self.hbp.match_fitting_data(count=500) # set to 500 for better representation of functional groups + self.hbp.match_fitting_data(count=None) # set to 500 for better representation of functional groups self.hbp.analyse_fitting_data() @@ -97,73 +91,9 @@ def calculate(self): print('Area under ROC curve: {} -- {}'.format(round(model.area_under_roc_curve, 3), model.advice_comment)) propensities = self.hbp.calculate_propensities() + groups = self.hbp.generate_hbond_groupings(min_donor_prob=0.1, min_acceptor_prob=0.1) - return propensities, self.hbp.donors, self.hbp.acceptors - - -def cm2inch(*tupl): - inch = 2.54 - if isinstance(tupl[0], tuple): - return tuple(i / inch for i in tupl[0]) - else: - return tuple(i / inch for i in tupl) - - -def launch_word_processor(output_file): - '''This function launches the platform specific word processor - when not running under continuous integration''' - if 'TEAMCITY_VERSION' in os.environ: - return - if sys.platform == 'win32': - os.startfile(output_file) - elif sys.platform.startswith('linux'): - subprocess.Popen(['xdg-open', output_file]) - else: - subprocess.Popen(['open', output_file]) - - -def make_diagram(mol, directory): - # Generates a diagram from a given structure - molecule_diagram_generator = DiagramGenerator() - molecule_diagram_generator.settings.line_width = 1.6 - molecule_diagram_generator.settings.font_size = 12 - molecule_diagram_generator.settings.image_height = 300 - img = molecule_diagram_generator.image(mol) - fname = str(os.path.join(directory, '%s_diagram.png' % mol.identifier)) - if img: - img.save(fname) - return fname - - -def make_mc_chart(dictionary, directory, mol): - results = [value[0] for key, value in dictionary if isinstance(value[0], float)] - ymin = min(results) - ymax = max(results) - indices = range(1, len(results) + 1) - - fig = plt.figure(figsize=cm2inch(22, 18)) - ax = fig.add_subplot(1, 1, 1) - color = 'cornflowerblue' - color1 = 'royalblue' - ls = '' - plt.plot(indices, results, marker='D', markersize=10, color=color, ls=ls, markeredgecolor=color1, alpha=0.7) - plt.axhline(y=0, color='gray') - plt.xlabel('Co-Former Rank', fontweight='bold', fontsize='12') - plt.ylabel('Multi-Component Score', fontweight='bold', fontsize='12') - plt.title('MCHBP screening results', fontweight='bold', fontsize='12') - ax.axhspan(0.025, ymax + 0.1, facecolor='lightgreen', alpha=0.5) - ax.axhspan(-0.025, ymin - 0.1, facecolor='lightpink', alpha=0.5) - ax.axhspan(0.025, -0.025, facecolor='grey', alpha=0.5) - plt.ylim(ymin - 0.025, ymax + 0.025) - fname = str(os.path.join(directory, '%s_MC_HBP_plot.png' % mol.identifier)) - plt.savefig(fname, format='png', dpi=600) - return fname - - -def add_picture_subdoc(picture_location, docx_template, wd=7): - # This function adds a picture to the .docx file - return docxtpl.InlineImage( - docx_template, image_descriptor=picture_location, width=Cm(wd)) + return propensities, groups, self.hbp.donors, self.hbp.acceptors def coordination_scores_calc(crystal, directory): @@ -224,6 +154,20 @@ def get_mc_scores(propensities, identifier): identifier] +def chart_output(groups, directory, mol): + # Write out the data points of the HBP chart to a file + with open(os.path.join(directory, f'{mol.identifier}_chart_data.csv'), 'w', newline='') as outfile: + csv_writer = csv.writer(outfile) + csv_writer.writerow(['Mean Propensity', 'Mean Coordination Score', 'Hydrogen Bonds']) + # Write the remaining putative networks + for group in groups: + csv_writer.writerow( + [group.hbond_score, + group.coordination_score, + '; '.join(['%s - %s' % (g.donor.label, g.acceptor.label) for g in group.hbonds])] + ) + + def make_pair_file(api_molecule, tempdir, f): # Creates a file for the api/coformer pair with io.MoleculeReader(f) as reader: @@ -257,73 +201,7 @@ def make_molecule_pair(api_molecule, coformer_molecule): return molecule_pair -def pair_output(identifier, propensities, donors, acceptors, coordination_scores, directory): - # Writes out the output from a single HBP calculation for the multi-component pair - - # This looks for the .docx template that is used to generate the report from - if os.path.isfile(PAIR_TEMPLATE_FILE): - docx_template = docxtpl.DocxTemplate(PAIR_TEMPLATE_FILE) - else: - print('Error! {} not found!'.format(PAIR_TEMPLATE_FILENAME)) - quit() - - dscores = {} - ascores = {} - - for d in donors: - coord_cols = [i for i in range(len(coordination_scores.predictions_for_label(d.label, 'd')[1]))] - dscores[d.label] = [round((coordination_scores.predictions_for_label(d.label, 'd')[1])[j], 3) - for j in coord_cols] - - for a in acceptors: - ascores[a.label] = [round((coordination_scores.predictions_for_label(a.label, 'a')[1])[k], 3) - for k in coord_cols] - - context = { - 'identifier': identifier, - 'propensities': propensities, - 'coord_cols': coord_cols, - 'donors': donors, - 'dscores': dscores, - 'acceptors': acceptors, - 'ascores': ascores - } - docx_template.render(context) - output_file = os.path.join(directory, '%s_pair_output.docx' % identifier) - docx_template.save(output_file) - - -def make_mc_report(identifier, results, directory, diagram_file, chart_file): - # Write the MC-HBP report from the results - - # This looks for the .docx template that is used to generate the report from - if os.path.isfile(TEMPLATE_FILE): - docx_template = docxtpl.DocxTemplate(TEMPLATE_FILE) - else: - print('Error! {} not found!'.format(TEMPLATE_FILENAME)) - quit() - - # Generate content for the report - diagram = add_picture_subdoc(diagram_file, docx_template) - chart = add_picture_subdoc(chart_file, docx_template, wd=18) - - # The context is the information that is given to the template to allow it to be populated - context = { - 'identifier': str(identifier).split('.')[0], - 'diagram': diagram, - 'chart': chart, - 'results': results - } - - # Send all the information to the template file then open up the final report - docx_template.render(context) - output_file = os.path.join(directory, '%s_MC_HBP_report.docx' % str(identifier).split('.')[0]) - docx_template.save(output_file) - - launch_word_processor(output_file) - - -def main(structure, work_directory, failure_directory, library, csdrefcode, force_run): +def main(structure, work_directory, library, csdrefcode): # This loads up the CSD if a refcode is requested, otherwise loads the structural file supplied if csdrefcode: try: @@ -331,10 +209,7 @@ def main(structure, work_directory, failure_directory, library, csdrefcode, forc except RuntimeError: print('Error! %s is not in the database!' % structure) quit() - if io.CrystalReader('CSD').entry(structure).has_disorder and not force_run: - raise RuntimeError("Disorder can cause undefined behaviour. It is not advisable to run this " - "script on disordered entries.\n To force this script to run on disordered entries" - " use the flag --force_run_disordered.") + else: crystal = io.CrystalReader(structure)[0] @@ -347,8 +222,6 @@ def main(structure, work_directory, failure_directory, library, csdrefcode, forc # find the coformers and set up the calculations coformer_files = glob.glob(os.path.join(library, '*.mol2')) tempdir = tempfile.mkdtemp() - mc_dictionary = {} - failures = [] hbp_calculator = PropensityCalc() @@ -363,33 +236,19 @@ def main(structure, work_directory, failure_directory, library, csdrefcode, forc if os.path.exists(os.path.join(directory, "success.json")): with open(os.path.join(directory, "success.json"), "r") as file: tloaded = json.load(file) - mc_dictionary[coformer_name] = tloaded + else: try: hbp_calculator.crystal = crystal hbp_calculator.directory = directory - propensities, donors, acceptors = hbp_calculator.calculate() - coordination_scores = coordination_scores_calc(crystal, directory) - pair_output(crystal.identifier, propensities, donors, acceptors, coordination_scores, directory) - with open(os.path.join(directory, "success.json"), "w") as file: - tdata = get_mc_scores(propensities, crystal.identifier) - json.dump(tdata, file) - mc_dictionary[coformer_name] = get_mc_scores(propensities, crystal.identifier) - except Exception as error_message: + propensities, groups, donors, acceptors = hbp_calculator.calculate() + + chart_output(groups, directory, crystal) + + + except (RuntimeError, TypeError): print("Propensity calculation failure for %s!" % coformer_name) - error_string = f"{coformer_name}: {error_message}" - warnings.warn(error_string) - mc_dictionary[coformer_name] = ["N/A", "N/A", "N/A", "N/A", "N/A", crystal.identifier] - failures.append(error_string) - # Make sense of the outputs of all the calculations - mc_hbp_screen = sorted(mc_dictionary.items(), key=lambda e: 0 if e[1][0] == 'N/A' else e[1][0], reverse=True) - diagram_file = make_diagram(api_molecule, work_directory) - chart_file = make_mc_chart(mc_hbp_screen, work_directory, api_molecule) - make_mc_report(structure, mc_hbp_screen, work_directory, diagram_file, chart_file) - if failure_directory is not None: - with open(os.path.join(failure_directory, 'failures.txt'), 'w', encoding='utf-8', newline='') as file: - file.write('\n'.join(map(str, failures))) if __name__ == '__main__': @@ -427,18 +286,13 @@ def main(structure, work_directory, failure_directory, library, csdrefcode, forc parser.add_argument('-c', '--coformer_library', type=str, help='the directory of the desired coformer library', default=ccdc_coformers_dir) - parser.add_argument('-f', '--failure_directory', type=str, - help='The location where the failures file should be generated') - - parser.add_argument('--force_run_disordered', action="store_true", - help='Forces running the script on disordered entries. (NOT RECOMMENDED)', default=False) args = parser.parse_args() - refcode = False + csdrefcode = False args.directory = os.path.abspath(args.directory) if not os.path.isfile(args.input_structure): if len(str(args.input_structure).split('.')) == 1: - refcode = True + csdefcode = True else: parser.error('%s - file not found.' % args.input_structure) if not os.path.isdir(args.directory): @@ -446,5 +300,4 @@ def main(structure, work_directory, failure_directory, library, csdrefcode, forc if not os.path.isdir(args.coformer_library): parser.error('%s - library not found.' % args.coformer_library) - main(args.input_structure, args.directory, args.failure_directory, args.coformer_library, refcode, - args.force_run_disordered) + main(args.input_structure, args.directory, args.coformer_library, csdrefcode)