UNRES-Dock example output
[unres.git] / examples / UNRES-Dock / output / prepare_all.csh
1 #!/bin/csh -f
2
3 # Provide names of the pdb files
4 set protein1 = 1UDI_A.pdb
5 set protein2 = 1UDI_B.pdb
6
7 # Provide number of steps in milions
8 set nstep = 0.05
9
10 # Provide number of clusters
11 set nclust = 10
12
13 # Provide temperature of clustering
14 set temper = 300
15
16 # Provide UNRES instalation pathway
17 set UNRES = /users2/vetinari/unres_hom_m/unres
18
19
20
21 ################################################################################
22 ################################################################################
23 ################################################################################
24 sed -e "s=UUU=$UNRES=" start_unres_example.mat > start_unres.mat
25 sed -e "s=UUU=$UNRES=" start_wham_example.mat > start_wham.mat
26 sed -e "s=UUU=$UNRES=" start_cluster_example.mat > start_cluster.mat
27
28
29 # Generation of initial orientations
30 echo "Generation of initial orientations"
31 echo "TER" > TER
32 #gfortran -O3 generator_v13b.f -o generator_v13b
33 #./generator_v13b $protein1 $protein2 0 > log
34 grep -v "TER" $protein1 > prot1
35 grep -v "TER" $protein2 > prot2
36 cat prot1 "TER" prot2 > reference.pdb
37 rm TER prot1 prot2
38
39
40 ################################################################################
41 # Setting up and running MREMD
42 echo "Setting up and running MREMD"
43 set nstep2 = `echo $nstep \* 1000000 | bc | awk -F "." '{print $1}'`
44 set box = `grep -A1 "Boxsize (suggested)" log | tail -1 | cut -c 3-8`
45 sed "s/boxxx/$box/g" unres_example.inp > unres_dock.inp
46 sed "s/NNN/$nstep2/g" unres_dock.inp > unres_dock.in
47 mv unres_dock.in unres_dock.inp
48
49 #qsub start_unres.mat
50
51
52
53 ################################################################################
54 # Setting up and running WHAM analysis
55 echo "Setting up and running WHAM analysis"
56 #set resnum = `grep -A 1 "PDB data will be read from file" dimer.out_GB000 | tail -1 | awk '{print $2}'`
57 #set resnum = `echo $resnum + 5 | bc`
58
59 #grep -A $resnum "PDB data will be read from file" dimer.out_GB000 | grep -B 9999 "nsup=" | egrep -v "PDB|Nres:|nsup=" | awk '{print $3}' > seq_raw
60
61 egrep "CA|TER" reference.pdb > temp.pdb
62 set res1 = `grep -B 9999 "TER" temp.pdb | grep -v TER | grep -c CA`
63 set res2 = `grep -A 9999 "TER" temp.pdb | grep -v TER | grep -c CA`
64 set start1 = `head -1 temp.pdb | awk '{print $4}'`
65 set end2 = `tail -1 temp.pdb | awk '{print $4}'`
66 set start2 = `grep -A 1 TER temp.pdb | tail -1 | awk '{print $4}'`
67 set end1 = `grep -B 1 TER temp.pdb | head -1 | awk '{print $4}'`
68
69 echo First protein has $res1 residues, while second molecule has $res2 residues.
70 echo $start1 $end1 $start2 $end2
71
72 rm seq_raw
73
74 if ($start1 != GLY) then
75 echo D > seq_raw
76 endif
77
78 head -$res1 temp.pdb | awk '{print $4}' >> seq_raw
79
80 if ($end1 != GLY) then
81 echo D >> seq_raw
82 endif
83
84 if ($start2 != GLY) then
85 echo D >> seq_raw
86 endif
87
88 tail -$res2 temp.pdb | awk '{print $4}' >> seq_raw
89
90 if ($end2 != GLY) then
91 echo D >> seq_raw
92 endif
93
94
95 set resnum = `grep -c "." seq_raw`
96
97 echo System contains $resnum residues including dummy atoms.
98
99 set i = 1
100 set j = 1
101 set seqq = ""
102 rm seq
103 while ($i <= $resnum)
104
105 set aa = `head -$i seq_raw | tail -1`
106 if ($aa == D) then
107 set seqq = ${seqq}X
108 else if ($aa == ALA) then
109 set seqq = ${seqq}A
110 else if ($aa == ARG) then
111 set seqq = ${seqq}R
112 else if ($aa == ASN) then
113 set seqq = ${seqq}N
114 else if ($aa == ASP) then
115 set seqq = ${seqq}D
116 else if ($aa == CYS) then
117 set seqq = ${seqq}C
118 else if ($aa == GLU) then
119 set seqq = ${seqq}E
120 else if ($aa == GLN) then
121 set seqq = ${seqq}Q
122 else if ($aa == GLY) then
123 set seqq = ${seqq}G
124 else if ($aa == HIS) then
125 set seqq = ${seqq}H
126 else if ($aa == ILE) then
127 set seqq = ${seqq}I
128 else if ($aa == LEU) then
129 set seqq = ${seqq}L
130 else if ($aa == LYS) then
131 set seqq = ${seqq}K
132 else if ($aa == MET) then
133 set seqq = ${seqq}M
134 else if ($aa == PHE) then
135 set seqq = ${seqq}F
136 else if ($aa == PRO) then
137 set seqq = ${seqq}P
138 else if ($aa == SER) then
139 set seqq = ${seqq}S
140 else if ($aa == THR) then
141 set seqq = ${seqq}T
142 else if ($aa == TRP) then
143 set seqq = ${seqq}W
144 else if ($aa == TYR) then
145 set seqq = ${seqq}Y
146 else if ($aa == VAL) then
147 set seqq = ${seqq}V
148 endif
149
150 if ($j == 80) then
151 echo $seqq >> seq
152 set seqq = ""
153 set j = 0
154 endif
155
156 set i = `echo $i + 1 | bc`
157 set j = `echo $j + 1 | bc`
158 end
159 echo $seqq >> seq
160 rm seq_raw
161
162 set range1 = `echo $nstep2 / 2500 | bc | awk -F "." '{print $1}'`
163 set range1 = `echo $range1 | bc`
164 set div = `echo 12000 \* 10 / 24 | bc | awk -F "." '{print $1}'`
165 set range0 = `echo $range1 - $div | bc | awk -F "." '{print $1}'`
166 if ($range0 < 0 ) then
167 set range0 = 10
168 endif
169
170 #echo $nstep2 $div
171 echo WHAM will analyze range of snapshots from $range0 to $range1.
172
173 sed "s/xxxxx/$box/g" wham_example.inp > wham.inp
174 sed "s/XXX/$resnum/g" wham.inp > wham
175 head -4 wham > wham.in
176 grep -B 1 -A 999 "HOMOL_DIST" wham_example.inp > wham
177 cat wham.in seq wham > wham.inp
178 sed "s/BBB/$range0/g" wham.inp > wham
179 sed "s/CCC/$range1/g" wham > wham.inp
180 rm wham wham.in
181 #qsub start_wham.mat
182
183
184 ################################################################################
185 # Setting up and running cluster analysis
186 echo "Setting up and running cluster analysis"
187
188 sed "s/XXX/$resnum/g" cluster_example.inp > clust
189 head -8 clust > clust.in
190 grep -B 1 -A 999 "HOMOL_DIST" cluster_example.inp > clus
191 cat clust.in seq clus > cluster.inp
192 sed "s/nclust=ZZZ/nclust=$nclust/g" cluster.inp > clut
193 sed "s/temper=AAA/temper=$temper/g" clut > cluster.inp
194 rm clust clust.in clus
195 #qsub start_cluster.mat