OSCR

A Disorder-Aware Computational Framework to Identify Structurally Tractable Targets in Proliferative Vitreoretinopathy.

Code ↔ Paper

4 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 4 matches
  1. [1] § Methods › Evaluation of a Target Protein–SNAIL1 › RFdiffusion-Based Binder Design ↔ rf/examples/diffusion.ipynb, lines 475–514 · score 0.67 · initial_guess, rm aa, ProteinMPNN, soluble, folded, binder
  2. [2] § Methods › Evaluation of a Target Protein–SNAIL1 › RFdiffusion-Based Binder Design ↔ af/examples/RSO.ipynb, lines 257–320 · score 0.60 · initial guess, rm aa, Cysteine, soluble, backbone, sequences
  3. [3] § Methods › Evaluation of a Target Protein–SNAIL1 › RFdiffusion-Based Binder Design ↔ af/examples/RSO.ipynb, lines 767–801 · score 0.55 · multimer model, generated sequences, recycles, RMSD, template, pLDDT
  4. [4] § Results › RFdiffusion-Based Binder Design for SNAIL1 ↔ af/examples/RSO.ipynb, lines 121–231 · score 0.53 · designed sequence, ProteinMPNN sequence, AlphaFold, confidence, pLDDT, score

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Jupyter notebook · 1,098 lines · 38 KB · other · 3 matches

  1. # %% [markdown]
  2. # <a href="https://colab.research.google.com/github/sokrypton/ColabDesign/blob/main/af/examples/RSO.ipynb" target="_parent"><img src="https://colab.research.google.com/assets/colab-badge.svg" alt="Open In Colab"/></a>
  4. # %% [markdown]
  5. #
  6. #
  7. #
  8. # #Protein Design using Relaxed Sequence Optimization
  9. #
  10. #
  11. # **Scalable protein design using optimization in a relaxed sequence space**
  12. #
  13. #
  14. #
  15. #
  16. # Christopher Frank, Ali Khoshouei, Lara Fuß, Lara Weber Dominik Schiewitz,Zhixuan Zhao, Motoyuki Hattori, Yosta de Stigter, Shihao Feng, Sergey Ovchinnikov and Hendrik Dietz
  17. #
  18. #
  19. # This notebook contains code to run relaxed sequence optimisation for de novo protein design as described in the manuscript. There are additional options to modify the pipeline according to ones needs
  20. #
  21. # We recommend using at least an L4 GPU to run this notebook, as the free T4 GPU struggles with larger proteins
  22. #
  23. # Alternativly a local installation of ColabDesign is strongly recommendet, especially for the design of larger proteins.
  24. #
  25. # For questions feel free to reach out to the authors
  27. # %%
  28. #@title setup
  29. %%time
  30. import os
  31. if not os.path.isdir("params"):
  32. # get code
  33. os.system("pip -q install pyppeteer nest_asyncio")
  34. os.system("pip -q install git+https://github.com/sokrypton/ColabDesign.git")
  35. # for debugging
  36. os.system("ln -s /usr/local/lib/python3.*/dist-packages/colabdesign colabdesign")
  37. # download params
  38. os.system("mkdir params")
  39. os.system("apt-get install aria2 -qq")
  40. os.system("aria2c -q -x 16 https://storage.googleapis.com/alphafold/alphafold_params_2022-12-06.tar")
  41. os.system("tar -xf alphafold_params_2022-12-06.tar -C params")
  43. import warnings
  44. warnings.simplefilter(action='ignore', category=FutureWarning)
  46. import os
  47. from colabdesign import mk_afdesign_model, clear_mem
  48. from colabdesign.mpnn import mk_mpnn_model
  50. from IPython.display import HTML
  51. from google.colab import files
  52. import numpy as np
  54. import requests, time
  55. if not os.path.isfile("TMscore"):
  56. os.system("wget -qnc https://zhanggroup.org/TM-score/TMscore.cpp")
  57. os.system("g++ -static -O3 -ffast-math -lm -o TMscore TMscore.cpp")
  58. def tmscore(x,y):
  59. # pass to TMscore
  60. output = os.popen(f'./TMscore {x} {y}')
  61. # parse outputs
  62. parse_float = lambda x: float(x.split("=")[1].split()[0])
  63. o = {}
  64. for line in output:
  65. line = line.rstrip()
  66. if line.startswith("RMSD"): o["rms"] = parse_float(line)
  67. if line.startswith("TM-score"): o["tms"] = parse_float(line)
  68. if line.startswith("GDT-TS-score"): o["gdt"] = parse_float(line)
  69. return o
  71. import asyncio
  72. import nest_asyncio
  73. from pyppeteer import launch
  74. import base64
  76. # Apply nest_asyncio to enable nested event loops
  77. nest_asyncio.apply()
  79. async def fetch_blob_content(page, blob_url):
  80. blob_to_base64 = """
  81. async (blobUrl) => {
  82. const blob = await fetch(blobUrl).then(r => r.blob());
  83. return new Promise((resolve) => {
  84. const reader = new FileReader();
  85. reader.onloadend = () => resolve(reader.result);
  86. reader.readAsDataURL(blob);
  87. });
  88. }
  89. """
  90. base64_data = await page.evaluate(blob_to_base64, blob_url)
  91. _, encoded = base64_data.split(',', 1)
  92. return base64.b64decode(encoded)
  94. async def extract_pdb_file_download_link_and_content(url):
  95. browser = await launch(headless=True, args=['--no-sandbox', '--disable-setuid-sandbox'])
  96. page = await browser.newPage()
  97. await page.goto(url, {'waitUntil': 'networkidle0'})
  98. elements = await page.querySelectorAll('a.btn.bg-purple')
  99. for element in elements:
  100. href = await page.evaluate('(element) => element.getAttribute("href")', element)
  101. if 'blob:https://esmatlas.com/' in href:
  102. content = await fetch_blob_content(page, href)
  103. await browser.close()
  104. return href, content
  105. await browser.close()
  106. return "No PDB file link found.", None
  108. def esmfold_api(sequence):
  109. url = f'https://esmatlas.com/resources/fold/result?fasta_header=%3Eunnamed&sequence={sequence}'
  110. result = asyncio.get_event_loop().run_until_complete(extract_pdb_file_download_link_and_content(url))
  111. if result[1]:
  112. pdb_str = result[1].decode('utf-8')
  113. return pdb_str
  114. else:
  115. return "Failed to retrieve PDB content."
  117. import jax
  118. import jax.numpy as jnp
  119. from colabdesign.af.alphafold.common import residue_constants
  121. # %%
  122. #@title # Unconditional Generation (Custom)
  123. #@markdown For a given length, generate/hallucinate a protein sequence that AlphaFold thinks folds into a well structured protein (high plddt, low pae, many contacts).
  124. LENGTH = 100 #@param {type:"integer"}
  125. #@markdown With copies you can specify the number of identical sequences design, resulting in homo oligomers. Copies = 1 is the standard, resulting in a monomer
  127. COPIES = 1 #@param ["1", "2", "3", "4", "5", "6", "7", "8"] {type:"raw"}
  128. MODE = "manuscript"
  130. #@markdown Select the losses you want to use. For unconditional generation as reported in the manuscript use all the losses. To increase the diversity of designes remove confidence losses and/or increase the weight of the helix loss.
  132. use_rg_loss = True #@param {type:"boolean"}
  133. #@markdown A too strong rg loss can lead to problems and clashes. Use 0.1 for backbones smaller then 600 AA and 0.01 for larger proteins (0.001 for 1000 AA).
  134. rg_weight = 0.1 #@param {type:"raw"}
  135. use_helix_loss = True #@param {type:"boolean"}
  136. use_con_loss = True #@param {type:"boolean"}
  137. use_confidence_loss = True #@param {type:"boolean"}
  138. #@markdown How many halluicnation iteration you want to perform. The standard in the manuscript is 100.
  140. iters = 50 #@param ["100", "50", "30"] {type:"raw"}
  143. #@markdown Select if you want to use the 'standard" ProteinMPNN weights or the soluble ones. The soluble ones usually result in higher in silico as well as experimental sucess, but will increase the negative net charge of the protein which sould potentially interfer with certain protein design problems. The manuscript settings are soluble MPNN
  145. use_solubleMPNN = True #@param {type:"boolean"}
  146. #@markdown Select this to use an experimental ProteinMPNN loss, also backpropagating through ProteinMPNN. This was not used in the manuscript
  148. use_mpnn_loss = False #@param {type:"boolean"}
  149. #@markdown
  151. def add_rg_loss(self, weight=0.1):
  152. '''add radius of gyration loss'''
  153. def loss_fn(inputs, outputs):
  154. xyz = outputs["structure_module"]
  155. ca = xyz["final_atom_positions"][:,residue_constants.atom_order["CA"]]
  156. if self.protocol == "binder":
  157. ca = ca[-self._binder_len:]
  158. #This uses a scaled version of the rg loss, only looking at every 5th residue
  159. if MODE == "manuscript":
  160. ca = ca[::5]
  161. rg = jnp.sqrt(jnp.square(ca - ca.mean(0)).sum(-1).mean() + 1e-8)
  163. if MODE == "original":
  164. rg_th = 2.38 * ca.shape[0] ** 0.365
  165. rg = jax.nn.elu(rg - rg_th)
  166. return {"rg":rg}
  167. self._callbacks["model"]["loss"].append(loss_fn)
  168. self.opt["weights"]["rg"] = weight
  170. def add_mpnn_loss(self, mpnn=0.1, mpnn_seq=0.0):
  171. '''
  172. add mpnn loss
  173. mpnn = maximize confidence of proteinmpnn
  174. mpnn_seq = push designed sequence to match proteinmpnn logits
  175. '''
  177. self._mpnn = mk_mpnn_model(weights = "soluble" if use_solubleMPNN else "original")
  178. def loss_fn(inputs, outputs, aux, key):
  180. # get structure
  181. atom_idx = tuple(residue_constants.atom_order[k] for k in ["N","CA","C","O"])
  182. I = {"S": inputs["aatype"],
  183. "residue_idx": inputs["residue_index"],
  184. "chain_idx": inputs["asym_id"],
  185. "X": outputs["structure_module"]["final_atom_positions"][:,atom_idx],
  186. "mask": outputs["structure_module"]["final_atom_mask"][:,1],
  187. "lengths": self._lengths,
  188. "key": key}
  190. if "offset" in inputs:
  191. I["offset"] = inputs["offset"]
  193. # set autoregressive mask
  194. L = sum(self._lengths)
  195. if self.protocol == "binder":
  196. I["ar_mask"] = 1 - np.eye(L)
  197. I["ar_mask"][-self._len:,-self._len:] = 0
  198. else:
  199. I["ar_mask"] = np.zeros((L,L))
  201. # get logits
  202. logits = self._mpnn._score(**I)["logits"][:,:20]
  203. if self.protocol == "binder":
  204. logits = logits[-self._len:]
  205. else:
  206. logits = logits[:self._len]
  207. aux["mpnn_logits"] = logits
  209. # compute loss
  210. log_q = jax.nn.log_softmax(logits)
  211. p = inputs["seq"]["hard"]
  212. q = jax.nn.softmax(logits)
  213. losses = {}
  214. losses["mpnn"] = -log_q.max(-1).mean()
  215. losses["mpnn_seq"] = -(p * jax.lax.stop_gradient(log_q)).sum(-1).mean()
  216. return losses
  218. self._callbacks["model"]["loss"].append(loss_fn)
  219. self.opt["weights"]["mpnn"] = mpnn
  220. self.opt["weights"]["mpnn_seq"] = mpnn_seq
  222. clear_mem()
  223. af_model = mk_afdesign_model(protocol="hallucination")
  224. af_model.prep_inputs(length=LENGTH, copies=COPIES)
  226. # add extra losses
  228. if use_mpnn_loss: add_mpnn_loss(af_model)
  230. print("length",af_model._lengths)
  231. print("weights",af_model.opt["weights"])
  233. # %%
  234. #This cell runs the design loop. Run this in a for loop for design of multiple proteins
  236. af_model.restart()
  237. af_model.set_seq(mode=["gumbel","soft"])
  238. if use_rg_loss: add_rg_loss(af_model,rg_weight)
  239. if use_helix_loss : af_model.set_weights(helix=-0.2)
  240. if use_con_loss : af_model.set_weights(con=1.0)
  241. if use_confidence_loss : af_model.set_weights(plddt=0.5, pae=0.5)
  242. print("weights",af_model.opt["weights"])
  243. af_model.design_logits(iters-10)
  244. af_model.design_logits(10, save_best=True)
  246. # %%
  247. #This cell plots and saves the results as a pdb file
  248. af_model.save_pdb(f"{af_model.protocol}.pdb")
  249. af_model.plot_pdb()
  251. # %%
  252. HTML(af_model.animate())
  254. # %%
  255. af_model.get_seqs()
  257. # %%
  258. import pandas as pd
  259. #@title # Designability test
  260. #@markdown Test the designability of the backbone, taking in the backbone, generating sequences with solubleMPNN and predicting the sequence with AF2 in single sequence mode.
  261. #@markdown Use Initial Guess (IG) and All Atom Initialisation (AA) for larger proteins
  263. AA = False #@param {type:"boolean"}
  264. IG = False #@param {type:"boolean"}
  265. #@markdown NOTE: we remove cysteines from all designed proteins. Additionally for large proteins we also exclude methions to reduce the number of internal start codons
  267. def designability_test(af_model_test, mpnn_model_test,
  268. num_seqs=8, sampling_temp=0.1, num_recycles=3,
  269. model_num=4, best_metric="rmsd",
  270. in_pdb="init.pdb", out_pdb="final.pdb",
  271. verbose=False):
  272. alphafold_model = f"model_{model_num}_ptm"
  274. af_model_test.prep_inputs(in_pdb)
  275. af_model_test.restart(rm_aa="C,M")
  276. af_model_test._args["best_metric"] = best_metric
  277. L = sum(af_model_test._lengths)
  278. mpnn_model_test.get_af_inputs(af_model_test)
  279. out = mpnn_model_test.sample(num=num_seqs // 8, batch=8,
  280. temperature=sampling_temp)
  282. af_terms = ["plddt", "ptm", "pae", "rmsd", "dgram_cce"]
  283. for k in af_terms: out[k] = []
  285. for n in range(num_seqs):
  286. seq = out["seq"][n]
  287. af_model_test.predict(seq=seq,
  288. num_recycles=num_recycles,
  289. num_models=1,
  290. verbose=False,
  291. models=alphafold_model)
  293. for k in af_terms: out[k].append(af_model_test.aux["log"][k])
  294. out["pae"][-1] = out["pae"][-1] * 31
  295. af_model_test._save_results(save_best=True, verbose=verbose)
  296. af_model_test._k += 1
  298. af_model_test.save_pdb(out_pdb)
  299. labels = ["score"] + af_terms + ["seq"]
  300. data = [[out[k][n] for k in labels] for n in range(num_seqs)]
  301. labels[0] = "mpnn"
  302. df = pd.DataFrame(data, columns=labels)
  303. return df
  305. af_model_test = mk_afdesign_model(protocol="fixbb",best_metric="rmsd",use_initial_guess=IG,use_initial_atom_pos=AA,use_templates=False)
  306. mpnn_model_test = mk_mpnn_model(weights="soluble")
  309. lowest_rmsd = float('inf')
  310. lowest_rmsd_data = None
  312. in_pdb = f"{af_model.protocol}.pdb"
  313. out_pdb = f"{af_model.protocol}_out.pdb"
  316. out = designability_test(af_model_test, mpnn_model_test,
  317. num_seqs=8, sampling_temp=0.1, num_recycles=3,
  318. model_num=4, best_metric="rmsd",
  319. in_pdb=in_pdb, out_pdb=out_pdb,
  320. verbose=True)
  324. # %%
  325. from colabdesign import mk_afdesign_model, clear_mem
  326. from colabdesign.af.alphafold.common import residue_constants
  327. import jax
  328. import jax.numpy as jnp
  329. #@title # OPTIONAL Unconditional Generation (Manuscript Code)
  331. #@markdown This code generates a sample of 10 unconditional proteins for lengths between 100 and 800 AA exactly as in the manuscript. For larger proteins CUDA_UNIFIED_MEMORY is needed. This can be done by localy running the code on a CUDA capeable GPU with sufficient memory (A100 80GB e.g.) and running the code with the environment variables XLA_PYTHON_CLIENT_MEM_FRACTION=100.0 TF_FORCE_UNIFIED_MEMORY=1
  332. def rg_loss(inputs, outputs):
  333. positions = outputs["structure_module"]["final_atom_positions"]
  334. ca = positions[::5,residue_constants.atom_order["CA"]]
  335. center = ca.mean(0)
  336. rg = jnp.sqrt(jnp.square(ca - center).sum(-1).mean() + 1e-8)
  337. rg_th = 2.38 * ca.shape[0] ** 0.365
  338. rg = jax.nn.elu(rg - rg_th)
  339. return {"rg":rg}
  344. for length in [100,200,300,400,500,600,700,800]:
  345. model = mk_afdesign_model(protocol="hallucination",loss_callback=rg_loss)
  346. model.prep_inputs(length=length)
  347. print("weights",model.opt["weights"])
  348. print('Starting up and compiling JAX model....')
  350. for i in range(10):
  351. model.restart(mode=["gumbel", "soft"],rm_aa="C")
  352. model.opt["weights"]["rg"] = 0.1
  353. if length > 600:
  354. model.opt["weights"]["rg"] = 0.01
  355. #model.opt["weights"]['helix'] = -0.1
  356. model.opt["weights"]['plddt'] = 1.0
  357. model.opt["weights"]['pae'] = 1.0
  358. model.opt["weights"]['helix'] = -0.1
  359. print("weights", model.opt["weights"])
  360. model.design_logits(100)
  362. #change the output path for local execution
  363. model.save_pdb(f"Hallo_{i}.pdb")
  365. # %%
  366. #@markdown #Redesign with ProteinMPNN for ESMFold prediction
  367. #@markdown The standard manuscript settings were 8 sequences, 0.1 sampling temperature and the removal of cysteines
  368. import pickle
  369. num_seqs = 8 #@param ["8", "16", "32", "64"] {type:"raw"}
  370. mpnn_sampling_temp = 0.1 #@param ["0.0001", "0.1", "0.15", "0.2", "0.25", "0.3", "0.5", "1.0"] {type:"raw"}
  371. rm_aa = "C" #@param {type:"string"}
  372. use_solubleMPNN = False #@param {type:"boolean"}
  373. #@markdown - `mpnn_sampling_temp` - control diversity of sampled sequences. (higher = more diverse).
  374. #@markdown - `rm_aa='C'` - do not use [C]ysteines.
  375. #@markdown - `use_solubleMPNN` - use weights trained only on soluble proteins.
  376. #@markdown
  378. from colabdesign.shared.protein import alphabet_list as chain_list
  379. mpnn_model = mk_mpnn_model()
  380. mpnn_model.prep_inputs(pdb_filename=f"{af_model.protocol}.pdb",
  381. chain=",".join(chain_list[:COPIES]),
  382. homooligmer=COPIES>1,
  383. rm_aa=rm_aa,
  384. weights = "soluble" if use_solubleMPNN else"original")
  385. out = mpnn_model.sample(num=num_seqs//8,
  386. batch=8,
  387. temperature=mpnn_sampling_temp)
  388. for seq,score in zip(out["seq"],out["score"]):
  389. print(score,seq.split("/")[0])
  390. df = pd.DataFrame(out["seq"])
  392. # Define the output path for saving the sequences as a .pkl file
  393. output_pkl_file = "redesigned_sequences.pkl"
  395. # Save the DataFrame to a .pkl file
  396. with open(output_pkl_file, 'wb') as f:
  397. pickle.dump(df, f)
  399. # %%
  400. #@markdown #Run ESMFold to test designability
  401. #@markdown This cells runs ESMFold from huggingface and automatically calculates the RMSD to the designed backbone
  402. #@markdown NOTE: GPU memory can be a big problem here. If you get memory errors please restart the runtime and run this cell again. It should be self contained. Additionally, after finish the ESMFold prediction rerun the setup cell
  404. import os
  405. import pandas as pd
  406. from Bio.PDB import PDBParser, Superimposer
  407. import pickle
  408. import torch
  409. import numpy as np
  410. from transformers import AutoTokenizer, EsmForProteinFolding
  411. from transformers.models.esm.openfold_utils.protein import to_pdb, Protein as OFProtein
  412. from transformers.models.esm.openfold_utils.feats import atom14_to_atom37
  413. output_pkl_file = "redesigned_sequences.pkl"
  415. with open(output_pkl_file, 'rb') as f:
  416. seq = pickle.load(f)
  417. seq_list = []
  418. for i in np.asarray(seq):
  419. seq_list.append(i[0])
  421. pdb_file = "hallucination.pdb"
  422. print(seq_list)
  425. tokenizer = AutoTokenizer.from_pretrained("facebook/esmfold_v1")
  426. model = EsmForProteinFolding.from_pretrained("facebook/esmfold_v1", low_cpu_mem_usage=True)
  428. device = 'cuda:0'
  429. model = model.cuda(device)
  430. model.esm = model.esm.half()
  431. model.trunk.set_chunk_size(64)
  432. torch.backends.cuda.matmul.allow_tf32 = True
  434. def convert_outputs_to_pdb(outputs):
  435. final_atom_positions = atom14_to_atom37(outputs["positions"][-1], outputs)
  436. outputs = {k: v.to("cpu").numpy() for k, v in outputs.items()}
  437. final_atom_positions = final_atom_positions.cpu().numpy()
  438. final_atom_mask = outputs["atom37_atom_exists"]
  439. pdbs = []
  440. for i in range(outputs["aatype"].shape[0]):
  441. aa = outputs["aatype"][i]
  442. pred_pos = final_atom_positions[i]
  443. mask = final_atom_mask[i]
  444. resid = outputs["residue_index"][i] + 1
  445. pred = OFProtein(
  446. aatype=aa,
  447. atom_positions=pred_pos,
  448. atom_mask=mask,
  449. residue_index=resid,
  450. b_factors=outputs["plddt"][i],
  451. chain_index=outputs["chain_index"][i] if "chain_index" in outputs else None,
  452. )
  453. pdbs.append(to_pdb(pred))
  454. return pdbs
  456. def calculate_ca_rmsd(pdb_file1, pdb_file2):
  457. parser = PDBParser(QUIET=True)
  459. structure1 = parser.get_structure("Protein1", pdb_file1)
  460. structure2 = parser.get_structure("Protein2", pdb_file2)
  462. ca_atoms1 = [atom for atom in structure1.get_atoms() if atom.get_name() == "CA"]
  463. ca_atoms2 = [atom for atom in structure2.get_atoms() if atom.get_name() == "CA"]
  465. super_imposer = Superimposer()
  466. super_imposer.set_atoms(ca_atoms1, ca_atoms2)
  467. super_imposer.apply(structure2.get_atoms())
  468. rmsd = super_imposer.rms
  469. return rmsd
  471. def process_sequences(seq_list, pdb_file):
  472. lowest_rmsd = float('inf')
  473. lowest_rmsd_data = None
  474. out_ss_path = "./output"
  476. if not os.path.exists(out_ss_path):
  477. os.mkdir(out_ss_path)
  479. for test_protein in seq_list:
  480. data = {}
  481. tokenized_input = tokenizer([test_protein], return_tensors="pt", add_special_tokens=False)['input_ids']
  482. tokenized_input = tokenized_input.cuda(device)
  484. with torch.no_grad():
  485. output = model(tokenized_input)
  487. data['out'] = output
  488. data["plddt"] = torch.mean(output['plddt']).item()
  489. data['pae'] = torch.mean(output['predicted_aligned_error']).item()
  491. pdb_data = convert_outputs_to_pdb(output)
  492. tmp_pdb_file = os.path.join(out_ss_path, "TMP.pdb")
  494. with open(tmp_pdb_file, 'w') as file:
  495. for line in pdb_data:
  496. file.write(line)
  498. data['rmsd'] = calculate_ca_rmsd(tmp_pdb_file, pdb_file)
  499. print(f'Sequence: {test_protein}, plddt: {data["plddt"]}, PAE: {data["pae"]}, RMSD: {data["rmsd"]}')
  501. if data['rmsd'] < lowest_rmsd:
  502. lowest_rmsd = data['rmsd']
  503. lowest_rmsd_data = data
  505. if lowest_rmsd_data is not None:
  506. print(f'Lowest RMSD: {lowest_rmsd}')
  507. best_pdb_data = convert_outputs_to_pdb(lowest_rmsd_data['out'])
  508. best_pdb_file = os.path.join(out_ss_path, "best_structure.pdb")
  510. with open(best_pdb_file, 'w') as file:
  511. for line in best_pdb_data:
  512. file.write(line)
  514. original_dict = lowest_rmsd_data
  515. key_to_exclude = 'out'
  516. data_out = {k: v for k, v in original_dict.items() if k != key_to_exclude}
  518. with open(os.path.join(out_ss_path, "best_structure_data.pkl"), 'wb') as f:
  519. pickle.dump(data_out, f)
  521. return lowest_rmsd, best_pdb_file, data_out
  523. return None, None, None
  527. lowest_rmsd, best_pdb_file, best_data = process_sequences(seq_list, pdb_file)
  528. if lowest_rmsd is not None:
  529. print(f"Lowest RMSD: {lowest_rmsd}, Best PDB file: {best_pdb_file}")
  530. else:
  531. print("No valid result found.")
  534. # %%
  535. #@title # Heterodimer Design Prep
  536. #@markdown Design a set of heterodimeric proteins with two chains making a complex. The settings are excatly the ones used in the manuscript to design the heterodimer binders.
  537. LENGTH1 = 100 #@param {type:"integer"}
  538. LENGTH2 = 100 #@param {type:"integer"}
  540. #@markdown ProteinMPNN Settings
  541. use_solubleMPNN = True #@param {type:"boolean"}
  542. #@markdown
  545. from colabdesign.af.alphafold.common import residue_constants
  546. import jax
  547. import jax.numpy as jnp
  549. def hd_loss(inputs, outputs):
  550. positions = outputs["structure_module"]["final_atom_positions"]
  551. ca1 = positions[:LENGTH1, residue_constants.atom_order["CA"]]
  552. center1 = ca1.mean(0)
  553. rg1 = jnp.sqrt(jnp.square(ca1 - center1).sum(-1).mean() + 1e-8)
  554. rg_th = 2.38 * ca1.shape[0] ** 0.365
  555. rg1 = jax.nn.elu(rg1 - rg_th)
  558. ca2 = positions[LENGTH2:, residue_constants.atom_order["CA"]]
  559. center2 = ca2.mean(0)
  560. rg2 = jnp.sqrt(jnp.square(ca2 - center2).sum(-1).mean() + 1e-8)
  561. rg_th = 2.38 * ca2.shape[0] ** 0.365
  562. rg2 = jax.nn.elu(rg2 - rg_th)
  566. return {"hd":rg1+rg2}
  568. total_length = LENGTH1 + LENGTH2
  569. clear_mem()
  570. af_model = mk_afdesign_model(protocol="hallucination", loss_callback=hd_loss)
  571. af_model.prep_inputs(length=total_length)
  572. af_model._inputs['residue_index'][LENGTH1:] = np.arange(LENGTH2) + 50 + LENGTH1
  573. # add extra losses
  574. af_model.restart(mode=["gumbel", "soft"])
  575. af_model.opt["weights"]["hd"] = 0.1
  576. af_model.opt["weights"]['plddt'] = 1.0
  577. af_model.opt["weights"]['pae'] = 1.0
  578. af_model.opt["weights"]['helix'] = -0.5
  579. print("weights", af_model.opt["weights"])
  580. print('Starting up and compiling JAX model....')
  583. # %%
  584. #@title # Run Design
  585. af_model.design_logits(100)
  586. af_model.save_pdb("Heterodimer.pdb")
  588. # %%
  589. af_model.save_pdb("Heterodimer.pdb")
  590. af_model.plot_pdb()
  592. # %%
  593. #@title # Design Sequence using Homooligomer Filter
  594. #@markdown We first test if the two protomers are predictd to fold into a high confidence protein on their own, removing proteins that are not likely to be expressed on their own. Then we predict the heterodimer using the AF multimer model. Generally the AF multimer model has a hard time predicting de novo designed proteins. This is why we use templates and remove any interchain information. Finally we predict each individual protomer with a copy of itself, testing for homooligomerisation.
  597. file_path ="Heterodimer.pdb"
  599. folder_path = "/content/"
  601. ######## make A - B chain file
  603. from Bio.PDB import PDBParser, PDBIO, Chain
  605. # Set the input and output PDB file names
  606. input_pdb_file = file_path
  607. if not os.path.exists(os.path.join(folder_path, 'AB')):
  608. os.mkdir(os.path.join(folder_path, 'AB'))
  609. output_pdb_file = os.path.join(folder_path, 'AB',"Heterodimer.pdb")
  611. # Create a PDB parser and read the input PDB file
  612. parser = PDBParser()
  613. structure = parser.get_structure("input_structure", input_pdb_file)
  615. # Find the initial chain id
  616. initial_chain_id = None
  617. for chain in structure[0]:
  618. initial_chain_id = chain.get_id()
  619. break
  621. # Create new chains A and B
  622. chain_A = Chain.Chain("A")
  623. chain_B = Chain.Chain("B")
  625. # Iterate over the residues in the original chain
  626. for residue in structure[0][initial_chain_id]:
  627. res_id = residue.get_id()[1]
  629. # Add residues 1-200 to chain A
  630. if 1 <= res_id <= 100:
  631. chain_A.add(residue.copy())
  633. # Add residues 201-400 to chain B
  634. elif 151 <= res_id <= 450:
  635. chain_B.add(residue.copy())
  637. # Remove the existing chain
  638. for model in structure:
  639. model.detach_child(initial_chain_id)
  641. # Add the new chains to the model
  642. structure[0].add(chain_A)
  643. structure[0].add(chain_B)
  645. # Save the modified structure to a new PDB file
  646. io = PDBIO()
  647. io.set_structure(structure)
  648. io.save(output_pdb_file)
  651. clear_mem()
  652. he_model = mk_afdesign_model(protocol="fixbb", use_templates=True, use_multimer=True)
  655. ho_model = mk_afdesign_model(protocol="hallucination")
  656. ho_model.prep_inputs(length=LENGTH1, copies=2)
  658. ho_model.set_weights(i_pae=1.0)
  659. s_model = mk_afdesign_model(protocol="hallucination")
  660. s_model.prep_inputs(length=LENGTH2)
  661. mpnn_model = mk_mpnn_model(weights="soluble")
  664. mpnn_model.prep_inputs(pdb_filename=output_pdb_file, chain='A,B',rm_aa="C")
  665. samples = mpnn_model.sample_parallel(8)
  667. he_model.prep_inputs(pdb_filename=output_pdb_file, chain='A,B',rm_template_ic=True)
  668. he_model._inputs['residue_index'][LENGTH1:] = np.arange(LENGTH2) + 50 + LENGTH1
  670. k = 0
  671. for seq in samples['seq']:
  672. print('Predicting Protomer 1...')
  673. s_model.predict(seq=seq[:LENGTH1], num_recycles=3)
  674. plddt1 = s_model.aux['losses']['plddt']
  675. print('Predicting Protomer 2...')
  676. s_model.predict(seq=seq[LENGTH1+1:], num_recycles=3)
  677. plddt2 = s_model.aux['losses']['plddt']
  678. k = k + 1
  679. if plddt1 < 0.20 and plddt2 < 0.20:
  680. print('Passed Protomer Check! Predicting Heterodimer...')
  681. he_model.predict(seq=''.join([seq[:LENGTH1], seq[LENGTH1+1:]]), num_recycles=3)
  683. if he_model.aux['losses']['plddt'] < 0.15 and he_model.aux['losses']['rmsd'] < 2.0:
  684. print('Passed Heterodimer Check! Predicting Homodimer 1...')
  685. ho_model.predict(seq=seq[:LENGTH1],num_recycles=3)
  686. print('Predicting Homodimer 2...')
  687. ipae1 = ho_model.aux['losses']['i_pae']
  688. ho_model.predict(seq=seq[LENGTH1+1:],num_recycles=3)
  689. ipae2 = ho_model.aux['losses']['i_pae']
  690. if ipae1 > 0.8 and ipae2 > 0.8:
  691. print('Passed Homodimer check!')
  692. he_model.save_pdb(f'Heterodimer_seq_{k}.pdb')
  696. # %%
  697. def get_pdb(pdb_code=""):
  698. if pdb_code is None or pdb_code == "":
  699. upload_dict = files.upload()
  700. pdb_string = upload_dict[list(upload_dict.keys())[0]]
  701. with open("tmp.pdb","wb") as out: out.write(pdb_string)
  702. return "tmp.pdb"
  703. elif os.path.isfile(pdb_code):
  704. return pdb_code
  705. elif len(pdb_code) == 4:
  706. os.system(f"wget -qnc https://files.rcsb.org/view/{pdb_code}.pdb")
  707. return f"{pdb_code}.pdb"
  708. else:
  709. os.system(f"wget -qnc https://alphafold.ebi.ac.uk/files/AF-{pdb_code}-F1-model_v3.pdb")
  710. return f"AF-{pdb_code}-F1-model_v3.pdb"
  712. def add_rg_loss(self, weight=0.1):
  713. '''add radius of gyration loss'''
  714. def loss_fn(inputs, outputs):
  715. xyz = outputs["structure_module"]
  716. ca = xyz["final_atom_positions"][:,residue_constants.atom_order["CA"]]
  718. ca = ca[-self._binder_len:]
  720. rg = jnp.sqrt(jnp.square(ca - ca.mean(0)).sum(-1).mean() + 1e-8)
  721. rg_th = 2.38 * ca.shape[0] ** 0.365
  722. rg = jax.nn.elu(rg - rg_th)
  723. return {"rg":rg}
  724. self._callbacks["model"]["loss"].append(loss_fn)
  725. self.opt["weights"]["rg"] = weight
  729. #@title # Binder Design
  730. #@markdown For a given length, generate/hallucinate a protein sequence that AlphaFold thinks folds into a well structured protein (high plddt, low pae, many contacts).
  731. LENGTH = 100 #@param {type:"integer"}
  732. binder_pdb = '5NGV' #@param {type:"string"}
  733. binder_chain ='A' #@param {type:"string"}
  734. hotspot ='' #@param {type:"string"}
  735. if hotspot == "": hotspot = None
  736. #@markdown ProteinMPNN Settings
  737. use_solubleMPNN = True #@param {type:"boolean"}
  738. #@markdown
  740. clear_mem()
  741. af_model = mk_afdesign_model(protocol="binder")
  742. add_rg_loss(af_model)
  743. af_model.prep_inputs(pdb_filename=get_pdb(binder_pdb), chain=binder_chain,hotspot=hotspot, binder_len=LENGTH)
  746. af_model.restart(mode=["gumbel", "soft"])
  748. af_model.opt["weights"]["rg"] = 0.5
  750. af_model.opt["weights"]['helix'] = -0.2
  751. af_model.opt["weights"]['plddt'] = 0.1
  752. af_model.opt["weights"]['pae'] = 0.1
  753. af_model.opt["weights"]['i_pae'] = 0.1
  754. af_model.opt["weights"]['i_con'] = 2.0
  756. print("weights", af_model.opt["weights"])
  757. print('Starting up and compiling JAX model....')
  760. # %%
  761. af_model.design_logits(100)
  762. af_model.save_pdb("Binder.pdb")
  764. # %%
  765. af_model.plot_pdb()
  767. # %%
  769. #@title # Binder Sequence Design with AF Multimer filtering
  770. #@markdown Use this to generate sequences for the binder candidate generated in the previous step
  772. #@markdown First we use the AF2 PTM model to predict the binder without receptor, acting as a fast pre filter. Then we use the AF Multimer model to predict the Receptor Binder complex. Again we use a template for the binder to help AF Multimer predicting the de novo designed protein
  774. binder_model = mk_afdesign_model(protocol="binder",use_multimer=True,use_initial_guess=True)
  775. hall_model = mk_afdesign_model(protocol="fixbb")
  778. binder_model.set_weights(i_pae=1.0)
  781. mpnn_model = mk_mpnn_model(weights="soluble")
  782. mpnn_model.prep_inputs(pdb_filename="Binder.pdb", chain='A,B', fix_pos='A',rm_aa="C")
  784. samples = mpnn_model.sample_parallel(8,temperature=0.01)
  785. hall_model.prep_inputs(pdb_filename="Binder.pdb", chain='B')
  786. binder_model.prep_inputs(pdb_filename="Binder.pdb", chain='A', binder_chain='B',use_binder_template=True,rm_template_ic=True)
  787. k=0
  788. for seq in samples['seq']:
  789. print("Predicting binder only")
  790. hall_model.predict(seq=seq[-LENGTH:], num_recycles=3)
  791. if hall_model.aux['losses']['rmsd'] < 2.0 :
  792. print("Passed! Predicting binder with receptor using AF Multimer")
  793. binder_model.predict(seq=seq[-LENGTH:], num_recycles=3)
  794. plddt1 = binder_model.aux['losses']['plddt']
  795. i_pae = binder_model.aux['losses']['i_pae']
  796. if plddt1 < 0.15 and i_pae < 0.4:
  797. print(f"Passed! Final I_PAE is {i_pae*31}")
  798. binder_model.save_pdb(f'Binder_seq_{k}.pdb')
  799. binder_model.plot_pdb()
  801. k = k + 1
  803. # %%
  804. #@title # Site scaffolding example
  805. #@markdown This cell provides the code to perform the site scaffolding in bulk.
  806. #@markdown Just go to the commented section with names, contigs and length to insert the desired PDB identifier, contigs and final size and start designing.
  807. #@markdown Num_designs controls how many backbones one designes per PDB file
  809. num_designs = 1 #@param {type:"integer"}
  812. def get_pdb(pdb_code=""):
  813. if pdb_code is None or pdb_code == "":
  814. upload_dict = files.upload()
  815. pdb_string = upload_dict[list(upload_dict.keys())[0]]
  816. with open("tmp.pdb","wb") as out: out.write(pdb_string)
  817. return "tmp.pdb"
  818. elif os.path.isfile(pdb_code):
  819. return pdb_code
  820. elif len(pdb_code) == 4:
  821. os.system(f"wget -qnc https://files.rcsb.org/view/{pdb_code}.pdb")
  822. return f"{pdb_code}.pdb"
  823. else:
  824. os.system(f"wget -qnc https://alphafold.ebi.ac.uk/files/AF-{pdb_code}-F1-model_v3.pdb")
  825. return f"AF-{pdb_code}-F1-model_v3.pdb"
  829. from colabdesign import mk_afdesign_model, clear_mem
  830. import contextlib
  832. from colabdesign.af.alphafold.common import residue_constants
  833. import jax
  834. import jax.numpy as jnp
  835. import pickle
  836. from colabdesign.mpnn import mk_mpnn_model
  838. import re
  839. import os
  842. #Add the names of the PDB files for the scaffolding problem here
  843. names = [
  844. "1PRW"
  845. ]
  846. print(len(names))
  847. #Add the design contigs here
  849. inputs = [
  850. "5-20,A16-35,10-25,A52-71,5-20"
  851. ]
  852. #Add the total length here. We only use the maximum length specified
  853. lengths = [
  854. "60-105"
  855. ]
  858. def rg_loss(inputs, outputs):
  859. positions = outputs["structure_module"]["final_atom_positions"]
  860. ca = positions[::5, residue_constants.atom_order["CA"]]
  861. center = ca.mean(0)
  862. rg = jnp.sqrt(jnp.square(ca - center).sum(-1).mean() + 1e-8)
  863. rg_th = 2.38 * ca.shape[0] ** 0.365
  864. rg = jax.nn.elu(rg - rg_th)
  865. return {"rg": rg}
  868. clear_mem()
  870. for _name, _input, _length in zip(
  871. names, inputs, lengths
  872. ):
  873. print(f"Starting on {_name}")
  874. _input = _input.replace(" ", "")
  875. __name = _name.split("_")[0]
  876. model = mk_afdesign_model(
  877. protocol="partial"
  878. )
  879. wire_loop_repr = ["l" if re.search("[A-Z]", x) else "w" for x in _input.split(",")]
  881. _lengths = []
  882. for _id, rep in zip(wire_loop_repr, _input.split(",")):
  883. if "-" in rep: # loop or range
  884. if _id == "l": # loop
  885. rep = rep[1:]
  886. _len = int(rep.split("-")[1]) - int(rep.split("-")[0]) + 1
  887. else: # range
  888. _len = int(rep.split("-")[1])
  889. else:
  890. if _id == "l":
  891. rep = rep[1:]
  892. _len = int(1)
  893. _lengths.append(_len)
  895. overall_length = sum(_lengths)
  896. print(overall_length)
  897. old_pos = list(filter(lambda x: re.search("[A-Z]", x), _input.split(",")))
  898. order = list(range(len(old_pos)))
  899. old_pos = ",".join(old_pos)
  900. wires = list(filter(lambda x: not re.search("[A-Z]", x), _input.split(",")))
  901. wires = [
  902. int(wire) if "-" not in wire else int(wire.split("-")[1]) for wire in wires
  903. ]
  904. offset = wires[0] if not wire_loop_repr[0] == "l" else 0
  905. if wire_loop_repr[0] == "w":
  906. wires = wires[1:]
  907. if wire_loop_repr[-1] == "w":
  908. wires = wires[:-1]
  910. chain = re.findall("[A-Z]", _input)
  911. chain = list(set(chain))
  912. assert len(chain) == 1
  913. chain = chain[0]
  914. if "-" in _length:
  915. _length = _length.split("-")[1]
  916. _length = int(_length)
  917. if _length < overall_length:
  918. _length = overall_length
  920. debug = False
  921. if debug:
  922. print("chain " + str(chain))
  923. print("old_pos " + str(old_pos))
  924. print("wires " + str(wires))
  925. print("offset " + str(offset))
  926. print("_length " + str(_length))
  927. print("order " + str(order))
  928. print(_name)
  929. pdb_file = get_pdb(_name)
  930. model.prep_inputs(
  931. pdb_file,
  932. chain=chain,
  933. pos=old_pos,
  934. length=_length,
  935. fix_seq=True,
  936. )
  939. model.rewire(
  940. order=order, # set order of segments
  941. loops=wires, # change loop length inbetween segments
  942. offset=offset,
  943. ) # essentially loop length at the N term
  945. print(" Starting up and compiling JAX model....")
  947. for i in range(num_designs):
  948. print(f" Iteration {i} of 100")
  949. model.restart(mode=["gumbel", "soft"], rm_aa="C")
  950. model.opt["weights"]["rg"] = 0.1
  951. model.opt["weights"]["dgram_cce"] = 2.0
  952. model.opt["weights"]["plddt"] = 0.1
  953. model.opt["weights"]["pae"] = 0.1
  954. model.opt["weights"]["rmsd"] = 1.0
  955. model.opt["weights"]['sc_rmsd'] = 1.0
  956. # model.opt["weights"]['fape'] = 1.0
  958. model.design_logits(190)
  959. model.design_logits(10, save_best=True)
  960. outfile = f"out_sc/{_name}_resesigned/{_name}_redesigned_{i}.pdb"
  961. os.makedirs(os.path.dirname(outfile), exist_ok=True)
  962. model.save_pdb(outfile)
  963. mpnn_model = mk_mpnn_model()
  965. p = (
  966. []
  967. ) # [homo if not n in _interfaceFixturesIndexSecChain else hetero for n, (homo, hetero) in enumerate(zip(list(ho2), list(he[-len(ho2):])))]
  968. for k in model.opt["pos"]:
  969. p.append(str(k + 1)) # Might be wrong
  970. p.append(",")
  971. posf = "".join(p[:-1])
  973. repredictionModel = mk_afdesign_model(
  974. protocol="fixbb", use_templates=False
  975. )
  976. os.makedirs(os.path.dirname('out_sc_Redesigned/'), exist_ok=True)
  978. for j in range(num_designs):
  979. print(f" Reprediction Iteration {j} of 100")
  980. repredictionModel.prep_inputs(
  981. f"out_sc/{_name}_resesigned/{_name}_redesigned_{j}.pdb"
  982. )
  984. mpnn_model.prep_inputs(
  985. pdb_filename=f"out_sc/{_name}_resesigned/{_name}_redesigned_{j}.pdb",
  986. chain="A",
  987. fix_pos=posf,
  988. rm_aa="C",
  989. )
  990. out = mpnn_model.sample(num=1, batch=8, temperature=0.1)
  992. for n, i in enumerate(out["seq"]):
  994. repredictionModel.predict(seq=i, num_recycles=3)
  995. if (
  996. repredictionModel.aux["log"]["rmsd"] < 2.0
  997. and repredictionModel.aux["log"]["plddt"] > 0.85
  998. ):
  999. filename = f'out_sc_Redesigned/{_name}_resesigned/{_name}_redesigned-{j}_num-{n}_rmsd-{int(repredictionModel.aux["log"]["rmsd"]*100)}.pdb'
  1000. os.makedirs(os.path.dirname(filename), exist_ok=True)
  1001. repredictionModel.save_pdb(filename)
  1003. for _name, _input, _length in zip(
  1004. names, inputs, lengths
  1005. ):
  1006. print(f"Starting on {_name}")
  1007. clear_mem()
  1009. _input = _input.replace(" ", "")
  1010. __name = _name.split("_")[0]
  1012. test_model = mk_afdesign_model(protocol='fixbb')
  1013. model = mk_afdesign_model(
  1014. protocol="partial", use_templates=False
  1015. ) # set True to constrain positions using template input
  1016. # define positions we want to constrain (input PDB numbering)
  1018. wire_loop_repr = ["l" if re.search("[A-Z]", x) else "w" for x in _input.split(",")]
  1020. _lengths = []
  1021. for _id, rep in zip(wire_loop_repr, _input.split(",")):
  1022. if "-" in rep: # loop or range
  1023. if _id == "l": # loop
  1024. rep = rep[1:]
  1025. _len = int(rep.split("-")[1]) - int(rep.split("-")[0]) + 1
  1026. else: # range
  1027. _len = int(rep.split("-")[1])
  1028. else:
  1029. if _id == "l":
  1030. rep = rep[1:]
  1031. _len = 1
  1032. _lengths.append(_len)
  1034. overall_length = sum(_lengths)
  1036. old_pos = list(filter(lambda x: re.search("[A-Z]", x), _input.split(",")))
  1037. order = list(range(len(old_pos)))
  1038. old_pos = ",".join(old_pos)
  1039. wires = list(filter(lambda x: not re.search("[A-Z]", x), _input.split(",")))
  1040. wires = [
  1041. int(wire) if "-" not in wire else int(wire.split("-")[1]) for wire in wires
  1042. ]
  1043. offset = wires[0] if not wire_loop_repr[0] == "l" else 0
  1044. if wire_loop_repr[0] == "w":
  1045. wires = wires[1:]
  1046. if wire_loop_repr[-1] == "w":
  1047. wires = wires[:-1]
  1049. chain = re.findall("[A-Z]", _input)
  1050. chain = list(set(chain))
  1051. assert len(chain) == 1
  1052. chain = chain[0]
  1053. if "-" in _length:
  1054. _length = _length.split("-")[1]
  1055. _length = int(_length)
  1056. if _length < overall_length:
  1057. _length = overall_length
  1060. print(_name)
  1061. pdb_file = get_pdb(_name)
  1065. model.prep_inputs(
  1066. pdb_file,
  1067. chain=chain,
  1068. pos=old_pos, # define positions to contrain
  1069. length=_length, # define if the desired length is different from input PDB
  1070. fix_seq=True,
  1071. ) # set True to constrain the sequence
  1073. # set positions (if different from PDB)
  1074. # reorder the segments,
  1075. model.rewire(
  1076. order=order, # set order of segments
  1077. loops=wires, # change loop length inbetween segments
  1078. offset=offset,
  1079. ) # essentially loop length at the N term
  1081. in_files = os.listdir(f'out_sc_Redesigned/{_name}_resesigned/')
  1082. if not os.path.exists(f'out_sc_Redesigned/{_name}_resesigned/out/'):
  1083. os.mkdir(f'out_sc_Redesigned/{_name}_resesigned/out/')
  1084. for ii in in_files:
  1085. if ii[-1] == 'b':
  1086. test_model.prep_inputs(pdb_filename=f'out_sc_Redesigned/{_name}_resesigned/{ii}')
  1087. seq = test_model._inputs['batch']["aatype"]
  1088. #print(seq)
  1089. model.predict(seq=seq, num_recycles=3)
  1090. if model.aux["losses"]["rmsd"] < 1.0:
  1091. model.save_pdb(f'out_sc_Redesigned/{_name}_resesigned/out/{ii}')
  1092. with open(f'out_sc_Redesigned/{_name}_resesigned/out/{ii[:-4]}_data.pkl', 'wb') as f:
  1093. pickle.dump(model.aux["losses"]["rmsd"], f)
  1097. # %%

