-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmakecst.py
executable file
·117 lines (102 loc) · 5.55 KB
/
makecst.py
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
115
116
117
#!/usr/bin/env python3
'''
ver 1.0.1
This script works on Python 3.6 or more.
'''
import argparse
parser = argparse.ArgumentParser(description='Generate cst file for Rosetta AbinitioRelax from a fasta and a map file.')
parser.add_argument("-f","--fasta", dest="fasta",
help="Required. Input your fasta file to generate cst file for Rosetta AbinitioRelax", required=True)
parser.add_argument("-m","--map", dest="map",
help="Contact map file. Default is 'hmmrenumbered.map'.",
required=False, default="hmmrenumbered.map")
parser.add_argument("-c","--coarse", dest="coarse",
help="Output file for coarse cst file (-cst_file). Default is 'coarse.cst'.",
required=False, default="coarse.cst")
parser.add_argument("-o","--full", dest="full",
help="Output file for full-atom cst file (-cst_fa_file). Default is 'full.cst'.",
required=False, default="full.cst")
args = parser.parse_args()
aa = "G A S V C T P D N I L E Q M H K F Y R W".split()
rawx0 = """4.467,5.201,5.510,5.671,5.777,5.619,6.140,6.135,6.321,6.413,6.554,7.036,7.297,7.383,7.472,8.216,7.966,9.098,9.166,8.966
,5.381,5.829,5.854,6.057,5.982,6.412,6.388,6.766,6.587,6.707,7.124,7.583,7.605,7.591,8.327,8.162,9.121,9.365,9.252
,,6.190,6.567,6.590,6.450,6.937,6.760,7.081,7.142,7.394,7.483,7.807,8.010,8.051,8.792,8.694,9.594,9.753,9.770
,,,6.759,6.941,6.791,7.063,6.972,7.219,7.441,7.633,7.404,8.008,8.335,8.179,8.077,9.057,9.442,9.513,10.021
,,,,6.426,6.801,7.157,6.985,7.205,7.476,7.685,7.449,7.962,8.265,8.422,8.494,9.026,9.362,9.460,9.752
,,,,,6.676,7.062,6.971,7.159,7.442,7.642,7.628,8.055,8.397,8.221,8.715,9.030,9.813,9.764,9.980
,,,,,,7.288,7.321,7.497,7.554,7.751,7.938,8.308,8.247,8.537,9.198,8.895,9.965,10.266,9.719
,,,,,,,8.001,7.672,7.472,7.696,8.945,8.601,8.401,8.634,9.306,9.111,9.979,10.123,9.867
,,,,,,,,7.682,7.631,7.889,8.485,8.502,8.550,8.672,9.319,9.168,10.039,10.135,9.976
,,,,,,,,,8.096,8.342,7.949,8.302,8.874,8.523,8.329,9.602,9.719,9.746,10.470
,,,,,,,,,,8.522,8.077,8.480,9.122,8.676,8.479,9.900,9.889,9.852,10.707
,,,,,,,,,,,9.863,9.328,8.870,9.454,9.842,9.403,10.544,10.713,10.303
,,,,,,,,,,,,9.074,9.102,9.391,9.667,9.506,10.534,10.610,10.429
,,,,,,,,,,,,,9.530,9.396,9.096,10.253,10.400,10.250,11.110
,,,,,,,,,,,,,,10.606,9.582,9.602,10.843,10.879,10.661
,,,,,,,,,,,,,,,10.662,9.344,10.627,11.322,10.136
,,,,,,,,,,,,,,,,10.903,10.999,10.577,11.758
,,,,,,,,,,,,,,,,,11.536,11.615,11.807
,,,,,,,,,,,,,,,,,,12.050,11.355
,,,,,,,,,,,,,,,,,,,12.806""".split("\n")
rawm = """0.017,,,,,,,,,,,,,,,,,,,
0.269,0.262,,,,,,,,,,,,,,,,,,
0.153,0.291,0.292,,,,,,,,,,,,,,,,,
0.107,0.312,0.205,0.145,,,,,,,,,,,,,,,,
0.129,0.394,0.240,0.173,0.178,,,,,,,,,,,,,,,
0.120,0.378,0.214,0.138,0.181,0.188,,,,,,,,,,,,,,
0.245,0.399,0.321,0.298,0.259,0.320,0.339,,,,,,,,,,,,,
0.193,0.289,0.323,0.287,0.299,0.307,0.416,0.392,,,,,,,,,,,,
0.169,0.349,0.305,0.232,0.240,0.262,0.334,0.337,0.249,,,,,,,,,,,
0.179,0.214,0.342,0.242,0.295,0.259,0.336,0.341,0.341,0.321,,,,,,,,,,
0.125,0.250,0.287,0.179,0.206,0.190,0.317,0.348,0.279,0.261,0.198,,,,,,,,,
0.249,0.340,0.446,0.510,0.538,0.409,0.475,0.354,0.423,0.453,0.475,0.389,,,,,,,,
0.216,0.356,0.408,0.359,0.347,0.378,0.410,0.357,0.373,0.406,0.411,0.450,0.436,,,,,,,
0.255,0.394,0.369,0.295,0.439,0.292,0.388,0.361,0.310,0.327,0.318,0.511,0.498,0.457,,,,,,
0.206,0.380,0.435,0.383,0.203,0.417,0.457,0.325,0.289,0.379,0.401,0.443,0.401,0.342,0.333,,,,,
0.358,0.550,0.445,0.634,0.521,0.464,0.550,0.343,0.398,0.582,0.591,0.434,0.521,0.611,0.714,0.738,,,,
0.219,0.260,0.394,0.246,0.286,0.264,0.425,0.351,0.393,0.347,0.260,0.512,0.451,0.377,0.542,0.441,0.460,,,
0.267,0.443,0.467,0.535,0.585,0.430,0.506,0.676,0.586,0.589,0.611,0.469,0.547,0.661,0.554,0.704,0.767,0.855,,
0.334,0.485,0.483,0.514,0.491,0.477,0.506,0.327,0.372,0.557,0.578,0.363,0.535,0.641,0.595,0.648,0.738,0.822,0.704,
0.239,0.290,0.497,0.271,0.417,0.315,0.462,0.475,0.458,0.397,0.331,0.493,0.490,0.397,0.458,0.470,0.447,0.684,0.889,0.473""".split("\n")
x0 = [[0.0]*20 for _ in range(20)]
for i in range(20):
ln = rawx0[i].split(",")
for j in range(20):
if len(ln[j]) > 0:
x0[i][j] = x0[j][i] = float(ln[j])
m = [[0.0]*20 for _ in range(20)]
for i in range(20):
ln = rawm[i].split(",")
for j in range(20):
if len(ln[j]) > 0:
m[i][j] = m[j][i] = 1.0/float(ln[j])
def readfasta(fastafile):
fasta = ''
with open(fastafile) as f:
for line in f:
if line.startswith('>'):
continue
else:
fasta += line.strip()
return fasta
def main():
#ここ二行の入力があまりにも雑
# ファイル書き出し
with open(args.coarse, mode='w') as c:
with open(args.full, mode='w') as fu:
fasta = readfasta(args.fasta)
for ln in open(args.map):
lnt = ln.split()
a = int(lnt[1]) - 1
b = int(lnt[2]) - 1
if a*b == 0 : continue
weight = float(lnt[3])
aa_a = fasta[a]; aa_b = fasta[b]
out_x0 = x0[ aa.index(aa_a) ][ aa.index(aa_b) ]
out_m = m[ aa.index(aa_a) ][ aa.index(aa_b) ]
cb_a = "CB" if not aa_a == 'G' else "CA"
cb_b = "CB" if not aa_b == 'G' else "CA"
fu.write(f"AtomPair {cb_a} {a+1} {cb_b} {b+1} SCALARWEIGHTEDFUNC {weight} SIGMOID {out_x0} {out_m:.4}\n")
c.write(f"AtomPair {cb_a} {a+1} {cb_b} {b+1} SCALARWEIGHTEDFUNC {weight} BOUNDED 0 {out_x0} 1 0.5 #\n")
if __name__ == '__main__':
main()