#!/bin/csh
#
#This software is Copyright (C) 2005 by Andrew Firth, University of Otago, 
#Dunedin, New Zealand. All rights reserved.
#
#Script to redo the plots part of run_mlrgd script for new window and clipping 
#  parameters.
limit stacksize unlimited
setenv PATH ${PATH}:${HOME}/bin

@ argc = `echo $argv | awk '{print NF}'`
if ($argc != 7) then
  echo "Usage: redo_plots prefix window1 window2 circular clip1 clip2"\
       "newprefix"
  exit
endif

set prefix = $1
set refseq = `head -1 $prefix.seqs`
set window1 = $2 
set window2 = $3
set circular = $4
set clip1 = $5
set clip2 = $6
set newprefix = $7

cp TIDYUP/seq.001.aln .
cp TIDYUP/orfs.001 .
cp TIDYUP/mlrgd.R .

echo -n "Preparing plot data..."

set npairs = `wc $prefix.pairs | awk '{print $1}'`
set np2 = `wc $prefix.pairs | awk '{print $1*0.5}'`

foreach name (nc fc)
  set mname = `echo $name | awk '{printf "%sm\n",$1}'`
  set cname = `echo $name | awk '{printf "%sc\n",$1}'`
  awk '{if ((0.+$4)>(0.5*"'"$npairs"'")) print $1,$2-$3,2.*$9,$10}' \
    $prefix.$name.plot | \
    awk '{if (0.+$4>0) {print $1,$2/$4,$3/$4} else {print $1,0,0}}' | \
    sed 's/nan/0/g' > tghyv.$name.n0
  awk '{if ((0.+$5)>(0.5*"'"$npairs"'")) print $1,$2-$3,2.*$9,$10}' \
    $prefix.$name.plot | \
    awk '{if (0.+$4>0) {print $1,$2/$4,$3/$4} else {print $1,0,0}}' | \
    sed 's/nan/0/g' > tghyv.$name.n1
  awk '{if ((0.+$6)>(0.5*"'"$npairs"'")) print $1,$2-$3,2.*$9,$10}' \
    $prefix.$name.plot | \
    awk '{if (0.+$4>0) {print $1,$2/$4,$3/$4} else {print $1,0,0}}' | \
    sed 's/nan/0/g' > tghyv.$name.n2
  awk '{if ((0.+$7)>(0.5*"'"$npairs"'")) print $1,$2-$3,2.*$9,$10}' \
    $prefix.$name.plot | \
    awk '{if (0.+$4>0) {print $1,$2/$4,$3/$4} else {print $1,0,0}}' | \
    sed 's/nan/0/g' > tghyv.$name.n3
  awk '{if (($4+$5+$6+$7)>0.) print $1,$2-$3,2.*$9,$10}' \
    $prefix.$name.plot | \
    awk '{if ($4>0) {print $1,$2/$4,$3/$4} else {print $1,0,0}}' | \
    sed 's/nan/0/g' > tghyv.$name.all
  awk '{if ((0.+$16)>=(0.5*"'"$npairs"'")) print $1,$13-$14,2.*$9,$15}' \
    $prefix.$name.plot | \
    awk '{if (0.+$4>0) {print $1,$2/$4,$3/$4} else {print $1,0,0}}' | \
    sed 's/nan/0/g' > tghyv.$name.syn
  awk '{print $1,0.25*$10}' $prefix.$name.plot | sed 's/nan/0/g' \
    > tghyv.$name.lambda
  runmean tghyv.$name.n0  $window2 $circular > tghyv.$mname.n0  || exit
  runmean tghyv.$name.n1  $window1 $circular > tghyv.$mname.n1  || exit
  runmean tghyv.$name.n2  $window1 $circular > tghyv.$mname.n2  || exit
  runmean tghyv.$name.n3  $window1 $circular > tghyv.$mname.n3  || exit
  runmean tghyv.$name.all $window2 $circular > tghyv.$mname.all || exit
  runmean tghyv.$name.syn $window1 $circular > tghyv.$mname.syn || exit
  runclipmean tghyv.$name.n0  $window2 $clip1 $clip2 $circular \
    > tghyv.$cname.n0  || exit
  runclipmean tghyv.$name.n1  $window1 $clip1 $clip2 $circular \
    > tghyv.$cname.n1  || exit
  runclipmean tghyv.$name.n2  $window1 $clip1 $clip2 $circular \
    > tghyv.$cname.n2  || exit
  runclipmean tghyv.$name.n3  $window1 $clip1 $clip2 $circular \
    > tghyv.$cname.n3  || exit
  runclipmean tghyv.$name.all $window2 $clip1 $clip2 $circular \
    > tghyv.$cname.all || exit
  runclipmean tghyv.$name.syn $window1 $clip1 $clip2 $circular \
    > tghyv.$cname.syn || exit
end

#Add gaps back in.
foreach i (nc ncm ncc fc fcm fcc)
  foreach k (n0 n1 n2 n3 all syn) 
    mv tghyv.$i.$k tghyv.temp
    addgaps tghyv.temp > tghyv.$i.$k || exit
  end
end

#If file is empty, write '0 0 0' to it so that R doesn't barf.
foreach i (ncm ncc fcm fcc)
  foreach k (n0 n1 n2 n3 all syn) 
    set t = `wc -l tghyv.$i.$k | awk '{print $1}'`
    if (0 == $t) then
      echo 0 0 0 >> tghyv.$i.$k
    endif
  end
end