RSO.ipynb at commit e31a56f, under other · at the source

Overview

Authors: Mak B Djulbegovic1, Nedym Hadzijahic2, David J Taylor Gonzalez3, Michael Antonietti4, Sidra Zafar1,5, Ajay E Kuriyan1,2
ORCID iDs: Mak B Djulbegovic, Nedym Hadzijahic, Ajay E Kuriyan
  1. Wills Eye Hospital, Thomas Jefferson University Hospital, Philadelphia, Pennsylvania
  2. University of Miami, Miami, Florida
  3. Department of Ophthalmology, Broward Health North, Pompano Beach, Florida
  4. Department of Ophthalmology, Massachusetts Eye and Ear Infirmary, Harvard Medical School, Boston, Massachusetts
  5. Mid Atlantic Retina at Wills Eye Hospital, Philadelphia, Pennsylvania
Institutions: Thomas Jefferson University Hospital (United States); Wills Eye Hospital (United States); University of Miami (United States); Broward Health (United States); Massachusetts Eye and Ear Infirmary (United States); Harvard University (United States); Mid Atlantic Retina (United States)
Journal: Ophthalmology science, volume 6, issue 8, article 101249
Dates: received 2 February 2026; accepted 18 May 2026; published online 25 May 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1016/j.xops.2026.101249 · PMID 42421755 · PMCID PMC13343151 · OpenAlex W7162299899
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: human (organism)
Keywords: Proliferative vitreoretinopathy, Epithelial–mesenchymal transition, Intrinsic disorder, SNAIL1, Artificial intelligence
Topic: Retinal Imaging and Analysis (Radiology, Nuclear Medicine and Imaging, Medicine), according to OpenAlex
Funding: VitreoRetinal Surgery Foundation
Citations: not cited yet (Europe PMC); 62 references in the paper

