bshor commited on
Commit
f9da277
1 Parent(s): bfb5fa1

add dataset generation script

Browse files
.gitignore ADDED
@@ -0,0 +1,6 @@
 
 
 
 
 
 
 
1
+ # Ignore __pycache__ folders
2
+ __pycache__/
3
+ .idea/
4
+
5
+ # Ignore .DS_Store files
6
+ .DS_Store
README.md CHANGED
@@ -1,11 +1,11 @@
1
  ---
2
- title: Pinder Inference Template
3
- emoji: 👁
4
  colorFrom: indigo
5
  colorTo: indigo
6
  sdk: docker
7
  pinned: false
8
- license: mit
9
  ---
10
 
11
  Check out the configuration reference at https://huggingface.co/docs/hub/spaces-config-reference
 
1
  ---
2
+ title: DockFormerPP
3
+ emoji:
4
  colorFrom: indigo
5
  colorTo: indigo
6
  sdk: docker
7
  pinned: false
8
+ license: apache-2.0
9
  ---
10
 
11
  Check out the configuration reference at https://huggingface.co/docs/hub/spaces-config-reference
env_consts.py CHANGED
@@ -3,7 +3,7 @@ import os
3
  TEST_INPUT_DIR = None
4
  TEST_OUTPUT_DIR = None
5
  THIS_FILE_DIR = os.path.dirname(os.path.abspath(__file__))
6
- CKPT_PATH = os.path.join(THIS_FILE_DIR, "resources", "77-182500_only_weights.ckpt")
7
  RUN_CONFIG_PATH = os.path.join(THIS_FILE_DIR, "resources", "run_config.json")
8
 
9
  OUTPUT_PATH = os.path.join(THIS_FILE_DIR, "predicted_out.pdb")
 
3
  TEST_INPUT_DIR = None
4
  TEST_OUTPUT_DIR = None
5
  THIS_FILE_DIR = os.path.dirname(os.path.abspath(__file__))
6
+ CKPT_PATH = os.path.join(THIS_FILE_DIR, "resources", "only_weights_102-240750.ckpt")
7
  RUN_CONFIG_PATH = os.path.join(THIS_FILE_DIR, "resources", "run_config.json")
8
 
9
  OUTPUT_PATH = os.path.join(THIS_FILE_DIR, "predicted_out.pdb")
prepare_pinder_dataset.py ADDED
@@ -0,0 +1,622 @@
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
 
