md-analysis/dump_PCA_CU_code.py at main · EngAsal/md-analysis · 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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
base_dir = ("/home/pghw87/Documents/md-sim/5ue6/MDAniA/data/")
os.chdir(base_dir)
apotrimer = mda.Universe("trimer/trimer.gro", "trimer/trimer.xtc")
holotrimer3_T1 = mda.Universe("trimer3-T1/bonded/trimer.gro", "trimer3-T1/bonded/trimer.xtc")
holotrimer3_T2 = mda.Universe("trimer3-T2/bonded/trimer.gro", "trimer3-T2/bonded/trimer.xtc")
holotrimer6 = mda.Universe("trimer6/bonded/trimer.gro", "trimer6/bonded/trimer.xtc")
# %%
CU_chainA_apotri = apotrimer.select_atoms("resid 86 91 126 127 135 140 591 and not name H*")
CU_chainB_apotri = apotrimer.select_atoms("resid 396 401 436 437 445 450 901 and not name H*")
# CU_chainC_apotri = apotrimer.select_atoms("resid 706 711 746 747 755 760 281 and not name H*")
CU_chainA_apotri.write("CU_chainA_apotri.pdb")
CU_chainB_apotri.write("CU_chainB_apotri.pdb")
# CU_chainC_apotri.write("CU_chainC_apotri.pdb")
# Chain C: first take the 706..760 set, then append resid 281
CU_chainC_main_apotri = apotrimer.select_atoms(
    "resid 706 711 746 747 755 760 and not name H*"
)
CU_chainC_281_apotri = apotrimer.select_atoms(
    "resid 281 and not name H*"
)

CU_chainC_apotri = CU_chainC_main_apotri + CU_chainC_281_apotri
CU_chainC_apotri.write("CU_chainC_apotri.pdb")

with mda.Writer("CU_chainA_apotri.xtc", CU_chainA_apotri.n_atoms) as W:
    for ts in apotrimer.trajectory:
        W.write(CU_chainA_apotri)
with mda.Writer("CU_chainB_apotri.xtc", CU_chainB_apotri.n_atoms) as W:
    for ts in apotrimer.trajectory:
        W.write(CU_chainB_apotri)
with mda.Writer("CU_chainC_apotri.xtc", CU_chainC_apotri.n_atoms) as W:
    for ts in apotrimer.trajectory:
        W.write(CU_chainC_apotri)

CU_chainA_apotri_universe = mda.Universe("CU_chainA_apotri.pdb", "CU_chainA_apotri.xtc")
CU_chainB_apotri_universe = mda.Universe("CU_chainB_apotri.pdb", "CU_chainB_apotri.xtc")
CU_chainC_apotri_universe = mda.Universe("CU_chainC_apotri.pdb", "CU_chainC_apotri.xtc")
#%%
ref = CU_chainA_apotri_universe
ref.trajectory[0]
align.AlignTraj(CU_chainA_apotri_universe, ref, select="name CA", filename="CU_chainA_apotri_aligned.xtc").run()
align.AlignTraj(CU_chainB_apotri_universe, ref, select="name CA", filename="CU_chainB_apotri_aligned.xtc").run()
align.AlignTraj(CU_chainC_apotri_universe, ref, select="name CA", filename="CU_chainC_apotri_aligned.xtc").run()

#%%
CU_chainA_holotri3 = holotrimer3.select_atoms("resid 86 91 126 127 135 140 592 and not name H*")
CU_chainB_holotri3 = holotrimer3.select_atoms("resid 397 402 437 438 446 451 903 and not name H*")
# CU_chainC_holotri3 = holotrimer3.select_atoms("resid 708 713 748 749 757 762 281 and not name H*")
# Chain C: first take the 706..760 set, then append resid 281
CU_chainC_main_holotri3 = holotrimer3.select_atoms(
    "resid 708 713 748 749 757 762 and not name H*"
)
CU_chainC_281_holotri3 = holotrimer3.select_atoms(
    "resid 281 and not name H*"
)

CU_chainC_holotri3 = CU_chainC_main_holotri3 + CU_chainC_281_holotri3

CU_chainA_holotri3.write("CU_chainA_holotri3.pdb")
CU_chainB_holotri3.write("CU_chainB_holotri3.pdb")
CU_chainC_holotri3.write("CU_chainC_holotri3.pdb")

with mda.Writer("CU_chainA_holotri3.xtc", CU_chainA_holotri3.n_atoms) as W:
    for ts in holotrimer3.trajectory:
        W.write(CU_chainA_holotri3)
with mda.Writer("CU_chainB_holotri3.xtc", CU_chainB_holotri3.n_atoms) as W:
    for ts in holotrimer3.trajectory:
        W.write(CU_chainB_holotri3)
