Adds physicochemical properties to output

.gitignore

@@ -141,4 +141,5 @@ Benchmarking/
 Ertl_benchmarking/
 
 # Ignore any graphics folders
-*.ai
+*.ai
+repotrectinib_hcie_results_june_2025/

hcie/molecule.py

@@ -4,7 +4,7 @@
 from collections import defaultdict
 import rdkit
 from rdkit import Chem
-from rdkit.Chem import AllChem, rdDetermineBonds
+from rdkit.Chem import AllChem, rdDetermineBonds, Descriptors, Crippen, rdMolDescriptors
 from rdkit.Geometry import Point3D
 
 from hcie.constants import distance_bins, angle_bins
@@ -37,6 +37,7 @@ def __init__(
         self.exit_vectors = []
         self.exit_vector_properties = []
         self.exit_vector_properties_by_hash = []
+        self.phys_chem_props = {}
 
         self.shape_scores = {}
         self.esp_scores = {}
@@ -81,6 +82,9 @@ def initialize(self):
         # Generate dictionary of exit vectors ordered by hash
         self.exit_vector_properties_by_hash = self.get_exit_vector_properties_by_hash()
 
+        # Generate the dictionary of physicochemical properties
+        self.phys_chem_props = self.calculate_physicochemical_properties()
+
         return None
 
     @classmethod
@@ -686,6 +690,35 @@ def get_atom_ids_of_ring_plane(self, functionalisable_bond: tuple) -> tuple:
 
         return functionalisable_bond[0], neighbours[0], neighbours[1]
 
+    def calculate_physicochemical_properties(self):
+        """
+        Calculates the physicochemical properties of the parent molecule.
+        Returns
+        -------
+        dict[str, float]: the physicochemical properties of the parent molecule
+        {
+        "MW": molecular weight,
+        "cLogP": Crippen lipophilicity,
+        "TPSA": Total polar surface area,
+        "HBD": number of hydrogen bond donors,
+        "HBA": number of hydrogen bond acceptors,
+        "HeavyAtoms": number of heavy atoms,
+        "Heteroatoms": number of heteroatoms
+        }
+
+        """
+        clean_smiles = self.smiles.replace("[*]", "[H]").replace("*", "[H]")
+        mol = Chem.MolFromSmiles(clean_smiles)
+        return {
+            "MW": Descriptors.MolWt(mol),
+            "cLogP": Crippen.MolLogP(mol),
+            "TPSA": rdMolDescriptors.CalcTPSA(mol),
+            "HBD": rdMolDescriptors.CalcNumHBD(mol),
+            "HBA": rdMolDescriptors.CalcNumHBA(mol),
+            "HeavyAtoms": rdMolDescriptors.CalcNumHeavyAtoms(mol),
+            "Heteroatoms": rdMolDescriptors.CalcNumHeteroatoms(mol),
+        }
+
 
 class MoleculeError(Exception):
     """Base class for exceptions in the Molecule class."""

hcie/outputs.py

@@ -3,6 +3,7 @@
 """
 
 import os
+import csv
 from datetime import datetime
 
 from rdkit import Chem
@@ -37,49 +38,52 @@ def wrapper(*args, **kwargs):
 
 
 @new_directory
-def print_results(results: list, query_smiles: str, query_name: str) -> None:
+def print_results(mol_dict: dict, results: list, query_smiles: str, query_name: str) -> None:
     """
     Prints out the results of a search against the VEHICLe database.
+    :param mol_dict:
     :param results:
     :param query_smiles:
     :param query_name:
     :return:
     """
     current_datetime = datetime.now()
-    with open(f"{query_name}_results.txt", "w") as output_file:
-        # Write the title line
-        output_file.write(
-            f'Generated by HCIE at {current_datetime.strftime("%H:%M, %d/%m/%Y")}'
-            + "\n"
-        )
-        # Write the query file out
-        output_file.write(f"Query molecule: {query_name}" + "\n")
-        output_file.write(f"Query SMILES: {query_smiles}" + "\n")
-
-        output_headers = [
+    filename = f"{query_name}_results.csv"
+    query_props = mol_dict["query"].phys_chem_props
+
+    with open(filename, "w", newline="") as csvfile:
+        writer = csv.writer(csvfile)
+
+        # Metadata rows
+        writer.writerow(["Generated by HCIE:", current_datetime.strftime("%H:%M, %d/%m/%Y")])
+        writer.writerow(["Query molecule:", query_name])
+        writer.writerow(["Query SMILES:", query_smiles])
+
+        writer.writerow(["MW", "cLogP", "total polar surface area (TPSA)", "hydrogen bond donors (HBD)", "hydrogen bond acceptors (HBA)", "heavy atoms", "heteroatoms"])
+        prop_keys = ["MW", "cLogP", "TPSA", "HBD", "HBA", "HeavyAtoms", "Heteroatoms"]
+        query_props_list = [query_props[k] for k in prop_keys]
+
+        writer.writerow(query_props_list)
+        writer.writerow([]) # blank line before table starts
+
+        # Header row
+        writer.writerow([
             "Rank",
             "RegID",
             "SMILES",
             "Score",
             "ESP Score",
             "Shape Score",
+            "Delta_MW",
+            "Delta_cLogP",
+            "Delta_TPSA",
+            "Delta_HBD",
+            "Delta_HBA",
+            "Delta_Heavy atoms",
+            "Delta_Heteroatoms",
             "Conformer ID",
-        ]
-
-        output_file.write(
-            "-" * 120
-            + "\n"
-            + f"{output_headers[0]:10}\t"
-            + f"{output_headers[1]:6}\t"
-            + f"{output_headers[2]:35}\t"
-            + f"{output_headers[3]:3}\t"
-            + f"{output_headers[4]:3}\t"
-            + f"{output_headers[5]:3}\t"
-            + f"{output_headers[6]}"
-            + "\n"
-            + "-" * 120
-            + "\n"
-        )
+        ])
+
         for rank, result in enumerate(results):
             regid, score, conf_id, esp_score, shape_score, smiles = (
                 str(result[0]),
@@ -93,16 +97,28 @@ def print_results(results: list, query_smiles: str, query_name: str) -> None:
             # RDKit insists on putting dummy atoms in SMILES that aren't recognised by ChemDraw, so these need replacing
             smiles = smiles.replace("[*:1]", "[R1]").replace("[*:2]", "[R2]")
 
-            row_line = (
-                f"{rank:<10}\t"
-                + f"{regid:6}\t"
-                + f"{smiles:35}\t"
-                + f"{score:<5.2f}\t"
-                + f"{esp_score:<9.2f}\t"
-                + f"{shape_score:<11.2f}\t"
-                + f"{conf_id}"
-            )
-            output_file.write(row_line + "\n")
+            # Calculate the deltas (difference between query and probe) for the physicochemical properties
+            deltas = {
+                prop: mol_dict[regid].phys_chem_props[prop] - query_props[prop]
+                for prop in prop_keys
+            }
+
+            writer.writerow([
+                rank,
+                regid,
+                smiles,
+                f"{score:.2f}",
+                f"{esp_score:.2f}",
+                f"{shape_score:.2f}",
+                f"{deltas['MW']:.1f}",
+                f"{deltas['cLogP']:.2f}",
+                f"{deltas['TPSA']:.2f}",
+                int(deltas['HBD']),
+                int(deltas['HBA']),
+                int(deltas['HeavyAtoms']),
+                int(deltas['Heteroatoms']),
+                conf_id
+            ])
 
     return None
 

hcie/vehicle_search.py

@@ -130,7 +130,7 @@ def results_to_file(self, results: list, mols: dict) -> None:
         mols["query"] = self.query
         query_label = self.query.smiles if self.query.smiles else f"<XYZ:{self.query.name or 'query'}>"
         print_results(
-            results, query_smiles=query_label, query_name=self.query.name
+            mols, results, query_smiles=query_label, query_name=self.query.name
         )
         alignments_to_sdf(
             results=results, mol_alignments=mols, query_name=self.query.name