Tools

PDB class

class gromacs_py.gromacs.tools.pdb_manip.Coor

Topologie base on coordinates like pdb or gro.

The coor object containt a dictionnary of atoms indexed on the atom num and the crystal packing info.

Parameters:

• atom_dict (dict) – dictionnary of atom
• crystal_pack (str) – crystal packing

Atom dictionnary parameters

Parameters:

• field (str) – pdb field
• num (int) – atom number
• name (str) – atom name
• alter_loc (str) – atom number
• res_name (str) – residue name (3 letters)
• chain (str) – chain ID
• res_num (int) – residue number (based on pdb file)
• uniq_resid (int) – unique residue number
• insert_res (str) – atom number
• xyz – coordinate
• occ (float) – occupation
• beta (float) – beta flactor

Note

The atom num index in the dictionnary, is not the same as the atom_num field of the dictionnary.

Note

Files necessary for testing : ../test/input/1y0m.pdb, ../test/input/1rxz.pdb and ../test/input/4n1m.pdb. To do the unitary test, execute pdb_mani.py (-v for verbose mode)

Todo

Add an atom class ?

add_zinc_finger(ZN_pdb, cutoff=3.2)

Change protonation state of cysteins and histidine coordinating Zinc atoms. To do after correct_his_name and correct_cys_name, in order that protonation is recognize by pdb2gmx.

Example:

>>> try: #doctest: +ELLIPSIS
...   print("Start import")
...   from . import pdb2pqr
... except ImportError:
...   import pdb2pqr
Start import...
>>> TEST_OUT = str(getfixture('tmpdir'))
>>>
>>> # Read the pdb 1jd4 and keep only chain A
>>> input_pdb = Coor()
>>> input_pdb.read_pdb(os.path.join(TEST_PATH, '1jd4.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1jd4.pdb ,  1586 atoms found
>>> chain_A = input_pdb.select_part_dict(selec_dict = {'chain' : ['A']})
>>> chain_A.write_pdb(os.path.join(TEST_OUT, '1jd4_A.pdb')) #doctest: +ELLIPSIS
Succeed to save file .../1jd4_A.pdb
>>>
>>> # Compute protonation with pdb2pqr:
>>> pdb2pqr.compute_pdb2pqr(os.path.join(TEST_OUT, '1jd4_A.pdb'), os.path.join(TEST_OUT, '1jd4.pqr')) #doctest: +ELLIPSIS
Succeed to read file .../1jd4_A.pdb ,  793 atoms found
Succeed to save file .../tmp_pdb2pqr.pdb
pdb2pqr.py --ff CHARMM --ffout CHARMM --chain .../tmp_pdb2pqr.pdb .../1jd4.pqr
0
>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_OUT, '1jd4.pqr'), pqr_format = True)
Succeed to read file .../1jd4.pqr ,  1548 atoms found
>>> prot_coor.correct_cys_name() #doctest: +ELLIPSIS
<...Coor object at 0x...
>>> prot_coor.correct_his_name() #doctest: +ELLIPSIS
<...Coor object at 0x...
>>> prot_coor.correct_chain() #doctest: +ELLIPSIS
Chain: A  Residue: 0 to 95
<...Coor object at 0x...
>>> ZN_index = prot_coor.get_index_selection({'name':['ZN']})
>>> print(len(ZN_index))
0
>>> prot_coor.add_zinc_finger(os.path.join(TEST_OUT, '1jd4_A.pdb')) #doctest: +ELLIPSIS
Succeed to read file .../1jd4_A.pdb ,  793 atoms found
Presence of 1 Zinc detected
change cystein residue(s) : [48, 51, 75]
change histidine residue(s) : [68]
True
>>> ZN_index = prot_coor.get_index_selection({'name':['ZN']})
>>> print(len(ZN_index))
1

Note

This function seems useless. Since last version of pdb2pqr residue name seems correct.

static atom_dist(atom_1, atom_2)

Compute the distance between 2 atoms.

Parameters:

• atom_1 (dict) – atom dictionnary
• atom_2 (dict) – atom dictionnary