Abstract

Objective: Proliferative vitreoretinopathy (PVR) remains a major cause of failure after rhegmatogenous retinal detachment repair and lacks effective pharmacologic therapies. Although epithelial–mesenchymal transition (EMT) is central to PVR pathogenesis, the structural determinants governing the tractability of EMT regulators, particularly those involving intrinsic disorder, remain poorly defined. We developed a disorder-aware, artificial intelligence–enabled computational framework to evaluate EMT-associated proteins in PVR and prioritize structurally tractable regulators for structure-based targeting.

Design: A computational, hypothesis-generating study employing an in silico screening and structural modeling pipeline.

Subjects: No human subjects or biological specimens were included. The dataset comprised 25 EMT-associated proteins implicated in PVR, curated through a narrative review of peer-reviewed literature.

Methods: Candidate proteins were evaluated using a multistage pipeline integrating intrinsic disorder profiling (Rapid Intrinsic Disorder Analysis Online), redox-sensitive disorder-to-order transition (DOT) analysis (AIUPred), and protein–protein interaction network assessment (Search Tool for the Retrieval of Interacting Genes/Proteins [STRING]). Structure-based modeling and generative binder design were then applied to the top-ranked candidate using RFdiffusion for de novo backbone generation, protein message passing neural network for sequence design, and AlphaFold2 for structural validation.

