#!/bin/csh

################ estano.csh ###################################
# Version 5/7/05

# Description: Annotate ESTs. 
#              Create a .csv file with a table listing the following,
#              for each fasta (.fsa) nucleotide sequence file in
#              the current working directory:
#
#              ESTname - name of the EST
#              GInum - GI number of highest-scoring BLASTX hit
#                 from the GenBank non-redundant protein database
#              Species - Genus and species associated with GInum
#
#

###############################################################


#---------------------------------------------
# Global parameters
#---------------------------------------------

set OUTFILE = $1
set JOBID = $$
set SEP = '^' # character to separate fields in output. ^ by default.
              # It would be preferrable to use a TAB character here,
	      # but the C-shell appears to replace TABs with blank spaces
	      # As a workaround, we temporarily set $SEP to some other
	      # character that is not likely to appear in the output,
	      # When the output has been completed, we then translate
	      # '^' to a TAB.
	      
set DEBUG = 1 # if DEBUG == 1, turn on debugging code, which could
              # be inserted anywhere in the program      

# RM_CMD - command to be used for removing files and directories
# Should work on any Unix system.
if (-e /usr/bin/rm) then
   set RM_CMD = /usr/bin/rm
else
   if (-e /bin/rm) then
      set RM_CMD = /bin/rm
   else
      set RM_CMD = rm
   endif
endif

#---------------------------------------------
# Create a temporary working directory
#---------------------------------------------
set TEMPDIR = estano.tempdir.$JOBID
mkdir $TEMPDIR
cd $TEMPDIR


#################################################
#       MAIN LOOP
#################################################

# The first line of outfile contains column headings.
echo 'NAME'$SEP'GI'$SEP'TAXNAME'$SEP'PROTEINNAME'$SEP'3D_IDS'$SEP'E-VALUE'> outfile.tmp

foreach file (../*.fsa)

    # - - - - - - - - - - - - - - - - - - - -
    # extract name from first line of fasta-format file, 
    # which is in the form '>name description'
    # - - - - - - - - - - - - - - - - - - - -
echo Extracting EST name...
    set ESTNAME = `head -1 $file | cut -f1 -d " " | sed -e 's/^>//'`

    # - - - - - - - - - - - - - - - - - - - -
    # Run blastx against non-redundant protein translations
    # For debugging purposes, we can read blast table files from
    # the parent directory, instead of running blast every time.
    # - - - - - - - - - - - - - - - - - - - -

    if ($DEBUG == 0) then
echo Running blast...
       blastcl3 -i $file -p blastx -d nr -e 0.00001 -o blastx.table -m 8
       # These lines create table files for debugging purposes
       #set basename = $file:r
       #cp blastx.table $basename.blastx.table        
    else
       set basename = $file:r
       cp $basename.blastx.table blastx.table       
    endif

    # If no hits were found, set GINUM to -1. Otherwise, set GINUM to the
    # GI number of the highest scoring hit, which will be on
    # the first line of the blastx.table file. 
    # Also extract EVALUE
    if ( -z blastx.table) then
       set GINUM = '-1'
       set EVALUE = ' '
    else
echo Extracting GI number...
       head -1 blastx.table > firstline
       set GINUM = `cat firstline | cut -f2 -d "|" `
       echo 'GINUM =' $GINUM
    endif
    
    # - - - - - - - - - - - - - - - - - - - -
    # Get the name of the species the protein comes from. This
    # has to be done in two steps. First, get the Taxonomy ID from the GI,
    # and then get the species name from the Taxonomy ID.
    # If no hits were found above, or if there is no Taxonomy ID,
    # TAXNAME = 'n/a' 
    # - - - - - - - - - - - - - - - - - - - -
echo Getting Species name...
    if ($GINUM == '-1') then
       set TAXNAME = 'n/a'
    else
       leash -mn SHoundTaxIDFromGi -mpi $GINUM -of taxid
       set TAXID = `cat taxid`
       if ($TAXID == '-1') then
	  set TAXNAME = 'n/a'
       else
	  leash -mn SHoundGetTaxNameFromTaxID -mpi $TAXID -of taxname
	  set TAXNAME = `cat taxname`
       endif
    endif

# ==============================================Ben Added Codes Start===============================
    
    # (1) added- - - - - - - - - - - - - - - 
    # Get protein name from GenBank description line
    # First, check if GI exists. If it does, go to get protein name; otherwise, set protein name to n/a
    # Second,If GI exists,get the GenBank description line using SHoundGetDefLine based on GI
    # At last, get the protein name form the Genbank description line using 2 'cut' commands
    # - - - - - - - - - - - - - - - - - - - -   
echo Getting Protein Name ...
    if ($GINUM == '-1') then
    	set ProteinName = 'n/a'
    else
    	leash -mn SHoundGetDefLine -mpi $GINUM -of defline
	
    	cat defline|cut -f5 -d"|">item
    	cat item|cut -f1 -d"[">finalitem
    	set ProteinName = `cat finalitem`
    endif


    # (2) added- - - - - - - - - - - - - - - 
    # Get 3D structure IDs of the retrieved protein
    # First, check if GI exists, If it does, go to get 3D ids; otherwise, set 3D ids to n/a
    # Second, If GI exists, retrieve 3D ids using SHound3DNeighboursFromGi based on GI and E-value
    # and export ids into the file 3did
    # At last, retrieve each id from file 3did using foreach loop and append to variable PIDS
    # one-by-on with semicolon-seperated
    # - - - - - - - - - - - - - - - - - - - -   
echo Getting 3D Structure IDs...
    if ($GINUM == '-1') then
    	set PIDS = 'n/a'
    else
    	echo $GINUM start
	set PIDS = ''
    	leash -mn SHound3DNeighboursFromGi -mpi $GINUM -mpi 1e-10 -of 3did
	foreach newline (`cat 3did`)
		set PIDS = $PIDS$newline';'
	end
	echo $PIDS
    endif
    
    
    # (3) added- - - - - - - - - - - - - - - 
    # Get E-value of the highest blastx hit from the first line of blastx.table
    # First, check if GI exists. If it does, go to get E-value; otherwise, set E-value to n/a
    # If GI exists, get E-value from the first line of blastx.table using one 'cut' command
    # - - - - - - - - - - - - - - - - - - - - 
echo Getting E-Value ...
    if ($GINUM == '-1') then 
        set EVALUE = 'n/a'
    else
    	set EVALUE = `cat firstline | cut -f11`  
    endif
    
# ===========================================Ben Added Codes Stop===================================


    # - - - - - - - - - - - - - - - - - - - -
    # Append the  output line to outfile.Each output line is a set of tab-separated fields:
    # name<TAB>gi<TAB>taxname
    # - - - - - - - - - - - - - - - - - - - -
    
echo Printing output line...
    echo $ESTNAME$SEP$GINUM$SEP$TAXNAME$SEP$ProteinName$SEP$PIDS$SEP$EVALUE >> outfile.tmp

end

# Translate the seperator characters into tab characters
tr $SEP ',' <outfile.tmp > outfile.csv

#################################################


#---------------------------------------------
# Move the output file to the parent directory,
# and delete the temporary directory.
#---------------------------------------------
mv outfile.csv ../$OUTFILE
cd ..
$RM_CMD -r $TEMPDIR