Returns:

distance

Return type:

float

Example:

>>> atom_1 = {'xyz': np.array([0.0, 0.0, 0.0])}
>>> atom_2 = {'xyz': np.array([0.0, 1.0, 0.0])}
>>> atom_3 = {'xyz': np.array([1.0, 1.0, 1.0])}
>>> Coor.atom_dist(atom_1, atom_2)
1.0
>>> Coor.atom_dist(atom_1, atom_3)
1.7320508075688772

center_of_mass(selec_dict={})

Compute the center of mass of a selection Avoid using atoms with 2 letters atom name like NA Cl …

Parameters:selec_dict (dict, default={}) – selection dictionnary Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> com_1y0m = prot_coor.center_of_mass()
>>> print("x:{:.2f} y:{:.2f} z:{:.2f}".format(*com_1y0m))
x:16.01 y:0.45 z:8.57

Warning

Atom name must start with its type letter (H, C, N, O, P, S).

change_index_pdb_field(index_list, change_dict)

Change all atom field of a part of coor object defined by index, the change is based on the change_dict dictionnary.

Parameters:

• index_list (list) – list of atom index to change
• change_dict (dict) – change ditionnay eg. {“chain” : “A”}

Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> res_826_852 = prot_coor.get_index_selection({'res_num' : range(826,852)})
>>> prot_coor.change_index_pdb_field(index_list = res_826_852, change_dict = {"chain" : "B"}) #doctest: +ELLIPSIS
<...Coor object at ...>
>>> prot_seq = prot_coor.get_aa_seq()
>>> prot_seq == {'A': 'TFKSAVKALFDYKAQREDELTFTKSAIIQNVEKQD', 'B': 'GGWWRGDYGGKKQLWFPSNYVEEMIN'}
True

change_pdb_field(change_dict)

Change all atom field of a coor object, the change is based on the change_dict dictionnary.

Parameters:change_dict (dict) – change ditionnay eg. {“chain” : “A”} Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> prot_coor.get_aa_seq()
{'A': 'TFKSAVKALFDYKAQREDELTFTKSAIIQNVEKQDGGWWRGDYGGKKQLWFPSNYVEEMIN'}
>>> prot_coor.change_pdb_field(change_dict = {"chain" : "B"}) #doctest: +ELLIPSIS
<...Coor object at ...>
>>> prot_coor.get_aa_seq()
{'B': 'TFKSAVKALFDYKAQREDELTFTKSAIIQNVEKQDGGWWRGDYGGKKQLWFPSNYVEEMIN'}

static concat_pdb(*pdb_in_files, pdb_out)

Concat a list of pdb files in one.

Parameters:

• pdb_in_files (list) – list of pdb files
• pdb_out (dict) – atom dictionnary

Example:

>>> TEST_OUT = str(getfixture('tmpdir'))
>>> Coor.concat_pdb(os.path.join(TEST_PATH, '1y0m.pdb'),
...                 os.path.join(TEST_PATH, '1rxz.pdb'),
...                 pdb_out = os.path.join(TEST_OUT, 'tmp_2.pdb')) #doctest: +ELLIPSIS
Succeed to save concat file: .../tmp_2.pdb

correct_chain(Ca_cutoff=4.5)

Correct the chain ID’s of a coor object, by checking consecutive Calphas atoms distance. If the distance is higher than Ca_cutoff, the former atoms are considered as in a different chain.

Parameters:Ca_cutoff (float, default=4.5) – cutoff for distances between Calphas atoms (X) Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> res_810 = prot_coor.get_index_selection({'res_num' : [810]})
>>> prot_coor = prot_coor.del_atom_index(index_list = res_810)
>>> prot_coor.get_aa_seq()
{'A': 'TFKSAVKALFDYKAQREDETFTKSAIIQNVEKQDGGWWRGDYGGKKQLWFPSNYVEEMIN'}
>>> prot_coor.correct_chain() #doctest: +ELLIPSIS
Chain: A  Residue: 0 to 18
Chain: B  Residue: 20 to 60
<...Coor object at ...>
>>> # As a residue is missing, Calphas after residue 18 is no more consecutive

