forked from gcorso/DiffDock
/
evaluate_files.py
180 lines (161 loc) · 11.8 KB
/
evaluate_files.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
# small script to extract the ligand and save it in a separate file because GNINA will use the ligand position as initial pose
import os
import time
from argparse import FileType, ArgumentParser
import numpy as np
from biopandas.pdb import PandasPdb
from rdkit import Chem
from tqdm import tqdm
from datasets.pdbbind import read_mol
from datasets.process_mols import read_molecule
from utils.utils import read_strings_from_txt, get_symmetry_rmsd
parser = ArgumentParser()
parser.add_argument('--config', type=FileType(mode='r'), default=None)
parser.add_argument('--data_dir', type=str, default='data/PDBBind_processed', help='')
parser.add_argument('--results_path', type=str, default='results/user_predictions_testset', help='Path to folder with trained model and hyperparameters')
parser.add_argument('--file_suffix', type=str, default='_baseline_ligand.pdb', help='Path to folder with trained model and hyperparameters')
parser.add_argument('--project', type=str, default='ligbind_inf', help='')
parser.add_argument('--file_to_exclude', type=str, default=None, help='')
parser.add_argument('--all_dirs_in_results', action='store_true', default=True, help='Evaluate all directories in the results path instead of using directly looking for the names')
parser.add_argument('--num_predictions', type=int, default=10, help='')
parser.add_argument('--no_id_in_filename', action='store_true', default=False, help='')
parser.add_argument('--test_names_path', type=str, default='data/splits/timesplit_test', help='Path to text file with the folder names in the test set')
parser.add_argument('--no_overlap_names_path', type=str, default='data/splits/timesplit_test_no_rec_overlap', help='Path text file with the folder names in the test set that have no receptor overlap with the train set')
args = parser.parse_args()
print('Reading paths and names.')
names = read_strings_from_txt(args.test_names_path)
names_no_rec_overlap = read_strings_from_txt(args.no_overlap_names_path)
results_path_containments = os.listdir(args.results_path)
all_times = []
successful_names_list = []
rmsds_list = []
centroid_distances_list = []
min_cross_distances_list = []
min_self_distances_list = []
without_rec_overlap_list = []
start_time = time.time()
for i, name in enumerate(tqdm(names)):
mol = read_mol(args.data_dir, name, remove_hs=True)
mol = Chem.RemoveAllHs(mol)
orig_ligand_pos = np.array(mol.GetConformer().GetPositions())
if args.all_dirs_in_results:
directory_with_name_list = [directory for directory in results_path_containments if name in directory]
if directory_with_name_list == []:
print('Did not find a directory for ', name, '. We are skipping that complex')
continue
else:
directory_with_name = directory_with_name_list[0]
ligand_pos = []
debug_paths = []
for i in range(args.num_predictions):
file_paths = sorted(os.listdir(os.path.join(args.results_path, directory_with_name)))
if args.file_to_exclude is not None:
file_paths = [path for path in file_paths if not args.file_to_exclude in path]
file_path = [path for path in file_paths if f'rank{i+1}_' in path][0]
mol_pred = read_molecule(os.path.join(args.results_path, directory_with_name, file_path),remove_hs=True, sanitize=True)
mol_pred = Chem.RemoveAllHs(mol_pred)
ligand_pos.append(mol_pred.GetConformer().GetPositions())
debug_paths.append(file_path)
ligand_pos = np.asarray(ligand_pos)
else:
if not os.path.exists(os.path.join(args.results_path, name, f'{"" if args.no_id_in_filename else name}{args.file_suffix}')): raise Exception('path did not exists:', os.path.join(args.results_path, name, f'{"" if args.no_id_in_filename else name}{args.file_suffix}'))
mol_pred = read_molecule(os.path.join(args.results_path, name, f'{"" if args.no_id_in_filename else name}{args.file_suffix}'), remove_hs=True, sanitize=True)
if mol_pred == None:
print("Skipping ", name, ' because RDKIT could not read it.')
continue
mol_pred = Chem.RemoveAllHs(mol_pred)
ligand_pos = np.asarray([np.array(mol_pred.GetConformer(i).GetPositions()) for i in range(args.num_predictions)])
try:
rmsd = get_symmetry_rmsd(mol, orig_ligand_pos, [l for l in ligand_pos], mol_pred)
except Exception as e:
print("Using non corrected RMSD because of the error:", e)
rmsd = np.sqrt(((ligand_pos - orig_ligand_pos) ** 2).sum(axis=2).mean(axis=1))
rmsds_list.append(rmsd)
centroid_distances_list.append(np.linalg.norm(ligand_pos.mean(axis=1) - orig_ligand_pos[None,:].mean(axis=1), axis=1))
rec_path = os.path.join(args.data_dir, name, f'{name}_protein_processed.pdb')
if not os.path.exists(rec_path):
rec_path = os.path.join(args.data_dir, name,f'{name}_protein_obabel_reduce.pdb')
rec = PandasPdb().read_pdb(rec_path)
rec_df = rec.df['ATOM']
receptor_pos = rec_df[['x_coord', 'y_coord', 'z_coord']].to_numpy().squeeze().astype(np.float32)
receptor_pos = np.tile(receptor_pos, (args.num_predictions, 1, 1))
cross_distances = np.linalg.norm(receptor_pos[:, :, None, :] - ligand_pos[:, None, :, :], axis=-1)
self_distances = np.linalg.norm(ligand_pos[:, :, None, :] - ligand_pos[:, None, :, :], axis=-1)
self_distances = np.where(np.eye(self_distances.shape[2]), np.inf, self_distances)
min_cross_distances_list.append(np.min(cross_distances, axis=(1,2)))
min_self_distances_list.append(np.min(self_distances, axis=(1, 2)))
successful_names_list.append(name)
without_rec_overlap_list.append(1 if name in names_no_rec_overlap else 0)
performance_metrics = {}
for overlap in ['', 'no_overlap_']:
if 'no_overlap_' == overlap:
without_rec_overlap = np.array(without_rec_overlap_list, dtype=bool)
rmsds = np.array(rmsds_list)[without_rec_overlap]
centroid_distances = np.array(centroid_distances_list)[without_rec_overlap]
min_cross_distances = np.array(min_cross_distances_list)[without_rec_overlap]
min_self_distances = np.array(min_self_distances_list)[without_rec_overlap]
successful_names = np.array(successful_names_list)[without_rec_overlap]
else:
rmsds = np.array(rmsds_list)
centroid_distances = np.array(centroid_distances_list)
min_cross_distances = np.array(min_cross_distances_list)
min_self_distances = np.array(min_self_distances_list)
successful_names = np.array(successful_names_list)
np.save(os.path.join(args.results_path, f'{overlap}rmsds.npy'), rmsds)
np.save(os.path.join(args.results_path, f'{overlap}names.npy'), successful_names)
np.save(os.path.join(args.results_path, f'{overlap}min_cross_distances.npy'), np.array(min_cross_distances))
np.save(os.path.join(args.results_path, f'{overlap}min_self_distances.npy'), np.array(min_self_distances))
performance_metrics.update({
f'{overlap}steric_clash_fraction': (100 * (min_cross_distances < 0.4).sum() / len(min_cross_distances) / args.num_predictions).__round__(2),
f'{overlap}self_intersect_fraction': (100 * (min_self_distances < 0.4).sum() / len(min_self_distances) / args.num_predictions).__round__(2),
f'{overlap}mean_rmsd': rmsds[:,0].mean(),
f'{overlap}rmsds_below_2': (100 * (rmsds[:,0] < 2).sum() / len(rmsds[:,0])),
f'{overlap}rmsds_below_5': (100 * (rmsds[:,0] < 5).sum() / len(rmsds[:,0])),
f'{overlap}rmsds_percentile_25': np.percentile(rmsds[:,0], 25).round(2),
f'{overlap}rmsds_percentile_50': np.percentile(rmsds[:,0], 50).round(2),
f'{overlap}rmsds_percentile_75': np.percentile(rmsds[:,0], 75).round(2),
f'{overlap}mean_centroid': centroid_distances[:,0].mean().__round__(2),
f'{overlap}centroid_below_2': (100 * (centroid_distances[:,0] < 2).sum() / len(centroid_distances[:,0])).__round__(2),
f'{overlap}centroid_below_5': (100 * (centroid_distances[:,0] < 5).sum() / len(centroid_distances[:,0])).__round__(2),
f'{overlap}centroid_percentile_25': np.percentile(centroid_distances[:,0], 25).round(2),
f'{overlap}centroid_percentile_50': np.percentile(centroid_distances[:,0], 50).round(2),
f'{overlap}centroid_percentile_75': np.percentile(centroid_distances[:,0], 75).round(2),
})
top5_rmsds = np.min(rmsds[:, :5], axis=1)
top5_centroid_distances = centroid_distances[np.arange(rmsds.shape[0])[:,None],np.argsort(rmsds[:, :5], axis=1)][:,0]
top5_min_cross_distances = min_cross_distances[np.arange(rmsds.shape[0])[:,None],np.argsort(rmsds[:, :5], axis=1)][:,0]
top5_min_self_distances = min_self_distances[np.arange(rmsds.shape[0])[:,None],np.argsort(rmsds[:, :5], axis=1)][:,0]
performance_metrics.update({
f'{overlap}top5_steric_clash_fraction': (100 * (top5_min_cross_distances < 0.4).sum() / len(top5_min_cross_distances)).__round__(2),
f'{overlap}top5_self_intersect_fraction': (100 * (top5_min_self_distances < 0.4).sum() / len(top5_min_self_distances)).__round__(2),
f'{overlap}top5_rmsds_below_2': (100 * (top5_rmsds < 2).sum() / len(top5_rmsds)).__round__(2),
f'{overlap}top5_rmsds_below_5': (100 * (top5_rmsds < 5).sum() / len(top5_rmsds)).__round__(2),
f'{overlap}top5_rmsds_percentile_25': np.percentile(top5_rmsds, 25).round(2),
f'{overlap}top5_rmsds_percentile_50': np.percentile(top5_rmsds, 50).round(2),
f'{overlap}top5_rmsds_percentile_75': np.percentile(top5_rmsds, 75).round(2),
f'{overlap}top5_centroid_below_2': (100 * (top5_centroid_distances < 2).sum() / len(top5_centroid_distances)).__round__(2),
f'{overlap}top5_centroid_below_5': (100 * (top5_centroid_distances < 5).sum() / len(top5_centroid_distances)).__round__(2),
f'{overlap}top5_centroid_percentile_25': np.percentile(top5_centroid_distances, 25).round(2),
f'{overlap}top5_centroid_percentile_50': np.percentile(top5_centroid_distances, 50).round(2),
f'{overlap}top5_centroid_percentile_75': np.percentile(top5_centroid_distances, 75).round(2),
})
top10_rmsds = np.min(rmsds[:, :10], axis=1)
top10_centroid_distances = centroid_distances[np.arange(rmsds.shape[0])[:,None],np.argsort(rmsds[:, :10], axis=1)][:,0]
top10_min_cross_distances = min_cross_distances[np.arange(rmsds.shape[0])[:,None],np.argsort(rmsds[:, :10], axis=1)][:,0]
top10_min_self_distances = min_self_distances[np.arange(rmsds.shape[0])[:,None],np.argsort(rmsds[:, :10], axis=1)][:,0]
performance_metrics.update({
f'{overlap}top10_self_intersect_fraction': (100 * (top10_min_self_distances < 0.4).sum() / len(top10_min_self_distances)).__round__(2),
f'{overlap}top10_steric_clash_fraction': ( 100 * (top10_min_cross_distances < 0.4).sum() / len(top10_min_cross_distances)).__round__(2),
f'{overlap}top10_rmsds_below_2': (100 * (top10_rmsds < 2).sum() / len(top10_rmsds)).__round__(2),
f'{overlap}top10_rmsds_below_5': (100 * (top10_rmsds < 5).sum() / len(top10_rmsds)).__round__(2),
f'{overlap}top10_rmsds_percentile_25': np.percentile(top10_rmsds, 25).round(2),
f'{overlap}top10_rmsds_percentile_50': np.percentile(top10_rmsds, 50).round(2),
f'{overlap}top10_rmsds_percentile_75': np.percentile(top10_rmsds, 75).round(2),
f'{overlap}top10_centroid_below_2': (100 * (top10_centroid_distances < 2).sum() / len(top10_centroid_distances)).__round__(2),
f'{overlap}top10_centroid_below_5': (100 * (top10_centroid_distances < 5).sum() / len(top10_centroid_distances)).__round__(2),
f'{overlap}top10_centroid_percentile_25': np.percentile(top10_centroid_distances, 25).round(2),
f'{overlap}top10_centroid_percentile_50': np.percentile(top10_centroid_distances, 50).round(2),
f'{overlap}top10_centroid_percentile_75': np.percentile(top10_centroid_distances, 75).round(2),
})
for k in performance_metrics:
print(k, performance_metrics[k])