#!/usr/bin/perl -w

=pod

=head1  NAME

final_count.pl
 
=head1  SYNOPSIS

final.pl [options] 
-arg0 <Path to genome bowtie2bd fasta file (/bank/bowtie2db/genome.fa)> 
-arg1 <Path to your merged.gtf file> 
-arg2 <Result directory> 
-arg3 <Your raw transcripts tsv files> 

=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 <Result directory> 
-arg3 <Your raw transcripts tsv files>


Example : ./final.pl /bank/bowtie2db/genome.fa All/merged.gtf RESULTS tests/1.tsv tests/2.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;
use File::Basename;


#######################################################################################################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############            DECLARATION DES VARIABLES                            ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
#######################################################################################################
#Ligne de commande pour tester l'outil :
#./final_count_file.pl /bank/bowtie2db/ensembl_sus_scrofa_genome All/merged.gtf

#my $list = `\\ls *_L*/raw_transcripts.tsv`;
#Suppression (faux pour Galaxy)
#my $list = $ARGV[3];
#Ajout du nombre de fichiers tsv entrants
my $nbfile=$#ARGV+1;


my $listsom ="";
my $list ="";
my $listnom ="";
my $somname ="";
#separer liste des tsv et noms des tsv
for (my $l=3; $l<$nbfile; $l += 2)
{
	my $listfile=$ARGV[$l];
	$listsom .= $listfile."\t";
	$list = $listsom;
	#print "liste des fichiers : $listfile \n $listsom\t";
	my $listname=$ARGV[$l+1];
	#print "liste des noms : $listname \n\n";
	$listnom .= $listname."\t";
	$somname = $listnom;
}

print STDERR "Liste des fichiers tsv entrants : $list \n\n";
print STDERR "Liste des noms des fichiers tsv entrants : $somname \n\n";


my $cmd = '';
my $i = '';
my $j='';
my $file_list = '' ;
my $c_header = 'transcript	gene	';
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time); 
#my $ALL =  "RESULTATS".$mday."-".($mon+1)."_".$hour."h".$min."mn".$sec."sec"; 
my $ALL  = $ARGV[2];

#######################################################################################################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############         REPERTOIRE RESULSTATS                                   ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
#######################################################################################################


#pour regrouper tous les fichiers sortants dans un repertoire RESULTATS
$cmd = "mkdir $ALL; chmod -R 777  $ALL;";
print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 1 / 6 - Results directory creation : \n";
print STDERR "OK / $ALL \n";
system $cmd;


#######################################################################################################
##############                                                                 ########################
##############                                                                 ########################
##############                            p                                    ########################
##############                                                                 ########################
##############         FICHIER DE COMPTAGE DES GENES ET DES TRANSCRITS         ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
#######################################################################################################

#GENERATION DU FICHIER DES COMPTAGE GENE ET TRANSCRITS
print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 2 / 6 - Generate count file for genes and transcripts : \n";
print STDERR "OK / Sort of count file.\n";

#TODO : Mis en commentaires pour les tests
#foreach $i ( split /\s+/,$list )
#    {
#    $cmd = "sort -S 3G -k2,2 $i | cut -f3 > $i.cut3";
#    system $cmd;
#    $c_header .= $i."\t";
#    $file_list .= $i.".cut3 ";
#    $j = $i;
#    }

my $k ='';
#Ajout pour separer la liste de nom de fichiers de la liste des noms des fichiers
foreach $i ( split /\s+/,$list )
    {
    $cmd = "sort -S 3G -k2,2 $i | cut -f3 > $i.cut3";
    system $cmd;
    #$c_header .= $i."\t";
    $file_list .= $i.".cut3 ";
    $k = $i;
    }

foreach $i ( split /\s+/,$somname )
    {
    $cmd = "sort -S 3G -k2,2 $i | cut -f3 > $i.cut3";
    system $cmd;
    $c_header .= $i."\t";
   # $file_list .= $i.".cut3 ";
    $j = $i;
    }








print STDERR "OK / counts.\n";
$c_header .= "\n";
#print $j."\n";
$cmd = "sort -S 3G -k2 $k | cut -f1,2 > $ALL/header";
system $cmd;
print STDERR "OK / header : $cmd \n";


open(FH,">$ALL/c_header") or die "Unable to open file $ALL/c_header";
#print $c_header."\n";
print FH $c_header;
close(FH);

print STDERR "OK / Grouping of count files in one global count file : $ALL/global_count.csv \n";

#Voir le header pour le modifier avec les noms de fichiers entrants
#print "HEADER : $ALL/header\n\n";
#print "c HEADER : $c_header\n\n";


$cmd = "paste $ALL/header $file_list  > $ALL/global_count.csv ";
#print "$cmd"."\n";
system $cmd;

$cmd = "cat $ALL/c_header $ALL/global_count.csv > $ALL/global_count_with_header.csv ";
system $cmd;
print STDERR "OK / Generation of genes and transcrits global count file. \n";


