search

forlilab / Meeko

Interfacing RDKit and AutoDock
GNU Lesser General Public License v2.1
182 stars 43 forks source link

AttributeError: 'PDBQTMolecule' object has no attribute 'export_rdkit_mol' #27

Closed BJWiley233 closed 1 year ago

BJWiley233 commented 2 years ago

When I read in the results from Vina I get this error???

pdbqt_mol = PDBQTMolecule.from_file("/Users/brian/test/lab9/KRAS_pocket/Cpd1_scaffold/ZINC000000293947/results2/out.pdbqt", 
                                    skip_typing=True)

for pose in pdbqt_mol:
    print(pose.pose_id)
    rdkit_mol = pose.export_rdkit_mol()  
Traceback (most recent call last):

  Input In [23] in <cell line: 1>
    rdkit_mol = pose.export_rdkit_mol()

AttributeError: 'PDBQTMolecule' object has no attribute 'export_rdkit_mol'
BJWiley233 commented 2 years ago

Additionally does the commandline handle flexible residues in results file from Vina? Because only the rotatable bond atoms are include in the PDBQT results file while the backbone atoms are still in the rigid receptor file? Also since the results PDBQT have the flexible residues in addition to small-molecule, does the template need to include these flex residue atoms in additional the ligand?

I am getting this error:

$ mk_copy_coords.py -i /Users/brian/test/lab9/KRAS_pocket/Cpd1_scaffold/ZINC000000293947/results2/out.best.mol2 -s /Users/brian/test/lab9/KRAS_pocket/Cpd1_scaffold/ZINC000000293947/results2/meeko.out /Users/brian/test/lab9/KRAS_pocket/Cpd1_scaffold/ZINC000000293947/results2/out.pdbqt
Traceback (most recent call last):
  File "/Users/brian/anaconda3/envs/vina/bin/mk_copy_coords.py", line 76, in <module>
    pose.copy_coordinates_to_obmol(copy_obmol)
  File "/Users/brian/tools/Meeko/meeko/molecule_pdbqt.py", line 644, in copy_coordinates_to_obmol
    raise RuntimeError('Heavy atom in OBMol is missing, only hydrogens can be missing')
RuntimeError: Heavy atom in OBMol is missing, only hydrogens can be missing
BJWiley233 commented 2 years ago

I see what the issue is. Here is it looking for an index_map

if ob_index in index_map:
    pdbqt_index = index_map[ob_index] - 1
    x, y, z = self._positions[self._current_pose][pdbqt_index, :]
    atom.SetVector(x, y, z)
    n_matched_atoms += 1
elif atom.GetAtomicNum() != 1:
    raise RuntimeError('Heavy atom in OBMol is missing, only hydrogens can be missing')

but each of my poses has no index_map:

pose._pose_data['index_map']
Out[32]: {}
pose._pose_data['n_poses']
Out[45]: 2

Also see below I have one Hydrogen Donor, atom number 13 (HD in last column) in my PDBQT results but none are showing in pose._atom_annotations but instead showing under pose._pose_data['smiles_h_parent']

pose._atom_annotations['hb_don']
Out[47]: []
pose._pose_data['smiles_h_parent']
Out[50]: [13, 13]

This is example of my Vina output. I tested removing most of the poses except 2. Why don't you ask Dr. Sorli if I can collaborate on fixing these issues since I will eventually be using the covalent docking for cryptic cysteine pockets and I would like to use Meeko to parse my results and I have more faith in Dr. Sorli lab's new code than Dr. Sanner's ADRF/ADT older code for Vina.

