I read a PDB with a CSO3 group into pybel (The example was edited in Maestro Schrodinger to obscure the original chemistry). Note that since I use a lot of non natural amino acids and other molecule building blocks the amino acid name is inaccurate but easily read by Pymol and Schrodinger. Schrodinger also adds implicit hydrogens without the error! here is an original example:
TITLE FullModel_2569-std-man
REMARK 4 COMPLIES WITH FORMAT V. 3.0, 1-DEC-2006
REMARK 888
REMARK 888 WRITTEN BY MAESTRO (A PRODUCT OF SCHRODINGER, LLC)
ATOM 1 N ALA X 1 17.799 25.136 52.626 1.00 16.34 N1+
ATOM 2 CA ALA X 1 16.591 25.195 51.764 1.00 16.34 C
ATOM 3 C ALA X 1 16.975 24.592 50.433 1.00 16.34 C
ATOM 4 O ALA X 1 17.063 23.378 50.330 1.00 16.34 O
ATOM 5 N ALA X 2 17.423 25.444 49.499 1.00 6.89 N
ATOM 6 CA ALA X 2 18.667 25.160 48.776 1.00 6.89 C
ATOM 7 C ALA X 2 19.807 25.237 49.806 1.00 6.89 C
ATOM 8 O ALA X 2 19.691 26.012 50.771 1.00 6.89 O
ATOM 9 N ALA X 3 20.783 24.323 49.735 1.00 10.36 N
ATOM 10 CA ALA X 3 21.710 24.047 50.849 1.00 10.36 C
ATOM 11 C ALA X 3 20.940 23.401 52.024 1.00 10.36 C
ATOM 12 O ALA X 3 19.722 23.200 51.950 1.00 10.36 O
ATOM 13 CB ALA X 3 22.871 23.144 50.393 1.00 10.36 C
ATOM 14 SG ALA X 3 24.414 23.586 51.244 1.00 20.00 S
ATOM 15 OD1 ALA X 3 24.591 24.986 50.847 1.00 20.00 O1-
ATOM 16 OD2 ALA X 3 24.087 23.405 52.666 1.00 20.00 O
ATOM 17 OD3 ALA X 3 25.366 22.624 50.685 1.00 20.00 O
ATOM 18 N ASP X 4 21.584 23.207 53.175 1.00 21.55 N
ATOM 19 CA ASP X 4 20.963 22.850 54.454 1.00 21.55 C
ATOM 20 C ASP X 4 19.873 23.860 54.847 1.00 21.55 C
ATOM 21 O ASP X 4 18.792 23.497 55.298 1.00 21.55 O
ATOM 22 N ARG X 5 20.063 25.132 54.464 1.00 29.63 N
ATOM 23 CA ARG X 5 19.068 26.199 54.538 1.00 29.63 C
ATOM 24 C ARG X 5 17.759 25.922 53.771 1.00 29.63 C
ATOM 25 O ARG X 5 16.712 26.441 54.157 1.00 29.63 O
TER 26 ARG X 5
CONECT 13 14
CONECT 14 13 15 16 17
CONECT 14 16 17
CONECT 15 14
CONECT 16 14
CONECT 16 14
CONECT 17 14
CONECT 17 14
END
I ingest the PDB using the function below [it has to be a PDB and not mol2 or sdf because the original is PDB]
def pep_to_formats(data_path, mol_name, add_hydrogen=True):
base = f"{data_path}/data/{mol_name}/{mol_name}"
pep_pdb = f"{base}.pdb"
std_pdb = f"{base}-std.pdb"
babel_pdbh = f"{base}-babelh.pdb"
babel_sdfh = f"{base}-babelh.sdf"
babel_mol2h = f"{base}-babelh.mol2"
babel_smileh = f"{base}-babelh.smile"
babel_smilesh = f"{base}-babelh.smiles"
##input pepticom chemical expansion pdb file output and output standard hydrogenated pdb file
babel_mol = pb.readfile(format="pdb", filename=std_pdb).__next__()
if add_hydrogen:
babel_mol.addh() # add Hs for 3D
# dd.make3D()
babel_mol.write(format='pdb', filename=babel_pdbh, overwrite=True)
babel_mol.write("sdf", babel_sdfh, True)
babel_mol.write("mol2", babel_mol2h, True)
babel_mol.write("smi", babel_smileh, True)
babel_mol.write("smiles", babel_smilesh, True)
return babel_mol, base
When I add hydrogens the output SDF adds hydrogens to one of the oxygens and the sulfur. clearly a mistake.
Here is the resulting SDF
/Users/gideonbar/devel/pep-chem/md_python/data/FullModel_2569-std-man/FullModel_2569-std-man-std.pdb
OpenBabel02022411333D
43 43 0 0 1 0 0 0 0 0999 V2000
17.7990 25.1360 52.6260 N 0 0 0 0 0 0 0 0 0 0 0 0
16.5910 25.1950 51.7640 C 0 0 0 0 0 0 0 0 0 0 0 0
16.9750 24.5920 50.4330 C 0 0 0 0 0 0 0 0 0 0 0 0
17.0630 23.3780 50.3300 O 0 0 0 0 0 0 0 0 0 0 0 0
17.4230 25.4440 49.4990 N 0 0 0 0 0 0 0 0 0 0 0 0
18.6670 25.1600 48.7760 C 0 0 0 0 0 0 0 0 0 0 0 0
19.8070 25.2370 49.8060 C 0 0 0 0 0 0 0 0 0 0 0 0
19.6910 26.0120 50.7710 O 0 0 0 0 0 0 0 0 0 0 0 0
20.7830 24.3230 49.7350 N 0 0 0 0 0 0 0 0 0 0 0 0
21.7100 24.0470 50.8490 C 0 0 1 0 0 0 0 0 0 0 0 0
20.9400 23.4010 52.0240 C 0 0 0 0 0 0 0 0 0 0 0 0
19.7220 23.2000 51.9500 O 0 0 0 0 0 0 0 0 0 0 0 0
22.8710 23.1440 50.3930 C 0 0 0 0 0 0 0 0 0 0 0 0
24.4140 23.5860 51.2440 S 0 0 0 0 0 0 0 0 0 0 0 0
24.5910 24.9860 50.8470 O 0 0 0 0 0 0 0 0 0 0 0 0
24.0870 23.4050 52.6660 O 0 5 0 0 0 0 0 0 0 0 0 0
25.3660 22.6240 50.6850 O 0 0 0 0 0 0 0 0 0 0 0 0
21.5840 23.2070 53.1750 N 0 0 0 0 0 0 0 0 0 0 0 0
20.9630 22.8500 54.4540 C 0 0 0 0 0 0 0 0 0 0 0 0
19.8730 23.8600 54.8470 C 0 0 0 0 0 0 0 0 0 0 0 0
18.7920 23.4970 55.2980 O 0 0 0 0 0 0 0 0 0 0 0 0
20.0630 25.1320 54.4640 N 0 0 0 0 0 0 0 0 0 0 0 0
19.0680 26.1990 54.5380 C 0 0 0 0 0 0 0 0 0 0 0 0
17.7590 25.9220 53.7710 C 0 0 0 0 0 0 0 0 0 0 0 0
16.7120 26.4410 54.1570 O 0 0 0 0 0 0 0 0 0 0 0 0
18.5772 24.5766 52.4008 H 0 0 0 0 0 0 0 0 0 0 0 0
16.2814 26.2108 51.6327 H 0 0 0 0 0 0 0 0 0 0 0 0
15.7766 24.6608 52.2070 H 0 0 0 0 0 0 0 0 0 0 0 0
16.9149 26.2630 49.2980 H 0 0 0 0 0 0 0 0 0 0 0 0
18.8183 25.8854 48.0041 H 0 0 0 0 0 0 0 0 0 0 0 0
18.6312 24.1954 48.3143 H 0 0 0 0 0 0 0 0 0 0 0 0
20.8868 23.8113 48.9004 H 0 0 0 0 0 0 0 0 0 0 0 0
22.1321 24.9731 51.1792 H 0 0 0 0 0 0 0 0 0 0 0 0
23.0086 23.2573 49.3380 H 0 0 0 0 0 0 0 0 0 0 0 0
22.6282 22.1299 50.6330 H 0 0 0 0 0 0 0 0 0 0 0 0
24.6600 24.5520 50.3188 H 0 0 0 0 0 0 0 0 0 0 0 0
25.3568 21.8133 51.2175 H 0 0 0 0 0 0 0 0 0 0 0 0
22.5627 23.3128 53.1652 H 0 0 0 0 0 0 0 0 0 0 0 0
21.7153 22.8355 55.2147 H 0 0 0 0 0 0 0 0 0 0 0 0
20.5093 21.8863 54.3517 H 0 0 0 0 0 0 0 0 0 0 0 0
20.9483 25.3660 54.1024 H 0 0 0 0 0 0 0 0 0 0 0 0
19.5055 27.0896 54.1376 H 0 0 0 0 0 0 0 0 0 0 0 0
18.7981 26.2859 55.5697 H 0 0 0 0 0 0 0 0 0 0 0 0
1 2 1 0 0 0 0
1 24 1 0 0 0 0
1 26 1 0 0 0 0
2 3 1 0 0 0 0
2 27 1 0 0 0 0
2 28 1 0 0 0 0
3 4 2 0 0 0 0
3 5 1 0 0 0 0
5 6 1 0 0 0 0
5 29 1 0 0 0 0
6 7 1 0 0 0 0
6 30 1 0 0 0 0
6 31 1 0 0 0 0
7 8 2 0 0 0 0
7 9 1 0 0 0 0
9 10 1 0 0 0 0
9 32 1 0 0 0 0
10 11 1 0 0 0 0
10 13 1 0 0 0 0
10 33 1 1 0 0 0
11 12 2 0 0 0 0
11 18 1 0 0 0 0
13 14 1 0 0 0 0
13 34 1 0 0 0 0
13 35 1 0 0 0 0
14 15 2 0 0 0 0
14 16 1 0 0 0 0
14 36 1 0 0 0 0
17 14 1 0 0 0 0
17 37 1 0 0 0 0
18 19 1 0 0 0 0
18 38 1 0 0 0 0
19 20 1 0 0 0 0
19 39 1 0 0 0 0
19 40 1 0 0 0 0
20 21 2 0 0 0 0
20 22 1 0 0 0 0
22 23 1 0 0 0 0
22 41 1 0 0 0 0
23 24 1 0 0 0 0
23 42 1 0 0 0 0
23 43 1 0 0 0 0
24 25 2 0 0 0 0
M CHG 1 16 -1
M END
$$$$
In my code I artificially add double bonds to two of the sulfur oxygen bonds in the PDB, because PDB doesn't capture resonance.
How can I prevent this or just delete the hydrogens. I tried
for neighbour_atom in ob.OBAtomAtomIter(at):
atom_type = neighbour_atom.GetType()
bond = at.GetBond(neighbour_atom)
neighbour_index = bond.GetNbrAtomIdx(at)
print(bond.GetLength())
# print(bond.GetBondOrder())
print(neighbour_index)
# print(bond.IsAromatic())
# print(bond.IsInRing())
# print(bond.IsRotor())
# print(bond.IsAmide())
print(atom_type)
if atom_type == 'H':
pass
# at.GetParent().DeleteAtom(neighbour_atom)
mol.OBMol.DeleteAtom(neighbour_atom)
The code freezes
I converted it to a script and received this error:
Process finished with exit code 139 (interrupted by signal 11:SIGSEGV)
I found a way to remove hydrogen bonds after adding hydrogens:
I will be glad to see if there is a way to set ip the input PDB more correctly or add the hydrogens with better precision.