Note

This is specially usefull for pdb2gmx which cut the protein chains based on the chain ID’s.

correct_cys_name()

Correct the CYS resname from pdb2pqr

Example:

>>> try: #doctest: +ELLIPSIS
...   print("Start import")
...   from . import pdb2pqr
... except ImportError:
...   import pdb2pqr
Start import...
>>> TEST_OUT = str(getfixture('tmpdir'))
>>> # Compute protonation with pdb2pqr:
>>> pdb2pqr.compute_pdb2pqr(os.path.join(TEST_PATH, '1dpx.pdb'), os.path.join(TEST_OUT, '1dpx.pqr')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1dpx.pdb ,  1192 atoms found
Succeed to save file .../tmp_pdb2pqr.pdb
pdb2pqr.py --ff CHARMM --ffout CHARMM --chain .../tmp_pdb2pqr.pdb .../1dpx.pqr
0
>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_OUT, '1dpx.pqr'), pqr_format = True)
Succeed to read file .../1dpx.pqr ,  1960 atoms found
>>> Isu_index = prot_coor.get_index_selection({'res_name' : ['DISU']})
>>> print(len(Isu_index))
16
>>> prot_coor.correct_cys_name() #doctest: +ELLIPSIS
<...Coor object at 0x...
>>> Isu_index = prot_coor.get_index_selection({'res_name' : ['DISU']})
>>> print(len(Isu_index))
0

correct_his_name()

Get his protonation state from pdb2pqr and replace HIS resname with HSE, HSD, HSP resname. To do after pdb2pqr, in order that protonation is recognize by pdb2gmx.

Example:

>>> try: #doctest: +ELLIPSIS
...   print("Start import")
...   from . import pdb2pqr
... except ImportError:
...   import pdb2pqr
Start import...
>>> TEST_OUT = str(getfixture('tmpdir'))
>>> # Compute protonation with pdb2pqr:
>>> pdb2pqr.compute_pdb2pqr(os.path.join(TEST_PATH, '4n1m.pdb'), os.path.join(TEST_OUT, '4n1m.pqr')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/4n1m.pdb ,  2530 atoms found
Succeed to save file .../tmp_pdb2pqr.pdb
pdb2pqr.py --ff CHARMM --ffout CHARMM --chain .../tmp_pdb2pqr.pdb .../4n1m.pqr
0
>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_OUT, '4n1m.pqr'), pqr_format = True) #doctest: +ELLIPSIS
Succeed to read file .../4n1m.pqr ,  2548 atoms found
>>> HSD_index = prot_coor.get_index_selection({'res_name' : ['HSD'], 'name':['CA']})
>>> print(len(HSD_index))
5
>>> HSE_index = prot_coor.get_index_selection({'res_name' : ['HSE'], 'name':['CA']})
>>> print(len(HSE_index))
0
>>> HSP_index = prot_coor.get_index_selection({'res_name' : ['HSP'], 'name':['CA']})
>>> print(len(HSP_index))
0
>>> prot_coor.correct_his_name() #doctest: +ELLIPSIS
<...Coor object at 0x...
>>> HIS_index = prot_coor.get_index_selection({'res_name' : ['HIS'], 'name':['CA']})
>>> print(len(HIS_index))
0

Note

This function seems useless. Since last version of pdb2pqr residue name seems correct.

correct_water_name()

Correct the water resname from pdb2pqr

Example:

