Including saccharides

molsysmt/_private/exceptions/argument_error.py

@@ -29,7 +29,7 @@ class ArgumentError(Exception):
     .. admonition:: See Also
        :class: attention
 
-        :ref:`Developer Guide \> Exceptions \> BadCallError <developer:exceptions:BadCallError>`
+        :ref:`Developer Guide > Exceptions > BadCallError <developer:exceptions:BadCallError>`
 
     """
 

molsysmt/build/get_missing_bonds.py

@@ -42,6 +42,7 @@ def get_missing_bonds(molecular_system, threshold='2 angstroms', selection='all'
         from molsysmt.element.group.amino_acid import get_bonded_atom_pairs as _bonds_in_amino_acid
         from molsysmt.element.group.terminal_capping import get_bonded_atom_pairs as _bonds_in_terminal_capping
         from molsysmt.element.group.small_molecule import get_bonded_atom_pairs as _bonds_in_small_molecule
+        from molsysmt.element.group.saccharide import get_bonded_atom_pairs as _bonds_in_saccharide
         from molsysmt.element.group.terminal_capping import is_n_terminal_capping, is_c_terminal_capping
 
         old_bonds = get(molecular_system, selection=selection, inner_bonded_atom_pairs=True)
@@ -121,7 +122,12 @@ def get_missing_bonds(molecular_system, threshold='2 angstroms', selection='all'
 
                 elif group_type=='saccharide':
 
-                    raise NotImplementedError('Group type "saccharide" not implemented')
+                    aux_bonds = _bonds_in_saccharide(group_name, atom_names, atom_indices, sorted=False)
+                    if aux_bonds is None:
+                        aux_bonds = _bonds_in_unknown_group(molecular_system, atom_indices, atom_names,
+                                                            structure_index=structure_index, threshold=threshold,
+                                                            sorted=False)
+                    bonds += aux_bonds
 
                 elif group_type=='oligosaccharide':
 

molsysmt/data/databases/saccharides/make_saccharides_db.py

@@ -0,0 +1,67 @@
+# wget https://files.wwpdb.org/pub/pdb/data/monomers/components.cif
+
+from molsysmt.native import CIFFileHandler
+import pickle
+
+handler = CIFFileHandler('../components.cif')
+cif_entries = handler.parse()
+
+types_dict = {}
+
+for ii in cif_entries:
+    jj = cif_entries[ii]['_chem_comp']['type']
+    if jj in types_dict:
+        types_dict[jj].append(ii)
+    else:
+        types_dict[jj] = [ii]
+
+# Adding only non-polymer and other
+types = ['saccharide', 'D-saccharide', '"L-saccharide, beta linking"', '"L-saccharide, alpha linking"',
+        'D-SACCHARIDE', 'L-saccharide', '"D-saccharide, beta linking"', '"D-saccharide, alpha linking"',
+        'L-SACCHARIDE', 'SACCHARIDE']
+
+data = []
+for ii in types:
+    data += types_dict[ii]
+
+output = {}
+
+for value in data:
+
+    key = cif_entries[value]
+
+    tmp_dict={}
+
+    tmp_dict['name']=key['_chem_comp']['name']
+    tmp_dict['topology']=[]
+    tmp_dict['topology'].append({'atoms':key['_chem_comp_atom']['atom_id'], 'bonds':[]})
+    tmp_dict['topology'].append({'atoms':key['_chem_comp_atom']['alt_atom_id'], 'bonds':[]})   
+    if '_chem_comp_bond' in key:
+        if not isinstance(key['_chem_comp_bond']['atom_id_1'],list):
+            key['_chem_comp_bond']['atom_id_1'] = [key['_chem_comp_bond']['atom_id_1']]
+            key['_chem_comp_bond']['atom_id_2'] = [key['_chem_comp_bond']['atom_id_2']]
+        for atom1,atom2 in zip(key['_chem_comp_bond']['atom_id_1'], key['_chem_comp_bond']['atom_id_2']):
+            tmp_dict['topology'][0]['bonds'].append([atom1,atom2])
+            ii = tmp_dict['topology'][0]['atoms'].index(atom1)
+            jj = tmp_dict['topology'][0]['atoms'].index(atom2)
+            atom1 = tmp_dict['topology'][1]['atoms'][ii]
+            atom2 = tmp_dict['topology'][1]['atoms'][jj]
+            tmp_dict['topology'][1]['bonds'].append([atom1,atom2])
+
+    output[value]=tmp_dict
+
+split_output = {}
+for name,value in output.items():
+    if name[0] not in split_output:
+        split_output[name[0]]={}
+    split_output[name[0]][name]=value
+
+import gzip
+
+for file_name, aux_output in split_output.items():
+    with gzip.open(file_name+'.pkl.gz', 'wb', compresslevel=9) as fff:
+        pickle.dump(aux_output, fff)
+
+with gzip.open('group_names.pkl.gz', 'wb', compresslevel=9) as fff:
+        pickle.dump(sorted(list(output.keys())), fff)
+

molsysmt/data/databases/small_molecules/make_small_molecules_db.py

@@ -15,13 +15,14 @@
     else:
         types_dict[jj] = [ii]
 
+# Adding only non-polymer and other
 types = ['NON-POLYMER', 'non-polymer', 'other']
 data = []
 for ii in types:
     data += types_dict[ii]
 data.remove('UNL')
 
-# quito iones
+# removing ions
 ions = []
 for value in data:
     key = cif_entries[value]
@@ -31,6 +32,13 @@
 for ion in ions:
     data.remove(ion)
 
+# removing saccharides
+types = ['saccharide', 'D-saccharide', '"L-saccharide, beta linking"', '"L-saccharide, alpha linking"',
+        'D-SACCHARIDE', 'L-saccharide', '"D-saccharide, beta linking"', '"D-saccharide, alpha linking"',
+        'L-SACCHARIDE', 'SACCHARIDE', '']
+
+
+
 output = {}
 
 for value in data:

molsysmt/element/atom/names.py

@@ -482,6 +482,9 @@
 'O2D'  : 'O',
 'O31'  : 'O',
 'O32'  : 'O',
+'O1P'  : 'O',
+'O2P'  : 'O',
+'O3P'  : 'O',
 'OT1'  : 'O',
 'OT2'  : 'O',
 'OC1'  : 'O',
@@ -512,6 +515,10 @@
 'CL'    :  'CL',
 'Cl'    :  'CL',
 'CLA'   :  'CL',
-'FE'   :  'FE'
+'FE'   :  'FE',
+
+## Unknown
+'UNK'   :  'X',
+
 }
 

molsysmt/element/component/get_component_type.py

@@ -93,7 +93,7 @@ def _get_component_type_from_group_names_and_types(group_names, group_types):
     last_group_type = group_types[-1]
     first_group_name = group_names[0]
 
-    if first_group_type in ['water', 'ion', 'small molecule', 'lipid']:
+    if first_group_type in ['water', 'ion', 'small molecule', 'lipid', 'saccharide']:
         tmp_type = first_group_type
     elif (first_group_type == 'amino acid') or (first_group_type == 'terminal capping'):
         if first_group_type == 'terminal capping':

molsysmt/element/group/get_bonded_atom_pairs.py

@@ -8,6 +8,7 @@ def get_bonded_atom_pairs(group_name, atom_names, atom_indices=None, sorted=True
     from .small_molecule import get_bonded_atom_pairs as get_bonded_atom_pairs_from_small_molecule
     from .terminal_capping import get_bonded_atom_pairs as get_bonded_atom_pairs_from_terminal_capping
     from .water import get_bonded_atom_pairs as get_bonded_atom_pairs_from_water
+    from .saccharide import get_bonded_atom_pairs as get_bonded_atom_pairs_from_saccharide
 
     group_type = get_group_type_from_group_name(group_name)
 
@@ -38,7 +39,7 @@ def get_bonded_atom_pairs(group_name, atom_names, atom_indices=None, sorted=True
 
         case 'saccharide':
 
-            pass
+            bonds = get_bonded_atom_pairs_from_saccharide(group_name, atom_names, atom_indices=atom_indices, sorted=sorted)
 
         case 'small molecule':
 

molsysmt/element/group/saccharide/__init__.py

@@ -1,2 +1,4 @@
-from .saccharide_names import saccharide_names
 from .is_saccharide import is_saccharide
+from .group_names import group_names
+from .get_group_db import get_group_db
+from .get_bonded_atom_pairs import get_bonded_atom_pairs

molsysmt/element/group/saccharide/get_bonded_atom_pairs.py

@@ -0,0 +1,45 @@
+import numpy as np
+
+_sorted=sorted
+
+def get_bonded_atom_pairs(group_name, atom_names, atom_indices=None, sorted=True):
+
+    from . import group_names, get_group_db
+
+    if group_name not in group_names:
+        raise ValueError
+
+    if atom_indices is None:
+        atom_indices = np.arange(len(atom_names), dtype=int).tolist()
+
+    db = get_group_db(group_name)
+
+    is_in = -1
+    for ii,jj in enumerate(db['topology']):
+        if np.all(np.isin(atom_names, jj['atoms'])):
+            is_in=ii
+            break
+
+    if is_in!=-1:
+
+        bonds = []
+        for ii,jj in db['topology'][is_in]['bonds']:
+            if ii in atom_names:
+                if jj in atom_names:
+                    iii = atom_indices[atom_names.index(ii)]
+                    jjj = atom_indices[atom_names.index(jj)]
+                    if iii<jjj:
+                        bonds.append([iii,jjj])
+                    else:
+                        bonds.append([jjj,iii])
+
+    else:
+
+        print(group_name, atom_names)
+
+        raise ValueError
+
+    if sorted:
+        return _sorted(bonds)
+    else:
+        return bonds

molsysmt/element/group/saccharide/get_group_db.py

@@ -0,0 +1,22 @@
+import pickle
+import sys
+import gzip
+import numpy as np
+from molsysmt.element.group.saccharide import group_names
+
+from importlib.resources import files
+def path(package, file):
+    return files(package).joinpath(file)
+
+
+def get_group_db(group_name):
+
+    if group_name not in group_names:
+        raise ValueError
+
+    with gzip.open(path('molsysmt.data.databases.saccharides',group_name[0]+'.pkl.gz'), 'rb') as fff:
+        dbs = pickle.load(fff)
+
+    db = dbs[group_name]
+
+    return db

molsysmt/element/group/saccharide/group_names.py

@@ -0,0 +1,15 @@
+import pickle
+import sys
+import gzip
+
+from importlib.resources import files
+def path(package, file):
+    return files(package).joinpath(file)
+
+try:
+    with gzip.open(path('molsysmt.data.databases.saccharides','group_names.pkl.gz'), 'rb') as fff:
+        group_names = pickle.load(fff)
+except:
+    group_names = None
+    print('The file molsysmt.data.databases.saccharides.group_names.pkl.gz was not loaded.')
+

molsysmt/element/group/saccharide/is_saccharide.py

@@ -1,6 +1,9 @@
-from .saccharide_names import saccharide_names
+from .group_names import group_names
 
 def is_saccharide(name):
+    """
+    To be written soon...
+    """
 
-    return (name in saccharide_names)
+    return (name in group_names)
 

molsysmt/element/molecule/__init__.py

@@ -25,6 +25,7 @@
         'dna',
         'rna',
         'lipid',
+        'saccharide',
         'oligosaccharide'
         ]
 
@@ -37,6 +38,7 @@
     'dna': 'dnas',
     'rna': 'rnas',
     'lipid': 'lipids',
+    'saccharide': 'saccharides',
     'oligosaccharide': 'oligosaccharides',
 }
 
@@ -49,6 +51,7 @@
     'dnas': 'dna',
     'rnas': 'rna',
     'lipids': 'lipid',
+    'saccharides': 'saccharide',
     'oligosaccharides': 'oligosaccharide',
 }