#!/usr/bin/perl -w

=pod

=head1  NAME

final_count.pl
 
=head1  SYNOPSIS

final_count.pl [options] 
-arg1 <Path to genome bowtie2bd fasta file (/bank/bowtie2db/genome.fa)> 
-arg2 <Path to your merged.gtf file> 
-arg3 <Path to your project> 
-arg4 <Path to your results>

=head1  OPTIONS

=over 8

=item B<-help>

-arg0    Path to genome bowtie2bd fasta file (/bank/bowtie2db/genome.fa)
-arg1    Path to your merged.gtf file
-arg2    Path to your project
-arg3    <Path to your results>
-arg4    Liste des fichiers issus de sigcufflinks ls *_L*/raw_transcripts.tsv

Example : ./final_count_file.pl /bank/bowtie2db/genome.fa All/merged.gtf /work/sigenae/Project_MONOPOLY.222 DIRNAMERESULTS/ ls *_L*/raw_transcripts.tsv

=item B<-man>

Prints the manual page and exits.

=back

=head1 DESCRIPTION

Prepare final count file for RNAseq analysis (header + count of raw_transcripts files).
Steps :
Grouping all output files in e RESULTS directory.
Generation of genes and transcripts global count files.
Production of transcripts fasta from merged.gtf of Cuffmerge.
Loading and processing data in R.

=head1 DATE

19/03/2013

=head1 AUTHORS

Sarah Maman avec l'aide de Christophe Klopp et Patrice Dehais.

=head1 CONTACT

Questions can be posted to the sigenae mailing list:
sigenaesupport@listes.inra.fr

=cut

use strict;
use warnings;
use List::Util qw(sum);
use Carp;
use Statistics::R;
use Pod::Usage;


my $cmd = '';


print "\n __________________________________________________________________________________________ \n";
print "\n - 9 - R statistics : dendogramme and graphic correlation. \n";

#./Rtest.pl /usr/local/bioinfo/src/galaxy/galaxy-dist/tools/sm_RNAseq testCOUNT3


#déclaration de l'espace de travail et du chemin de l'exécutable R :
#use Carp;
#use Statistics::R;
#!/usr/bin/perl
# Espace de travail, répertoire SPÉCIFIQUE pour l'établissement du pont entre R et Perl
my $home = $ARGV[0]; #dirname input SOIT /usr/local/bioinfo/src/galaxy/galaxy-dist/tools/sm_RNAseq
my $ALL =  $ARGV[1]; #rep de travail testCOUNT/ SOIT ?
my $log_dir = $home."/".$ALL;
# Rbin path
my $r_bin   = '/usr/local/bioinfo/src/R/current/bin/R';
#déclaration de l'objet $R et ouverture du pont :
# Déclaration de l'objet
my $R = Statistics::R->new(
    "r_bin"   => $r_bin,
    "log_dir" => $log_dir,
) or die "Problem with R : $!\n";


# Ouverture du pont  
$R->startR;


`cut -f2- $ALL/global_count_genes_transcripts.csv > $ALL/global_count_nogene.csv`;

## Chargement d'un fichier de données

$R->send(
    qq`
         MONOPOLY <- read.table("$ALL/global_count_nogene.csv", header=TRUE, row.names=1 ) 
         \n print(MONOPOLY)
    `
);
#my $MONOPOLY = $R->read;
#print $MONOPOLY,"\n\n";

  

# Chargement d'une bibliothèque ggplot2
$R->send(
    q`
        library(ggplot2) 
    `
);

# Chargement d'une bibliothèque reshape
$R->send(
    q`
        library(reshape) 
    `
);

# Création d'une image graphic correlation
my $graphic_correlation = "$log_dir/graphic_correlation.pdf";
$R->send(
    qq`
        pdf("$graphic_correlation")
        y <- cor(log10(MONOPOLY+1))
        print (y)
        y1.m <- melt(y)
        print (y1.m)
        q1 <- ggplot(y1.m, aes(X1, X2, fill = value)) + geom_tile() + scale_fill_gradient(low = \"blue\", high = \"yellow\")
	print (q1)
        q1 + theme(axis.text.x = element_text(angle =-90))
	print (q1 + theme(axis.text.x = element_text(angle =-90)))

        dev.off()
      `
);


# Fermeture du pont
$R->stopR();


print "OK / graphic correlation \n";


print "\n __________________________________________________________________________________________ \n";
print "THE END ...\n";

