Closed Spartanzhao closed 1 year ago
stefdoerr
commented
1 year ago The file you posted here works for me:
from moleculekit.molecule import Molecule
from moleculekit.tools.atomtyper import prepareProteinForAtomtyping
mol = Molecule("/tmp/test.pdb")
prepareProteinForAtomtyping(mol)
Spartanzhao
commented
1 year ago Thanks, I'll take a try.
Spartanzhao
commented
1 year ago I still have problems: here are my codes:
protein=sys.argv[1]
i=protein
name=i.split('/')[-1].split('.pdb')[0]
os.system('/gdata01/user/zhaoyue01/anaconda3/envs/pymol/bin/python /gdata01/user/zhaoyue01/test/learning/drug_response/PDB_descriptor/addH.py {} {}.addH.pdb'.format(i,name))
protein=Molecule('{}.addH.pdb'.format(name))
prepareProteinForAtomtyping(protein)
protein.get('resid', sel='resname CYS')
protein.filter("protein")
protein.guessBonds()
protein.bonds = protein._getBonds()
protein.bonds, protein.bondtype = calculateUniqueBonds(protein.bonds, protein.bondtype)
#print ('############',protein.bondtype)
#newprotein = autoSegment(protein)
#newprotein.bonds, newprotein.bondtype = calculateUniqueBonds(newprotein.bonds,newprtein.bondtype)
#newprotein.bonds = newprotein._getBonds()
#profeat, procenters,N = voxeldescriptors.getVoxelDescriptors(newprotein, buffer=8)
profeat, procenters,N = voxeldescriptors.getVoxelDescriptors(protein, buffer=8)
name=i.split('/')[-1].split('.pdb')[0]+'.npy'
np.save(name,profeat)
df=pd.DataFrame(procenters)
df.to_csv(name+'_procenters.xls',sep='\t',index=False,header=False)
files=open(name+'_dim.txt','w')
for i in files:
files.write(i+';')and bugs are:
Traceback (most recent call last):
File "/gdata01/user/zhaoyue01/test/PDB_Chem/PDB_download/decode/desicriptor2.py", line 31, in <module>
profeat, procenters,N = voxeldescriptors.getVoxelDescriptors(protein, buffer=8)
File "/gdata01/user/zhaoyue01/anaconda3/envs/htmd2/lib/python3.10/site-packages/moleculekit/tools/voxeldescriptors.py", line 250, in getVoxelDescriptors
channels, mol = getChannels(mol, aromaticNitrogen, version, validitychecks)
File "/gdata01/user/zhaoyue01/anaconda3/envs/htmd2/lib/python3.10/site-packages/moleculekit/tools/voxeldescriptors.py", line 120, in getChannels
mol.atomtype, mol.charge = getPDBQTAtomTypesAndCharges(
File "/gdata01/user/zhaoyue01/anaconda3/envs/htmd2/lib/python3.10/site-packages/moleculekit/tools/atomtyper.py", line 275, in getPDBQTAtomTypesAndCharges
atomtypingValidityChecks(mol)
File "/gdata01/user/zhaoyue01/anaconda3/envs/htmd2/lib/python3.10/site-packages/moleculekit/tools/atomtyper.py", line 233, in atomtypingValidityChecks
raise ValueError(
ValueError: The protein has less bonds than (number of atoms - 1). This seems incorrect. You can assign bonds with `mol.bonds = mol._getBonds()`But I had already used the command mol.bonds = mol._getBonds(). How do you successfully getVoxelDescriptors ?
Spartanzhao
commented
1 year ago If the command line is:
from htmd.ui import *
from moleculekit.tools import voxeldescriptors
from moleculekit.util import boundingBox
import numpy as np
from moleculekit.molecule import Molecule
from moleculekit.tools.atomtyper import prepareProteinForAtomtyping
from moleculekit.molecule import calculateUniqueBonds
import os
import sys
import pandas as pd
protein=sys.argv[1]
i=protein
name=i.split('/')[-1].split('.pdb')[0]
os.system('/gdata01/user/zhaoyue01/anaconda3/envs/pymol/bin/python /gdata01/user/zhaoyue01/test/learning/drug_response/PDB_descriptor/addH.py {} {}.addH.pdb'.format(i,name))
protein=Molecule('{}.addH.pdb'.format(name))
prepareProteinForAtomtyping(protein)
protein.get('resid', sel='resname CYS')
protein.filter("protein")
protein.guessBonds()
protein.bonds = protein._getBonds()
protein.bonds, protein.bondtype = calculateUniqueBonds(protein.bonds, protein.bondtype)
#print ('############',protein.bondtype)
newprotein = autoSegment(protein)
newprotein.bonds, newprotein.bondtype = calculateUniqueBonds(newprotein.bonds,newprotein.bondtype)
newprotein.bonds = newprotein._getBonds()
profeat, procenters,N = voxeldescriptors.getVoxelDescriptors(newprotein, buffer=8)
#profeat, procenters,N = voxeldescriptors.getVoxelDescriptors(protein, buffer=8)
name=i.split('/')[-1].split('.pdb')[0]+'.npy'
np.save(name,profeat)
df=pd.DataFrame(procenters)
df.to_csv(name+'_procenters.xls',sep='\t',index=False,header=False)
files=open(name+'_dim.txt','w')
for i in files:
files.write(i+';')Which is ._getBonds after calculateUniqueBonds, Then I got an error:
Traceback (most recent call last):
File "/gdata01/user/zhaoyue01/test/PDB_Chem/PDB_download/decode/desicriptor2.py", line 28, in <module>
profeat, procenters,N = voxeldescriptors.getVoxelDescriptors(newprotein, buffer=8)
File "/gdata01/user/zhaoyue01/anaconda3/envs/htmd2/lib/python3.10/site-packages/moleculekit/tools/voxeldescriptors.py", line 250, in getVoxelDescriptors
channels, mol = getChannels(mol, aromaticNitrogen, version, validitychecks)
File "/gdata01/user/zhaoyue01/anaconda3/envs/htmd2/lib/python3.10/site-packages/moleculekit/tools/voxeldescriptors.py", line 120, in getChannels
mol.atomtype, mol.charge = getPDBQTAtomTypesAndCharges(
File "/gdata01/user/zhaoyue01/anaconda3/envs/htmd2/lib/python3.10/site-packages/moleculekit/tools/atomtyper.py", line 275, in getPDBQTAtomTypesAndCharges
atomtypingValidityChecks(mol)
File "/gdata01/user/zhaoyue01/anaconda3/envs/htmd2/lib/python3.10/site-packages/moleculekit/tools/atomtyper.py", line 241, in atomtypingValidityChecks
raise RuntimeError(
RuntimeError: The protein has duplicate bond information. This will mess up atom typing. Please keep only unique bonds in the molecule. If you want you can use moleculekit.molecule.calculateUniqueBonds for this.
Sorry but this protein structure is really broken. It's not possible to reconstruct the bonds from such a broken structure.
Spartanzhao
commented
1 year ago Okay, thank you so much
I was trying to use getVoxelDescriptors to generate a matrix from pdb file for machine learning, the pdb file comes from PDBbind database, one of the protein's pocket structure which ligand can get insert into it. but for most of this kind of pdb file are failed to transform into a matrix, and the whole bug report is :
This PDB file is missing too many (83 out of 567, 0.146384) heavy atoms to accurately repair the file.The current repair limit is set at 0.1. You may also see this message if PDB2PQR does not have parameters for enough residues in your biomolecule.
This is my pdb file:
Hope you can help me to solve the issue!