>>> try: #doctest: +ELLIPSIS
...   print("Start import")
...   from . import pdb2pqr
... except ImportError:
...   import pdb2pqr
Start import...
>>> TEST_OUT = str(getfixture('tmpdir'))
>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1dpx.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1dpx.pdb ,  1192 atoms found
>>> prot_coor.water_to_ATOM() #doctest: +ELLIPSIS
<...Coor object at 0x...
>>> prot_coor.write_pdb(os.path.join(TEST_OUT, '1dpx_water.pdb')) #doctest: +ELLIPSIS
Succeed to save file .../1dpx_water.pdb
>>> # Compute protonation with pdb2pqr:
>>> pdb2pqr.compute_pdb2pqr(os.path.join(TEST_OUT, '1dpx_water.pdb'), os.path.join(TEST_OUT, '1dpx_water.pqr')) #doctest: +ELLIPSIS
Succeed to read file .../1dpx_water.pdb ,  1192 atoms found
Succeed to save file .../tmp_pdb2pqr.pdb
pdb2pqr.py --ff CHARMM --ffout CHARMM --chain .../tmp_pdb2pqr.pdb .../1dpx_water.pqr
0
>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_OUT, '1dpx_water.pqr'), pqr_format = True) #doctest: +ELLIPSIS
Succeed to read file .../1dpx_water.pqr ,  2491 atoms found
>>> water_index = prot_coor.get_index_selection({'res_name':['TP3M'], 'name':['OH2']})
>>> print(len(water_index))
177

del_atom_index(index_list)

Delete atoms of a coor object defined by their index.

Parameters:index_list (list) – list of atom index to delete Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> prot_coor.get_aa_seq()
{'A': 'TFKSAVKALFDYKAQREDELTFTKSAIIQNVEKQDGGWWRGDYGGKKQLWFPSNYVEEMIN'}
>>> res_810_852 = prot_coor.get_index_selection({'res_num' : range(810,852)})
>>> prot_coor.del_atom_index(index_list = res_810_852) #doctest: +ELLIPSIS
<...Coor object at ...>
>>> prot_coor.get_aa_seq()
{'A': 'TFKSAVKALFDYKAQREDE'}

dist_under_index(atom_sel_2, cutoff=10.0)

Check is distance between atoms of self.coor is under cutoff with atoms of group 1. Then return list of index of atoms of self.coor under ctuoff ditance.

Parameters:

• atom_sel_1 (dict) – atom dictionnary
• atom_sel_2 (dict) – atom dictionnary
• cutoff (float, default=10.0) – distance cutoff

get_aa_num()

Get the amino acid number of a coor object.

Returns:Number of residues Return type:int Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> prot_coor.get_aa_num()
61

Note

Only count Ca atoms, this may not be the best choice ?

get_aa_seq()

Get the amino acid sequence from a coor object.

Returns:dictionnary of sequence indexed by the chain ID Return type:dict Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> prot_coor.get_aa_seq()
{'A': 'TFKSAVKALFDYKAQREDELTFTKSAIIQNVEKQDGGWWRGDYGGKKQLWFPSNYVEEMIN'}

Warning

If atom chains are not arranged sequentialy (A,A,A,B,B,A,A,A …), the first atom seq will be overwritten by the last one.

get_attribute_selection(selec_dict={}, attribute='uniq_resid', index_list=None)

Select atom of a coor object based on the change_dict dictionnary. Return the list of unique attribute of the selected atoms.

Parameters:selec_dict (dict) – select ditionnay eg. {“chain” : [“A”,”G”]} Returns:list of atom index Return type:list of int Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> prot_coor.get_attribute_selection({'res_num' : [826,827]}, attribute='uniq_resid')
[35, 36]

get_index_dist_between(atom_sel_2, cutoff_min=0, cutoff_max=10)

Check is distance between atoms of self.atom_dict is under cutoff with the atoms of group 1. Then return list of index of atoms of self.coor under cutoff ditance.

Parameters:

• atom_sel_1 (dict) – atom dictionnary
• atom_sel_2 (dict) – atom dictionnary
• cutoff_min (float, default=0.0) – maximum distance cutoff
• cutoff_max (float, default=10.0) – minimum distance cutoff

Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> res_810 = prot_coor.select_part_dict({'res_num' : [810]})
>>> close_r810 = prot_coor.get_index_dist_between(res_810, cutoff_min = 3, cutoff_max = 5)
>>> print(len(close_r810))
65

get_index_selection(selec_dict)

Select atom of a coor object based on the change_dict dictionnary. Return the list of index of selected atoms.