MODEL 1
REMARK VINA RESULT:    -9.386      0.000      0.000
REMARK INTER + INTRA:         -10.881
REMARK INTER:                  -6.293
REMARK INTRA:                  -4.588
REMARK UNBOUND:                -0.398
REMARK SMILES c1ccc2c(Cn3ccnc3)c[nH]c2c1
REMARK SMILES IDX 6 1 7 2 9 3 8 4 10 5 11 6 5 7 12 8 4 9 13 10 14 11 3 12
REMARK SMILES IDX 15 14 2 15 1 16
REMARK H PARENT 13 13
REMARK Flexibility Score: inf
ROOT
ATOM      1  C   UNL     1      25.069   2.956 -10.736  1.00  0.00     0.160 C 
ENDROOT
BRANCH   1   2
ATOM      2  N   UNL     1      24.536   2.378  -9.523  1.00  0.00    -0.333 N 
ATOM      3  C   UNL     1      24.303   1.880  -7.422  1.00  0.00     0.130 A 
ATOM      4  C   UNL     1      25.244   2.254  -8.374  1.00  0.00     0.107 A 
ATOM      5  N   UNL     1      23.187   1.495  -8.086  1.00  0.00    -0.245 NA
ATOM      6  C   UNL     1      23.560   1.447  -9.389  1.00  0.00     0.199 A 
ENDBRANCH   1   2
BRANCH   1   7
ATOM      7  C   UNL     1      25.343   1.860 -11.721  1.00  0.00    -0.003 A 
ATOM      8  C   UNL     1      26.167   0.757 -11.534  1.00  0.00     0.089 A 
ATOM      9  C   UNL     1      25.321   2.009 -13.102  1.00  0.00     0.004 A 
ATOM     10  N   UNL     1      26.580   0.285 -12.735  1.00  0.00    -0.361 N 
ATOM     11  C   UNL     1      26.322   1.218 -13.654  1.00  0.00     0.046 A 
ATOM     12  C   UNL     1      24.510   2.777 -13.924  1.00  0.00     0.010 A 
ATOM     13  H   UNL     1      26.996  -0.592 -12.901  1.00  0.00     0.166 HD
ATOM     14  C   UNL     1      26.754   1.492 -14.948  1.00  0.00     0.026 A 
ATOM     15  C   UNL     1      24.741   2.750 -15.293  1.00  0.00     0.001 A 
ATOM     16  C   UNL     1      26.044   2.445 -15.669  1.00  0.00     0.002 A 
ENDBRANCH   1   7
TORSDOF 2
BEGIN_RES LYS A   5
REMARK  5 active torsions:
REMARK  status: ('A' for Active; 'I' for Inactive)
REMARK    1  A    between atoms: CA   and  CB  
REMARK    2  A    between atoms: CB   and  CG  
REMARK    3  A    between atoms: CG   and  CD  
REMARK    4  A    between atoms: CD   and  CE  
REMARK    5  A    between atoms: CE   and  NZ  
ROOT
ATOM      1  CA  LYS A   5      20.838  -3.457  -5.071  1.00  0.00     0.176 C 
ENDROOT
BRANCH   1   2
ATOM      2  CB  LYS A   5      22.139  -3.861  -5.749  1.00  0.00     0.035 C 
BRANCH   2   3
ATOM      3  CG  LYS A   5      22.949  -2.644  -6.171  1.00  0.00     0.004 C 
BRANCH   3   4
ATOM      4  CD  LYS A   5      22.981  -2.559  -7.701  1.00  0.00     0.027 C 
BRANCH   4   5
ATOM      5  CE  LYS A   5      21.719  -1.965  -8.290  1.00  0.00     0.229 C 
BRANCH   5   6
ATOM      6  NZ  LYS A   5      21.294  -0.733  -7.565  1.00  0.00    -0.079 N 
ATOM      7  HZ1 LYS A   5      21.017   0.012  -8.468  1.00  0.00     0.274 HD
ATOM      8  HZ2 LYS A   5      20.187  -0.915  -7.130  1.00  0.00     0.274 HD
ATOM      9  HZ3 LYS A   5      21.832  -0.143  -6.681  1.00  0.00     0.274 HD
ENDBRANCH   5   6
ENDBRANCH   4   5
ENDBRANCH   3   4
ENDBRANCH   2   3
ENDBRANCH   1   2
END_RES LYS A   5
BEGIN_RES GLU A  37
REMARK  3 active torsions:
REMARK  status: ('A' for Active; 'I' for Inactive)
REMARK    6  A    between atoms: CA   and  CB  
REMARK    7  A    between atoms: CB   and  CG  
REMARK    8  A    between atoms: CG   and  CD  
ROOT
ATOM     10  CA  GLU A  37      19.784   8.059 -14.581  1.00  0.00     0.177 C 
ENDROOT
BRANCH  10  11
ATOM     11  CB  GLU A  37      20.414   7.083 -15.556  1.00  0.00     0.045 C 
BRANCH  11  12
ATOM     12  CG  GLU A  37      20.159   5.633 -15.194  1.00  0.00     0.116 C 
BRANCH  12  13
ATOM     13  CD  GLU A  37      21.418   4.830 -15.394  1.00  0.00     0.172 C 
ATOM     14  OE1 GLU A  37      21.355   3.583 -15.321  1.00  0.00    -0.648 OA
ATOM     15  OE2 GLU A  37      22.473   5.439 -15.668  1.00  0.00    -0.648 OA
ENDBRANCH  12  13
ENDBRANCH  11  12
ENDBRANCH  10  11
END_RES GLU A  37
BEGIN_RES ASP A  54
REMARK  2 active torsions:
REMARK  status: ('A' for Active; 'I' for Inactive)
REMARK    9  A    between atoms: CA   and  CB  
REMARK   10  A    between atoms: CB   and  CG  
ROOT
ATOM     16  CA  ASP A  54      20.476   1.001  -3.137  1.00  0.00     0.180 C 
ENDROOT
BRANCH  16  17
ATOM     17  CB  ASP A  54      21.854   0.374  -3.098  1.00  0.00     0.147 C 
BRANCH  17  18
ATOM     18  CG  ASP A  54      22.912   1.205  -3.733  1.00  0.00     0.175 C 
ATOM     19  OD2 ASP A  54      22.744   2.452  -3.754  1.00  0.00    -0.648 OA
ATOM     20  OD1 ASP A  54      23.952   0.641  -4.105  1.00  0.00    -0.648 OA
ENDBRANCH  17  18
ENDBRANCH  16  17
END_RES ASP A  54
ENDMDL
MODEL 2
REMARK VINA RESULT:    -9.087      0.545      1.045
REMARK INTER + INTRA:         -10.547
REMARK INTER:                  -5.972
REMARK INTRA:                  -4.575
REMARK UNBOUND:                -0.398
REMARK SMILES c1ccc2c(Cn3ccnc3)c[nH]c2c1
REMARK SMILES IDX 6 1 7 2 9 3 8 4 10 5 11 6 5 7 12 8 4 9 13 10 14 11 3 12
REMARK SMILES IDX 15 14 2 15 1 16
REMARK H PARENT 13 13
REMARK Flexibility Score: inf
ROOT
ATOM      1  C   UNL     1      25.773   2.488 -10.317  1.00  0.00     0.160 C 
ENDROOT
BRANCH   1   2
ATOM      2  N   UNL     1      24.834   2.083  -9.296  1.00  0.00    -0.333 N 
ATOM      3  C   UNL     1      23.982   1.687  -7.337  1.00  0.00     0.130 A 
ATOM      4  C   UNL     1      25.188   1.791  -8.020  1.00  0.00     0.107 A 
ATOM      5  N   UNL     1      22.998   1.599  -8.265  1.00  0.00    -0.245 NA
ATOM      6  C   UNL     1      23.649   1.441  -9.444  1.00  0.00     0.199 A 
ENDBRANCH   1   2
BRANCH   1   7
ATOM      7  C   UNL     1      25.810   1.442 -11.390  1.00  0.00    -0.003 A 
ATOM      8  C   UNL     1      26.523   0.249 -11.367  1.00  0.00     0.089 A 
ATOM      9  C   UNL     1      25.650   1.689 -12.748  1.00  0.00     0.004 A 
ATOM     10  N   UNL     1      26.742  -0.177 -12.634  1.00  0.00    -0.361 N 
ATOM     11  C   UNL     1      26.489   0.840 -13.460  1.00  0.00     0.046 A 
ATOM     12  C   UNL     1      24.843   2.594 -13.420  1.00  0.00     0.010 A 
ATOM     13  H   UNL     1      27.037  -1.079 -12.899  1.00  0.00     0.166 HD
ATOM     14  C   UNL     1      26.799   1.158 -14.778  1.00  0.00     0.026 A 
ATOM     15  C   UNL     1      24.913   2.639 -14.807  1.00  0.00     0.001 A 
ATOM     16  C   UNL     1      26.124   2.228 -15.353  1.00  0.00     0.002 A 
ENDBRANCH   1   7
TORSDOF 2
BEGIN_RES LYS A   5
REMARK  5 active torsions:
REMARK  status: ('A' for Active; 'I' for Inactive)
REMARK    1  A    between atoms: CA   and  CB  
REMARK    2  A    between atoms: CB   and  CG  
REMARK    3  A    between atoms: CG   and  CD  
REMARK    4  A    between atoms: CD   and  CE  
REMARK    5  A    between atoms: CE   and  NZ  
ROOT
ATOM      1  CA  LYS A   5      20.838  -3.457  -5.071  1.00  0.00     0.176 C 
ENDROOT
BRANCH   1   2
ATOM      2  CB  LYS A   5      22.139  -3.861  -5.749  1.00  0.00     0.035 C 
BRANCH   2   3
ATOM      3  CG  LYS A   5      23.200  -2.780  -5.609  1.00  0.00     0.004 C 
BRANCH   3   4
ATOM      4  CD  LYS A   5      22.734  -1.512  -6.333  1.00  0.00     0.027 C 
BRANCH   4   5
ATOM      5  CE  LYS A   5      22.320  -1.764  -7.768  1.00  0.00     0.229 C 
BRANCH   5   6
ATOM      6  NZ  LYS A   5      21.121  -0.964  -8.150  1.00  0.00    -0.079 N 
ATOM      7  HZ1 LYS A   5      21.290  -0.798  -9.329  1.00  0.00     0.274 HD
ATOM      8  HZ2 LYS A   5      20.175  -1.706  -8.202  1.00  0.00     0.274 HD
ATOM      9  HZ3 LYS A   5      20.663   0.013  -7.646  1.00  0.00     0.274 HD
ENDBRANCH   5   6
ENDBRANCH   4   5
ENDBRANCH   3   4
ENDBRANCH   2   3
ENDBRANCH   1   2
END_RES LYS A   5
BEGIN_RES GLU A  37
REMARK  3 active torsions:
REMARK  status: ('A' for Active; 'I' for Inactive)
REMARK    6  A    between atoms: CA   and  CB  
REMARK    7  A    between atoms: CB   and  CG  
REMARK    8  A    between atoms: CG   and  CD  
ROOT
ATOM     10  CA  GLU A  37      19.784   8.059 -14.581  1.00  0.00     0.177 C 
ENDROOT
BRANCH  10  11
ATOM     11  CB  GLU A  37      20.414   7.083 -15.556  1.00  0.00     0.045 C 
BRANCH  11  12
ATOM     12  CG  GLU A  37      20.341   5.643 -15.086  1.00  0.00     0.116 C 
BRANCH  12  13
ATOM     13  CD  GLU A  37      21.644   4.946 -15.379  1.00  0.00     0.172 C 
ATOM     14  OE1 GLU A  37      22.658   5.638 -15.617  1.00  0.00    -0.648 OA
ATOM     15  OE2 GLU A  37      21.673   3.698 -15.332  1.00  0.00    -0.648 OA
ENDBRANCH  12  13
ENDBRANCH  11  12
ENDBRANCH  10  11
END_RES GLU A  37
BEGIN_RES ASP A  54
REMARK  2 active torsions:
REMARK  status: ('A' for Active; 'I' for Inactive)
REMARK    9  A    between atoms: CA   and  CB  
REMARK   10  A    between atoms: CB   and  CG  
ROOT
ATOM     16  CA  ASP A  54      20.476   1.001  -3.137  1.00  0.00     0.180 C 
ENDROOT
BRANCH  16  17
ATOM     17  CB  ASP A  54      21.854   0.374  -3.098  1.00  0.00     0.147 C 
BRANCH  17  18
ATOM     18  CG  ASP A  54      22.911   1.203  -3.737  1.00  0.00     0.175 C 
ATOM     19  OD2 ASP A  54      24.022   0.659  -3.970  1.00  0.00    -0.648 OA
ATOM     20  OD1 ASP A  54      22.610   2.346  -4.110  1.00  0.00    -0.648 OA
ENDBRANCH  17  18
ENDBRANCH  16  17
END_RES ASP A  54
ENDMDL
BJWiley233 commented 2 years ago