1
+ import json
2
+ import os
3
+ import shutil
4
+ import random
5
+ import sys
6
+ from typing import List, Tuple, Optional
7
+
8
+ import Bio.PDB
9
+ import Bio.SeqUtils
10
+ import pandas as pd
11
+ import numpy as np
12
+
13
+
14
+
15
+ OUTPUT_FOLDER = "/tmp/output"
16
+ PINDER_ANNOTATIONS = "/tmp/index.parquet"
17
+ GSUTIL_PATH = "/tmp/google-cloud-sdk/bin/gsutil"
18
+
19
+
20
+ MAX_SYSTEMS_FOR_CLUSTER = 2
21
+ MAX_LENGTH = 350
22
+ MAX_TRIES_OF_METHOD = 5
23
+
24
+
25
+ def do_robust_chain_object_renumber(chain: Bio.PDB.Chain.Chain, new_chain_id: str) -> Optional[Bio.PDB.Chain.Chain]:
26
+ all_residues = [res for res in chain.get_residues()
27
+ if "CA" in res and Bio.SeqUtils.seq1(res.get_resname()) not in ("X", "", " ")]
28
+ if not all_residues:
29
+ return None
30
+
31
+ res_and_res_id = [(res, res.get_id()[1]) for res in all_residues]
32
+
33
+ min_res_id = min([i[1] for i in res_and_res_id])
34
+ if min_res_id < 1:
35
+ print("Negative res id", chain, min_res_id)
36
+ factor = -1 * min_res_id + 1
37
+ res_and_res_id = [(res, res_id + factor) for res, res_id in res_and_res_id]
38
+
39
+ res_and_res_id_no_collisions = []
40
+ for res, res_id in res_and_res_id[::-1]:
41
+ if res_and_res_id_no_collisions and res_and_res_id_no_collisions[-1][1] == res_id:
42
+ # there is a collision, usually an insertion residue
43
+ res_and_res_id_no_collisions = [(i, j + 1) for i, j in res_and_res_id_no_collisions]
44
+ res_and_res_id_no_collisions.append((res, res_id))
45
+
46
+ first_res_id = min([i[1] for i in res_and_res_id_no_collisions])
47
+ factor = 1 - first_res_id # start from 1
48
+ new_chain = Bio.PDB.Chain.Chain(new_chain_id)
49
+
50
+ res_and_res_id_no_collisions.sort(key=lambda x: x[1])
51
+
52
+ for res, res_id in res_and_res_id_no_collisions:
53
+ chain.detach_child(res.id)
54
+ res.id = (" ", res_id + factor, " ")
55
+ new_chain.add(res)
56
+
57
+ return new_chain
58
+
59
+
60
+ def robust_renumber_protein(pdb_path: str, output_path: str):
61
+ if pdb_path.endswith(".pdb"):
62
+ pdb_parser = Bio.PDB.PDBParser(QUIET=True)
63
+ pdb_struct = pdb_parser.get_structure("original_pdb", pdb_path)
64
+ elif pdb_path.endswith(".cif"):
65
+ pdb_struct = Bio.PDB.MMCIFParser().get_structure("original_pdb", pdb_path)
66
+ else:
67
+ raise ValueError("Unknown file type", pdb_path)
68
+ assert len(list(pdb_struct)) == 1, "can't extract if more than one model"
69
+ model = next(iter(pdb_struct))
70
+ chains = list(model.get_chains())
71
+ new_model = Bio.PDB.Model.Model(0)
72
+ chain_ids = "ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789"
73
+ for chain, chain_id in zip(chains, chain_ids):
74
+ new_chain = do_robust_chain_object_renumber(chain, chain_id)
75
+ if new_chain is None:
76
+ continue
77
+ new_model.add(new_chain)
78
+ new_struct = Bio.PDB.Structure.Structure("renumbered_pdb")
79
+ new_struct.add(new_model)
80
+ io = Bio.PDB.PDBIO()
81
+ io.set_structure(new_struct)
82
+ io.save(output_path)
83
+
84
+
85
+ def get_chain_object_to_seq(chain: Bio.PDB.Chain.Chain) -> str:
86
+ res_id_to_res = {res.get_id()[1]: res for res in chain.get_residues() if "CA" in res}
87
+
88
+ if len(res_id_to_res) == 0:
89
+ print("skipping empty chain", chain.get_id())
90
+ return ""
91
+ seq = ""
92
+ for i in range(1, max(res_id_to_res) + 1):
93
+ if i in res_id_to_res:
94
+ seq += Bio.SeqUtils.seq1(res_id_to_res[i].get_resname())
95
+ else:
96
+ seq += "X"
97
+ return seq
98
+
99
+
100
+ def get_sequence_from_pdb(pdb_path: str) -> Tuple[str, List[int]]:
101
+ pdb_parser = Bio.PDB.PDBParser(QUIET=True)
102
+ pdb_struct = pdb_parser.get_structure("original_pdb", pdb_path)
103
+ # chain_to_seq = {chain.id: get_chain_object_to_seq(chain) for chain in pdb_struct.get_chains()}
104
+ all_chain_seqs = [get_chain_object_to_seq(chain) for chain in pdb_struct.get_chains()]
105
+ chain_lengths = [len(seq) for seq in all_chain_seqs]
106
+ return ("X" * 20).join(all_chain_seqs), chain_lengths
107
+
108
+
109
+ from Bio import PDB
110
+ from Bio import pairwise2
111
+
112
+
113
+ def extract_sequence(chain):
114
+ seq = ''
115
+ residues = []
116
+ for res in chain.get_residues():
117
+ seq_res = Bio.SeqUtils.seq1(res.get_resname())
118
+ if seq_res in ('X', "", " "):
119
+ continue
120
+ seq += seq_res
121
+ residues.append(res)
122
+ return seq, residues
123
+
124
+
125
+ def map_residues(alignment, residues_gt, residues_pred):
126
+ idx_gt = 0
127
+ idx_pred = 0
128
+ mapping = []
129
+ for i in range(len(alignment.seqA)):
130
+ aa_gt = alignment.seqA[i]
131
+ aa_pred = alignment.seqB[i]
132
+ res_gt = None
133
+ res_pred = None
134
+ if aa_gt != '-':
135
+ res_gt = residues_gt[idx_gt]
136
+ idx_gt += 1
137
+ if aa_pred != '-':
138
+ res_pred = residues_pred[idx_pred]
139
+ idx_pred += 1
140
+ if res_gt and res_pred:
141
+ mapping.append((res_gt, res_pred))
142
+ return mapping
143
+
144
+
145
+ class ResidueSelect(PDB.Select):
146
+ def __init__(self, residues_to_select):
147
+ self.residues_to_select = set(residues_to_select)
148
+
149
+ def accept_residue(self, residue):
150
+ return residue in self.residues_to_select
151
+
152
+
153
+ def count_gapped_single_aa(alignment):
154
+ count_non_gap = 0
155
+ count_fully_gapped = 0
156
+ for i in range(1, len(alignment.seqA) - 1):
157
+ if alignment.seqA[i] != '-':
158
+ count_non_gap += 1
159
+ if alignment.seqA[i - 1] == '-' and alignment.seqA[i + 1] == '-':
160
+ count_fully_gapped += 1
161
+ top_ratio = count_fully_gapped / count_non_gap
162
+
163
+ count_non_gap = 0
164
+ count_fully_gapped = 0
165
+ for i in range(1, len(alignment.seqB) - 1):
166
+ if alignment.seqA[i] != '-':
167
+ count_non_gap += 1
168
+ if alignment.seqA[i - 1] == '-' and alignment.seqA[i + 1] == '-':
169
+ count_fully_gapped += 1
170
+
171
+ if count_fully_gapped / count_non_gap > top_ratio:
172
+ top_ratio = count_fully_gapped / count_non_gap
173
+
174
+ return top_ratio
175
+
176
+
177
+ def copy_residue_numbering(gt_pdb_path, input_pdb_path):
178
+ parser = PDB.PDBParser(QUIET=True)
179
+ gt_structure = parser.get_structure('gt', gt_pdb_path)
180
+ input_structure = parser.get_structure('input', input_pdb_path)
181
+
182
+ for res in list(input_structure.get_residues()):
183
+ res.id = (' ', res.get_id()[1] + 10000, ' ')
184
+
185
+ for gt_res, input_res in zip(gt_structure.get_residues(), input_structure.get_residues()):
186
+ input_res.id = gt_res.id
187
+
188
+ io = PDB.PDBIO()
189
+ io.set_structure(input_structure)
190
+ io.save(input_pdb_path)
191
+
192
+
193
+ def align_gt_and_input(gt_pdb_path, input_pdb_path, output_gt_path, output_input_path):
194
+ # print("aligning", gt_pdb_path, input_pdb_path, output_gt_path, output_input_path)
195
+ parser = PDB.PDBParser(QUIET=True)
196
+ gt_structure = parser.get_structure('gt', gt_pdb_path)
197
+ pred_structure = parser.get_structure('pred', input_pdb_path)
198
+ matched_residues_gt = []
199
+ matched_residues_pred = []
200
+
201
+ total_gt_size = len([res for res in gt_structure.get_residues() if "CA" in res])
202
+
203
+ used_chain_pred = []
204
+ total_mapping_size = 0
205
+ for chain_gt in gt_structure.get_chains():
206
+ seq_gt, residues_gt = extract_sequence(chain_gt)
207
+ best_alignment = None
208
+ best_chain_pred = None
209
+ best_score = -1
210
+ best_residues_pred = None
211
+ # Find the best matching chain in pred
212
+ for chain_pred in pred_structure.get_chains():
213
+ # print("checking", chain_pred.get_id(), chain_gt.get_id())
214
+ if chain_pred in used_chain_pred:
215
+ continue
216
+ seq_pred, residues_pred = extract_sequence(chain_pred)
217
+ # print(seq_gt)
218
+ # print(seq_pred)
219
+ # alignments = pairwise2.align.globalxx(seq_gt, seq_pred, one_alignment_only=True)
220
+ alignments = pairwise2.align.globalms(seq_gt, seq_pred, 2, -10000, -1, 0, one_alignment_only=True)
221
+ if not alignments:
222
+ continue
223
+ # print("checking2", chain_pred.get_id(), chain_gt.get_id())
224
+
225
+ alignment = alignments[0]
226
+ score = alignment.score
227
+ if score > best_score:
228
+ best_score = score
229
+ best_alignment = alignment
230
+ best_chain_pred = chain_pred
231
+ best_residues_pred = residues_pred
232
+ if best_alignment and count_gapped_single_aa(best_alignment) < 0.2:
233
+ mapping = map_residues(best_alignment, residues_gt, best_residues_pred)
234
+ total_mapping_size += len(mapping)
235
+ used_chain_pred.append(best_chain_pred)
236
+ for res_gt, res_pred in mapping:
237
+ matched_residues_gt.append(res_gt)
238
+ matched_residues_pred.append(res_pred)
239
+ else:
240
+ print(f"No matching chain found for chain {chain_gt.get_id()}")
241
+ assert total_mapping_size / total_gt_size > 0.8, \
242
+ f"Mapping size too low ({total_mapping_size}/{total_gt_size}), skipping"
243
+ print(f"Total mapping size: {total_mapping_size}")
244
+
245
+ # Write new PDB files with only matched residues
246
+ io = PDB.PDBIO()
247
+ io.set_structure(gt_structure)
248
+ io.save(output_gt_path, ResidueSelect(matched_residues_gt))
249
+ io = PDB.PDBIO()
250
+ io.set_structure(pred_structure)
251
+ io.save(output_input_path, ResidueSelect(matched_residues_pred))
252
+
253
+ copy_residue_numbering(output_gt_path, output_input_path)
254
+
255
+
256
+ def validate_matching_input_gt(gt_pdb_path, input_pdb_path):
257
+ gt_residues = [res for res in PDB.PDBParser().get_structure('gt', gt_pdb_path).get_residues()]
258
+ input_residues = [res for res in PDB.PDBParser().get_structure('input', input_pdb_path).get_residues()]
259
+
260
+ if len(gt_residues) != len(input_residues):
261
+ print(f"Residue count mismatch: {len(gt_residues)} vs {len(input_residues)}")
262
+ return -1
263
+
264
+ for res_gt, res_input in zip(gt_residues, input_residues):
265
+ if res_gt.get_resname() != res_input.get_resname():
266
+ print(f"Residue name mismatch: {res_gt.get_resname()} vs {res_input.get_resname()}")
267
+ return -1
268
+ return len(input_residues)
269
+
270
+
271
+ def download_pdb(pdb_name, output_folder):
272
+ output_path = os.path.join(output_folder, pdb_name)
273
+ if os.path.exists(output_path):
274
+ return output_path
275
+ print("downloading", pdb_name)
276
+ os.system(f'{GSUTIL_PATH} -m -q cp "gs://pinder/2024-02/pdbs/{pdb_name}" {output_path}')
277
+ return output_path
278
+
279
+
280
+ INTERFACE_MIN_ATOM_DIST = 5
281
+
282
+
283
+ def get_filtered_res(gt_r_res, gt_l_res, max_length: int):
284
+ gt_r_ca = np.array([res["CA"].coord for res in gt_r_res])
285
+ gt_l_ca = np.array([res["CA"].coord for res in gt_l_res])
286
+
287
+ if len(gt_r_res) + len(gt_l_res) < max_length:
288
+ # continue without cropping
289
+ print("no cropping needed", len(gt_r_res), len(gt_l_res))
290
+ return gt_r_res, gt_l_res
291
+
292
+ # close_residues = np.argwhere(scipy.spatial.distance.cdist(gt_r_ca, gt_l_ca) < INTERFACE_MIN_ATOM_DIST)
293
+ # gt_r_interface, gt_l_interface = set(), set()
294
+ # for i, j in close_residues:
295
+ # gt_r_interface.add(gt_r_res[i].id[1])
296
+ # gt_l_interface.add(gt_l_res[j].id[1])
297
+
298
+ inter_dists = gt_r_ca[:, np.newaxis, :] - gt_l_ca[np.newaxis, :, :]
299
+ inter_dists = np.sqrt((inter_dists ** 2).sum(-1))
300
+ min_inter_dist_per_gt_l_res = inter_dists.min(axis=0)
301
+ min_inter_dist_per_gt_r_res = inter_dists.min(axis=1)
302
+
303
+ assert min_inter_dist_per_gt_l_res.shape[0] == len(gt_l_res)
304
+ assert min_inter_dist_per_gt_r_res.shape[0] == len(gt_r_res)
305
+
306
+ min_r_res, max_r_res = min(min_inter_dist_per_gt_r_res), max(min_inter_dist_per_gt_r_res)
307
+ min_l_res, max_l_res = min(min_inter_dist_per_gt_l_res), max(min_inter_dist_per_gt_l_res)
308
+
309
+ r_pocket = [res for res in gt_r_res if min_r_res <= res.id[1] <= max_r_res]
310
+ l_pocket = [res for res in gt_l_res if min_l_res <= res.id[1] <= max_l_res]
311
+
312
+ if len(r_pocket) + len(l_pocket) < max_length:
313
+ # add extra residues to both chains to get a total of max_length
314
+ res_r_before = [res for res in gt_r_res if res.id[1] < min_r_res]
315
+ res_r_after = [res for res in gt_r_res if res.id[1] > max_r_res]
316
+ res_l_before = [res for res in gt_l_res if res.id[1] < min_l_res]
317
+ res_l_after = [res for res in gt_l_res if res.id[1] > max_l_res]
318
+
319
+ extra_to_add = max_length - len(r_pocket) - len(l_pocket)
320
+
321
+ actions = []
322
+ if len(res_r_before) > 0:
323
+ actions.append("add_r_before")
324
+ if len(res_r_after) > 0:
325
+ actions.append("add_r_after")
326
+ if len(res_l_before) > 0:
327
+ actions.append("add_l_before")
328
+ if len(res_l_after) > 0:
329
+ actions.append("add_l_after")
330
+ while extra_to_add > 0 and actions:
331
+ action = random.choice(actions)
332
+
333
+ if action == "add_r_before":
334
+ r_pocket.insert(0, res_r_before.pop())
335
+ if not len(res_r_before):
336
+ actions.remove("add_r_before")
337
+ elif action == "add_r_after":
338
+ r_pocket.append(res_r_after.pop())
339
+ if not len(res_r_after):
340
+ actions.remove("add_r_after")
341
+ elif action == "add_l_before":
342
+ l_pocket.insert(0, res_l_before.pop())
343
+ if not len(res_l_before):
344
+ actions.remove("add_l_before")
345
+ elif action == "add_l_after":
346
+ l_pocket.append(res_l_after.pop())
347
+ if not len(res_l_after):
348
+ actions.remove("add_l_after")
349
+ extra_to_add -= 1
350
+ print("Extended pocket sizes", len(r_pocket), len(l_pocket), "extra_to_add", extra_to_add)
351
+ return r_pocket, l_pocket
352
+
353
+ print("cropping simply")
354
+ # remove residues that are farthest from the interface
355
+ res_and_dist_r = [(res, min_inter_dist_per_gt_r_res[res_idx]) for res_idx, res in enumerate(gt_r_res)]
356
+ res_and_dist_l = [(res, min_inter_dist_per_gt_l_res[res_idx]) for res_idx, res in enumerate(gt_l_res)]
357
+
358
+ res_and_dist_r = [(res, dist) for res, dist in res_and_dist_r if res in r_pocket]
359
+ res_and_dist_l = [(res, dist) for res, dist in res_and_dist_l if res in l_pocket]
360
+
361
+ res_and_dist_r = sorted(res_and_dist_r, key=lambda x: x[1], reverse=True)
362
+ res_and_dist_l = sorted(res_and_dist_l, key=lambda x: x[1], reverse=True)
363
+
364
+ while len(res_and_dist_r) + len(res_and_dist_l) > max_length:
365
+ if res_and_dist_r[0][1] > res_and_dist_l[0][1]:
366
+ res_and_dist_r.pop(0)
367
+ else:
368
+ res_and_dist_l.pop(0)
369
+
370
+ return [res for res, _ in res_and_dist_r], [res for res, _ in res_and_dist_l]
371
+
372
+
373
+ def prepare_holo(row, tmp_dir_path, max_length: int):
374
+ tmp_gt_r_pdb = os.path.join(tmp_dir_path, f"tmp_{row.id}_gt_r.pdb")
375
+ tmp_gt_l_pdb = os.path.join(tmp_dir_path, f"tmp_{row.id}_gt_l.pdb")
376
+
377
+ if os.path.exists(tmp_gt_r_pdb) and os.path.exists(tmp_gt_l_pdb):
378
+ return tmp_gt_r_pdb, tmp_gt_l_pdb
379
+
380
+ holo_r_pdb = download_pdb(row.holo_R_pdb, tmp_dir_path)
381
+ holo_l_pdb = download_pdb(row.holo_L_pdb, tmp_dir_path)
382
+
383
+ # make gt and apo that matches
384
+ robust_renumber_protein(holo_r_pdb, tmp_gt_r_pdb)
385
+ robust_renumber_protein(holo_l_pdb, tmp_gt_l_pdb)
386
+
387
+ parser = PDB.PDBParser(QUIET=True)
388
+ gt_r_prot = parser.get_structure('r', tmp_gt_r_pdb)
389
+ gt_l_prot = parser.get_structure('l', tmp_gt_l_pdb)
390
+
391
+ assert len(list(gt_r_prot.get_chains())) == 1, "can't extract if more than one chain"
392
+ assert len(list(gt_l_prot.get_chains())) == 1, "can't extract if more than one chain"
393
+
394
+ gt_r_res = [res for res in gt_r_prot.get_residues() if "CA" in res]
395
+ gt_l_res = [res for res in gt_l_prot.get_residues() if "CA" in res]
396
+
397
+ to_keep_r, to_keep_l = get_filtered_res(gt_r_res, gt_l_res, max_length)
398
+
399
+ io = PDB.PDBIO()
400
+ io.set_structure(gt_r_prot)
401
+ io.save(tmp_gt_r_pdb, ResidueSelect(to_keep_r))
402
+ io = PDB.PDBIO()
403
+ io.set_structure(gt_l_prot)
404
+ io.save(tmp_gt_l_pdb, ResidueSelect(to_keep_l))
405
+
406
+ return tmp_gt_r_pdb, tmp_gt_l_pdb
407
+
408
+
409
+ def generate_input_pdbs(tmp_input_r_pdb, tmp_input_l_pdb, tmp_gt_r_pdb, tmp_gt_l_pdb,
410
+ input_r_output_pdb, input_l_output_pdb, gt_r_output_pdb, gt_l_output_pdb):
411
+ # print("preparing input pdbs", gt_r_output_pdb)
412
+ if not os.path.exists(tmp_input_r_pdb) or not os.path.exists(tmp_input_l_pdb):
413
+ raise False
414
+
415
+ try:
416
+ align_gt_and_input(tmp_gt_r_pdb, tmp_input_r_pdb, gt_r_output_pdb, input_r_output_pdb)
417
+ protein_size_r = validate_matching_input_gt(gt_r_output_pdb, input_r_output_pdb)
418
+ assert protein_size_r > -1, "Failed to validate matching input and gt"
419
+
420
+ align_gt_and_input(tmp_gt_l_pdb, tmp_input_l_pdb, gt_l_output_pdb, input_l_output_pdb)
421
+ protein_size_l = validate_matching_input_gt(gt_l_output_pdb, input_l_output_pdb)
422
+ assert protein_size_l > -1, "Failed to validate matching input and gt"
423
+ except Exception as e:
424
+ print("Failed to align", e)
425
+ if os.path.exists(gt_r_output_pdb):
426
+ os.remove(gt_r_output_pdb)
427
+ if os.path.exists(gt_l_output_pdb):
428
+ os.remove(gt_l_output_pdb)
429
+ if os.path.exists(input_r_output_pdb):
430
+ os.remove(input_r_output_pdb)
431
+ if os.path.exists(input_l_output_pdb):
432
+ os.remove(input_l_output_pdb)
433
+ return False
434
+
435
+ return True
436
+
437
+
438
+ def _get_rel_path(abs_path):
439
+ return os.path.join(os.path.basename(os.path.dirname(abs_path)), os.path.basename(abs_path))
440
+
441
+
442
+ def main(start_ind: Optional[int] = None, end_ind: Optional[int] = None):
443
+ print("running with", start_ind, end_ind)
444
+
445
+ os.makedirs(OUTPUT_FOLDER, exist_ok=True)
446
+ output_models_folder = os.path.join(OUTPUT_FOLDER, "pinder_models")
447
+ output_train_jsons_folder = os.path.join(OUTPUT_FOLDER, "pinder_jsons_train")
448
+ output_val_jsons_folder = os.path.join(OUTPUT_FOLDER, "pinder_jsons_val")
449
+ output_test_jsons_folder = os.path.join(OUTPUT_FOLDER, "pinder_jsons_test")
450
+ output_info = os.path.join(OUTPUT_FOLDER, "pinder_generation_info.csv")
451
+
452
+ os.makedirs(output_models_folder, exist_ok=True)
453
+ os.makedirs(output_train_jsons_folder, exist_ok=True)
454
+ os.makedirs(output_val_jsons_folder, exist_ok=True)
455
+ os.makedirs(output_test_jsons_folder, exist_ok=True)
456
+
457
+ split_to_folder = {
458
+ "train": output_train_jsons_folder,
459
+ "val": output_val_jsons_folder,
460
+ "test": output_test_jsons_folder
461
+ }
462
+
463
+ # output_info_file = open(output_info, "a+")
464
+ systems = pd.read_parquet(PINDER_ANNOTATIONS)
465
+ systems = systems[systems.split.isin(['train', 'val', 'test'])]
466
+
467
+ cluster_ids = systems["cluster_id"].value_counts()
468
+ cluster_ids = cluster_ids[cluster_ids >= 1]
469
+ print("There are", len(cluster_ids), "clusters")
470
+
471
+ # clusters_with_data = 0
472
+ # for cluster_id in cluster_ids.index:
473
+ # cluster_systems = systems[systems["cluster_id"] == cluster_id]
474
+ # with_apo = cluster_systems[cluster_systems.apo_R & cluster_systems.apo_L]
475
+ # if len(with_apo) > 0:
476
+ # print("Cluster", cluster_id, "has", len(with_apo), "systems with apo")
477
+ # clusters_with_data += 1
478
+ # continue
479
+ # with_pred = cluster_systems[cluster_systems.predicted_R & cluster_systems.predicted_L]
480
+ # if len(with_pred) > 0:
481
+ # print("Cluster", cluster_id, "has", len(with_pred), "systems with pred")
482
+ # clusters_with_data += 1
483
+ # continue
484
+ # print("There are", clusters_with_data, "clusters with data out of", len(cluster_ids))
485
+
486
+ for cluster_ind, cluster_id in enumerate(sorted(cluster_ids.index)):
487
+ if (start_ind is not None and cluster_ind < start_ind) or (end_ind is not None and cluster_ind >= end_ind):
488
+ continue
489
+ # if cluster_id != "cluster_10004_p":
490
+ # continue
491
+ tmp_dir_path = os.path.join(OUTPUT_FOLDER, "tmp_" + cluster_id)
492
+ os.makedirs(tmp_dir_path, exist_ok=True)
493
+ system_id_to_method = {}
494
+
495
+ cluster_systems = systems[systems["cluster_id"] == cluster_id]
496
+ print("--- Starting cluster", cluster_ind, cluster_id, "size", cluster_systems.shape)
497
+
498
+ with_apo = cluster_systems[cluster_systems.apo_R & cluster_systems.apo_L]
499
+ print("*** APO *** Cluster", cluster_id, "has", len(with_apo), "systems with apo")
500
+ for try_counter, row in enumerate(with_apo.itertuples()):
501
+ if row.split not in ("test", "val") \
502
+ and (try_counter >= MAX_TRIES_OF_METHOD or len(system_id_to_method) >= MAX_SYSTEMS_FOR_CLUSTER):
503
+ continue
504
+ print("-- Trying to prepare apo", row.id, row.split)
505
+ try:
506
+ tmp_gt_r_pdb, tmp_gt_l_pdb = prepare_holo(row, tmp_dir_path, MAX_LENGTH)
507
+
508
+ gt_r_output_path = os.path.join(output_models_folder, f"{row.id}_gt_r.pdb")
509
+ gt_l_output_path = os.path.join(output_models_folder, f"{row.id}_gt_l.pdb")
510
+
511
+ input_r_output_path = os.path.join(output_models_folder, f"{row.id}_input_r.pdb")
512
+ input_l_output_path = os.path.join(output_models_folder, f"{row.id}_input_l.pdb")
513
+
514
+ input_r_pdb_path = download_pdb(row.apo_R_pdb, tmp_dir_path)
515
+ input_l_pdb_path = download_pdb(row.apo_L_pdb, tmp_dir_path)
516
+
517
+ if generate_input_pdbs(input_r_pdb_path, input_l_pdb_path, tmp_gt_r_pdb, tmp_gt_l_pdb,
518
+ input_r_output_path, input_l_output_path, gt_r_output_path, gt_l_output_path):
519
+ system_id_to_method[row.id] = "apo"
520
+
521
+ except Exception as e:
522
+ print("Failed to prepare apo", row.id, e)
523
+ continue
524
+
525
+ with_pred = cluster_systems[cluster_systems.predicted_R & cluster_systems.predicted_L]
526
+ print("*** Pred *** Cluster", cluster_id, "has", len(with_pred), "systems with pred")
527
+ for try_counter, row in enumerate(with_pred.itertuples()):
528
+ if row.id in system_id_to_method:
529
+ continue
530
+ if row.split not in ("test", "val") \
531
+ and (try_counter >= MAX_TRIES_OF_METHOD or len(system_id_to_method) >= MAX_SYSTEMS_FOR_CLUSTER):
532
+ continue
533
+ print("-- Trying to prepare pred", row.id, row.split)
534
+ try:
535
+ tmp_gt_r_pdb, tmp_gt_l_pdb = prepare_holo(row, tmp_dir_path, MAX_LENGTH)
536
+
537
+ gt_r_output_path = os.path.join(output_models_folder, f"{row.id}_gt_r.pdb")
538
+ gt_l_output_path = os.path.join(output_models_folder, f"{row.id}_gt_l.pdb")
539
+
540
+ input_r_output_path = os.path.join(output_models_folder, f"{row.id}_input_r.pdb")
541
+ input_l_output_path = os.path.join(output_models_folder, f"{row.id}_input_l.pdb")
542
+
543
+ input_r_pdb_path = download_pdb(row.predicted_R_pdb, tmp_dir_path)
544
+ input_l_pdb_path = download_pdb(row.predicted_L_pdb, tmp_dir_path)
545
+
546
+ if generate_input_pdbs(input_r_pdb_path, input_l_pdb_path, tmp_gt_r_pdb, tmp_gt_l_pdb,
547
+ input_r_output_path, input_l_output_path, gt_r_output_path, gt_l_output_path):
548
+ system_id_to_method[row.id] = "pred"
549
+
550
+ except Exception as e:
551
+ print("Failed to prepare pred", row.id, e)
552
+
553
+ # default - use gt
554
+ print("*** GT *** ")
555
+ for row in cluster_systems.itertuples():
556
+ if row.id in system_id_to_method:
557
+ continue
558
+ if row.split not in ("test", "val") and len(system_id_to_method) >= MAX_SYSTEMS_FOR_CLUSTER:
559
+ continue
560
+ try:
561
+ tmp_gt_r_pdb, tmp_gt_l_pdb = prepare_holo(row, tmp_dir_path, MAX_LENGTH)
562
+
563
+ gt_r_output_path = os.path.join(output_models_folder, f"{row.id}_gt_r.pdb")
564
+ gt_l_output_path = os.path.join(output_models_folder, f"{row.id}_gt_l.pdb")
565
+
566
+ input_r_output_path = os.path.join(output_models_folder, f"{row.id}_input_r.pdb")
567
+ input_l_output_path = os.path.join(output_models_folder, f"{row.id}_input_l.pdb")
568
+
569
+ shutil.copyfile(tmp_gt_r_pdb, gt_r_output_path)
570
+ shutil.copyfile(tmp_gt_r_pdb, input_r_output_path)
571
+
572
+ shutil.copyfile(tmp_gt_l_pdb, gt_l_output_path)
573
+ shutil.copyfile(tmp_gt_l_pdb, input_l_output_path)
574
+
575
+ system_id_to_method[row.id] = "gt"
576
+
577
+ except Exception as e:
578
+ print("Failed to prepare gt", row.id, e)
579
+
580
+ # save jsons
581
+ for row in cluster_systems.itertuples():
582
+ if row.id not in system_id_to_method:
583
+ continue
584
+
585
+ output_json_path = os.path.join(split_to_folder[row.split], f"{row.id}.json")
586
+
587
+ gt_r_output_path = os.path.join(output_models_folder, f"{row.id}_gt_r.pdb")
588
+ gt_l_output_path = os.path.join(output_models_folder, f"{row.id}_gt_l.pdb")
589
+
590
+ input_r_output_path = os.path.join(output_models_folder, f"{row.id}_input_r.pdb")
591
+ input_l_output_path = os.path.join(output_models_folder, f"{row.id}_input_l.pdb")
592
+
593
+ protein_r_seq_len = validate_matching_input_gt(gt_r_output_path, input_r_output_path)
594
+ protein_l_seq_len = validate_matching_input_gt(gt_l_output_path, input_l_output_path)
595
+
596
+ json_data = {
597
+ "input_r_structure": _get_rel_path(input_r_output_path),
598
+ "input_l_structure": _get_rel_path(input_l_output_path),
599
+ "gt_r_structure": _get_rel_path(gt_r_output_path),
600
+ "gt_l_structure": _get_rel_path(gt_l_output_path),
601
+ "resolution": 1.0,
602
+ "protein_r_seq_len": protein_r_seq_len,
603
+ "protein_l_seq_len": protein_l_seq_len,
604
+ "uniprot_r": row.uniprot_R,
605
+ "uniprot_l": row.uniprot_L,
606
+ "cluster": row.cluster_id,
607
+ "input_protein_source": system_id_to_method[row.id],
608
+ "pdb_id": row.id,
609
+ }
610
+ open(output_json_path, "w").write(json.dumps(json_data, indent=4))
611
+
612
+ print("******* saved", row.id, system_id_to_method[row.id], flush=True)
613
+ shutil.rmtree(tmp_dir_path)
614
+
615
+ print("total systems", len(systems))
616
+
617
+
618
+ if __name__ == '__main__':
619
+ if len(sys.argv) == 3:
620
+ main(int(sys.argv[1]), int(sys.argv[2]))
621
+ else:
622
+ main()
resources/{77-182500_only_weights.ckpt → only_weights_102-240750.ckpt} RENAMED
@@ -1,3 +1,3 @@
1
  version https://git-lfs.github.com/spec/v1
2
- oid sha256:3412e505a8fce15011c589bb008f0f9c92d0ce9703bc75d70a085db733c2b16d
3
- size 52579155
 
1
  version https://git-lfs.github.com/spec/v1
2
+ oid sha256:3e32cd27b63f684ed813e65db4369f5af28b1fceb3c81df66bdd4952f6b78853
3
+ size 52579856