Parameters:selec_dict (dict) – select ditionnay eg. {“chain” : [“A”,”G”]} Returns:list of atom index Return type:list of int Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> prot_coor.get_index_selection({'res_num' : [826,827]})
[297, 298, 299, 300, 301, 302, 303, 304]

insert_mol(pdb_out, out_folder, mol_chain, check_file_out=True)

Insert molecules defined by chain ID mol_chain in a water solvant. Check which water molecules are within cutoff_prot_off=12.0 X and cutoff_prot_in=15.0 X of protein and peptide C alpha atoms. Move the molecules to be inserted at the position of water molecules. Then delete all water molecules within cutoff_water_clash=1.2 X of the inserted molecule atoms.

Parameters:

• pdb_out (str) – name of output pdb file
• out_folder (str) – path of the ouput directory
• mol_chain (str) – chain ID of the molecule to be inserted,
• check_file_out (bool, optional, default=True) – flag to check or not if file has already been created. If the file is present then the command break.

Warning

self.atom_dict file must contain alredy a concatenated system with a ligand (chain: mol_chain) and a solvated system.

make_peptide(sequence, pdb_out, check_file_out=True)

Create a linear peptide structure.

Parameters:

• sequence (str) – peptide sequence
• pdb_out (str) – name of output pdb file
• check_file_out (bool, optional, default=True) – flag to check or not if file has already been created. If the file is present then the command break.

read_pdb(pdb_in, pqr_format=False)

Read a pdb file and return atom informations as a dictionnary indexed on the atom num. The fonction can also read pqr files if specified with pqr_format = True, it will only change the column format of beta and occ factors.

Parameters:

• pdb_in (str) – path of the pdb file to read
• pqr_format (bool, default=False) – Flag for .pqr file format reading.

Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found

select_part_dict(selec_dict)

Select atom of a coor object defined, the selection is based on the change_dict dictionnary. Return a new coor object.

Parameters:selec_dict (dict) – change ditionnay eg. {“chain” : [“A”,”G”]} Returns:a new coor object Return type:coor Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> prot_coor.get_aa_num()
61
>>> prot_20_coor = prot_coor.select_part_dict(selec_dict = {'res_num' : list(range(791,800))})
>>> prot_20_coor.get_aa_seq()
{'A': 'TFKSAVKAL'}
>>> prot_20_coor.get_aa_num()
9
>>> prot_N_atom = prot_coor.select_part_dict(selec_dict = {'name' : ['ZN']})
>>> # WARNING using selec_dict = {'name' : 'ZN'} will give you 61 residues !!
>>> print(len(prot_N_atom.atom_dict))
0

translate(vector)

Translate all atoms of a coor object by a given vector

Parameters:vector (list) – 3d translation vector Example:

>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> com_1y0m = prot_coor.center_of_mass()
>>> print("x:{:.2f} y:{:.2f} z:{:.2f}".format(*com_1y0m))
x:16.01 y:0.45 z:8.57
>>> prot_coor.translate(-com_1y0m)
>>> com_1y0m = prot_coor.center_of_mass()
>>> print("x:{:.2f} y:{:.2f} z:{:.2f}".format(*com_1y0m))
x:-0.00 y:0.00 z:0.00

water_to_ATOM()

Change HETATM field of water to ATOM, as pdb2pqr only use ATOM field.

Example:

>>> try: #doctest: +ELLIPSIS
...   print("Start import")
...   from . import pdb2pqr
... except ImportError:
...   import pdb2pqr
Start import...
>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1dpx.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1dpx.pdb ,  1192 atoms found
>>> hetatm_index = prot_coor.get_index_selection({'field':['HETATM']})
>>> print(len(hetatm_index))
179
>>> prot_coor.water_to_ATOM() #doctest: +ELLIPSIS
<...Coor object at 0x...
>>> hetatm_index = prot_coor.get_index_selection({'field':['HETATM']})
>>> print(len(hetatm_index))
2
>>> water_index = prot_coor.get_index_selection({'res_name':['HOH']})
>>> print(len(water_index))
177

