5


rdplot
#!/bin/sh

#==================================== Description =====================================#
### (Leave help function at beginning of script, as it's also the doc string.)       ###
function help { cat <<'HELP_'
    
================================================
Read-depth plot from .xls file
================================================
                                                                                      
Primitive script with little error checking. A gnuplot script is created in a
here-document and passed directly to gnuplot, avoiding unnecessary file I/O.
Input data file is of type .xls generated by "sol" workflow.

Usage:
    $ rdplot [-h] <.xls input file> <max. range> <bin size> <plot color> [plot title]

 Options:

     -h             help

 Arguments (strict order)

     $in = .xls data file
         (e.g., "sty.xls")

     $maxrange = maximum range value in units of bin size
          (e.g., 5000 for sty chromosome)

     $bin = bin size
          (e.g., 1000 for sty chromosome)

     $color = plot color
          (e.g., "blue")

     $title = plot title (optional; defaults to data file name)
         (e.g., "TT26528")

 Author:

    Eric Kofoid, 4/2/2012
    eckofoid@ucdavis.edu
    Last change: 8/30/12

HELP_
}
#======================================================================================#

## Check arguments and process when possible

if [ $# = 0 ] || [ $1 = "-h" ] || [ $1 = "--help" ]
then
    help
    exit

elif [ $# -lt 4 ]
then
    echo "Insufficient parameters"
    exit
else
    in=$1
    maxrange=$2
    bin=$3
    color=$4
    outf=`basename ${in%.*} 2>/dev/null`.ps
fi

if [ $# = 5 ]
then
    title=$5
else
    title=$in
fi

if [ $bin = 1000 ]
then
    xunit="KB"
else
    xunit=$bin" Bases"
fi

###### Adjusts ymax to a convenient multiple of 10 #####
ymax=$(awk -v mxval=0 '{if ($1>mxval) mxval=$1};END {printf "%d", mxval/'$bin'}' $in)
ynorm=1
foo=$ymax
while [[ $foo -ne 0 ]]
do
    (( ynorm=10*ynorm ))
    (( foo=foo/10 ))
done
(( ynorm=ynorm/10 ))
(( ymax = ( ( ynorm + ymax - 1 ) / ynorm ) * ynorm ))
if (( (ymax/ynorm)<3 )) # too few labelled tics
then
    (( ynorm=ynorm/2 ))
fi
    
########################################################

gnuplot <<SCRIPT
unset key
set xrange [0:$maxrange]
set yrange [0:$ymax]
set xlabel "Chromosomal Position ($xunit)"
set ylabel "Read Depth"
set xtics nomirror out
set ytics 0,$ynorm out
set mytics 5
set title "$title"
set size ratio 0.1
# uncomment following line for ssh imaging using X11:
# set terminal x11 persist
plot "$in" using (\$1/$bin) w boxes lt rgb "$color"
set t push
set t postscript
set o "$outf"
replot
set o
set t pop
SCRIPT