-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathmd_multianalyze.py
167 lines (142 loc) · 5.81 KB
/
md_multianalyze.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
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
# %%
import os
import shutil
import configparser
import subprocess
import Bio.PDB as pdb
def make_dirs_for_results(id_list: list) -> list:
"""測定結果を格納するディレクトリを作成する
Args:
id_list (list): config.iniなどから読み込んだIDのリスト
Returns:
list: 指定したIDの解析結果を格納する用のディレクトリリスト
"""
pwd = os.path.abspath(os.path.dirname(__file__))
pdb_dir = os.path.join(pwd, "init_pdb")
dihed_dir = os.path.join(pwd, "dihed")
dists_dir = os.path.join(pwd, "dists")
os.makedirs(name=pdb_dir, exist_ok=True)
os.makedirs(name=dihed_dir, exist_ok=True)
os.makedirs(name=dists_dir, exist_ok=True)
return [pdb_dir, dihed_dir, dists_dir]
def dihed_to_in_file(dihed_dict: dict,
resi: str,
id: str):
"""辞書を読み取り、二面角測定用のinputファイルを作成する
Args:
dists_dict (dict): 測定したい距離の辞書 (key: 結果名, value: 距離を構成する2原子のリスト)
resi (str): 二面角を測定したい残基
id (str): 対象となるタンパク質のID
"""
dihed_script = open(f"{id}/amber/pr/dihed.in", "w")
dihed_script.write(f"""trajin ./001/mdcrd 1 last 50
trajin ./002/mdcrd 1 last 50
""")
dihed_script.close()
for name, atoms in dihed_dict.items():
dihed_script = open(f"{id}/amber/pr/dihed.in", "a")
dihed_script.write(f"""dihedral {name} :{resi}@{atoms[0]} :{resi}@{atoms[1]} :{resi}@{atoms[2]} :{resi}@{atoms[3]} out dihed.txt range360
strip :SOD,CLA,WAT,TIP3,Na,Cl-
""")
dihed_script.close()
def dists_to_in_file(dists_dict: dict,
resi: str,
id: str):
"""辞書を読み取り、距離測定用のinputファイルを作成する
Args:
dists_dict (dict): 測定したい距離の辞書 (key: 結果名, value: 距離を構成する2原子のリスト)
resi (str): 距離を測定したい残基
id (str): 対象となるタンパク質のID
"""
dists_script = open(f"{id}/amber/pr/dists.in", "w")
dists_script.write(f"""trajin ./001/mdcrd 1 last 50
trajin ./002/mdcrd 1 last 50
""")
dists_script.close()
for name, atoms in dists_dict.items():
dists_script = open(f"{id}/amber/pr/dists.in", "a")
dists_script.write(f"""distance {name} :{resi}@{atoms[0]} :{resi}@{atoms[1]} out dists.txt
""")
dists_script.close()
def make_in_files(id_list: list,
dihed_dict: dict,
dists_dict: dict,
res_name: str) -> None:
"""初期構造pdbファイルから対象残基の番号を取得し、二面角と距離測定用のinputファイルを作成する
Args:
id_list (list): config.iniなどから読み込んだIDのリスト
dists_dict (dict): 測定したい距離の辞書 (key: 結果名, value: 距離を構成する2原子のリスト)
dists_dict (dict): 測定したい距離の辞書 (key: 結果名, value: 距離を構成する2原子のリスト)
res_name (str): 測定したい残基名
"""
for id in id_list:
pdbfile = os.path.join("./", id, "amber", "pr", "init.pdb")
pdb_parser = pdb.PDBParser()
try:
struct = pdb_parser.get_structure("all", pdbfile)
except FileNotFoundError:
continue
chain = struct[0].get_list()[0]
res = chain.get_list()
for r in res:
if r.get_resname() == res_name:
resi = str(r.get_id()[1])
dihed_to_in_file(dihed_dict=dihed_dict,
resi=resi,
id=id)
dists_to_in_file(dists_dict=dists_dict,
resi=resi,
id=id)
# config.iniの読み込みと結果を格納するディレクトリの作成
config = configparser.ConfigParser(strict=False, allow_no_value=True)
config.read("config.ini")
id = [k.upper() for k, v in config.items("ID")]
pdb_dir, dihed_dir, dists_dir = make_dirs_for_results(id_list=id)
# pdbの初期構造とトラジェクトリを取得
for i in id:
get_init_cmd = f"""
cd ./{i}/amber/pr
cpptraj -i ./trajfix.in -p ../../top/leap.parm7
cp ./init.pdb {pdb_dir}/{i}.pdb
cp ./traj.trr {pdb_dir}/{i}.trr
cd ../../../"""
subprocess.run(get_init_cmd, shell=True)
# 解析対象を定義
dihed_dict = {"O1~C3": ["O1", "C1", "C2", "C3"],
"C1~C5": ["C1", "C2", "C3", "C5"],
"C2~C6": ["C2", "C3", "C5", "C6"],
"C3~C7": ["C3", "C5", "C6", "C7"],
"C5~C8": ["C5", "C6", "C7", "C8"],
"C6~C9": ["C6", "C7", "C8", "C9"],
"C7~C11": ["C7", "C8", "C9", "C11"],
"C8~C12": ["C8", "C9", "C11", "C12"],
"C9~C13": ["C9", "C11", "C12", "C13"],
"fin1_1": ["C1", "C2", "C3", "C4"],
"fin1_2": ["C4", "C3", "C5", "C6"],
"fin2_1": ["C6", "C7", "C8", "C10"],
"fin2_2": ["C10", "C8", "C9", "C11"],
"fin3_1": ["C11", "C12", "C13", "C14"],
"fin3_2": ["C14", "C13", "C12", "C11"], }
dists_dict = {"C1-C5": ["C1", "C7"],
"C1-C6": ["C1", "C8"],
"C1-C10": ["C1", "C12"],
"C1-C11": ["C1", "C13"]}
make_in_files(id_list=id,
res_name="FPP",
dihed_dict=dihed_dict,
dists_dict=dists_dict)
# 解析の実行
for i in id:
get_dihed_cmd = f"""
cd ./{i}/amber/pr
cpptraj -i ./dihed.in -p ../../top/leap.parm7
cp ./dihed.txt {dihed_dir}/{i}.txt
cd ../../../"""
subprocess.run(get_dihed_cmd, shell=True)
get_dists_cmd = f"""
cd ./{i}/amber/pr
cpptraj -i ./dists.in -p ../../top/leap.parm7
cp ./dists.txt {dists_dir}/{i}.txt
cd ../../../"""
subprocess.run(get_dists_cmd, shell=True)
# %%