Main Outcome Measures: Primary measures were the proportion of intrinsically disordered residues, redox-sensitive disorder change, STRING network coherence within EMT-related pathways, and the structural consistency of the designed binder–target complex, assessed by root mean square deviation (RMSD) and mean per-residue confidence (predicted local distance difference test [pLDDT]).

Results: Of the 25 EMT-associated proteins screened, several exhibited intermediate intrinsic disorder profiles and measurable DOT potential. Snail Family Transcriptional Repressor 1 (SNAIL1) emerged as the highest-priority candidate, demonstrating an intermediate intrinsic disorder profile (∼35%), a pronounced redox-sensitive DOT region, and selective connectivity within EMT-related signaling networks. Functional mapping of the SNAIL1 C-terminal DOT segment identified 6 basic residues with literature-supported or motif-based regulatory significance (K187, R191, R224, K234, K253, and R264). Following sequence design and structural validation, the top-ranked binder exhibited the lowest structural deviation within the generated ensemble (RMSD 18.5 Å) and high per-residue confidence (mean pLDDT 0.84).

Conclusions: Our study introduces a disorder-informed computational framework for prioritizing structurally tractable EMT regulators in PVR. As a proof-of-concept, the pipeline nominates SNAIL1 and generates a structure-aware de novo binder targeting its C-terminal DOT region, providing a foundation for disorder-based therapeutic discovery in fibrotic retinal disease.

