#! /usr/local/ARCH/alpha5/bin/bash
#Lancer amiga a partir d'une soumission web en passant par LSF

###########
#Fonctions#
###########

#Usage du programme
function help {
        echo "Usage: $pg  <fna file> <directory of your matrix(OPTIONNAL)>" 
        exit 1
}

#Erreur
function Error {
        touch $dirAmiga"/_ERROR_"
	echo $2 >  $dirAmiga"/_ERROR_"
        exit 2
}

#################################################
#Verification des  aguements passes au programme#
#################################################

#Nom du programme
pg=`basename $0`

#Tester le nombre d'arguments passe au programme
if [ $# -eq 0 ]
	then help $pg
fi

if [ $# -ge 3 ] 
        then help $pg 
fi

if [ ! -f $1 ]
        then echo "$1" " not found, exit " ; exit 2
fi


#nom du repertoire des resultats amiga pour l'execution en cours
bank=$1

dirAmiga=`ls $bank| sed 's/[a-zA-Z0-9.]*$//'`

#Home amiga
AMIGAHOME="/env/cns/home/jlesaux/bbb/"

echo "Your Annotation file : " $bank

#Tester si le repertoire existe des resultats amiga
if [ ! -d $dirAmiga ] 
	then help $pg 
fi

#Tester si le repertoire amiga existe
if [ ! -d $AMIGAHOME ] 
	 then Error $dirAmiga 1
fi


##########################
#Variables d'environnment#
##########################

source /etc/profile

PATH=$AMIGAHOME:$PATH

#################
#Variable amiga #
#################

cat $bank | gawk -f $AMIGAHOME/micheckFilter.awk | gawk '{if ($0 ~/^ID|^AC|^OS|^FT/) { exit 1; } else { if ($0 ~ /^LOCUS|^DEFINITION|^REFERENCE|^     gene/) { exit 0; } } } '

if [ $? -eq 0 ] 
    then banktype=`echo "gbk"`
    else banktype=`echo "embl"`
fi 

if [ "$banktype" = "gbk" ]
	then labelseq=$bank.gbk
	else labelseq=$bank.embl
fi


seq="$labelseq.fna"

banktab="$labelseq.tab"

dirSpecies=$dirAmiga

amigares="$seq.cds.lst"

micheckres="$labelseq.res"

StringStart="/ATG/GTG/TTG/CTG/"

StringStop="/TAA/TAG/TGA/"

if [ $# -eq 2 ]
	then if [ -f $2 ];
		then dirSpecies=` ls $2 | sed 's/[a-zA-Z0-9.]*$//'`
		speciesFile=$2 
	     fi
	else rm $seq.* 2> /dev/null && rm -f *.bin
fi


########################################
#se deplacer dans le reperoire resultat#
########################################

cd $dirAmiga


###############################################
#Parser le fichier de banque $bank            #
#creation fu chier fna et tab $seq et $banktab#
###############################################

echo "run MICHECK parse  $labelseq $banktype"

if [ -f "$seq" ]
    then rm $seq
fi

# recuperation de la sequence et test des features 
cat $bank | gawk -f $AMIGAHOME/micheckFilter.awk | gawk -v filename=$seq -f $AMIGAHOME/fastaParser.awk

Lseq=`tail +2 $seq | tr -d "\012" | wc -c | sed 's/ //g'`

echo "$Lseq"

if [ "$Lseq" -eq 0 ]
	then exit 2
fi


echo "Get features and write it into  $banktab"

if [ "$banktype" = "gbk" ]
	then cat $bank | gawk -f $AMIGAHOME/micheckFilter.awk | gawk -v Lseq=$Lseq -f  $AMIGAHOME/genbankParser.awk | sort -n i| uniq | sed 's/^[0-9]* //' > $banktab
	else cat $bank | gawk -f $AMIGAHOME/micheckFilter.awk | gawk -v Lseq=$Lseq -f  $AMIGAHOME/emblParser.awk | sort -n | uniq | sed 's/^[0-9]* //' > $banktab
fi

if [ `wc -l $banktab | awk '{ print $1 }' ` -eq 0 ] 
	then exit 2
fi 

echo "Test number of features "

FeatureMean=`cat $banktab | awk -v Lseq=$Lseq 'BEGIN{nb=0;} {if ($0 ~ /^$/ ){nb++;}} END{if (nb/Lseq < 0.0005) {print "1"; } else { print "0";}}'`

#echo "Density of features : $FeatureMean"

if [ "$FeatureMean" -eq  0 ]
	then echo "Tow few CDS" #2> $dirAmiga"/_ERROR_"
		exit 2
fi
##############################################################
#Formatage de la sequence et calcul de sa longueur et du GC %#
##############################################################

translateCode=`cat translTable.txt 2> /dev/null`



cp $seq tmp.seq
cat tmp.seq | gawk 'BEGIN {getline ;print $O}
{
gsub("[^A-Z]","")
gsub("[^ATGC]","N")
if ( $0 != "") print $O
}' > $seq

rm tmp.seq

tail +2 $seq | tr -d "\012" | wc -c > $labelseq".length"

#longueur de la sequence
lSEQ=` cat  $labelseq.length | sed 's/ //g'`

#longueur de la sequence sans N
lSEQnoN=`tail +2 $seq | tr -d "\012" | tr -d "N" | wc -c | sed 's/ //g'`

echo "Length of the sequence = $lSEQ "
echo "Length of the  sequence without N = $lSEQnoN "

#% GC
GC=`tail +2 $seq | awk '$0 !~ /^$/ {
    long+=length($0)
    nbG+=gsub("G","X")
    nbC+=gsub("C","X")
}
END {
    printf( "%s", ((nbG+nbC)/long)*100 )
}' `

echo "GC pct = $GC"

GCcontent=`echo $GC | awk '{  print sprintf("%2.0d",$0) }'`


if [ "$GCcontent" -lt 45 ]
	then source $AMIGAHOME/Low.seuils
			echo "Low GC"
	else 
		if [ $GCcontent -ge 45 -a $GCcontent -lt 55 ]
			then source $AMIGAHOME/Medium.seuils
					echo "Medium GC"
			else 
				if [ $GCcontent -ge 55 ]
					then source $AMIGAHOME/High.seuils
						echo "High GC"
				fi
		fi
fi

##################################################################
#Si l'utilisateur a choisi de generer automatiquement une matrice#
#    AutoMatrix == "Automatic"                                   #
##################################################################


if [ `ls $dirAmiga/*.bin 2> /dev/null | wc -l ` -eq 0  ] 
	then AutoMatrix="Automatic"
	else AutoMatrix="Generated"
fi

if [ "$AutoMatrix" = "Automatic" ]; then

    echo "Matrix is being generated..."


    echo "$seq.codant is being generated ..."

    #generer du codant
    cat $banktab | awk -v FS="$" ' { if ( $1 !~ /^FRAG/  ) { print $1 " " $9 " " $2 " " $4  } }'  | $AMIGAHOME/Go2Nucleome $seq > $seq".codant" 

    echo "Test of the number of Ac Nuc in coding sequence"

    res=`grep '^[A-Z]' $seq.codant |tr -d '\n'|wc -c`

    echo "CDS total length :$res"
  
    echo "Determining order for CoDing Sequences..."

    if [ $res -lt 25600 ]
       then 
           echo "Not enough CoDing Sequences to build up matrix !!!" 
           echo "Bye....." 
       exit
    fi

    if [ $res -ge 25600 -a $res -lt 102400 ]
       then
       ordreMat=4
    fi

    if [ $res -ge 102400 -a $res -lt 409600 ]
       then
       ordreMat=5
    fi
 
    if [ $res -gt 409600 ]
       then
       ordreMat=6
    fi
    echo "k for CoDing Sequences : $ordreMat" 1>&2

    #Verif ordre de la matrice si 0 erreur
    if [ "$ordreMat" = "0" ]
	then  Error $dirAmiga 4
    fi

    
    #Construire la matrice
    prokov_learn -k $ordreMat -K 1 -I -o $seq".mat1.bin" $seq".codant"
    
    #Pattern du debut du nom de fichier de la matrice avec le nom de repertoire
    begMat=$seq".mat"
    
    #Pattern de fin du nom de fichier de la matrice
    endMat=".bin"
    
fi
 

##############################################################
#Si l'utilisateur a choisi d'utiliser une matrice specifique #
#d'un organisme $species et $AutoMatrix == "Generated"       #
##############################################################
if [ "$AutoMatrix" = "Generated" ]; then
    #Repertoire ou se trouve la matrice
    if [ -f $speciesFile ]; then
    	#tester si il existe au moins une matrice
    	begMat=`ls $speciesFile | sed 's/[0-9]*.bin//'`
    	#Pattern du debut du nom de fichier des matrices avec le nom de repertoire
   	endMat=".bin"
    	#Pattern de fin du nom de fichier des matrices
    	endMat=".bin"
    	else
	    Error $dirAmiga 2 
	fi
else
	begMat=$seq".mat"
	endMat=".bin"	
fi

##############
#Lancer amiga#
##############
#gestion du cas des mycoplasme


if [ "$translateCode" = "4" ]; then
	StringStop="/TAA/TAG/"
fi

#diviser les parametres par 100
Sure_ss_I=`echo $Sure_ss_I | awk '{print $1 / 100}' `
Sure_os_I=`echo $Sure_os_I | awk '{print $1 / 100}' `
Sure_prob_O=`echo $Sure_prob_O | awk '{print $1 / 100}' `
Prob_glob_IO=`echo $Prob_glob_IO | awk '{print $1 / 100}' `


echo "run Amigene"

$AMIGAHOME/amiga_run.sh $seq $begMat $endMat $StringStart $StringStop $Prob_Pc $Sure_Pc $Prob_LMin $Sure_LMin $Sure_ss_I $Sure_os_I  $Sure_prob_O  $Prob_glob_IO $AMIGAHOME  1> /dev/null


###############################################
#Si l'utilisateur a choisi d'executer PROFED  #
#$profed == "CHECKED"                        #
###############################################

##############################################
# reformatage des fichiers pour lancer profed#
# sur le genome annote                       #
##############################################

echo "Changing format"

cat $banktab | awk '{ gsub("[$]"," ",$0); print ">" $1 " " $2 " " $4 " " $10 }'  > $seq.pks.lst

$AMIGAHOME/bkprkS_micheck.sh $seq $seq.curve $seq.pks.lst $AMIGAHOME 2> /dev/null

cat $seq.bk.comb.lstort | awk '{if ($7 > 0) { strand="D" ; } else { strand="R" ; }  print $1 " " $2 " " $3 " " $4 " " ($4 - $2) " " ($4 - $2) " " $5 " " $6 " " $7 " " strand " " $8  }' > $seq.cds.user.lst

i=0

cat  $seq.cds.user.lst | awk '{ print $1 " " $2 " " $3 " " $4 " " $7 " " $8 " " $9 " " $11}' > $seq.cds.user.comb.lst


# on prend uniquement la matrice 1 pour eviter des artefacts avec le non codant contenu dans la matrice 3 ou 4


########################################
#reformatage des fichiers pour obtenir #
#les pcods pour les cds USER           #
########################################


cat $seq.cds.user.lst | awk '{ sub(/^>/,"",$0) ; gsub(" ","$",$0) ; print $0 ; } ' > $seq.user.pcod.lst

cat $banktab | awk 'BEGIN { FS="$" }{ print $1 "$"  $2 "$" $4 "$" $12 }' > $seq.user.product.lst

# reassociement des objets product avec leur information 
perl  -e  'open(P,$ARGV[1]) or die "error $ARGV[1]"; while (<P>) { if (m/(\w+)\$(\d+)\$(\d+)\$(.*)/) {$name=$1; $begin=$2; $end=$3; $product=$4;  open(L,$ARGV[0]) or die "error $ARGV[0]" ; while (<L>) { if (m/(\w+)\$(\d+)\$\d+\$(\d+)/) { if ($begin==$2 && $end==$3 && $name eq $1) {chomp $_; print  $_ ."\$".$product ."\n" ;}} } close L;  } }close P ; '  $seq".user.pcod.lst"  $seq".user.product.lst"   > $labelseq".pcod.tab"




#########################################
# lancement de compannot et des parsers #
#########################################

echo "Comparison of the annotations"

$AMIGAHOME/compannot.pl $amigares $labelseq.pcod.tab 2> /dev/null | sort -n | sed 's/^[0-9]* //' > $micheckres 

if [ $? -ne 0 ]
    then exit 2
fi

##########################################################################
# ajout des uniques amigene dans la liste qui va etre utilise par profed #
##########################################################################

cat $micheckres | grep "AMIGene" | awk -v FS="$" '{print $1 " " $2 " " $3 " " $4 " " $7 " " $8 " " $9 " " $11 }' > $seq.cds.user.comb.addon.lst

cat $micheckres | grep "AMIGene" | awk -v FS="$" '{print  " " $2 " " $3" " $4 " " $5 " " $6 " " $7 " " $8 " " $9 " " $10 " " $11 }' > $seq.cds.user.addon.lst

cat $seq.cds.user.comb.addon.lst  $seq.cds.user.comb.lst | awk '{print $2 " " $0}' | sort -n | sed 's/^[0-9]* //' > $seq.cds.user.comb.tot.lst

cat $seq.cds.user.addon.lst $seq.cds.user.lst | awk '{print $2 " " $0}' | sort -n | sed 's/^[0-9]* //' > $seq.cds.user.tot.lst

echo "run profed"
$AMIGAHOME/profed -r -c $seq.curve1 -f  $seq.cds.user.tot.lst -l $seq.cds.user.comb.tot.lst 2> /dev/null | sort -n > $seq.fed.1 


cat $seq.fed.* |  $AMIGAHOME/profedParser.pl |  awk '{ if ($2~/STRONG/) { print "1 " $0 ;} else { if ($2~/MEDIUM/) { print "2 " $0 ; } else { print "3 " $0; } } } ' |sort -n -k1 -k2 | sed 's/^[0-9] //' | awk '{gsub(">",""); print;}' > $seq.fed



##################################################
# creation du fichier .fna pour capturer les cds #
##################################################
cat $micheckres | awk -v FS="$" ' { print $1 " " $9 " " $2 " " $4  } '  | Go2Nucleome $seq > $seq".cds.fna" 

cat $micheckres | awk -v FS="$" ' {if ( $NF == "UNIQUE_AMIGene" ) {  print $1 " " $9 " " $2 " " $4 }  } '  | Go2Nucleome $seq > $seq".UNIQUE_AMIGene.fna" 

if [ -e $AMIGAHOME/pwg_translate ]; then
    #Si code micloplasme
    if [ "$translateCode" = "4" ]; then
    	code=3
    #sinon code standard
    else
    	code=0
    fi
    
    file_cds_fna=$seq".UNIQUE_AMIGene.fna"
    file_translate=$seq".UNIQUE_AMIGene.faa"
    
    pwg_translate -c $code < $file_cds_fna > $file_translate  #2> $dirAmiga"/_ERROR_"
fi







##################################################
#Si l'utilisateur a choisi d'executer translate  #
#$translate == ""CHECKED"                        #
##################################################
if [  -e $AMIGAHOME/pwg_translate ]; then
    #Si code micloplasme
    if [  "$translateCode" = "4" ]; then
    	code=3
    #sinon code standard
    else
    	code=0
    fi
    
    file_cds_fna=$seq".cds.fna"
    file_translate=$seq".cds.faa"
    
    pwg_translate -c $code < $file_cds_fna > $file_translate  #2> $dirAmiga"/_ERROR_"
fi


#################################################
#creation du fichier de statistiques sur les cds#
#################################################

echo "create stat file "

nbuniqueusercds=`grep -c "UNIQUE_USER" $micheckres`
nbuniqueamigenecds=`grep -c "UNIQUE_AMIGene" $micheckres`
nbcommoncds=`grep -c "COMMON" $micheckres`
nbbankcds=`wc -l $banktab | awk '{ gsub (" "," ",$0) ; print $1 } '`
nbamigenecds=`wc -l $amigares | awk '{ gsub (" "," ",$0) ; print $1 } '`

echo "$nbuniqueusercds $nbuniqueamigenecds $nbcommoncds $nbbankcds $nbamigenecds" > $seq.stat.tmp


cat $seq.stat.tmp | awk '{cdstot = $1 + $2 + $3 ; print $4 " " $5 " " $3 "(" sprintf("%.2f",( $3 / cdstot )*100) "%) " $1 "(" sprintf("%.2f",( $1 / cdstot )*100) "%) " $2 "(" sprintf("%.2f",( $2 / cdstot )*100) "%)" ; }' > $seq.stat

echo
echo
echo
echo "********************************************************************************"
echo "*-------------------------------Summary-Result---------------------------------*"
echo
echo "CDS_bank"$'\t'"CDS_amigene"$'\t'"Common%"$'\t'$'\t'"Unique_Bank%"$'\t'"Unique_amigene%"
cat $seq.stat | awk '{print $1 "\t\t" $2 "\t\t" $3 "\t" $4 "\t" $5 }'
echo "********************************************************************************"

rm -f $seq.stat.tmp
rm -f $seq.lst*
rm -f $seq.cds.user.*
rm -f start.pat*
rm -f $seq.user.*

#########################################################################
# creation des fichiers contenant les uniques amiga et les uniques user #
# trier par ordre decroissant de longueur de cds                        #
#########################################################################


cat $micheckres | awk ' BEGIN { FS="$" ; } { if ( $13 ~ /UNIQUE_USER/ ) { print $5 " " $0 ;} }' | sort -rn | sed 's/^[0-9]* //' > $micheckres.uniqueuser 

cat $micheckres | awk ' BEGIN { FS="$" ; } { if ( $13 ~ /UNIQUE_AMIGene/ ) { print $5 " " $0 "$" int( (($4 - $2 +1)/3)-1) ;} }' | sort -rn | sed 's/^[0-9]* //' > $micheckres.uniqueamigene

echo
echo "--------------------------------------------------------------------------------------"
echo
echo "You can find your results in :"
echo "- "$micheckres.uniqueuser" (list of unique CDS from file $bank)"
echo "- "$micheckres.uniqueamigene" (list of unique Amigene CDS)"
echo "- "$seq".UNIQUE_AMIGene.fna" "(unique Amigene CDS DNA sequences)"
echo "- "$seq".UNIQUE_AMIGene.faa" "(unique Amigene CDS protein sequences)"
echo "- "$seq.fed" (Profed results)"
echo
echo "--------------------------------------------------------------------------------------"
echo





#################################
#Fin Creation d'un fichier _END_#
#################################

touch $dirAmiga"/_END_"


