Prestin_Analysis/classifier/classifier.py at main · Becksteinlab/Prestin_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
import numpy as np
import pandas as pd
import MDAnalysis as mda
from MDAnalysis.analysis.rms import rmsd, RMSD
from MDAnalysis.analysis.align import alignto, rotation_matrix
from MDAnalysis.lib.transformations import rotation_from_matrix
from tqdm import tqdm

from glob import glob


pdb_path = '/nfs/homes3/chenou/Project/Anton-prestin/charmm/charmm-production-0mV-296/step6.6_equilibration.gro'
xtc_path = '/nfs/homes5/Projects/SLC26/chenou/processed/compact2-eq/run-1/processed.xtc'


ref_1 = mda.Universe('/nfs/homes5/Projects/SLC26/chenou/data/clustering/CS/prestin-excompact-protein.pdb')
ref_2 = mda.Universe('/nfs/homes3/chenou/Project/Anton-prestin/charmm/charmm-production-0mV-296/step6.6_equilibration.gro')
ref_3 = mda.Universe('/nfs/homes5/Projects/SLC26/chenou/data/clustering/CS/prestin-excompact-protein.pdb')

u = mda.Universe(pdb_path, xtc_path)
size = len(u.select_atoms('protein'))

size_1 = len(ref_1.select_atoms('protein'))
size_2 = len(ref_2.select_atoms('protein'))
size_3 = len(ref_3.select_atoms('protein'))


TM_A_protomer   = 'name CA and bynum 1-' + str(int(size/2))
TM_B_protomer   = 'name CA and bynum ' + str(int(size/2)) + '-' + str(int(size))

# selections all CA

TM_CA_A = TM_A_protomer + ' and resid 76-155 167-504'
TM_CA_B = TM_B_protomer + ' and resid 76-155 167-504'

# selections gate CA 'protein and resid 207-307 437-504'

TM_CA_A_gate = TM_A_protomer + ' and resid 207-307 437-504'
TM_CA_B_gate = TM_B_protomer + ' and resid 207-307 437-504'

# selections core CA 'protein and resid 76-105 109-125 138-151 170-196 339-365 371-388 397-403 411-426'

TM_CA_A_core = TM_A_protomer + ' and resid 76-155 167-196 339-426'
TM_CA_B_core = TM_B_protomer + ' and resid 76-155 167-196 339-426'

# add a9, a6 core

TM_CA_A_a9_core = ' and resid 67-146 161-190 333-420'
TM_CA_A_a9_gate = ' and resid 201-301 431-498'

###########################################################################################

core_A = u.select_atoms(TM_CA_A_core)
core_B = u.select_atoms(TM_CA_B_core)

gate_A = u.select_atoms(TM_CA_A_gate)
gate_B = u.select_atoms(TM_CA_B_gate)


fit_1 = RMSD(u, reference=ref_1, select=TM_CA_B_gate, groupselections=[TM_CA_B_core])
fit_1.run(verbose=True, stop=10)
class_1 = fit_1.results.rmsd.T[3]


fit_2 = RMSD(u, reference=ref_2, select=TM_CA_B_gate, groupselections=[TM_CA_B_core])
fit_2.run(verbose=True, stop=10)
class_2 = fit_2.results.rmsd.T[3]


fit_3 = RMSD(u, reference=ref_3, select=TM_CA_B_gate, groupselections=[TM_CA_B_core])
fit_3.run(verbose=True, stop=10)
class_3 = fit_3.results.rmsd.T[3]

np.save('rmsd_1.npy', fit_1.results.rmsd.T)
np.save('rmsd_2.npy', fit_2.results.rmsd.T)
np.save('rmsd_3.npy', fit_3.results.rmsd.T)