print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 3 / 6 - Header modification \n";
# modifier l'en-tête des ces fichiers pour ne plus avoir les "/all_transcripts.tsv"
my $LigneGlobal ='';
my $InFileGlobal = "$ALL/global_count_with_header.csv";
open (InFileGlobal,"$InFileGlobal") or die "Can't open $InFileGlobal \n";
open (GLOBAL, ">","$ALL/global_count_genes_transcripts.csv") or die "Cannot create $ALL/global_count_genes_transcripts.csv";
 
while ($LigneGlobal=<InFileGlobal>)
{
    if ($LigneGlobal =~ /raw_transcripts.tsv/)
        {
        ($LigneGlobal =~ s/\/raw_transcripts.tsv/ /g); 
        }
        print GLOBAL $LigneGlobal;
}
close( GLOBAL );
$cmd = "rm $ALL/global_count_with_header.csv";
system $cmd;
print STDERR "OK / Header modification in global_count_genes_transcripts.csv file.\n";
  
  
#######################################################################################################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############     FICHIER DE COMPTAGE DES GENES    + RENAME                   ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
#######################################################################################################
# mettre à jour les noms des gènes et des transcrits dans ces fichiers à partir du fichier de correspondance
#Remplacer les gene_id par les gene_name s'il existe
#Remplacer transcript_id par nearest_ref s'il existe
print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 4 / 6 - Update of genes and transcripts names \n";
my $GTFrename = "$ALL/merged.rename.gtf";
my $correspondance = "$ALL/correspondance.txt";
my $ref ="";
my $gene_id;
my $old_gene_id;
my $transcript_id;
my $old_transcript_id;
my %correspondance = ();
open (GTFrename, ">",$GTFrename) or die "Cannot create $ALL/merged.rename.gtf";
open (Correspondance, ">",$correspondance) or die "Cannot create $ALL/correspondance.txt";
my $InFile = $ARGV[1];
open (InFile,"$InFile") or die "Can't open merged.gtf \n";
    while (<InFile>) {
    if (/gene_id\s+"([^"]+)"/){$old_gene_id=$1;}
    if (/transcript_id\s+"([^"]+)"/){$old_transcript_id=$1;}
        if (/gene_name\s+"([^"]+)"/)
        {
            my $t=qq(gene_id "$1");
            s/gene_id\s+"XLOC_[^"]+"/$t/;
        }
    if (/oId\s+"([^"]+)"/)
    {
        if (substr($1,0,4) ne "CUFF"){
          	my $t=qq(transcript_id "$1");
             	s/transcript_id\s+"TCONS_[^"]+"/$t/;
        }
    #
    #if (/nearest_ref\s+"([^"]+)"/)
    #{  
    #    my $t=qq(transcript_id "$1");
    #    s/transcript_id\s+"TCONS_[^"]+"/$t/;
    }
    if (/gene_id\s+"([^"]+)"/){$gene_id=$1;}
    if (/transcript_id\s+"([^"]+)"/){$transcript_id=$1;}
    if ($gene_id."\t".$transcript_id ne $ref){
    	print Correspondance $old_gene_id."\t".$old_transcript_id."\t".$gene_id."\t".$transcript_id."\n";
    	$correspondance{$old_transcript_id} = $old_gene_id."\t".$old_transcript_id."\t".$gene_id."\t".$transcript_id;
        $ref = $gene_id."\t".$transcript_id;
    }
    print GTFrename ;
}
close (GTFrename) ;
close (Correspondance) ;
print STDERR "OK / Names update from merge.gtf file in merge.rename.gtf file.\n";

#########################################################################################################
################                                                                 ########################
################                                                                 ########################
################                                                                 ########################
################                                                                 ########################
################                     Rename count files                          ########################
################                                                                 ########################
################                                                                 ########################
################                                                                 ########################
#########################################################################################################
#
my $fileGenesTranscrits = "$ALL/global_count_genes_transcripts.csv";
open (COUNT,"$fileGenesTranscrits") or die "Can't open $ALL/global_count_genes_transcripts.csv \n";
open (RENAME, ">", "$ALL/global_count_genes_transcripts_rename.csv") or die "Can't open $ALL/global_count_genes_transcripts_rename.csv'";
while (my $l=<COUNT>){
	chomp $l;
	my @F=split(/\t/,$l);
	if (exists($correspondance{$F[1]})){
            my @tab=split(/\t/,$correspondance{$F[1]});
	    $F[0]=$tab[3];
	    $F[1]=$tab[2];
	}
	printf RENAME "%s\n",join("\t",@F);
}
close RENAME;
close COUNT;
print STDERR "OK / Names update in $ALL/global_count_genes_transcripts_rename.csv file.\n";

my $filegenecount = "$ALL/global_count_genes_transcripts_rename.csv";
my %genes_and_sums;
open ("filegenecount","$filegenecount") or die "Can't open $filegenecount \n";
$_=<filegenecount>;

