Read non-standard residues from PDB

calvados/build.py

@@ -310,7 +310,47 @@ def geometry_from_pdb(pdb,use_com=False):
     else:
         cas = u.select_atoms('name CA')
         pos = cas.positions / 10.
-    return pos, u.dimensions
+    box = np.append(u.dimensions[:3]/10.,u.dimensions[3:])
+    return pos, box
+
+def geometry_from_pdb_rna(pdb,use_com=False):
+    """ positions in nm"""
+    backbone_atoms_name = [ "1H2'", "1H5'", "2H5'", "2HO'", "C1'",
+                            "C2'",   "C3'",  "C4'",  "C5'", "H1'",
+                            "H3'",   "H4'",  "O2'",  "O3'", "O5'",
+                            "O4'",   "OP1",  "OP2",  "P"  ]
+
+    with catch_warnings():
+        simplefilter("ignore")
+        u = Universe(pdb)
+    ag = u.atoms
+    ag.translate(-ag.center_of_mass())
+    if use_com:
+        pos = []
+        for res in u.residues:
+            backbone = res.atoms.select_atoms('name ' + ' '.join(backbone_atoms_name))
+            non_backbone = res.atoms.difference(backbone)
+            pos.append(backbone.center_of_mass())
+            pos.append(non_backbone.center_of_mass())
+        pos = np.array(pos) / 10.
+
+    else:
+        pos =  []
+        for res in u.residues:
+            ps = res.atoms.select_atoms('name P')
+            pos.append(ps.positions[0] / 10.)
+            if res.resname in ['U','C']:
+                ns = res.atoms.select_atoms('name N1')
+                pos.append(ns.positions[0] / 10.)
+            elif res.resname in ['A','G']:
+                ns = res.atoms.select_atoms('name N9')
+                pos.append(ns.positions[0] / 10.)
+            else:
+                raise ValueError(f"Invalid RNA resname, {res.resname}")
+        pos = np.array(pos)
+
+    box = np.append(u.dimensions[:3]/10.,u.dimensions[3:])
+    return pos, box
 
 def bfac_from_pdb(pdb,confidence=70.):
     """ get pLDDT encoded in pdb b-factor column """

calvados/data/default_component.yaml

@@ -13,15 +13,15 @@ pdb_folder : 'pdbs'
 restraint_type : 'harmonic'
 k_harmonic : 700.
 fdomains : 'domains.yaml'
-k_go : 10.
+k_go : 15.
 use_com : true
 periodic: false
 
 colabfold : 0
 bfac_shift : 0.8
 bfac_width : 50.
-pae_shift : 0.4
-pae_width : 8.
+pae_shift : 0.3
+pae_width : 15.
 
 rna_kb1 : 8033.0
 rna_kb2 : 8033.0

calvados/sequence.py

@@ -1,6 +1,8 @@
 import numpy as np
 
 from MDAnalysis import Universe
+from MDAnalysis.lib.util import convert_aa_code
+
 
 import random
 
@@ -10,6 +12,7 @@
 from Bio.Seq import Seq
 from Bio.SeqRecord import SeqRecord
 
+
 import os
 
 from localcider.sequenceParameters import SequenceParameters
@@ -52,7 +55,7 @@ def seq_from_pdb(pdb,selection='all',fmt='string'):
     res3 = ag.residues.resnames
     for res in res3:
         if len(res) == 3:
-            res1 = SeqUtils.seq1(res)
+            res1 = convert_aa_code(res)
         else:
             res1 = res
         if res1 == "":
@@ -136,6 +139,14 @@ def patch_terminal_qs(qs,n_termini,c_termini,loc='both'):
         qsnew[c_termini] -= qcoeff
     return qsnew
 
+def patch_terminal_mws(mws,n_termini,c_termini,loc='both'):
+    mwsnew = mws.copy()
+    if loc in ['N','both']:
+        mwsnew[n_termini] += 2
+    if loc in ['C','both']:
+        mwsnew[c_termini] += 16
+    return mwsnew
+
 def seq_dipole(seq):
     """ 1D charge dipole along seq """
     # print(seq)
@@ -747,7 +758,7 @@ def __init__(self,seq,residues=None,charge_termini=False,calc_dip=False,
 
             # q_intgrl_map = make_q_intgrl_map(residues)
             # self.q_ij = calc_q_ij(seq,q_intgrl_map)
-            
+
         if nu_file is not None:
             self.kappa = calc_kappa_manual(seq,residues=residues)
             if self.kappa == -1: # no charges
@@ -763,4 +774,4 @@ def __init__(self,seq,residues=None,charge_termini=False,calc_dip=False,
             feats_for_nu = [self.scd_noflex, self.shd, self.kappa, self.fcr, self.mean_lambda]
             model_nu = load(nu_file)
             X_nu = np.reshape(np.array(feats_for_nu),(1,-1))
-            self.nu_svr = model_nu.predict(X_nu)[0]
+            self.nu_svr = model_nu.predict(X_nu)[0]

examples/single_GO/README.md → examples/single_AF_CALVADOS/README.md

@@ -5,4 +5,4 @@ python prepare.py --name <protein name>
 python <protein name>/run.py --path <protein name>
 ```
 
-where `<protein name>` (hnRNPA1S or TIA1) is a protein for which a PDB file and the residues of the folded domains are provided in the `input` folder.
+where `<protein name>` (AF-P00698) is a protein for which a PDB file and a JSON PAE file are provided in the `input` folder.

examples/single_GO/prepare.py → examples/single_AF_CALVADOS/prepare.py

@@ -13,10 +13,10 @@
 sysname = f'{args.name:s}'
 
 # set the side length of the cubic box
-L = 25
+L = 50
 
 # set the saving interval (number of integration steps)
-N_save = 100
+N_save = 7000
 
 # set final number of frames to save
 N_frames = 1000
@@ -61,7 +61,7 @@
   restraint_type = 'go', # harmonic or go
   use_com = True, # apply on centers of mass instead of CA
   colabfold = 0, # PAE format (EBI AF=0, Colabfold=1&2)
-  k_go = 10., # Restraint force constant
+  k_go = 15., # Restraint force constant
 )
 components.add(name=args.name)
 

examples/single_dsRNA/input/domains.yaml

@@ -0,0 +1,3 @@
+dspolyR12:
+- [1,24]
+