#!/usr/bin/perl -w

# usage : perl sm_htseqcount.pl <BAM file> <output htseq-count_nonempty_nostranded file>
# Sarah Maman 2013
# Copyright (C) 2013 INRA
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
#

use strict;
use File::Basename;


my $input_accepted_hits_bam = $ARGV[0];
my $input_gtf               = $ARGV[1];
my $a                       = $ARGV[2];
my $t                       = $ARGV[3];
my $i                       = $ARGV[4];
my $input_s                 = $ARGV[5];
my $input_mode              = $ARGV[6];
my $bam_sorted_option       = $ARGV[7];
my $output_file             = $ARGV[8];

my $SAMTOOLS = '/usr/local/bioinfo/bin/samtools';
my $HTSEQ = '/usr/local/bioinfo/src/python/current/bin/htseq-count';
my $cmd = '';
my $cmd2 = '';
my $path = dirname($input_accepted_hits_bam);
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
my $ALL = $mday."-".($mon+1)."_".$hour."h".$min."mn".$sec."sec_htseq";
my ($nb) = ($output_file=~/galaxy_dataset_(\d+)\.\S+/);
#my $sam = $path."/tmp".$nb.".sam";
my $aoption ="",
my $toption ="";
my $ioption ="";
my $input_modeoption ="";

if ($a eq ""){$a = "0";}
if ($t eq "cds"){$toption = "-t CDS";}
if ($t eq "exon"){$toption = "-t exon";}
if ($t eq "5utr"){$toption = "-t 5UTR";}
if ($t eq "3utr"){$toption = "-t 3UTR";}


if ($i){$ioption = "-i $i";}else {$ioption = "";}

if ($input_mode eq "U"){$input_modeoption = "-m union ";}
if ($input_mode eq "IS"){$input_modeoption = "-m intersection-strict ";}
if ($input_mode eq "INE"){$input_modeoption = "-m intersection-nonempty ";}


$bam_sorted_option=~s/_//;


#commande : perl sm_htseqcount.pl tests/accepted_hits.bam tests/Danio_rerio_chr22.Zv9.62.gtf
#samtools view accepted_hits.bam | htseq-count --stranded=no -m intersection-nonempty - /work/.../Danio_rerio_chr22.Zv9.62.gtf -q > ERR022486.htseq-count_nonempty_nostranded &
#$cmd = "($SAMTOOLS view $input_accepted_hits_bam | $HTSEQ --stranded=no -m intersection-nonempty - $input_gtf -q > ${input_accepted_hits_bam}.htseq-count_nonempty_nostranded & ) >& ./htseqcount.log 2>&1";
#$cmd = "($SAMTOOLS view -@ 4 $input_accepted_hits_bam > $sam) >& /galaxydata/logs/samtoolsview_$nb.log 2>&1";
#Info pour les biologistes
#print STDOUT "Samtools : \n\n $cmd \n\n ";
#system $cmd;

#my $result_file = "${input_accepted_hits_bam}.htseq-count_nonempty_nostranded";
my $result_file = "/galaxydata/database/files/workspace/$nb.htseq-count_nonempty_nostranded";
$cmd2 = "($HTSEQ $input_s -a $a $toption $ioption $input_modeoption $bam_sorted_option $input_accepted_hits_bam $input_gtf  -f bam -q > $result_file ) >& ./htseqcount.log 2>&1";
#Info pour les biologistes
print STDOUT "HtSeq : \n\n $HTSEQ $input_s -a $a $toption $ioption $input_modeoption $bam_sorted_option $input_accepted_hits_bam $input_gtf  -f bam -q > $result_file \n\n ";
system $cmd2;

#`rm $sam;`;



if (! -e "$result_file")
{
print STDERR "htseq-count_nonempty_nostranded FILE NOT FOUND\n";
}
else
{
`cp -a $result_file $output_file`;
}