Financial Disclosure(s): Proprietary or commercial disclosure may be found in the Footnotes and Disclosures at the end of this article.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 4 matches between paragraphs and lines of code.

sokrypton/colabdesign

License: other
State: the link answers, verified on 28 September 2026
Evidence: files inventoried
Commit: e31a56fe1d9b4de25c8697f3a28b75892941cc72, 23 October 2025
Languages: Python (95), Jupyter (28), JavaScript (1)
Size: 155 files, 124 scripts
Software Heritage: not archived
Found in: the text, “RFdiffusion-Based Binder Design”
Holds: README, license file, environment (setup.py), continuous integration, 28 notebooks
Not found: CITATION.cff, tests, documentation
Tools: NumPy (77 files), JAX (71 files), Matplotlib (14 files), SciPy (10 files), pandas (7 files), Biopython (5 files), PyTorch (5 files), Plotly (2 files), Hugging Face Transformers (2 files), Keras (1 file), TensorFlow (1 file)
Availability: 1 check, the latest on 28 September 2026: the link answers
  • 28 September 2026: the link answers
126 files
At the source: colab.research.google.com/github/sokrypton/ColabDesign/blob/main/rf/examples/diffusion_ori.ipynb

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 124 scripts, each with its path and the digest of its content;
  • 4 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

