#!/bin/bash

explanation='
DESeq2 analysis
'
inputdef='
input_1::count table:*
input_2::sample list:*
input_3:optional:Trinotate added count table:*
'
optiondef='
opt_i:group 1 name (should be the same one in the sample list):
opt_j:group 2 name (should be the same one in the sample list):
opt_c:cpu threads:8
opt_m:memory limit (GB):64
opt_s:script name:run-DESeq2.R
opt_p:p:0.05
'
runcmd="$0 -c #opt_c# -m #opt_m# -s #opt_s# -i #opt_i# -j #opt_j# -p #opt_p# #input_1# #input_2# #input_3#"

export IM_R="c2997108/centos7:1-trinity_2.8.5-kallisto_0.46.0-blast_2.9.0-trinotate-3.1.1-R_4-kegg_4"

source $(dirname `readlink -f $0 || echo $0`)/common.sh

set -eux
set -o pipefail

script=$opt_s
inputsamplematrix=$input_2
inputmatrix=$input_1
g1=$opt_i
g2=$opt_j
if [ "$g1" = "" ]; then g1=`tail -n+2 $inputsamplematrix|cut -f 2|sort -V|uniq|head -n 1`; fi
if [ "$g2" = "" ]; then g2=`tail -n+2 $inputsamplematrix|cut -f 2|sort -V|uniq|head -n 2|tail -n 1`; fi
p=$opt_p
outputfile=DESeq2.`basename $inputmatrix`
samplex=DESeq2.`basename $inputsamplematrix`.`basename $inputmatrix`.x
sampley=DESeq2.`basename $inputsamplematrix`.`basename $inputmatrix`.y
outputfilexup=result.DESeq2.`basename $inputmatrix`.$g1.up.$g2.down.txt
outputfileyup=result.DESeq2.`basename $inputmatrix`.$g1.down.$g2.up.txt
if [ "$input_3" = "" ]; then
# touch dummy.annotation.txt
# annotationfile=dummy.annotation.txt
 annotationfile=$input_1
else
 annotationfile=${input_3}
fi


cat << EOF > $script
library(DESeq2)
colData <- read.table("$inputsamplematrix", header=T, row.names=1, sep="\\t")
countData=round(read.table("$inputmatrix",sep="\\t",header=T,row.names=1))
mydata=data.frame("group"=as.character(colData[,1]),stringsAsFactors=F)
dds <- DESeqDataSetFromMatrix(countData = countData, colData = mydata, design = ~ group)
dds <- DESeq(dds)
rld <- rlog(dds, blind=FALSE)
resMLE <- results(dds, alpha=$p, contrast=c("group","$g1","$g2"))
a=(resMLE\$padj<$p)
a[is.na(a)]=FALSE
write.table(as.data.frame(resMLE[a,]),file="$outputfile",quote=F,sep="\\t")
EOF

DO_R R --vanilla < $script || true
if [ ! -e "$outputfile" ]; then touch "$outputfile"; fi
tail -n+2 $outputfile|cut -f 1 > $outputfile.id
awk -F'\t' -v g1=$g1 '$2==g1{print $1}' $inputsamplematrix > $samplex
awk -F'\t' -v g2=$g2 '$2==g2{print $1}' $inputsamplematrix > $sampley
bash "$scriptdir"/statistics~scatter_with_variance -x $g1 -y $g2 $inputmatrix $outputfile.id $samplex $sampley
tail -n+2 $outputfile|awk '{if($3>=0){print $1"\tx"}else{print $1"\ty"}}' > $outputfile.xy
awk -F'\t' 'NR==1{print "id\t"$0} NR>1{print $0}' $outputfile| cut -f 1,7|
 awk -F'\t' 'NR==1{print -1"\t"$0} NR>1{print $2"\t"$0}'|sort -k1,1g|cut -f 2- > $outputfile.p
awk -F'\t' 'FILENAME==ARGV[1]{xy[$1]=$2} FILENAME==ARGV[2]{if(FNR==1){header=$0}; res[$1]=$0} FILENAME==ARGV[3]{if(FNR==1){print header"\t"$2}else{if(xy[$1]=="x"){print res[$1]"\t"$2}}}' $outputfile.xy $annotationfile $outputfile.p > $outputfilexup
awk -F'\t' 'FILENAME==ARGV[1]{xy[$1]=$2} FILENAME==ARGV[2]{if(FNR==1){header=$0}; res[$1]=$0} FILENAME==ARGV[3]{if(FNR==1){print header"\t"$2}else{if(xy[$1]=="y"){print res[$1]"\t"$2}}}' $outputfile.xy $annotationfile $outputfile.p > $outputfileyup
i=$outputfilexup; awk -F'\t' -v maxlen=30000 'FILENAME==ARGV[1]{for(i=1;i<=NF;i++){n=int((length($i)-1)/maxlen)+1; if(n>a[i]){a[i]=n}}} FILENAME==ARGV[2]{if(FNR==1){ORS=""; for(i=1;i<=NF;i++){if(a[i]<=1){print $i"\t"}else{for(j=1;j<=a[i];j++){print $i" (split "j")\t"}}}; print "\n"}else{for(i=1;i<=NF;i++){if(a[i]<=1){print $i"\t"}else{for(j=1;j<=a[i];j++){print substr($i,(j-1)*maxlen+1,maxlen)"\t" }}}; print "\n"}}' $i ./$i|sed 's/\t$//' > ${i}.temp
i=$outputfileyup; awk -F'\t' -v maxlen=30000 'FILENAME==ARGV[1]{for(i=1;i<=NF;i++){n=int((length($i)-1)/maxlen)+1; if(n>a[i]){a[i]=n}}} FILENAME==ARGV[2]{if(FNR==1){ORS=""; for(i=1;i<=NF;i++){if(a[i]<=1){print $i"\t"}else{for(j=1;j<=a[i];j++){print $i" (split "j")\t"}}}; print "\n"}else{for(i=1;i<=NF;i++){if(a[i]<=1){print $i"\t"}else{for(j=1;j<=a[i];j++){print substr($i,(j-1)*maxlen+1,maxlen)"\t" }}}; print "\n"}}' $i ./$i|sed 's/\t$//' > ${i}.temp
DO_R java -Xmx1G -jar /usr/local/bin/excel2.jar ${outputfilexup}.temp $outputfilexup.xlsx
DO_R java -Xmx1G -jar /usr/local/bin/excel2.jar ${outputfileyup}.temp $outputfileyup.xlsx

rm -f $outputfile.id $samplex $sampley $outputfile.xy $outputfile.p $outputfilexup.temp $outputfileyup.temp


post_processing


#<option detail>
#</option detail>