#!/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;

##------------------------------------------------------------
## options and parameters
##------------------------------------------------------------
## Options definitions (name,format,default,required)
#my @Options = ( ['help',undef,undef,0],
#		['man',undef,undef,0],
#	      );
## defaults values
#my %Values=();
#my %Required=();
#map {$Values{$_->[0]}=$_->[2];$Required{$_->[0]}=$_->[3];} @Options;
## build options list
#my %Options=();
#map {$Options{defined($_->[1])?$_->[0].'='.$_->[1]:$_->[0]}=\$Values{$_->[0]};} @Options;
## retrieve options
#GetOptions(%Options) || pod2usage(2);
#pod2usage(1) if ($Values{'help'});
#pod2usage(-exitstatus => 0, -verbose => 2) if ($Values{'man'});
#
## check required
#map {pod2usage("$_ is required") if ($Required{$_}&&!defined($Values{$_}));} keys %Values;


#Verification des parametres en entree ($#ARGV est le nombre d'arguments).
#$#ARGV==4 || pod2usage("missing one parameter");


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

#my $list = `\\ls *_L*/raw_transcripts.tsv`; Possible en executant le script dans un rep projet avec les fichiers issus du sigcufflings sous la main.
my $list = $ARGV[4];
my $cmd = '';
my $i = '';
my $j = '';
my $file_list = '' ;
my $c_header = 'gene	transcript	';
#my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time); 
my $ALL = $ARGV[3];
#my $ALL =  $dirnameresults."RESULTATS".$mday."-".($mon+1)."_".$hour."h".$min."mn".$sec."sec"; 


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


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


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

#GENERATION DU FICHIER DES COMPTAGE GENE ET TRANSCRITS
print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 2 - Generate count file for genes and transcripts : \n";
print STDERR "OK / Sort of count file.\n";
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;
    }
print STDERR "OK / counts.\n";
$c_header .= "\n";
#print $j."\n";
$cmd = "sort -S 3G -k2 $j | cut -f1,2 > $ALL/header";
print STDERR "OK / header. \n";
system $cmd;

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";

$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";

##Fichiers de comptages pour les gènes (en faisant la sommes des transcrits d'un gène, enlever la colonne transcrit)
#my $filecomptage1 = "$ALL/global_count_with_header.csv";
#$cmd = "cut -f1 $filecomptage1 > $ALL/genes_tmp.csv; sed '1d' $ALL/genes_tmp.csv > $ALL/genes.csv ";
#system $cmd;
#
#
##$cmd = "cut -f3- $filecomptage1 > $ALL/to_count_tmp.csv ; sed '1d' $ALL/to_count_tmp.csv > $ALL/to_count.csv";
#$cmd = "sed '1d' $filecomptage1 > $ALL/to_count_tmp.csv ; cut -f3- $ALL/to_count_tmp.csv > $ALL/to_count.csv";
#system $cmd;
#my $filecomptage2 = "$ALL/to_count.csv";
#my $Lignecomptage2 ='';
#my @ligne;
#open (filecomptage2,"$filecomptage2") or die "Can't open $filecomptage2 \n";
#open (GLOBALTranscripts, ">","$ALL/global_count_without_transcripts_somme.csv") or die "Cannot create $ALL/global_count_without_transcripts_somme.csv";
# 
#while ($Lignecomptage2=<filecomptage2>)
#{	
#	@ligne=split(/\t/,$Lignecomptage2);
#	my $somme = sum(@ligne)."\n";
#	print GLOBALTranscripts $somme;
#}
#close( GLOBALTranscripts );
#$cmd = "paste $ALL/genes.csv $ALL/global_count_without_transcripts_somme.csv > $ALL/global_count_genes.csv ;
#rm $ALL/global_count_without_transcripts_somme.csv;
#rm $ALL/to_count_tmp.csv;
#rm $ALL/genes_tmp.csv;
#rm $ALL/to_count.csv;
#rm $ALL/genes.csv;";
#system $cmd;
#print "OK / Remove transcripts column in the second count file. \n";



print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 3 - 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";
  