Datasets cited

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 2, 28 September 2026

  • Authors: added Nedym Hadzijahic (0009-0004-1454-2705); removed Nedym Hadzijahic

Version 1, 28 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 6 authors, 5 keywords, 1 funder, 62 references.

Cite

This paper

Djulbegovic, M. B., Hadzijahic, N., Taylor Gonzalez, D. J., Antonietti, M., Zafar, S., & Kuriyan, A. E. (2026). A Disorder-Aware Computational Framework to Identify Structurally Tractable Targets in Proliferative Vitreoretinopathy. Ophthalmology science, 6(8), 101249. https://doi.org/10.1016/j.xops.2026.101249

BibTeX

@article{djulbegovic2026disorder,
author = {Djulbegovic, Mak B and Hadzijahic, Nedym and Taylor Gonzalez, David J and Antonietti, Michael and Zafar, Sidra and Kuriyan, Ajay E},
title = {{A Disorder-Aware Computational Framework to Identify Structurally Tractable Targets in Proliferative Vitreoretinopathy}},
journal = {Ophthalmology science},
year = {2026},
month = may,
volume = {6},
number = {8},
pages = {101249},
publisher = {Elsevier},
issn = {2666-9145},
doi = {10.1016/j.xops.2026.101249},
url = {https://doi.org/10.1016/j.xops.2026.101249},
pmid = {42421755},
pmcid = {PMC13343151}
}

RIS

TY - JOUR
AU - Djulbegovic, Mak B
AU - Hadzijahic, Nedym
AU - Taylor Gonzalez, David J
AU - Antonietti, Michael
AU - Zafar, Sidra
AU - Kuriyan, Ajay E
TI - A Disorder-Aware Computational Framework to Identify Structurally Tractable Targets in Proliferative Vitreoretinopathy
T2 - Ophthalmology science
J2 - Ophthalmol Sci
PY - 2026
DA - 2026/05/25
VL - 6
IS - 8
SP - 101249
SN - 2666-9145
PB - Elsevier
DO - 10.1016/j.xops.2026.101249
UR - https://doi.org/10.1016/j.xops.2026.101249
LA - en
ER -

CSL-JSON

{
"id": "10.1016/j.xops.2026.101249",
"type": "article-journal",
"title": "A Disorder-Aware Computational Framework to Identify Structurally Tractable Targets in Proliferative Vitreoretinopathy",
"container-title": "Ophthalmology science",
"author": [
{
"family": "Djulbegovic",
"given": "Mak B"
},
{
"family": "Hadzijahic",
"given": "Nedym"
},
{
"family": "Taylor Gonzalez",
"given": "David J"
},
{
"family": "Antonietti",
"given": "Michael"
},
{
"family": "Zafar",
"given": "Sidra"
},
{
"family": "Kuriyan",
"given": "Ajay E"
}
],
"container-title-short": "Ophthalmol Sci",
"volume": "6",
"issue": "8",
"page": "101249",
"DOI": "10.1016/j.xops.2026.101249",
"PMID": "42421755",
"PMCID": "PMC13343151",
"ISSN": "2666-9145",
"publisher": "Elsevier",
"URL": "https://doi.org/10.1016/j.xops.2026.101249",
"language": "en",
"issued": {
"date-parts": [
[
2026,
5,
25
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: JAX, Biopython, Hugging Face Transformers, 8 other tools
[2] doi:10.1371/journal.pone.0346575 [code]
Statistically valid explainable black-box machine learning: applications in sex classification across species using brain imaging.
Journal: PloS one
In common: JAX, Hugging Face Transformers, TensorFlow, 6 other tools
[3] doi:10.1016/j.molcel.2026.07.006 [code]
DeorphaNN: Virtual screening of GPCR peptide agonists using AlphaFold-predicted active-state complexes and deep learning embeddings.
Journal: Molecular cell
In common: JAX, Biopython, TensorFlow, 5 other tools, 1 reference
[4] doi:10.1038/s41586-026-10658-6 [code]
An AI system to help scientists write expert-level empirical software.
Journal: Nature
In common: JAX, Hugging Face Transformers, TensorFlow, 5 other tools, 1 reference
[5] doi:10.1038/s41598-026-53415-5 [code]
Computational design and immunoinformatics validation of a T cell multi-epitope vaccine targeting glioblastoma stem cells.
Journal: Scientific reports
In common: JAX, Biopython, TensorFlow, 5 other tools
[6] doi:10.1038/s41586-026-10391-0 [code]
Cell-type-targeted mitochondrial transplantation rescues cell degeneration.
Journal: Nature
In common: JAX, Biopython, TensorFlow, 5 other tools
[7] doi:10.3390/ijms27156614 [code]
Candidalysin Inhibits &lt;i&gt;Porphyromonas gingivalis&lt;/i&gt; Lipoprotein-Induced IL-1β Production in BV-2 Microglia via Hydrophobic Microbial Interactions.
Journal: International journal of molecular sciences
In common: JAX, Biopython, TensorFlow, 5 other tools
[8] doi:10.1038/s41586-026-10670-w [code]
Zero-shot design of drug-binding proteins via neural iterative selection-expansion.
Journal: Nature
In common: Plotly, PyTorch, pandas, 3 other tools, 3 references
[9] doi:10.1038/s41592-026-03057-2 [code]
CREsted: modeling genomic and synthetic cell-type-specific enhancers across tissues and species.
Journal: Nature methods
In common: Biopython, Keras, TensorFlow, 5 other tools
[10] doi:10.1093/nar/gkag706 [code]
scDifformer: diffusion-based post-training for virtual cell modeling across large-scale single-cell data.
Journal: Nucleic acids research
In common: Keras, TensorFlow, PyTorch, 4 other tools, 1 reference

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.