#Reference sequence orfs.  Change from refseq to alignment coords.
set w = `wc -l orfs.001 | awk '{print $1}'`
if ($w > 0) then
  set norfs = `wc -l orfs.001 | awk '{print $1}'`
  cut --delimiter=" " -f 2- orfs.001 | sed 's/[^0-9]/\n/g' | grep "[0-9]" \
    > temp
  rm -f temp1; touch temp1
  foreach i (`cat temp`)
    echo `ntadjust seq.001.aln $i` >> temp1
  end
  set orfx1 = `awk '{if (NR%2==1) printf "%s,",$0}' temp1 | sed 's/,$//g'`
  set orfx2 = `awk '{if (NR%2==0) printf "%s,",$0}' temp1 | sed 's/,$//g'`
  awk '{print $1/2}' orfs.001 > temp1
  rm -f temp; touch temp
  @ b = 0
  foreach i (`cat temp1`)
    @ b ++
    @ a = 0
    while ($a < $i)
      echo $b $norfs | awk '{print ($1-0.5)/$2}' >> temp
      @ a ++
    end
  end
  set orfy = `awk '{printf "%s,",$0}' temp | sed 's/,$//g'`
else
  set orfx1 = "0"
  set orfx2 = "0"
  set orfy =  "0"
endif

#Known features.  Change from refseq to alignment coords.
if (-r $refseq.features) then
  set w = `wc -l $refseq.features | awk '{print $1}'`
  if ($w > 0) then
    sed 's/[^0-9]/ /g' $refseq.features | awk '{print NF,$0}' > temp2
    set nkfs = `wc -l temp2 | awk '{print $1}'`
    cut --delimiter=" " -f 2- temp2 | sed 's/[^0-9]/\n/g' | grep "[0-9]" > temp
    rm -f temp1; touch temp1
    foreach i (`cat temp`)
      echo `ntadjust seq.001.aln $i` >> temp1
    end
    set kfx1 = `awk '{if (NR%2==1) printf "%s,",$0}' temp1 | sed 's/,$//g'`
    set kfx2 = `awk '{if (NR%2==0) printf "%s,",$0}' temp1 | sed 's/,$//g'`
    awk '{print $1/2}' temp2 > temp1
    rm -f temp; touch temp
    @ b = 0
    foreach i (`cat temp1`)
      @ b ++
      @ a = 0
      while ($a < $i)
        echo $b $nkfs | awk '{print ($1-0.5)/$2}' >> temp
        @ a ++
      end
    end
    set kfy = `awk '{printf "%s,",$0}' temp | sed 's/,$//g'`
  else
    set kfx1 = "0"
    set kfx2 = "0"
    set kfy =  "0"
  endif
else
  set kfx1 = "0"
  set kfx2 = "0"
  set kfy =  "0"
endif

echo "Done"

echo -n "Making plots..."

set w1 = `echo $window1 | awk '{print 1+2*$1}'`
set w2 = `echo $window2 | awk '{print 1+2*$1}'`

#Conservation score plots.
foreach i (nc fc)
  set length = `head -2 $prefix.$i.info | tail -1 | awk '{print $2}'`
  set base1 = `head -2 $prefix.$i.info | tail -1 | awk '{print $3}'`
  set base2 = `head -2 $prefix.$i.info | tail -1 | awk '{print $4}'`
  set nmutsum = `head -3 $prefix.$i.info | tail -1 | awk '{print $2*0.25}'`
  set lambdasum = `head -4 $prefix.$i.info | tail -1 | awk '{printf "%.2f\n",$2*0.25}'`
  set lambda4sum = `tail -1 $prefix.$i.info | awk '{printf "%.2f\n",$2*0.25}'`
  set mname = `echo $i | awk '{printf "%sm\n",$1}'`
  set cname = `echo $i | awk '{printf "%sc\n",$1}'`
  foreach j ($mname $cname)
    cat tghyv.$j.all tghyv.$j.syn tghyv.$j.n[0123] | \
      awk '{printf "%14.6f\n", $2}' > temp1
    set ymin = `minmax temp1 | awk '{print $1}'`
    set ymax = `minmax temp1 | awk '{print $2}'`
    if ("fc" == $i) then
      set aaa = full
    else
      set aaa = non
    endif
    if ($mname == $j) then
      set bbb = ""
    else
      set bbb = "$clip1\/$clip2-clipped"
    endif
    cat mlrgd.R | sed 's/mmm/'$j'/g' | \
      sed 's/yyymax/'$ymax'/g' | sed 's/yyymin/'$ymin'/g' | \
      sed 's/lll/'$length'/g' | \
      sed 's/ddd1/'$base1'/g' | sed 's/ddd2/'$base2'/g' | \
      sed 's/qqq/'$np2'/g' | sed 's/nnn/'$nmutsum'/g' | \
      sed 's/rrr/'$lambdasum'/g' | sed 's/ccc/'$lambda4sum'/g' | \
      sed 's/aaa/'$aaa'-coding model/g' | sed 's/bbb/'$bbb' running mean/g' | \
      sed 's/  running mean/ running mean/' | sed 's/ppp/'$i'/' | \
      sed 's/www1/'$w1'/g' | sed 's/www2/'$w2'/g' | \
      sed 's/orfx1/'$orfx1'/g' | sed 's/orfx2/'$orfx2'/g' | \
      sed 's/orfy/'$orfy'/g' | sed 's/kfy/'$kfy'/g' | \
      sed 's/kfx1/'$kfx1'/g' | sed 's/kfx2/'$kfx2'/g' | \
      sed 's/xxx/'$newprefix'/g' > $newprefix.$j.R
    R --no-save --slave < $newprefix.$j.R
  end
end

rm -f seq.001.aln orfs.001 mlrgd.R temp temp[12] tghyv.*

echo "Done"