198 lines
8.2 KiB
Python
198 lines
8.2 KiB
Python
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() |