while (my $line =<filegenecount>){
    chomp $line,
    my @F = split(/\s+/, $line);
    for (my $i=2; $i<= $#F; $i++){
         $genes_and_sums{$F[1]}[$i-2]+=$F[$i];
    }
}


open (SYNTHESE2, ">","$ALL/global_count_genes_rename.csv") or die "Cannot create global_count_genes_rename.csv";
print SYNTHESE2 `head -n 1 $filegenecount | cut -f2,3-` ;

foreach my $key (sort {$a cmp $b} keys %genes_and_sums) {
    printf SYNTHESE2 "%s\t%s\n", $key, join("\t", @{$genes_and_sums{$key}});
      }
      print STDERR "OK / Grouping transcrits counts by genes and condition in $ALL/global_count_genes_conditions.csv avec les genes renommes \n";

print STDERR "OK / gene count file produced ion $ALL/global_count_genes_rename.csv file.\n";
#######################################################################################################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                     FASTA DES TRANSCRITS cuffmerge              ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
#######################################################################################################

print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 5 / 6 - Generation of fasta transcripts file from merged.rename.gtf file generated by Cuffmerge : \n";
# production du fichier fasta des transcrits avec dans l'en-tête :
#    >Ancien_nom_de_transcrits Ancien_nom_de_gene longueur_du_transcrit
my $merged_rename_file = "$ALL/merged.rename.gtf";
my $GTFTOFASTA ="/usr/local/bioinfo/bin/gtf_to_fasta";
my $bank = $ARGV[0];
#ligne de commande : gtf_to_fasta merged.rename.gtf /bank/bowtie2db/ensembl_sus_scrofa_genome results_from_merged_rename.fasta
$cmd = "$GTFTOFASTA $merged_rename_file $bank $ALL/results_from_merged_rename.fasta;";
print STDERR "OK / Production of fasta file. \n";
print STDERR "$cmd \n";
system $cmd;

my $LigneFasta ='';
my $gene_name_cuffmerge;
my $InFileFasta = "$ALL/results_from_merged_rename.fasta";
open (InFileFasta,"$InFileFasta") or die "Can't open $ALL/results_from_merged_rename.fasta \n";
#ouverture du fasta des transcrits Cuffmerge
open (FASTA, ">","$ALL/transcrits_Cuffmerge.fasta") or die "Cannot create $ALL/transcrits_Cuffmerge.fasta";
open (merged_rename_file,"$merged_rename_file") or die "Can't open $merged_rename_file \n"; 
while ($LigneFasta=<InFileFasta>)
{
    if ($LigneFasta =~ /^>.*/)
        {
            while (<merged_rename_file>) 
            { 
            if (/gene_id\s+"([^"]+)"/)
                {
                #$trans_name_cuffmerge = qq(nearest_ref "$1");
                $gene_name_cuffmerge = $1;                    
                }
            }
        ($LigneFasta =~ /^>(\S*)\s(\S*)\s(\S*)/);
        #~ /^>([^\ ]*)\ ([^\ ]*)\ ([^\ ]*)/);
        print FASTA ">".$3." ".$gene_name_cuffmerge." ".$2."\n";
        }
    else
        {
        print FASTA $LigneFasta;
        }
}
close( FASTA );
$cmd = "rm $ALL/results_from_merged_rename.fasta";
system $cmd;
print STDERR "OK / Modification of fasta file header : >Old_transcript_name Old_gene_name transcript_lenght \n";




#######################################################################################################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                     Statistiques R                              ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
#######################################################################################################

#print STDERR "\n __________________________________________________________________________________________ \n";
#print STDERR "\n - 6 / 6 - R statistics : dendogramme and graphic correlation. \n";
##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[2];
##use Cwd qw(abs_path);
##my $home = abs_path($0); #$path
#
##my $log_dir = $home."/".$ALL;
#my $log_dir = dirname("$ALL/global_count_genes_transcripts.csv");
## 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_pdf = "$log_dir/graphic_correlation.pdf";
#$R->send(
#    qq`
#        pdf("$graphic_correlation_pdf")
#        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 + opts(axis.text.x=theme_text(angle=-90))
#	print(q1 + opts(axis.text.x=theme_text(angle=-90)))
#        dev.off()
#    `
#);
#
## Création d'une image graphic correlation
#my $dendogramme_pdf = "$log_dir/dendogramme.pdf";
#$R->send(
#    qq`
#        pdf("$dendogramme_pdf")
#        heatmap(cor(log10(MONOPOLY+1)))
#        dev.off()
#    `
#);
#
## Fermeture du pont
#$R->stopR();
#
#print STDERR "OK / Graphic correlation \n";
#print STDERR "OK / Dendogramme \n";
#
#`rm $ALL/global_count_nogene.csv; rm $ALL/header; rm $ALL/c_header; rm $ALL/global_count.csv; rm $ALL/global_count_genes_transcripts.csv`;
#
#print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "THE END ...\n";