write_pdb(pdb_out, check_file_out=True)

Write a pdb file.

Parameters:pdb_out (str) – path of the pdb file to write Example:

>>> TEST_OUT = str(getfixture('tmpdir'))
>>> prot_coor = Coor()
>>> prot_coor.read_pdb(os.path.join(TEST_PATH, '1y0m.pdb')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/1y0m.pdb ,  648 atoms found
>>> prot_coor.write_pdb(os.path.join(TEST_OUT, 'tmp.pdb')) #doctest: +ELLIPSIS
Succeed to save file .../tmp.pdb

PDB2PQR function

######### PDB2PQR ##########

gromacs_py.gromacs.tools.pdb2pqr.compute_pdb2pqr(pdb_in, pdb_out, ff='CHARMM', check_file_out=True)

Use pdb2pqr to define protonation state of each residue of a protein.

Parameters:

• pdb_in (str) – path of input pdb file
• pdb_out (str) – path of output pdb file
• ff (str, optional, default="CHARMM") – forcefield nomenclature for atom names
• check_file_out (bool, optional, default=True) – flag to check or not if file has already been created. If the file is present then the command break.

Example:

>>> TEST_OUT = str(getfixture('tmpdir'))
>>> # Compute protonation with pdb2pqr:
>>> compute_pdb2pqr(os.path.join(TEST_PATH,'4n1m.pdb'), os.path.join(TEST_OUT, '4n1m.pqr')) #doctest: +ELLIPSIS
Succeed to read file ...test/input/4n1m.pdb ,  2530 atoms found
Succeed to save file .../tmp_pdb2pqr.pdb
pdb2pqr.py --ff CHARMM --ffout CHARMM --chain .../tmp_pdb2pqr.pdb .../4n1m.pqr
0
>>> prot_coor = pdb_manip.Coor()
>>> prot_coor.read_pdb(TEST_OUT+'/4n1m.pqr', pqr_format = True)
Succeed to read file .../4n1m.pqr ,  2548 atoms found
>>> HSD_index = prot_coor.get_index_selection({'res_name' : ['HSD'], 'name':['CA']})
>>> print(len(HSD_index))
5
>>> HSE_index = prot_coor.get_index_selection({'res_name' : ['HSE'], 'name':['CA']})
>>> print(len(HSE_index))
0
>>> HSP_index = prot_coor.get_index_selection({'res_name' : ['HSP'], 'name':['CA']})
>>> print(len(HSP_index))
0

Note

Idealy I would need a pdb file with 3 different histidine protonation. I couldn’t find one.

OS commands

Collection of function related to os and sys operations.

class gromacs_py.gromacs.tools.os_command.Command(list_cmd, my_env=None, **kwargs)

The Command class is a way to launch bash command and mainly gromacs

Parameters:

• cmd (list) – command list
• env (dict) – environment variable

Example:

>>> cmd_list = ['ls','-a',TEST_PATH]
>>> cmd_test = Command(list_cmd=cmd_list)
>>> cmd_test.display() #doctest: +ELLIPSIS
ls -a ...test/input
>>> return_code = cmd_test.run(out_data=True)
>>> print(return_code['stdout']) #doctest: +ELLIPSIS
.
..
1AWR-055_bestene1-mc.pdb
1dpx.pdb
1jd4.pdb
...
4n1m.pdb
...
<BLANKLINE>

define_env(my_env)

Define the environment of the Command object.

display()

Show Command object that will be launch. Show only the name of the command (eg. gmx) instead of the full path. Show relative path for files in the command.

display_raw()

Show Command object that will be launch. Show the full path of the command as well as the full path for files in the command.

run(com_input='', display=False, out_data=False)

Launch Command object that will be launch. return programm output is out_data is set to True

Parameters:

• com_input (str) – input for the command
• display (bool) – option to display output
• out_data (bool) – option to return output data

Returns:

Return Code or output dict

Return type:

bool/dict

run_background(func_input_dict, com_input='', display=False, out_data=False)

Launch Command object that will be launch. Will the command is running launch the function using func_input_dict as argument. return programm output is out_data is set to True

Parameters:

• function (function) – function to be launch
• func_input_dict (dict) – input for the function
• com_input (str) – input for the command
• display (bool) – option to display output
• out_data (bool) – option to return output data

Returns:

Return Code or output dict

Return type:

bool/dict

gromacs_py.gromacs.tools.os_command.check_directory_exist(directory)

Check is a directory exist.

Parameters:directory (str) – directory path Returns:if the file exist Return type:bool Example:

>>> test_exist = check_directory_exist(TEST_PATH)
>>> print("Directory {} exist: {}".format(TEST_PATH, test_exist)) #doctest: +ELLIPSIS
Directory ...test/input exist: True
>>> test_exist = check_directory_exist(os.path.join(TEST_PATH,'no_way'))
>>> print("Directory {} exist: {}".format(TEST_PATH+'/no_way', test_exist)) #doctest: +ELLIPSIS
Directory ...test/input/no_way exist: False

gromacs_py.gromacs.tools.os_command.check_file_and_create_path(file)

Check if file exist and create path if not available

Parameters:file (str) – file path Returns:File existance Return type:bool

gromacs_py.gromacs.tools.os_command.check_file_exist(file)

Check is a file exist.

Parameters:file (str) – file path Returns:if the file exist Return type:bool Example:

>>> test_exist = check_file_exist(os.path.join(TEST_PATH, "1y0m.pdb"))
>>> print("1y0m.pdb exist: ", test_exist)
1y0m.pdb exist:  True

gromacs_py.gromacs.tools.os_command.create_and_go_dir(dir_name)

Create the path to a directory and change path in it.

Parameters:dir_name – directorie name Example:

>>> TEST_OUT = str(getfixture('tmpdir'))
>>> start_dir = os.getcwd()
>>> create_and_go_dir(os.path.join(TEST_OUT, "tmp"))
>>> print("Path: ", os.getcwd()) #doctest: +ELLIPSIS
Path: .../tmp
>>> os.chdir(start_dir)

gromacs_py.gromacs.tools.os_command.create_dir(dir_name)

Create the path to a directory.

Parameters:dir_name – directorie name

gromacs_py.gromacs.tools.os_command.delete_directory(directory)

Delete a file.

Parameters:directory (str) – directory path Returns:operation sucess Return type:bool

gromacs_py.gromacs.tools.os_command.delete_file(file)

Delete a file.

Parameters:file (str) – file path Returns:operation sucess Return type:bool

gromacs_py.gromacs.tools.os_command.full_path_and_check(file)

Return the full path of a file

Parameters:file (str) – file path Returns:File path Return type:str

gromacs_py.gromacs.tools.os_command.get_directory(file)

Return the path of a file directory

Parameters:file (str) – file path Returns:File path Return type:str

gromacs_py.gromacs.tools.os_command.get_gmx_version()

Get gromacs version of mdrun.

Example:

>>> print('Version is {}'.format(get_gmx_version())) #doctest: +ELLIPSIS
Version is ...

gromacs_py.gromacs.tools.os_command.is_exe(fpath)

Check is a file path exist and if user has access to it

Parameters:fpath (str) – file path Returns:if the file is an executable Return type:bool Example:

>>> print(is_exe('/bin/ls'))
True

gromacs_py.gromacs.tools.os_command.which(*program_list)

find and return the path of a program Look for all combination within the $PATH env variable

Parameters:program_list (list of str) – list of program name Returns:path of the program Return type:str Example:

>>> ls_path = which('ls')
>>> print(ls_path)
/bin/ls

Monitor

Collection of function to monitor a simulation in real time.

gromacs_py.gromacs.tools.monitor.extract_log_dict(func_input_dict, tail_line_num=20)

Read the last lines of a gromacs .log file and return a dictionnary containing time, step and all energetic terms.

gromacs_py.gromacs.tools.monitor.isnotebook()

Return if the command is launch from a notebook or not Taken from: https://stackoverflow.com/questions/15411967/how-can-i-check-if-code-is-executed-in-the-ipython-notebook

gromacs_py.gromacs.tools.monitor.print_log_file(proc, func_input_dict, tail_line_num=20)

Monitor .log file information. The func_input_dict should contains several keys:

• terms: list of energetic terms to extract, eg. [‘Potential’, ‘Temperature’]
• log: path of the log file (Defined in os_command.run_background())
• refresh_time: time interval to refresh log extract (default=1.0 s)

Parameters:

• proc (subprocess object) – running subprocess
• func_input_dict (dict) – dictionnary containing parameters for log extract
• tail_line_num (int, default=20) – number of line to read at the end of .log file

Example:

>>> TEST_OUT = str(getfixture('tmpdir'))
>>> import sys
>>> #print(os.path.abspath(os.path.join(os.path.dirname(__file__), '../../..')))
>>> sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '../..')))
>>> import gromacs.gmx5 as gmx #doctest: +ELLIPSIS
Gromacs version is ...
FORCEFIELD_PATH = ...
>>> prot = gmx.GmxSys(name='1y0m', coor_file=TEST_PATH+'/1y0m.pdb')
>>> ###################################
>>> ####   Create the topologie:   ###
>>> ###################################
>>> prot.prepare_top(out_folder=os.path.join(TEST_OUT, 'top_SH3')) #doctest: +ELLIPSIS
Succeed to read file .../test/input/1y0m.pdb ,  648 atoms found
Succeed to save file tmp_pdb2pqr.pdb
pdb2pqr.py --ff CHARMM --ffout CHARMM --chain tmp_pdb2pqr.pdb 00_1y0m.pqr
Succeed to read file 00_1y0m.pqr ,  996 atoms found
Chain: A  Residue: 0 to 60
Succeed to save file 01_1y0m_good_his.pdb
-Create topologie
gmx pdb2gmx -f 01_1y0m_good_his.pdb -o 1y0m_pdb2gmx.pdb -p 1y0m_pdb2gmx.top -i 1y0m_posre.itp -water tip3p -ff charmm36-jul2017 -ignh -vsite hydrogens
Molecule topologie present in 1y0m_pdb2gmx.top , extract the topologie in a separate file: 1y0m_pdb2gmx.itp
Protein_chain_A
-ITP file: 1y0m_pdb2gmx.itp
-molecules defined in the itp file:
* Protein_chain_A
Rewrite topologie: 1y0m_pdb2gmx.top
>>> ######################################
>>> ### Monitor an energy minimisation ###
>>> ######################################
>>> monitor = {'function': print_log_file,           'terms':['Potential'],           'file_check_ext':'log'}
>>> prot.em(out_folder=os.path.join(TEST_OUT, 'em_SH3'), nsteps=100,    constraints='none', create_box_flag=True, monitor=monitor, nstlog=10)
-Create pbc box
gmx editconf -f .../top_SH3/1y0m_pdb2gmx.pdb -o .../top_SH3/1y0m_pdb2gmx_box.pdb -bt dodecahedron -d 1.0
-Create the tpr file  1y0m.tpr
gmx grompp -f 1y0m.mdp -c ../top_SH3/1y0m_pdb2gmx_box.pdb -r ../top_SH3/1y0m_pdb2gmx_box.pdb -p ../top_SH3/1y0m_pdb2gmx.top -po out_1y0m.mdp -o 1y0m.tpr -maxwarn 1
-Launch the simulation 1y0m.tpr
gmx mdrun -s 1y0m.tpr -deffnm 1y0m -nt 0 -ntmpi 0 -nsteps -2 -nocopyright

gromacs_py.gromacs.tools.monitor.simulation_plot(proc, func_input_dict, refresh_time=1.0)

This function is used for monitoring a simulation in real time. Function can be excecuted by the gromacs.tools.os_command.run_background() function. The function monitors a trajectory file, and launch the analysis if the file has been modified. It can plot as function of time an analysis of a simulation. Analysis is passed as input function.

Warning

Need to add the following lines to be run in jupyter notebook:

• %matplotlib notebook