B2AR_Dimer_Interface/RRCS.py at main · CompGenomeLab/B2AR_Dimer_Interface · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
# -*- coding: utf-8 -*-
import sys
def calc_contact(pdbbase):
    fi = open(pdbbase, 'r')
    all_lines = fi.readlines()
    fi.close()
    atom_lines = [l for l in all_lines if l[0:6] == 'ATOM  ']
    dict_coord = {} # dict to store coordinates. dict_coord[res][atom] = (x,y,z,occ)
    atomnum_2_name = {} # map atom number to atom name, in order to find N, CA, C, O
    contact_score = {} # dict to store final results. contact_score[ires][jres] = contact_score.
    for line in atom_lines:
        # retrive info from each atom line
        atom_num = int(line[6:11].strip())
        atom_name = line[12:16].replace(' ', '_')
        res_name = line[17:20]
        res_num = int(line[22:26].strip())
        chain_id = line[21:22]
        x = float(line[30:38].strip())
        y = float(line[38:46].strip())
        z = float(line[46:54].strip())
        occ = float(line[54:60].strip())
        res = chain_id + ':' + str(res_num) + '_' + res_name
        atomnum_2_name[atom_num] = atom_name
        if res_num <= 0:
            continue
        if res not in dict_coord:
            dict_coord[res] = {}
        dict_coord[res][atom_num] = (x, y, z, occ)
    #list shortening
    num_of_aa = len(dict_coord)
    res_list = dict_coord.keys()
    ires_contact = {}
    for ires in res_list:
        ires_contact[ires] = []
        ires_num = int(ires.split(':')[1].split('_')[0].strip())
        for jres in res_list:
            jres_num = int(jres.split(':')[1].split('_')[0].strip())
            if jres_num <= ires_num:
                continue
            jres_flag = 0
            atom_in_ires = dict_coord[ires].keys()
            atom_in_jres = dict_coord[jres].keys()
            for iatom in atom_in_ires:
                (ix, iy, iz, iocc) = dict_coord[ires][iatom]
                for jatom in atom_in_jres:
                    (jx, jy, jz, jocc) = dict_coord[jres][jatom]
                    dx = abs(ix-jx)
                    dy = abs(iy-jy)
                    dz = abs(iz-jz)
                    if dx < 4.63 and dy < 4.63 and dz < 4.63:
                        jres_flag = 1
                        ires_contact[ires].append(jres)
                        break
                if jres_flag:
                    break
        # loop over the shortened list
        contact_score[ires] = {}
        for kres in ires_contact[ires]:
            atom_in_kres = dict_coord[kres].keys()
            kres_num = int(kres.split(':')[1].split('_')[0].strip())
            contact_score[ires][kres] = 0
            total_score = 0
            if abs(ires_num - kres_num) < 5:
                for iatom in atom_in_ires:
                    iatom_name = atomnum_2_name[iatom]
                    if iatom_name in ['_N__', '_CA_', '_C__', '_O__']:
                        continue
                    (ix, iy, iz, iocc) = dict_coord[ires][iatom]
                    for katom in atom_in_kres:
                        katom_name = atomnum_2_name[katom]
                        if katom_name in ['_N__', '_CA_', '_C__', '_O__']:
                            continue
                        (kx, ky, kz, kocc) = dict_coord[kres][katom]
                        d2 = (ix-kx)**2 + (iy-ky)**2 + (iz-kz)**2
                        if d2 >= 21.4369:  # 4.63*4.63 = 21.4369
                            score = 0
                        elif d2 <= 10.4329:  # 3.23*3.23 = 10.4329
                            score = 1.0*iocc*kocc
                        else:
                            score = (1-(d2**0.5 - 3.23)/1.4)*iocc*kocc
                        total_score = total_score + score
            elif abs(ires_num - kres_num) > 4:
                for iatom in atom_in_ires:
                    (ix, iy, iz, iocc) = dict_coord[ires][iatom]
                    for katom in atom_in_kres:
                        (kx, ky, kz, kocc) = dict_coord[kres][katom]
                        d2 = (ix-kx)**2 + (iy-ky)**2 + (iz-kz)**2
                        if d2 >= 21.4369:  # 4.63*4.63 = 21.4369
                            score = 0
                        elif d2 <= 10.4329:  # 3.23*3.23 = 10.4329
                            score = 1.0*iocc*kocc
                        else:
                            score = (1-(d2**0.5 - 3.23)/1.4)*iocc*kocc
                        total_score = total_score + score
            contact_score[ires][kres] = total_score
    return contact_score

def main(pdbbase):
    contact = calc_contact(pdbbase)
    outf = pdbbase + '.cscore'
    fout = open(outf, 'w')
    for a_res in contact:
        b_res_list = contact[a_res].keys()
        for b_res in b_res_list:
            score = contact[a_res][b_res]
            if score > 0:
                fout.write('%-12s\t%-12s%10.6f\n' %(a_res, b_res, score))
    fout.close()

fin = sys.argv[1]
main(fin)