Nevermind, I am dumb as I forgot to add the --add_index_map to mk_prepare_ligand.py. We might want to add this as a note in the README.md but otherwise great job!

Last question, is it possible to make it optional to exclude the flexible residue atoms from the final sdf file of the results from mk_copy_coords.py. I think this is beneficial for let's say you want to reconstruct the flexible residues, rigid receptor, and best ligand pose together as a final PDB file. Else I would need to find a way to remove these atoms from the final sdf file because it's actually easier to convert the flexible atoms from Vina PDBQT to PDB and just write a bash script to recreate the receptor and then convert final sdf to PDB and add the HETATMs at the end. Thoughts?

diogomart commented 2 years ago

Since meeko v0.3 we have focused on using the smiles idx instead of the index_map because it does not require a template molecule to be passed with the "-i" option in mk_copy_coords. Handling the flexible residues with smiles idx is still being worked on: https://github.com/forlilab/Meeko/pull/21

We haven't tested the index_map for flexible residues, please keep an eye for suspicious behavior if this is what you are doing.

I think it is a good idea to have an option to export an rdkit molecule (or SDF) for the ligand and the flexible residues separately.

ale94mleon commented 1 year ago

Hi, I recently update to meeko-0.4.0. and export_rdkit_mol method is not available anymore for the class PDBQTMolecule. Why was it removed? Were there some errors? I think, as @diogomart mentioned before, it is a good idea to implement this method. For now I will just stick to 0.3.3 version.

diogomart commented 1 year ago

@ale94mleon the API changed, see the examples in the README

ale94mleon commented 1 year ago

Ok, got it. Thanks! Changing back to meeko-0.4.0 :-)

diogomart commented 1 year ago

Just reporting that #21 has been merged, flexres are now exported with mk_export.py (which was called mk_copy_coords.py before v0.4).