print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 4 - File : Global count genes per condition \n";
#Utiliser un fichier $ALL/global_count_with_header.csv 
#my $filegenecount = "RESULTATS4-4_11h48mn11sec/global_count_genes_transcripts.csv";
#my $filegenecount = "$ALL/global_count_with_header.csv";
my $filegenecount = "$ALL/global_count_genes_transcripts.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[0]}[$i-2]+=$F[$i];
    }
}
 
#Ecriture du fichier  
open (SYNTHESE2, ">","$ALL/global_count_gene_conditions.csv") or die "Cannot create global_count_gene_conditions.csv";
print SYNTHESE2 `head -n 1 $filegenecount | cut -f1,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";  


  
#######################################################################################################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############     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 - 5 - Update of genes and transcripts names \n";
open (GTFrename, ">","$ALL/merged.rename.gtf") or die "Cannot create $ALL/merged.rename.gft";
my $InFile = $ARGV[1];
open (InFile,"$InFile") or die "Can't open merged.gtf \n";
    while (<InFile>) {
        if (/gene_name\s+"([^"]+)"/)
        {
            my $t=qq(gene_id "$1");
            s/gene_id\s+"XLOC_[^"]+"/$t/;
        }
    
        if (/nearest_ref\s+"([^"]+)"/)
        {  
            my $t=qq(transcript_id "$1");
            s/transcript_id\s+"TCONS_[^"]+"/$t/;
        }
        print GTFrename ;
}
close (GTFrename) ;
print STDERR "OK / Names update from merge.gtf file in merge.rename.gtf file.\n";


#######################################################################################################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############       FICHIER DE CORRESPONDANCE DES GENES ET DES TRANSCRITS     ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
#######################################################################################################

print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 6 - Correspondance count file for genes \n";
#POUR LES GENES
# génération du fichier de correspondance entre les anciens et les nouveaux noms des gènes 
#et des transcrits (ne prendre le nouveau que quand ce n'est pas un CUFF*),
my $LigneTab ='';
my $Ligne ='';


open (InFile,"$InFile") or die "Can't open All/merged.gtf \n";

#ouverture d'un fichier de correspondance entre les anciens et les nouveaux noms
open (TXT, ">","$ALL/correspondances_genes.txt") or die "Cannot create $ALL/correspondances_genes.txt";
print TXT ("Gene id \t  Gene name \n");
 
while ($Ligne=<InFile>)
{
        if($Ligne=~ /(XLOC_\d\d\d\d\d\d)/)
        {
                print TXT ($1."\t");
        }  
        if ($Ligne=~ /gene_name\s+"([^"]+)"/)
        {
                print TXT ($1."\n");                
        }
}
close( TXT );

open (TXTsorted, ">","$ALL/correspondances_genes_sorted.txt") or die "Cannot create $ALL/correspondances_genes_sorted.txt";
    my $file = "$ALL/correspondances_genes.txt";
    open (TXT, "tail -n+2 $file | sort -S 3G -u |");
    while (<TXT>)
        {print TXTsorted $_;
    }
close (TXTsorted) ;


print STDERR "OK / Correspondance count file for GENES : $ALL/correspondances_genes_sorted.txt \n";

#POUR LES TRANSCRITS 
print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 7 - Correspondance count file for transcripts \n";
my $LigneT ='';
#ouverture d'un fichier de correspondance entre les anciens et les nouveaux noms
open (TXTT, ">","$ALL/correspondances_transcripts.txt") or die "Cannot create $ALL/correspondances_transcripts.txt";
    print TXTT ("Transcrit id \t nearest_ref \n");
open (InFile,"$InFile") or die "Can't All/open merged.gtf \n";

    while ($LigneT=<InFile>)
    {
        if($LigneT=~ /(TCONS_\d\d\d\d\d\d\d\d)/)
        {
                print TXTT ($1."\t");
        }  
         if ($LigneT=~ /nearest_ref\s+"([^"]+)"/)
        {  
              print TXTT ($1."\n");
        }     
    }
close( TXTT );
close (InFile);

open (TXTsortedT, ">","$ALL/correspondances_transcripts_sorted.txt") or die "Cannot create $ALL/correspondances_transcripts_sorted.txt";
    my $fileT = "$ALL/correspondances_transcripts.txt";
    open (TXTT, "tail -n+2 $fileT | sort -S 3G -u |");
    while (<TXTT>)
        {print TXTsortedT $_;
    }
close (TXTsortedT) ;
close (TXTT) ;

print STDERR "OK / Correspondance count file for transcripts : $ALL/correspondances_transcripts_sorted.txt \n";


#Renommage des genes  et des transcrits
my $fileGenesTranscrits = "$ALL/global_count_genes_transcripts.csv";

my $correstranscrits = "$ALL/correspondances_transcripts_sorted.txt";
open (CORRESTRANSCRITS,"$correstranscrits") or die "Can't open $ALL/correspondances_transcripts_sorted.txt \n";
my %corresp =();
while(my $l = <CORRESTRANSCRITS>){
	chomp $l;
	my @F=split(/\t/,$l);
    $corresp{$F[0]}= $F[1];
    }
close CORRESTRANSCRITS;

my $corresgenes = "$ALL/correspondances_genes_sorted.txt";
open (CORRESGENES,"$corresgenes") or die "Can't open $ALL/correspondances_genes_sorted.txt \n";
my %correspG =();
while(my $l = <CORRESGENES>){
	chomp $l;
	my @F=split(/\t/,$l);
    $correspG{$F[0]}= $F[1];
    }
close CORRESGENES;


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($corresp{$F[1]})){
	    $F[1]=$corresp{$F[1]};
	}
	if (exists($correspG{$F[0]})){
	    $F[0]=$correspG{$F[0]};
	}
	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";


#Renommage des genes:
my $fileGenes = "$ALL/global_count_gene_conditions.csv";

my $corresgenestab = "$ALL/correspondances_genes_sorted.txt";
open (CORRESGENESTAB,"$corresgenestab") or die "Can't open $ALL/correspondances_genes_sorted.txt \n";
my %corresp2 =();
while(my $l = <CORRESGENESTAB>){
	chomp $l;
	my @F=split(/\t/,$l);
    $corresp2{$F[0]}= $F[1];
    }
close CORRESGENESTAB;

open (COUNTGENES,"$fileGenes") or die "Can't open $ALL/global_count_gene_conditions.csv \n";
open (RENAMEGENES, ">", "$ALL/global_count_gene_conditions_rename.csv") or die "Can't open $ALL/global_count_gene_conditions_rename.csv'";
while (my $l=<COUNTGENES>){
	chomp $l;
	my @F=split(/\t/,$l);
	if (exists($corresp2{$F[0]})){
	    $F[0]=$corresp2{$F[0]};
	}
	printf RENAMEGENES "%s\n",join("\t",@F);
}
close RENAMEGENES;
close COUNTGENES;
print STDERR "OK / Names update in $ALL/global_count_gene_conditions_rename.csv file.\n";

#suppression des fichiers intermediaires inutiles
$cmd = "$ALL/rm correspondances_genes.txt; rm $ALL/correspondances_transcripts.txt; rm $ALL/global_count_gene_conditions.csv; ";
system $cmd;

#######################################################################################################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
##############                     FASTA DES TRANSCRITS cuffmerge              ########################
##############                                                                 ########################
##############                                                                 ########################
##############                                                                 ########################
#######################################################################################################

print STDERR "\n __________________________________________________________________________________________ \n";
print STDERR "\n - 8 - 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 $trans_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 (/transcript_id\s+"([^"]+)"/)
                {
                #$trans_name_cuffmerge = qq(nearest_ref "$1");
                $trans_name_cuffmerge = $1;                    
                }
            }
        ($LigneFasta =~ /^>(\S*)\s(\S*)\s(\S*)/);
        #~ /^>([^\ ]*)\ ([^\ ]*)\ ([^\ ]*)/);
        print FASTA ">".$trans_name_cuffmerge." ".$3." ".$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 - 9 - 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];
#my $log_dir = $home."/".$ALL;
my $log_dir = $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 $dendogramme_pdf = "$log_dir/dendogramme.pdf";
$R->send(
    qq`
        pdf("$dendogramme_pdf")
        heatmap(cor(log10(MONOPOLY+1))) 
        dev.off()
    `
);


## 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))
#        y1.m <- melt(y)
#        q1 <- ggplot(y1.m, aes(X1, X2, fill = value)) + geom_tile() + scale_fill_gradient(low = "blue", high = "yellow")
#        q1 + theme(axis.text.x = element_text(angle =-90))
#        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; rm $ALL/correspondances_genes.txt`;

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

