如何防止或减轻 pybel 错误地将氢添加到 CSO3 基团

问题描述 投票：0回答：1

我将带有 CSO3 基团的 PDB 读入 pybel（该示例在 Maestro Schrodinger 中进行了编辑，以掩盖原始化学成分）。请注意，由于我使用了大量非天然氨基酸和其他分子构建块，因此氨基酸名称不准确，但很容易被 Pymol 和 Schrodinger 读取。薛定谔还添加了隐式氢而没有错误！这是一个原始示例：

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

我使用下面的函数摄取 PDB [它必须是 PDB 而不是 mol2 或 sdf，因为原始是 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

当我添加氢气时，输出 SDF 会将氢气添加到氧和硫之一中。显然是一个错误。

这是生成的 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
$$$$

在我的代码中，我人为地向 PDB 中的两个硫氧键添加双键，因为 PDB 不会捕获共振。

我怎样才能防止这种情况发生，或者只是删除氢。我试过了

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)

代码冻结了

我将其转换为脚本并收到此错误：

Process finished with exit code 139 (interrupted by signal 11:SIGSEGV)
python jupyter openbabel pybel
1个回答
0
投票

我找到了加氢后去除氢键的方法：

# mol = next(pb.readfile("sdf", babel_sdfh)) # (readfile is described below)
mol = next(pb.readfile("pdb", std_pdb)) # (readfile is described below)
print(mol.molwt)
print(len(mol.atoms))
#%%
# dir(mol.OBMol)
#%%
print(f'num atoms: {mol.OBMol.NumAtoms()}')
print(f'mol.molwt: {mol.molwt}')
#%%
mol.OBMol.AddHydrogens()
#%%
print(f'num atoms: {mol.OBMol.NumAtoms()}')
print(f'mol.molwt: {mol.molwt}')
#%%
sulfur_atom = mol.OBMol.GetAtom(23)
#%%
sulfur_atom.GetType()
#%%
def remove_atom_hydrogens(mol_to_sheer, atom_to_sheer, adjacent_atoms=None, adjacent_atom_types=None):

    for neighbour_atom in ob.OBAtomAtomIter(atom_to_sheer):

        atom_type = neighbour_atom.GetType()
        bond = atom_to_sheer.GetBond(neighbour_atom)
        neighbour_index = bond.GetNbrAtomIdx(atom_to_sheer)

        bl = bond.GetLength()
        # print(bond.GetBondOrder())
        
        # print(bond.IsAromatic())
        # print(bond.IsInRing())
        # print(bond.IsRotor())
        # print(bond.IsAmide())

        print(f'neighbour_index: {neighbour_index} atom_type: {atom_type} bond length: {bl}')
        if atom_type[0] == 'H':

            # pass
            # at.GetParent().DeleteAtom(neighbour_atom)
            mol_to_sheer.OBMol.DeleteBond(bond)
            mol_to_sheer.OBMol.DeleteAtom(neighbour_atom)

        if adjacent_atoms is not None and adjacent_atom_types:
            if atom_type[0] in adjacent_atom_types:
                adjacent_atoms.append(neighbour_atom)


    return adjacent_atoms

adjacent_oxygens = remove_atom_hydrogens(mol, sulfur_atom, [], ['O'])
print(f'mol.OBMol.NumAtoms(): {mol.OBMol.NumAtoms()}')

for oxy in adjacent_oxygens:
    remove_atom_hydrogens(mol, oxy)

mol.OBMol.NumAtoms()
#%%
mol_name
#%%
mol
#%%
mol.OBMol.NumAtoms()
#%%
mol.write("sdf", babel_sdfh1, True)

我很高兴看到是否有一种方法可以更正确地设置输入 PDB 的 ip 或以更高的精度添加氢。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.