awsemmd_script/longRangeInteraction.py at master · luwei0917/awsemmd_script · 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
112
113
114
#!/usr/bin/python
# Modified from GetCACADistancesFile.py

import sys
from VectorAlgebra import *


def calc_dis(p1, p2):
    v = vector(p1, p2)
    return vabs(v)


def three2one(prot):
    """ translate a protein sequence from 3 to 1 letter code"""

    code = {"GLY": "G", "ALA": "A", "LEU": "L", "ILE": "I",
            "ARG": "R", "LYS": "K", "MET": "M", "CYS": "C",
            "TYR": "Y", "THR": "T", "PRO": "P", "SER": "S",
            "TRP": "W", "ASP": "D", "GLU": "E", "ASN": "N",
            "GLN": "Q", "PHE": "F", "HIS": "H", "VAL": "V",
            "M3L": "K", "MSE": "M", "CAS": "C"}

    newprot = ""
    for a in prot:
        newprot += code.get(a, "?")

    return newprot


def checkIfNative(ires, jres):
    xyz_CAi = ires['CA'].get_coord()
    xyz_CAj = jres['CA'].get_coord()
    v = vector(xyz_CAi, xyz_CAj)
    r = vabs(v)
    if r < 12.0:
        return True
    else:
        return False

if len(sys.argv) == 1:
    print "\nGetCACADistancesFile_multi.py PDB_Id Output_file \n"
    print "-s\tSplit into files for each chain"
#    sys.argv.append("1BG8")
    exit()

from Bio.PDB.PDBParser import PDBParser

p = PDBParser(PERMISSIVE=1)

struct_id = sys.argv[1]
filename = struct_id + ".pdb"

splite = False

output_fn = ""
if len(sys.argv) > 3:
    output_fn = sys.argv[2]
if output_fn[-4:] == ".dat":
    output_fn = output_fn[:-4]
SeqLen = int(sys.argv[3])

if output_fn != "" and not splite:
    out = open((output_fn + ".dat"), 'w')

k_bond = 100
k_angle = 20
k_dihedral = [1, 0.5]
epsilon = 1
epsilon2 = 1
sigma0 = 4

xyz_CA1 = []
xyz_CA2 = []
xyz_CA3 = []
xyz_CA4 = []

s = p.get_structure(struct_id, filename)
chains = s[0].get_list()
sequence = []
dis = [[0 for x in range(SeqLen)] for y in range(SeqLen)]
all_res = []
index = []

for chain in chains:
    for resid in chain:
        is_regular_res = resid.has_id('CA')
        res_id = resid.get_id()[0]
        if (res_id == ' ' or res_id == 'H_MSE' or res_id == 'H_M3L' or res_id == 'H_CAS') and is_regular_res:
            all_res.append(resid)
            sequence.append(resid.get_resname())
            index.append(int(resid.get_id()[1]))
print index

for i in range(0, len(all_res)):
    # dis.append([]);
    ires = all_res[i]
    xyz_CAi = ires['CA'].get_coord()
    for j in range(0, len(all_res)):
        jres = all_res[j]
        xyz_CAj = jres['CA'].get_coord()
        r = calc_dis(xyz_CAi, xyz_CAj)
        # print index[i]
        # print index[j]
        dis[index[i] - 1][index[j] - 1] = r

for ri in dis:
    for rij in ri:
        out.write(str(round(rij, 6)))
        out.write(' ')
    out.write('\n')
if not splite:
    out.write('\n')

out.close()