from pymol import cmd from pymol.cmd import util import sys import os import math def pi_pi_angle(x1,x2,x3,y1,y2,y3): import numpy as np #print(x1,x2,x3,y1,y2,y3) B1, B2, B3 = [x1[0] - x2[0], x1[1] - x2[1], x1[2] - x2[2]] C1, C2, C3 = [x1[0] - x3[0], x1[1] - x3[1], x1[2] - x3[2]] #print([x1[0] - x2[0], x1[1] - x2[1], x1[2] - x2[2]]) n1 = [B2 * C3 - C2 * B3, B3 * C1 - C3 * B1, B1 * C2 - C1 * B2] D1, D2, D3 = [y1[0] - y2[0], y1[1] - y2[1], y1[2] - y2[2]] E1, E2, E3 = [y1[0] - y3[0], y1[1] - y3[1], y1[2] - y3[2]] n2 = [D2 * E3 - E2 * D3, D3 * E1 - E3 * D1, D1 * E2 - E1 * D2] dot_product = np.dot(n1, n2) magnitude1 = np.linalg.norm(n1) # ???????????????? magnitude2 = np.linalg.norm(n2) #print(dot_product,magnitude1,magnitude2,dot_product / (magnitude1 * magnitude2)) xx = math.acos(dot_product / (magnitude1 * magnitude2)) degree = math.degrees(xx) if degree > 90: degree = 180 - degree return round(degree,2) def pi_pi(aro_ring_info,cutoff=5.0): aro_resn =['HIS','PHE','TRP','TYR'] aro_resn_list = '+'.join(aro_resn) chA_aro_res = f'inter_1 and resn {aro_resn_list}' chB_aro_res = f'inter_2 and resn {aro_resn_list}' chA_aro_res_list = list(set([f'resi {a.resi} and resn {a.resn}' for a in cmd.get_model(chA_aro_res).atom])) chB_aro_res_list = list(set([f'resi {a.resi} and resn {a.resn}' for a in cmd.get_model(chB_aro_res).atom])) for A in chA_aro_res_list: for B in chB_aro_res_list: A_type = A.split()[-1] B_type = B.split()[-1] #print(f'{A} and name {aro_ring_info[A_type]}') x = cmd.centerofmass(f'inter_1 and {A} and name {aro_ring_info[A_type]}') y = cmd.centerofmass(f'inter_2 and {B} and name {aro_ring_info[B_type]}') dist = math.sqrt((x[0] - y[0]) ** 2 + (x[1] - y[1]) ** 2 + (x[2] - y[2]) ** 2) if dist < cutoff: cmd.show('sticks',f'inter_1 and {A}') cmd.show('sticks', f'inter_2 and {B}') cmd.label(f'name CA and inter_1 and {A}', 'oneletter+resi') cmd.label(f'name CA and inter_2 and {B}', 'oneletter+resi') x1 = cmd.get_coords(f"{A} and name {aro_ring_info[A_type].split('+')[0]}")[0] x2 = cmd.get_coords(f"{A} and name {aro_ring_info[A_type].split('+')[1]}")[0] x3 = cmd.get_coords(f"{A} and name {aro_ring_info[A_type].split('+')[2]}")[0] y1 = cmd.get_coords(f"{B} and name {aro_ring_info[B_type].split('+')[0]}")[0] y2 = cmd.get_coords(f"{B} and name {aro_ring_info[B_type].split('+')[1]}")[0] y3 = cmd.get_coords(f"{B} and name {aro_ring_info[B_type].split('+')[2]}")[0] print(A,B) pi_angle = pi_pi_angle(x1, x2, x3, y1, y2, y3) pi_name = f'pi_{A.split()[1]}_{B.split()[1]}_{str(pi_angle)}' print(A, B,pi_name) cmd.distance(pi_name, f"inter_1 and resi {A.split()[1]} and name {aro_ring_info[A_type]}", f'inter_2 and resi {B.split()[1]} and name {aro_ring_info[B_type]}', mode=4) def Salt_bridge(ChA,ChB,cutoff=5.0): selection_1 = f'{ChB} and ((resn LYS and name NZ) or (resn ARG and name NE+NH*))' a_charge_res = set([f'chain {a.chain} and resi {a.resi} and resn {a.resn}' for a in cmd.get_model(f'inter_1 and resn LYS+ARG+ASP+GLU').atom]) b_charge_res = set([f'chain {a.chain} and resi {a.resi} and resn {a.resn}' for a in cmd.get_model(f'inter_2 and resn LYS+ARG+ASP+GLU').atom]) #print(a_charge_res) p_info = {'LYS':'NZ','ARG':'NE+NH*'} n_info = {'ASP':'OD*+OE*','GLU':'OD*+OE*'} for a in a_charge_res: a_resn = a.split()[-1] for b in b_charge_res: b_resn = b.split()[-1] name = f'SB_{a.split()[4]}_{b.split()[4]}' if (a_resn in ['LYS','ARG']) and (b_resn in ['ASP','GLU']): a_resn_atom = p_info[a_resn] b_resn_atom = n_info[b_resn] #print(a,b) #print(a_resn_atom,b_resn_atom) selection_1 = f'{a} and name {a_resn_atom}' selection_2 = f'{b} and name {b_resn_atom}' cal_saltbrige(selection_1, selection_2, name, a,b,cutoff) elif (a_resn in ['ASP','GLU'] ) and (b_resn in ['LYS','ARG']): a_resn_atom = n_info[a_resn] b_resn_atom = p_info[b_resn] selection_1 = f'{a} and name {a_resn_atom}' selection_2 = f'{b} and name {b_resn_atom}' cal_saltbrige(selection_1,selection_2,name, a,b,cutoff) def cal_saltbrige(selection1,selection2,name,a,b,cutoff): x = cmd.centerofmass(selection1) y = cmd.centerofmass(selection2) dist = math.sqrt((x[0] - y[0]) ** 2 + (x[1] - y[1]) ** 2 + (x[2] - y[2]) ** 2) if dist < cutoff: cmd.show('stick', a) cmd.show('stick', b) print(a,b) cmd.label(f'name CA and {a}', 'oneletter+resi') cmd.label(f'name CA and {b}', 'oneletter+resi') cmd.distance(name, selection1, selection2, mode=4) return name def DA_Hbond(DA_type,select1,select2,cutoff='3.5',angle='150'): hbond_pairs = [] pairs = cmd.find_pairs(select1, select2, mode=0, cutoff=cutoff) for pair in pairs: D_atom = f'index {pair[0][1]}' A_atom = f'index {pair[1][1]}' HD_atom = f'e. h and neighbor({D_atom})' index_1 = [a.index for a in cmd.get_model(A_atom).atom] index_2 = [a.index for a in cmd.get_model(HD_atom).atom] index_3 = [a.index for a in cmd.get_model(D_atom).atom] for h in index_2: h_angle = cmd.angle('angle', A_atom, f'index {h}', D_atom) #print(type, pair, h_angle) if h_angle > float(angle): D_atom_name = [f'{a.resi}-{a.name}' for a in cmd.get_model(f'index {pair[0][1]} ').atom][0] A_atom_name = [f'{a.resi}-{a.name}' for a in cmd.get_model(f'index {pair[1][1]} ').atom][0] #name = [f'HB-{a.resi}{a.resn}{a.chain}-{base_atom_name}-lig-{lig_atom_name}' for a in cmd.get_model(f'index {pair[0][1]} and polymer').atom][0] name = f'{DA_type}_{D_atom_name}_{A_atom_name}' cmd.show('stick', f'byres(index {pair[0][1]})') cmd.show('stick', f'byres(index {pair[1][1]})') cmd.label(f'name CA and byres(index {pair[0][1]})', 'oneletter+resi') cmd.label(f'name CA and byres(index {pair[1][1]})', 'oneletter+resi') cmd.distance(name, f'index {index_1[0]} ', f'index {h} ') break return hbond_pairs def run(pdb,ChA,ChB,dist): cmd.delete('all') cmd.load(pdb) pdb_name = pdb[:-4] #cmd.h_add('all') #cmd.remove('solvent and resn SO4, NA,CL,GOL') cmd.remove('not polymer') #util.cbc('all', first_color=7, quiet=1, legacy=0, _self=cmd) cmd.select('inter_1',f'{pdb_name} and chain {ChA} and byres(chain {ChB}) around {dist}') cmd.select('inter_2', f'{pdb_name} and chain {ChB} and byres(chain {ChA}) around {dist}') util.cba('lightblue',f'chain {ChA}') util.cba('gray90', f'chain {ChB}') else_part = cmd.select('else',f'{pdb_name} and not (chain {ChA}+{ChB})') if else_part > 0: util.cba('lime', 'else') chA_D = f'inter_1 and e. o+n and (neighbor e. h)' chB_A = f'inter_2 and e. o+n ' DA_Hbond('HDA',chA_D,chB_A) chB_D = f'inter_2 and e. o+n and (neighbor e. h)' chA_A = f'inter_1 and e. o+n ' DA_Hbond('HAD',chB_D,chA_A) aro_ring_name = {'HIS':'CG+ND1+CE1+NE2+CD2', 'TRP':'CH2+CZ3+CE3+CZ2+CE2+CD2', 'TYR':'CG+CD1+CE1+CZ+CE2+CD2', 'PHE':'CG+CD1+CE1+CZ+CE2+CD2'} pi_pi(aro_ring_name) Salt_bridge(ChA,ChB) cmd.do('set cartoon_transparency, 0.4') cmd.do('set dash_gap, 0.2') cmd.do('set dash_round_ends, 0') cmd.do('set label_size, 24') cmd.do('set dash_color, blue,pi*') cmd.do('set dash_color, violet,SB*') cmd.do('set dash_radius, 0.1') cmd.do('set stick_radius, 0.15') cmd.delete('angle') cmd.disable('angle') cmd.set('bg_rgb', 'white') cmd.save(f'{pdb_name}_interaction.pse') def main(): pdb_file = sys.argv[1] CHA = sys.argv[2] CHB = sys.argv[3] dist = sys.argv[4] run(pdb_file,CHA,CHB,dist) main()