4


sol
#!/bin/sh

#######################################################
########       Reusable Document String:       ########
#######################################################
help() { cat<<eof
 =======================
  Solexa data reduction
 =======================
 Usage:
    sol <infile 1> <infile 2> <outfile> <refseq> <bin size>

 Infiles are type "bfq" (infiles) and "bfa" (refseq).
 Outfile defines outpath and basename for outfiles.
 Script generates "map", "pileup" and binned "xls" files;
 same-named files will be overwritten without warning.

 Eric Kofoid, 3/29/2011

eof
}
#######################################################

## Check arguments and process when possible

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

elif [ $# -lt 5 ]
then
    echo "Insufficient number of parameters"
    exit
else
    in1=$1
    in2=$2
    outf=$3
    refseq=$4
    bins=$5
fi

if ! [ -f $in1 ]
then
    echo "Infile 1 does not exist!"
    exit

elif [ ${in1##*.} != "bfq" ]
then
    echo "Infiles must be of type 'bfq'"
    exit

elif ! [ -f $in2 ]
then
    echo "Infile 2 does not exist!"
    exit

elif [ ${in2##*.} != "bfq" ]
then
    echo "Infiles must be of type 'bfq'"
    exit

else ## Get base file and path names:
    f=`basename ${outf%.*} 2>/dev/null` # we don't care about extension
    p=${outf%/*}
fi

if [ -z $f ]
then
    echo "Outfile must include a file basename!"
    exit

elif ! [ -d $p ]
then
    mkdir -p $p
fi

if [ ${refseq##*.} != "bfa" ]
then
    echo "Refseq must be of type 'bfa'"
    exit

elif ! [ -f $refseq ]
then
    echo "Script requires valid reference sequence file!"
    exit

elif ! [ $bins -eq $bins 2> /dev/null ]
then
    echo "Bin size must be an integer!"
    exit

elif [ $bins -le 0 ]
then
    echo "Bin size must be a positive integer!"
    exit
fi

cat<<eof
==============  
SOL PARAMETERS
==============  
   in1 = $in1
   in2 = $in2
     f = $f
     p = $p
  outf = $outf
refseq = $refseq
  bins = $bins
==============  
eof

## Make 'map' file:
trap '' ABRT
maq map -a 320 $p/$f.map $refseq $in1 $in2 || {
    echo Map failed!
    exit
    }

echo "'Map' completed..."

## Make 'pileup' file:
trap '' ABRT
maq pileup $refseq $p/$f.map > $p/$f.pileup || {
    echo Pileup failed!
    exit
    }

echo "'Pileup' completed..."

## Make 'xls' binned file:
trap '' ABRT
gawk -v A=0 -v I=0 -v bin=$bins \
'{if ((I % bin == 0)&&(I != 0)) {print A;A=0};A=A+$4;I++}' \
$p/$f.pileup > $p/$f.xls || {
    echo Binning failed!
    exit
    }

echo "'Binning' completed..."

## Plot read depth
trap '' ABRT
size=`tail -n 1 $p/$f.pileup | awk '{print $2}'`
size=$(( (size+bins)/bins ))
rdplot $p/$f.xls $size $bins "purple" || {
    echo Plotting failed!
    exit
    }

echo "'Plotting' completed..."

echo "Sol run finished."