From c562e1f51ad66be8f2b36b7b83a8a2d31db51c9d Mon Sep 17 00:00:00 2001 From: Pablo Martinez-Bulit <59570168+pmbulit@users.noreply.github.com> Date: Thu, 13 Mar 2025 10:34:02 +0000 Subject: [PATCH 1/4] Fix for missing chart data file NO_JIRA --- ...mponent_hydrogen_bond_propensity_report.py | 31 +++++++++++++++---- 1 file changed, 25 insertions(+), 6 deletions(-) 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..0574e49 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,16 @@ - Performs a multi-component HBP calculation for a given library of co-formers """ -import sys import os +import sys +import csv +import time import glob +import json +import random import argparse import tempfile import subprocess -import json import matplotlib @@ -81,7 +84,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,8 +100,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_grouping(min_donor_prob=0.1, min_acceptor_prob=0.1) - return propensities, self.hbp.donors, self.hbp.acceptors + return propensities, groups, self.hbp.donors, self.hbp.acceptors def cm2inch(*tupl): @@ -224,6 +228,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: @@ -368,8 +386,9 @@ def main(structure, work_directory, failure_directory, library, csdrefcode, forc try: hbp_calculator.crystal = crystal hbp_calculator.directory = directory - propensities, donors, acceptors = hbp_calculator.calculate() + propensities, groups, donors, acceptors = hbp_calculator.calculate() coordination_scores = coordination_scores_calc(crystal, directory) + chart_output(groups, directory, crystal) 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) @@ -427,7 +446,7 @@ 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, + parser.add_argument('-f', '--failure_directory', type=str, default=os.getcwd(), help='The location where the failures file should be generated') parser.add_argument('--force_run_disordered', action="store_true", From dc195d700bdc94f48b5c4cc8d37dba82a972c86d Mon Sep 17 00:00:00 2001 From: Pablo Martinez-Bulit <59570168+pmbulit@users.noreply.github.com> Date: Thu, 13 Mar 2025 15:41:58 +0000 Subject: [PATCH 2/4] Removed unnecesary functions NO_JIRA --- ...mponent_hydrogen_bond_propensity_report.py | 222 ++++++------------ 1 file changed, 72 insertions(+), 150 deletions(-) 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 0574e49..e5eb988 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 @@ -105,71 +105,6 @@ def calculate(self): return propensities, groups, 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)) - - def coordination_scores_calc(crystal, directory): # Calculate coordination scores for the target structure @@ -275,70 +210,70 @@ 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 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): @@ -365,8 +300,7 @@ 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() @@ -387,19 +321,13 @@ def main(structure, work_directory, failure_directory, library, csdrefcode, forc hbp_calculator.crystal = crystal hbp_calculator.directory = directory propensities, groups, donors, acceptors = hbp_calculator.calculate() - coordination_scores = coordination_scores_calc(crystal, directory) + chart_output(groups, directory, crystal) - 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: + + + 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) @@ -446,14 +374,9 @@ 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, default=os.getcwd(), - 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: @@ -465,5 +388,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, refcode) From 064ead85a888396ae8663218827123206ae3fcfc Mon Sep 17 00:00:00 2001 From: Pablo Martinez-Bulit <59570168+pmbulit@users.noreply.github.com> Date: Thu, 13 Mar 2025 17:07:51 +0000 Subject: [PATCH 3/4] Removed more uneeded stuff, fixet typo NO_JIRA --- ...mponent_hydrogen_bond_propensity_report.py | 96 +------------------ 1 file changed, 4 insertions(+), 92 deletions(-) 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 e5eb988..8bcb047 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 @@ -18,13 +18,10 @@ import os import sys import csv -import time import glob import json -import random import argparse import tempfile -import subprocess import matplotlib @@ -50,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: @@ -100,7 +91,7 @@ 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_grouping(min_donor_prob=0.1, min_acceptor_prob=0.1) + groups = self.hbp.generate_hbond_groupings(min_donor_prob=0.1, min_acceptor_prob=0.1) return propensities, groups, self.hbp.donors, self.hbp.acceptors @@ -210,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: @@ -284,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] @@ -301,7 +223,6 @@ def main(structure, work_directory, failure_directory, library, csdrefcode, forc coformer_files = glob.glob(os.path.join(library, '*.mol2')) tempdir = tempfile.mkdtemp() - hbp_calculator = PropensityCalc() # for each coformer in the library, make a pair file for the api/coformer and run a HBP calculation @@ -315,7 +236,7 @@ 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 @@ -329,15 +250,6 @@ def main(structure, work_directory, failure_directory, library, csdrefcode, forc print("Propensity calculation failure for %s!" % coformer_name) - # 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__': # Set up the necessary arguments to run the script From 3aa1a472e499e105c0684e1c5c404adbc7aee623 Mon Sep 17 00:00:00 2001 From: Pablo Martinez-Bulit <59570168+pmbulit@users.noreply.github.com> Date: Thu, 13 Mar 2025 17:09:31 +0000 Subject: [PATCH 4/4] More typos NO_JIRA --- .../multi_component_hydrogen_bond_propensity_report.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) 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 8bcb047..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 @@ -292,7 +292,7 @@ def main(structure, work_directory, library, csdrefcode): 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): @@ -300,4 +300,4 @@ def main(structure, work_directory, library, csdrefcode): if not os.path.isdir(args.coformer_library): parser.error('%s - library not found.' % args.coformer_library) - main(args.input_structure, args.directory, args.coformer_library, refcode) + main(args.input_structure, args.directory, args.coformer_library, csdrefcode)