#!/bin/tcsh -f
#(ie run the tshell on this but don't read the .cshrc or .tcshrc)

echo version = 1.07 of run.curva 2003 Sep 2
# 2003 Sep  2, 1.07: upgrade documentation
# 2003 Sep  2, 1.06: needfile curvafeap
# 2003 Aug 31, 1.05: report more data to the curva outputfile
# 2003 Aug 29, 1.04: cleanup
# 2003 Aug 26, 1.03: generalize
# 2003 Aug 26, 1.02: renamed from mkcurva to run.curva, cleanup
# 2003 Aug 25, 1.01: functional
# 2003 Aug 25, 1.00: origin from run.mfold

# Compute curvature for each sequence in a Delila book.

# Input:  a delila book containing one or more DNA sequences

# Output: curvaxyin, output from the curvafea program.

# REFERENCE:

# @article{Shpigelman.Bolshoy1993,
# author = "E. S. Shpigelman
#  and E. N. Trifonov
#  and A. Bolshoy",
# title = "{CURVATURE: software for the analysis of curved DNA}",
# journal = "Comput Appl Biosci",
# volume = "9",
# pages = "435--440",
# year = "1993"}

# PROGRAM SOURCE:

# curva source (as of 2003 Aug 21):
# http://ext.haifa.ac.il/~leokol/curvature/index.html

# run.curva uses these programs:
# http://www.lecb.ncifcrf.gov/~toms/delila/curvafea.html
# http://www.lecb.ncifcrf.gov/~toms/delila/mfoldseq.html

# run.curva needs these data files:
# curvafeap

# OVERVIEW:

# The curva program can only handle one sequence at a time.  Therefore
# this script uses mfoldseq to extract one sequence at a time from the
# current delila book.  This is passed to curva which computes the
# curvature.  The results are collected in the curvaout file.  Then
# the curvafea program turns this data along with the original book
# into a set of marks for the lister program.

# SYNOPSIS:

# Take a book, generate DNA curvatures for each piece and put the
# results into file curvaout for use by the lister program.

# The input file is a delila book.  Each sequence in the book is used
# to find the curvature in turn, and the results are gathered in
# curvemarks

# After running this script, the results are in file mfoldfeatures. 
# The user can now:
#   cat mfoldfeatures >> features
# to attach the fold features to the end of the features file.  Then
# one can run the lister program to see the features.

# ** Preliminary setup *********************************************************

if !(-f angles.dat) then
  echo lk angles.dat $DELILA_DELILA/angles.dat
  lk angles.dat $DELILA_DELILA/angles.dat
  if !(-f angles.dat) then
     echo cannot locate $DELILA_DELILA/angles.dat
     echo This is provided with the curva program.
     echo Place a copy into the directory $DELILA_DELILA.
     exit
  endif
endif

set sourceseq = curvesequ

# ******************************************************************************

if !(-f book) then
   echo "usage: run.curva requires a delila book"
else
   echo 1 > columnp
   set sequences = `count | head -2 | tail -1 | column`
   echo $sequences sequences in the book to fold

   echo -n "" > curvaout

   set number = 1
   while ($number <= $sequences)
      echo "*********** sequence $number ***********"
      echo \
       "1.05  version of mfoldseq that this parameter file is designed for." \
       > mfoldseqp
      echo $number " which sequence to extract."  >> mfoldseqp

#echo ------------------ mfoldseqp:
#cat mfoldseqp
#echo ------------------

      mfoldseq

#echo $source -------
#cat $source
#echo -------

      # convert to the format needed by curva
      cat mfoldsequ |\
      grep -v ';' |\
      grep -v '#' |\
      sed -e 's/1//g' |\
      tr "[:upper:]" "[:lower:]" |\
      cat > $sourceseq

      cat $sourceseq

      echo '------------- curva:'
      curva -f $sourceseq
      # % curva
      # Infile? 
      # Correct syntax is
      #  curva -f  <infile> [-o <output file][ -p ] [ -a <starting position>]
      #  -b <length> [-m <path-mode>]
      # Defaults: p_opt=map; a_opt=1; b_opt=max; m_opt=3; o_opt=~
      echo '-------------'

      # clean up the files left behind by curva
      rm $sourceseq{ou2} work.0

      if !(-f $sourceseq.map) then
         printf "*****************************************************\n"
         printf "* The %10s file was not generated by mfold *\n" $sourceseq.map
         printf "*****************************************************\n"
         echo "* Something has gone wrong\!                 *"
         echo "* PROGRAM run.mfold ABORTED                  *"
         printf "*****************************************************\n"

         set debugging = 'not debugging'
         set debugging = 'debugging'

         if ('debugging' == $debugging) then
         endif

         exit
      endif

      echo "" >> curvaout
      echo "*********** sequence $number ***********" >> curvaout
      cat $sourceseq.map >> curvaout

      @ number++
   end

   echo use the curvaout to make the curvafeamarks and curvaxyin
   if !(-f curvafeap) then
      needfile curvafeap
   endif
   curvafea
endif