1
$\begingroup$

I have some PDB IDs and for every structure, I need all its ligands, so I want to automate the process.

The ligands are different for every chain:

$ cat ./7rdx_N_1N7.smi 
C1[C@@]2(C)[C@@H](C[C@@H](C1)O)[C][C]([C]1[C@H]3[C@]([C@H](C[C]21)O)([C@H](CC3)[C]([C])[C][C])C)O   1N7
$ cat ./7rdx_U_1N7.smi 
C1[C@@]2(C)[C@@H](C[C@@H](C1)O)C[C@H]([C@H]1[C@H]3[C@]([C@H](C[C@H]21)O)([C@H](CC3)[C@H](C)CCC(=O)NC[C][C][N+]([C])([C])[C])C)O 1N7

So, I cannot just download the ligands by its name (e.g. https://files.rcsb.org/ligands/download/ZN_ideal.sdf).

I've tried to find the structures in .pdb file for a given structure, but I see only their names:
pdb_file_part

Also, as you could see above, I'll convert the .sdf / .mol2 files to SMARTS which seems to be a simpler representation - maybe this information'll make the task easier.

I want to make a script for the process in Python, and the only way to access the RSCB PDB with a library I've found was using biopython but I found there an option to download just a .pdb file (Bio.PDB.PDBList.retrieve_pdb_file), not ligands.

I've also tried to make my own script, but still I don't know how to access ligands files automatically for every structure. The .pdb files I can access simply by ftp (e.g. https://ftp.wwpdb.org/pub/pdb/data/biounit/PDB/divided/rd/).

Improve this question
$\endgroup$
12

4 Answers 4

Reset to default
1
$\begingroup$

Ok, cannot try the API approach, no resources during vacation. Got the download and parse/split approach, added conversion to SMILES by Openbabel. I believe it is exactly the same question of : https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb from where I stole the juicy bit of the code:

## stealing from https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb

from Bio import __version__


print('BIOPYTHON Version : ' , __version__)

from Bio.PDB import MMCIFParser, PDBIO, Select, PDBList

import subprocess

filename = '7rdx.cif'


pdblist = PDBList()

pdblist.retrieve_pdb_file(filename.upper().split('.')[0] , pdir = '.')


pars = MMCIFParser(QUIET = True)

struct = pars.get_structure('7rdx',  filename)


print(struct, type(struct))


def is_het(residue):
    res = residue.id[0]
    
    if res != " " and res != "W":
        
        return True #1

    else:
    
        return False #0
 

class ResidueSelect(Select):
    def __init__(self, chain, residue):
        self.chain = chain
        self.residue = residue

    def accept_chain(self, chain):
    
        if chain.id == self.chain.id:
            
            return True #1
        
        else:
            
            return False #0
        
        
    def accept_residue(self, residue):
        """ Recognition of heteroatoms - Remove water molecules """
        
        if residue == self.residue and is_het(residue):
            
            return True #1
        
        else:
            
            return False #0


for model in struct :
    
    
    dizio = model.child_dict
    
    print('\n'+ str(dizio))
    
    

for model in struct:
    
    i = 1
    
    print(model, model.id)
    
    for chain in model:
        
        print(chain, chain.id)
        
        for res in chain:
            
            # print('\n'+'res : ', res, 'res.resname : ', res.resname)
            
            if res.id[0] != " " :
                
                
                print('\n'+'res : ', res, res.id, res.resname)
                
                io = PDBIO()
                
                io.set_structure(struct)
                
                io.save(f"lig_{chain.id}_{res.resname}_{i}.pdb", ResidueSelect(chain, res))
                
                try:
                    subprocess.run(["openbabel.obabel", "-ipdb", f"lig_{chain.id}_{res.resname}_{i}.pdb", "-osmi" ,"-O",
                                    
                                    f"lig_{chain.id}_{res.resname}_{i}.smi"])
                    

                    
                except Exception as excp:
                    
                    print('\n ', excp, '\n')
                    
                    
                
                i += 1

Problem is Biopython parses mmCIFs files; that are now the official standard for RCSB PDB Submission of PDBx/mmCIF format files for crystallographic depositions to the PDB will be mandatory from July 1st 2019 onward. PDB format files will no longer be accepted for deposition of structures solved by MX techniques. A nice review about the long history of the 'new' format :[ PDBx/mmCIF Ecosystem: Foundational Semantic Tools for Structural Biology.

Unfortunately the Bio.PDB.MMCIFParser module return a a structure that is:

<Structure id=7rdx> <class 'Bio.PDB.Structure.Structure'> so all the ligands are related to PDB chains classification and not PDBx entity.

I could be wrong here, I am a novice in both Biopython and Structural Bioinformatics but this is what I've noticed, please feel free to elaborate on this and add a solution using the RCSB APIes ... before curiosity kills the cat

PS:

NEEDs to be tested with sample containing water

Improve this answer
$\endgroup$
1
$\begingroup$

adding a different way to solve the question , still without using RCSB PDB database, but using Biopython.PDB.MMCIF2Dict module to start to get a feeling on how info are stored in .cif format files. Ligands are selected amongs the ones in the structure for their Mw.

Ligands are then saved using the same approach of first answer or a new one that takes advantage of Biopython.PDB.Dice module to extract and save portion of PDB. Again Openbabel is used to translate PDBs to SMILES.

Code :

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jan  5 11:43:07 2023

@author: bob

https://bioinformatics.stackexchange.com/questions/20301/how-to-download-ligands-for-pdb-structure?noredirect=1#comment29972_20301




NEEDs to be tested with sample containing water


"""

## stealing from https://stackoverflow.com/questions/61390035/how-to-save-each-ligand-from-a-pdb-file-separately-with-bio-pdb

from Bio import __version__


print('BIOPYTHON Version : ' , __version__)

from Bio.PDB import (
                    MMCIFParser, 
                    PDBIO, 
                    Select, 
                    PDBList, 
                    MMCIF2Dict,
                    Dice
                    )

import subprocess

filename = '7rdx.cif'

# filename = '7RDX.cif'  # this will  work , but biopython saves in lower case 


# pdblist = PDBList()

# pdblist.retrieve_pdb_file(filename.split('.')[0] , pdir = '.')


pars = MMCIFParser(QUIET = True)

struct = pars.get_structure('7rdx',  filename)


# print(struct, type(struct))


mmcifdict = MMCIF2Dict.MMCIF2Dict(filename)


# print(mmcifdict)

# for i in mmcifdict:
    
#     print('\n', i ,'   ---> ', mmcifdict[i])
    
print('##############################################################')
    
# print(mmcifdict['_entity'])

 
print(mmcifdict['_pdbx_nonpoly_scheme.entity_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.asym_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.entity_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.mon_id'])

ligands_map_mmcif = list(zip(mmcifdict['_pdbx_nonpoly_scheme.entity_id'],
                             
                             mmcifdict['_pdbx_nonpoly_scheme.asym_id'],
                             
                             mmcifdict['_pdbx_nonpoly_scheme.entity_id'],
                             
                             mmcifdict['_pdbx_nonpoly_scheme.mon_id']))


print('#########################################')

print(ligands_map_mmcif)

print(mmcifdict['_pdbx_nonpoly_scheme.entity_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.pdb_seq_num'])
print(mmcifdict['_pdbx_nonpoly_scheme.pdb_strand_id'])
print(mmcifdict['_pdbx_nonpoly_scheme.pdb_mon_id'])

ligands_map_pdb = list(zip(mmcifdict['_pdbx_nonpoly_scheme.entity_id'],
                           
                           mmcifdict['_pdbx_nonpoly_scheme.pdb_seq_num'],
                           
                           mmcifdict['_pdbx_nonpoly_scheme.pdb_strand_id'],
                           
                           mmcifdict['_pdbx_nonpoly_scheme.pdb_mon_id']))

print('#########################################')

print(ligands_map_pdb)


print('WWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWw')

print(mmcifdict['_entity.id'])
print(mmcifdict['_entity.formula_weight'])

entity_Mw = dict(zip(mmcifdict['_entity.id'],
                     
                     mmcifdict['_entity.formula_weight']))

print(entity_Mw)

print('-------------------------------------------------------------')


ligands_map_pdb_red = [i for i in ligands_map_pdb if float(entity_Mw[i[0]]) >= 200 ]

print(ligands_map_pdb_red)





#########################################################   first approach 


# class ResidueSelect(Select):
#     def __init__(self, chain, residue):
#         self.chain = chain
#         self.residue = residue

#     def accept_chain(self, chain):
    
#       if chain.id == self.chain.id:
            
#           return True #1
        
#       else:
            
#           return False #0
        
        
#     def accept_residue(self, residue):
#         """ Recognition of heteroatoms - Remove water molecules """
        
#         if residue == self.residue :
            
            
#             return True #1
        
#         else:
            
#             return False #0




# for model in struct:
    
#     print(model)
    
#     for chain in model:
        
#         if chain.id == 'A':
            
#             for res in chain:
                
#                 print(res, res.id, res.get_resname(), type(res))

# i = 1    

# for res in ligands_map_pdb_red :
    
#     print('\n'+'res : ', res[3], res[0], res[1], type(res[1]), res[2])
    
#     print(struct[0][res[2]] , type(struct[0][res[2]]) )
    
#     print(struct[0][res[2]][('H_'+res[3], int(res[1]), ' ')] , type(struct[0][res[2]][('H_'+res[3], int(res[1]), ' ')]) )
    
                
#     io = PDBIO()
                    
#     io.set_structure(struct)
                    
#     io.save(f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb", ResidueSelect(struct[0][res[2]],
                                                                          
#                                                     struct[0][res[2]][('H_'+res[3], int(res[1]), ' ')]))
                    
#     try:
#         subprocess.run(["openbabel.obabel", "-ipdb", f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb", "-osmi" ,"-O",
                                        
#         f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.smi"])
                        
    
                        
#     except Exception as excp:
                        
#         print('\n ', excp, '\n')
                        
                        
                    
#     i += 1
    
########################################################   second approach 


import warnings

from Bio import BiopythonWarning


def accept_residue_mod(self, residue):
        """Verify if a residue sequence is between the start and end sequence."""
        # residue - between start and end
        hetatm_flag, resseq, icode = residue.get_id()
        # if hetatm_flag != " ":
        #     # skip HETATMS
        #     return 0
        if icode != " ":
            warnings.warn(
                f"WARNING: Icode {icode} at position {resseq}", BiopythonWarning
            )
        if self.start <= resseq <= self.end:
            return 1
        return 0

Dice.ChainSelector.accept_residue = accept_residue_mod


i = 1

for res in ligands_map_pdb_red :
    
    print(struct, res[2], int(res[1]), int(res[1]))

    Dice.extract(struct, res[2], int(res[1]), int(res[1]), f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb")
    
    
    try:
        subprocess.run(["openbabel.obabel", "-ipdb", f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.pdb", "-osmi" ,"-O",
                                           
        f"lig_{res[3]}_{res[2]}_{str(i).zfill(5)}.smi"])
                            
        
    except Exception as excp:
                            
        print('\n ', excp, '\n')
                            
                            
                        
    i += 1

```
Improve this answer
$\endgroup$
1
$\begingroup$

OK, seems ligands and PDBs are not really a trending topic here.... (or my posts are really poor).

Nevertheless the question is literally eating my brain so I went into a second direction: "the API" already skimmed here Help me in understanding the PDB file.

So inspecting and copying directly from the PDB entry 7RDX webpage the download value: https://models.rcsb.org/v1/7rdx/ligand?auth_seq_id=1005&label_asym_id=M&encoding=sdf&filename=7rdx_M_1N7.sdfI got this bit of code:

import requests


def main():   
    risposta = requests.get(f"https://models.rcsb.org/v1/7rdx/ligand?auth_seq_id=1005&label_asym_id=M&encoding=sdf",                      
                      allow_redirects=True)   
    open(filename, 'wb').write(risposta.content)
    print(risposta.status_code)

if __name__ == '__main__':  
   filename = '7rdx_M_1N7.sdf'   
   main()

That is equivalent to feeding this:

{
  "atom_site": [
    {
      "auth_seq_id": 1005,
      "label_asym_id": "M"
    }
  ]
}

to the 'https://models.rcsb.org/#/General/ligand-post' test webpage Request body.

now need to figure out how to get the auth_seq_id and the label_asym_id parameters for all the 'non_polymer_entity_ids' of my targed PDB 7RDX from the https://data.rcsb.org/rest/v1/ DATA API.

Went till here:

import requests
import json

def main(entry_id: str):
    risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/entry/{entry_id}/")   
    parsed = json.loads(risposta.content.decode())    
    print(risposta.status_code)    
    print((parsed))    
    print(parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"])    
    print('\n', type(parsed))
    
    parsed_lig = parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"]    
    print('\n\n',parsed_lig, type(parsed_lig))    
    

if __name__ == '__main__':         
   main('7rdx')

output :

 ['7', '8', '9', '10', '11'] <class 'list'>

The list values are the ids of the 'non_polymer_entity_ids' of my targed PDB 7RDX.

Right now I am trying to recover the the auth_seq_id and the label_asym_id for every instance of the 'non_polymer_entity_ids' inside my mmCIF/PDBx using the DATA API.

Not sure this is the fastest and right way using the RCSB PDB API to accomplish the OP question (download all ligands of the model structure (not the ligands models)) from the PDB database. Please feel free to add more or point me towards the best solution.

EDITED:

It took a while, I am sure this is not the right way to make good use of all the databases, there should be a way to navigate the API schemas better and faster, and to get and manage the intermediate values needed to get the two keys for the download from ModelServer "Post /v1/{id}/ligand Coordinates of the first group satisfying the given criteria", anyway:

import requests
import json

parsed_lig_dict = {}

def main(entry_id: str):    
    risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/entry/{entry_id}/")   
    parsed = json.loads(risposta.content.decode())    
    # print(risposta.status_code)    
    # print((parsed))    
    # print(parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"])
    # print('\n', type(parsed))
    
    parsed_lig = parsed["rcsb_entry_container_identifiers"]["non_polymer_entity_ids"]    
    # print('\n\n',parsed_lig, type(parsed_lig))
    
    # for i in parsed:        
    #     print('\n', i  , parsed[i])
                  
    for lig in parsed_lig:    
        risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity/{entry_id}/{lig}")    
        parsed = json.loads(risposta.content.decode())    
        # print('\n',risposta.status_code, risposta.content)        
        
        # print(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"], 
              # parsed["pdbx_entity_nonpoly"]["comp_id"])
                
        parsed_lig_dict[lig] = [parsed["pdbx_entity_nonpoly"]["comp_id"] , 
                                dict.fromkeys(parsed["rcsb_nonpolymer_entity_container_identifiers"]["asym_ids"])]
    
    for lig in parsed_lig_dict:        
        for chain in parsed_lig_dict[lig][1]:    
            risposta = requests.get(f"https://data.rcsb.org/rest/v1/core/nonpolymer_entity_instance/{entry_id}/{chain}")        
            # print(risposta.status_code)
        
            parsed = json.loads(risposta.content.decode())        
            # print('\n',risposta.status_code, risposta.content)            
            # print(chain , '--> ',parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"])
            
            parsed_lig_dict[lig][1][chain] = parsed["rcsb_nonpolymer_entity_instance_container_identifiers"]["auth_seq_id"]        
               
    
    # print(parsed_lig_dict)
    
    # for i in parsed_lig_dict:
    #     print(parsed_lig_dict[i])
    
    # for i in parsed_lig_dict:
    #     for y in parsed_lig_dict[i][1]:
    #         print(i,parsed_lig_dict[i][0], y , parsed_lig_dict[i][1][y])
                        
    cnt = 1     
    for i in parsed_lig_dict:        
        for y in parsed_lig_dict[i][1]:
            seq_id = parsed_lig_dict[i][1][y]            
            chain = y            
            comp_id = parsed_lig_dict[i][0]            
            print(seq_id , y)

            risposta = requests.get(f"https://models.rcsb.org/v1/7rdx/ligand?auth_seq_id={seq_id}&label_asym_id={chain}&encoding=sdf" ,                              
                              allow_redirects=True)             
            print(risposta.status_code)
                        
            filename = f"lig_{comp_id}_{chain}_{seq_id}_{cnt}.sdf"            
            open(filename, 'wb').write(risposta.content)            
            cnt += 1

    
if __name__ == '__main__':         
   main('7rdx')

This approach is very slow on my machine, the other two answers that parse the downloaded full mmCIF file via Biopython Structure Object are faster even adding the translation .sdf files to SMILES via Openbabel. I think this is because I translate the DATA API JSON formatted schema into a dictionary to extract the data I need, there should be another way to make use of that data. This way doesn't need installed Biopython though.

Improve this answer
$\endgroup$
1
$\begingroup$

Here can use GraphiQL API provided by RCSB. Here is an example

import requests
from rdkit import Chem
pdbid = '7RDX'
ls_smi = []

query = f'''{{
  entry(entry_id: "{pdbid}") {{
    nonpolymer_entities {{
      nonpolymer_comp {{
        rcsb_id
        pdbx_chem_comp_descriptor {{
          type
          comp_id
          program
          descriptor
        }}
      }}
    }}
  }}
}}'''
url = "https://data.rcsb.org/graphql"
r = requests.post(url, json={'query': query})
for i in r.json()['data']['entry']['nonpolymer_entities']:
    if i['nonpolymer_comp']['pdbx_chem_comp_descriptor'][0]['type'] == 'SMILES':
        ls_smi.append(Chem.CanonSmiles(i['nonpolymer_comp']['pdbx_chem_comp_descriptor'][0]['descriptor']))
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks @AustinApple and welcome to the site. A good answer $\endgroup$
    – M__ ♦
    Mar 9 at 22:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.