口袋双共价残基体系筛选

#双共价，要求1.有CYS；2.CYS周围存在有含-OH,-SH，-NH2，-COOH的氨基酸
#用于检测可形成双共价键的cys, ser, thr, asp, glu, lys, gln, asn，同时要求这些氨基酸处于溶剂区，能与小分子接触
#USAGE: python3 covalent_site_detection.py
#输入文件：蛋白质pdb文件和配体sdf文件


from pymol import cmd
import pandas as pd
import numpy as np


def calculate_surface_area(sele):
    cmd.set('dot_solvent', 1)
    cmd.set('dot_density', 4)
    area = cmd.get_area(sele)
    return area


def get_chain_resi_resn_name(sele):
    space = {"res_info":[]}
    cmd.iterate(sele,"res_info.append([chain,resi,resn,name])", space=space)
    res_info = space['res_info']
    #remove duplicates
    res_info = [list(i) for i in set(tuple(_) for _ in res_info)]
    return res_info


df = pd.read_csv("/home/zjxu/ylshi/RAFT/Docking/test/D3CARP_protein.csv", header=0)


for protein in df['PDB code']:
    print(protein)
    cmd.delete("all")
    find_suitable_atm = 0
    cmd.load("/home/zjxu/ylshi/RAFT/Docking/Ligand/sdf/" + protein + "_ligand.sdf",'ligand')
    cmd.load("/home/zjxu/ylshi/RAFT/Docking/Protein/pdb/" + protein + "_protein.pdb",'protein')
    cmd.select("poc_surf_cys", "resn cys within 5 of " + 'ligand')
    poc_surf_cys = [i[:2] for i in get_chain_resi_resn_name("poc_surf_cys")]
    #去重
    poc_surf_cys = [list(i) for i in set(tuple(_) for _ in poc_surf_cys)]
    print(poc_surf_cys)
    if len(poc_surf_cys) != 0:
        for cys in poc_surf_cys:
            cys_s = "(chain "+cys[0]+" and resi "+str(cys[1])+" and resn CYS and name SG)"
            cys_s_sa = calculate_surface_area(cys_s)
            #判断半胱氨酸上的S是否暴露在外
            if cys_s_sa > 5:
                #选择半胱氨酸上的S周围5埃符合条件的原子
                cmd.select("cys_neigh_atm","("+cys_s+" around 5) and ((resn CYS and name SG) or (resn SER and name OG) or (resn THR and name OG1) or (resn ASP and (name OD1 or name OD2)) or (resn GLU and (name OE1 or name OE2)) or (resn LYS and name NZ) or (resn GLN and name NE2) or (resn ASN and name ND2) or (resn ARG and name NH2) or (resn TYR and name OH))")
                cys_neigh_atm = get_chain_resi_resn_name("cys_neigh_atm")
                print(cys_neigh_atm)
                if len(cys_neigh_atm) > 0:
                    for atm in cys_neigh_atm:
                        atm_ = "chain "+atm[0]+" and resi "+str(atm[1])+" and resn "+atm[2]+" and name "+atm[3]
                        atm_sa = calculate_surface_area(atm_)
                        #判断周围符合距离条件的原子是否暴露在外
                        if atm_sa > 5:
                            cmd.distance(str(find_suitable_atm),cys_s, atm_)
                            cmd.show("sticks", "byres "+atm_)
                            cmd.show("sticks","byres "+cys_s)
                            find_suitable_atm += 1
    if find_suitable_atm > 0:
        cmd.save("/home/zjxu/ylshi/RAFT/Docking/test/result/" + protein + ".pse")