#!/bin/bash
#Runs hole and generates plot of channel radius, sphpdb file and dotted surface. sphpdb file and dotsurface.pdb
#can be used together to generate a dotted surface colored by channel radius in Chimera.
if [ $# -lt 3 ]; then
    echo "Function: runs HOLE and converts output to format readable by Chimera."
    echo " "
    echo "Usage: hole_chimera_script [pdbfile] [end radius in Angstroems] [dot density (5-30)]"
    echo " "
    echo "The program qpt_conv will be run interactively as part of the script."
    echo "Select output format 'A' (ASCII) for the conversion and leave other options as defaults."
    echo " "
    echo "Expected output: dotsurf.pdb (dotted surface), sphpdb.pdb (channel axis pdb),"
    echo "channel_coord.txt (tabulated values of pore radius), and channel_coord.svg (graph of same)."
    exit 1
fi
hole > hole_out.txt <<eof
coord $1        ! Co-ordinates in pdb format
radius /Applications/hole2/rad/nussinov.rad       ! Use simple AMBER vdw radii
                                ! n.b. can use ~ in hole
!
! now optional cards
sphpdb hole_out.sph             ! pdb format output of hole sphere centre info
                                ! (for use in sph_process program)
endrad $2.                      ! This is the pore radius that is taken
                                ! as where channel ends. 5.0 Angstroms is good
                                ! for a narrow channel
eof
wait
grep 'ATOM' hole_out.sph | grep -v ' S-888 ' | sort -k2.22,2.25 -t'A' -n > sphpdb.pdb
sph_process -dotden $3 -nocen hole_out.sph dotsurface.qpt
wait
qpt_conv
wait
grep '^   4' dotsurface.asc | awk '{printf "ATOM %6.0f QSS SPH S-1 %13.3f %7.3f %7.3f 10.00 10.00 \n", NR, $2,$3,$4}' > dotsurface.pdb
egrep "mid-point|sampled" hole_out.txt | awk '{print $1 "  " $2}' > channel_coord.txt
gnuplot<<eof
set ylabel "Pore radius (Å)"
set xlabel "Channel coordinate (Å)"
set output "channel_coord.svg"
set term svg fname "Arial" fsize 12 round dashed
set xtics nomirror out scale 0.5
set ytics nomirror out scale 0.5
set border 3
set size ratio 0.3
plot "channel_coord.txt" using 1:2 with lines lc "blue" lw 3 notitle
eof