with mda.Writer("CU_chainC_holotri3.xtc", CU_chainC_holotri3.n_atoms) as W:
    for ts in holotrimer3.trajectory:
        W.write(CU_chainC_holotri3)
#%%
CU_chainA_holotri3_universe = mda.Universe("CU_chainA_holotri3.pdb", "CU_chainA_holotri3.xtc")
CU_chainB_holotri3_universe = mda.Universe("CU_chainB_holotri3.pdb", "CU_chainB_holotri3.xtc")
CU_chainC_holotri3_universe = mda.Universe("CU_chainC_holotri3.pdb", "CU_chainC_holotri3.xtc")

align.AlignTraj(CU_chainA_holotri3_universe, ref, select="name CA", filename="CU_chainA_holotri3_aligned.xtc").run()
align.AlignTraj(CU_chainB_holotri3_universe, ref, select="name CA", filename="CU_chainB_holotri3_aligned.xtc").run()
align.AlignTraj(CU_chainC_holotri3_universe, ref, select="name CA", filename="CU_chainC_holotri3_aligned.xtc").run()

# %%
CU_chainA_holotri6 = holotrimer6.select_atoms("resid 86 91 126 127 135 140 591 and not name H*")
CU_chainB_holotri6 = holotrimer6.select_atoms("resid 396 401 436 437 445 450 901 and not name H*")
# CU_chainC_holotri6 = holotrimer6.select_atoms("resid 706 711 746 747 755 760 281 and not name H*")
# Chain C: first take the 706..760 set, then append resid 281
CU_chainC_main = holotrimer6.select_atoms(
    "resid 706 711 746 747 755 760 and not name H*"
)
CU_chainC_281 = holotrimer6.select_atoms(
    "resid 281 and not name H*"
)

CU_chainC_holotri6 = CU_chainC_main + CU_chainC_281

CU_chainA_holotri6.write("CU_chainA_holotri6.pdb")
CU_chainB_holotri6.write("CU_chainB_holotri6.pdb")
CU_chainC_holotri6.write("CU_chainC_holotri6.pdb")

with mda.Writer("CU_chainA_holotri6.xtc", CU_chainA_holotri6.n_atoms) as W:
    for ts in holotrimer6.trajectory:
        W.write(CU_chainA_holotri6)
with mda.Writer("CU_chainB_holotri6.xtc", CU_chainB_holotri6.n_atoms) as W:
    for ts in holotrimer6.trajectory:
        W.write(CU_chainB_holotri6)
with mda.Writer("CU_chainC_holotri6.xtc", CU_chainC_holotri6.n_atoms) as W:
    for ts in holotrimer6.trajectory:
        W.write(CU_chainC_holotri6)
# %%
CU_chainA_holotri6_universe = mda.Universe("CU_chainA_holotri6.pdb", "CU_chainA_holotri6.xtc")
CU_chainB_holotri6_universe = mda.Universe("CU_chainB_holotri6.pdb", "CU_chainB_holotri6.xtc")
CU_chainC_holotri6_universe = mda.Universe("CU_chainC_holotri6.pdb", "CU_chainC_holotri6.xtc")

align.AlignTraj(CU_chainA_holotri6_universe, ref, select="name CA", filename="CU_chainA_holotri6_aligned.xtc").run()
align.AlignTraj(CU_chainB_holotri6_universe, ref, select="name CA", filename="CU_chainB_holotri6_aligned.xtc").run()
align.AlignTraj(CU_chainC_holotri6_universe, ref, select="name CA", filename="CU_chainC_holotri6_aligned.xtc").run()
#%%
cu_dir = ("/home/pghw87/Documents/md-sim/5ue6/MDAniA/data/pca/CU_site")
os.chdir(cu_dir)
# Combine trajectories into multi_pdb
traj_files = [
    "CU_chainA_apotri_aligned.xtc",
    "CU_chainB_apotri_aligned.xtc",
    "CU_chainC_apotri_aligned.xtc",
    "CU_chainA_holotri3_aligned.xtc",
    "CU_chainB_holotri3_aligned.xtc",
    "CU_chainC_holotri3_aligned.xtc",
    "CU_chainA_holotri6_aligned.xtc",
    "CU_chainB_holotri6_aligned.xtc",
    "CU_chainC_holotri6_aligned.xtc"
]
multi_pdb = mda.Universe('CU_chainA_apotri.pdb', *traj_files)
with mda.Writer('multipdb_CU_site.pdb', multi_pdb.atoms.n_atoms) as W:
    for ts in multi_pdb.trajectory:
        W.write(multi_pdb)