#!/bin/bash
echo "start run evolution module..."

binDir=$(cd `dirname $0`; pwd)
export PATH=$binDir/../../software/kentUtils:$PATH
transDir=$1
species=$2
transBedFile=$transDir/exon.bed
dataDir=$3/$species/evolution
outDir=$transDir/evolution

if [ ! -s $transBedFile ]; then
  echo "[Exit 1] $transBedFile is empty or not existed."
  echo -e  "evolution\tfailed" >> $transDir/status
  exit 1
fi
promoterBedFile=$transDir/promoter.bed
if [ ! -f $promoterBedFile ]; then
	python $binDir/../utils/get_relative_region.py -e ${transBedFile} -t promoter -o $promoterBedFile
fi

case $species in
	hg19)
		clade_list=("All" "Placental" "Primates")
		score_list=("phyloP" "phastCons" "DAF")
		;;
	hg38)
		clade_list=("All" "Mammals" "Primates")
		score_list=("phyloP" "phastCons" "DAF")
		;;
	mm10)
		clade_list=("All" "Euarchontoglire" "Glire" "Placental")
		score_list=("phyloP" "phastCons")
		;;
	*)
		echo "Unsupported species: $species"
		exit 1
esac

summaryFile=$transDir/evolution/conservation_score.summary
echo -e "region\tscore_type\tclade\tmin\tQ1\tmedian\tmean\tQ3\tmax" > $summaryFile
for region in exon promoter; do
	for score in ${score_list[@]}; do
		if [ $score == "DAF" ]; then
			OutFile=$outDir/$species.$region.$score.wig
			touch $OutFile
			Rscript $binDir/get_DAF.R $transDir/$region.bed $dataDir/DAF/ $OutFile
			if [ ! -s $OutFile ]; then
				echo "No DAF value found in $region"
			fi
			Rscript $binDir/stat_summary.R $OutFile $region $score "NA" >> $summaryFile
		else
			cd $dataDir/$score
			for clade in ${clade_list[@]}; do
				OutFile=$outDir/$species.$region.$score.$clade.wig
				touch $OutFile
				hgWiggle -bedFile=$transDir/$region.bed $species.$score.$clade | grep -P "^[1-9]"  > $OutFile
				if [ ! -s $OutFile ]; then
					echo "No $clade $score value found in $region"
				fi
				Rscript $binDir/stat_summary.R $OutFile $region $score $clade >> $summaryFile
			done
		fi
	done
done

### remove tmp files
rm -f $outDir/*.wig
tag=`cat $transDir/status |grep evolution`
if [[ $tag == ""  ]];then
echo -e "evolution\tfinished" >> $transDir/status
fi

echo "evolution module is done!"