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