﻿#!/usr/bin/perl

#Author
# Nicolas Cerveau
# University of Göttingen
# Geobiology department
# Evolution of the Metazoa research group
# Germany

# nicolas.cerveau1@gmail.com


# The aim of this script is to make serialy all the step needed in order to concatenate transcriptome outputs from several software and remove the coding redundancy.
# They are 6 steps that are explained in the figure 1 of our paper from step 4 to 9.
# The output in the folder 9-concatenated_assembly is composed of 3 files containing the transcripts, the CDS and peptides associated. The output file format is compatible
# with Trinotate software version 3 to perform the transcriptome annotation. There is a fourth file that contain informations about the CDS clustering step and the transcript
# selection.

#requirements
#	- cd-hit-est: http://weizhongli-lab.org/cd-hit/download.php
#	- TransDecoder V3: https://transdecoder.github.io/
#	- Warning: This version of our script is not working with TransDecoder V2, if you want a version working with TransDecoder V2 please use the email adress above.

# Mandatory options
#	- path of the folder containing the raw assembly outputs (-i)
#	- path for Transdecoder LongOrf and predict binary files (-Tr)
#	- path for cd-hit-est binary file (-Cd)

#Warning: This script was design and tested on unix system but is probably not compatible with Windows system

#example of command line:
#perl Supplementary_file2.pl -i raw_assemblies/ -w Concatenated_assembly -m 100 -aS 1.00 -Cd /usr/users/ncervea/bin -Tr /usr/users/ncervea/bin

#File names in the raw_assemblies folder for Drosophila simulans:
#	- Dsim_assembly_Trinity.fasta
#	- Dsim_assembly_CLC.fasta
#	- Dsim_assembly_idba.fasta
#	- Warning: The files name have to strickly follow this pattern: "organisme name or shortcut without space or underscore"_assembly_"assembler name without space or underscore".fasta
#	- Organism name or shortcut will be used to name files created during the process.
#	- The three first letter of the assembler name will be used to name the transcript in the concatenated assembly to follow their origin.

#to use a time point

use Time::localtime;

####################################################################################################################################################################################################################################################################
# Beginning of the work 
####################################################################################################################################################################################################################################################################

my $start_time=time();

# Variables
###########
my $a;
my $b = 0;
my $table;
my $c = 0;
my $table1;
my $d = "1.00";
my $assembler_shortcut;
my $liste;
my $liste_assembler;
my $working = "concatenated_assembly";
my $size = 100;
my $aS = 1.00;
my $div = 0.98;
my $threads = 5;

# Argument recovery
###################
for my $args (0 .. $#ARGV){

 SWITCH: {

    # path of the folder containing the raw assembly outputs. All the file names have to have this format: "studied species name"_assembly_"assembler name".fasta
	# example: Dsim_assembly_CLC.fasta.
	# Warning: do not use any other underscore ("_") in the assembly output file names than the one before and after "assembly" as those characters are used to 
		# split the file name and get the sample name and the assembler name for each file.
    if( $ARGV[$args] =~ /^-i/)
      { $args++; $input = $ARGV[$args]; last SWITCH; }

    # peptide size threshold for transdecoder (default value is 100 amino acides)
    if( $ARGV[$args] =~ /^-m/)
      { $args++; $size = $ARGV[$args]; last SWITCH; }

    # ratio of the smallet sequence for cds clustering (default value is 1.00, 100% of the smallest size)
    if( $ARGV[$args] =~ /^-aS/)
      { $args++; $aS = $ARGV[$args]; last SWITCH; }

    # nucleotide divergence for cds clustering (default value is 0.98, 98% nucleotide identity)
    if( $ARGV[$args] =~ /^-d/)
      { $args++; $div = $ARGV[$args]; last SWITCH; }

    # thread number for cds clustering (default value is 5)
    if( $ARGV[$args] =~ /^-Th/)
      { $args++; $threads = $ARGV[$args]; last SWITCH; }

    # path of the working folder in order to create all the files and folders (default is creating a concatenated_assembly directory in the current directory)
    if( $ARGV[$args] =~ /^-w/)
      { $args++; $working = $ARGV[$args]; last SWITCH; }

    # path for Transdecoder LongOrf and predict binary files
    if( $ARGV[$args] =~ /^-Tr/)
      { $args++; $TransDecoder = $ARGV[$args]; last SWITCH; }

    # path for cd-hit-est binary file
    if( $ARGV[$args] =~ /^-Cd/)
      { $args++; $cd_hit_est = $ARGV[$args]; last SWITCH; }

    # use 1 to clean the folders from intermediate step, use 0 or no option if you want to keep all the intermediate folders
    if( $ARGV[$args] =~ /^-Cl/)
      { $args++; $clean = $ARGV[$args]; last SWITCH; }
 }
}

# Recover the sample and assembler names 
####################################################################################################################################################################################################################################################################

#open the input folder
opendir FOLDER,"$input" or die "The input folder was not found";

#reading the foler
while ($file = readdir FOLDER){

	#extract the first character of the file name
	$a = substr $file,0,1;

	#if $a variable is not a dot
	if ($a ne "."){

		#Split the file name with "." character
		@decompose1 = split (/\./,$file);

		#Split the file name with "_" character
		@decompose = split (/\_/,$decompose1[0]);

		#create the assembler name shortcut
		$assembler_shortcut = substr $decompose[2],0,3;

		#store sample name and assembler names
		######################################
		#create the table name
		$table = "assembler".$b;

		#print table name to control
		#print $table,"\n";

		#store the info into a table
		push @$table, {	sample_name =>$decompose[0],
				assembler_name =>$decompose[2],
				assembler_shortcut =>$assembler_shortcut,
					};

		#implement the $b variable
		$b ++;
	}
}

closedir FOLDER;

# intra-assembly analyses 
####################################################################################################################################################################################################################################################################

#create a folder to contain the data of size verification and stats
###################################################################
system "mkdir $working";
system "mkdir $working/2-size_verification";

#treat single assembler output one by one
while ($c < $b){

	#create table name
	$table1 = "assembler".$c;

	#Statistic on single assembler output

	#open single assembler output file
	open FICHIER,"<$input/$$table1[0]{'sample_name'}_assembly_$$table1[0]{'assembler_name'}.fasta" or die "The single assembler output file for $$table1[0]{'assembler_name'} assembler was not found!";

	#Open log file
	open WRITER,">$working/2-size_verification/Log_$$table1[0]{'sample_name'}_assembly_$$table1[0]{'assembler_name'}_statistic.txt" or die "The single assembler static log file for $$table1[0]{'assembler_name'} assembler was not created! ";

	#variables
	my $x;
	my $nom = "";
	my $size;
	my $seq;
	my $sum = 0;
	my $contig_number = 0;
	my $min = "";
	my $max = "";
	my $average;
	my $contig_number_sup_300 = 0;
	my $sum_sup_300 = 0;
	my $average_sup_300;

	#reading of single assembler output file
	while ($line=<FICHIER>){

		#extraction of first line character
		$x = substr $line,0,1;

		#get contig name
		if ($x eq ">"){

			#if the variable is not empty
			if ($nom ne ""){

				#get sequence information
				$size = length($seq);
				$sum = $sum + $size;
				$contig_number ++;

				#store the size into a table
				push(@lengths, $size);

				#empty $seq variable
				$seq = "";

				#store min and max variable
				#if it is the first transcript
				if($min eq ""){

					#store the data
					$min = $size;
					$max = $size;
				}

				#if the transcript size is lower than the actual minimum
				elsif($size < $min){

					#store the data
					$min = $size;
				}

				#if the transcript size is higher than the actual maximum
				elsif($size > $max){

					#store the data
					$max = $size;
				}
			}

			#Store the new transcript name
			$line =~ s/\n//g;
			$line =~ s/\r//g;
			$line =~ s/\>//g;
			$nom = $line;
		}

		#get transcript sequence
		elsif ($x ne ">"){

			#remove useless character
			$line =~ s/\n//g;
			$line =~ s/\r//g;
			$line =~ s/\ //g;

			#store sequence into a variable
			$seq = $seq.$line;
		}
	}

	#get sequence information for the last contig
	$size = length($seq);
	$sum = $sum + $size;
	$contig_number ++;

	#store the size into a table
	push(@lengths, $size);

	#if the transcript size is lower than the actual minimum
	if($size < $min){

		#store the data
		$min = $size;
	}

	#if the transcript size is higher than the actual maximum
	elsif($size > $max){

		#store the data
		$max = $size;
	}

	#calculate average
	$average = $sum/$contig_number;

	#to do in order to get the N50 of the whole transcriptome
	@lengths = reverse sort {$y <=> $z} @lengths;

	my $rsum = 0;
	my $i = -1;
	while ($rsum < $sum/2) {
	  $i++;
	  $rsum += $lengths[$i];
	}
	my $N50 = $lengths[$i];

	#print info into log file
	print WRITER "The total size of the assembly is $sum bp and there are $contig_number contigs\n";
	print WRITER "The minimum contig size is $min bp, the maximum is $max bp, the N50 is $N50 and the average is $average bp\n\n";

	#Fermer le fichier de sortie des données
	close WRITER;

	#implement $c variable
	$c ++;

}

#reset $c variable to 0
$c = 0;

#modify the transcript sequence file to put the transcript in only one line
###########################################################################
system "mkdir $working/3-modifed_sequence";

#treat single assembler output one by one
while ($c < $b){

	#create table name
	$table1 = "assembler".$c;

	#open single assembler output file
	open FICHIER,"<$input/$$table1[0]{'sample_name'}_assembly_$$table1[0]{'assembler_name'}.fasta" or die "le fichier n'existe pas!";

	#create new file name
	my $new_nom = "$$table1[0]{'sample_name'}_assembly_$$table1[0]{'assembler_name'}.fasta";
	$new_nom =~ s/\.fasta//g;
	open WRITER,">$working/3-modifed_sequence/".$new_nom."_modif.fasta";;

	#read single assembler output file
	while ($line=<FICHIER>){

		#extract first line character
		$a = substr $line, 0, 1;

		#if the line contain sequence name
		if ($a eq ">"){
			$line =~ s/\n//g;
			$line =~ s/\r//g;
			print WRITER "\n$line\n";
		}

		#if the line contain no sequence name
		else{
			$line =~ s/\n//g;
			$line =~ s/\r//g;
			print WRITER $line;
		}
	}

	#close single assembler output file
	close FICHIER;

	#close output file
	close WRITER;

	#implement $c variable
	$c ++;
}

#reset $c variable to 0
$c = 0;

#4-rename the contig name
#######################
system "mkdir $working/4-rename";

#treat single assembler output one by one
while ($c < $b){

	#create table name
	$table1 = "assembler".$c;

	#Variable
	my $a;
	my $read_number = 0;

	#open single assembler output file
	open INPUT,"<$working/3-modifed_sequence/$$table1[0]{'sample_name'}_assembly_$$table1[0]{'assembler_name'}_modif.fasta" or die "The raw file was not found";

	#open output file with modified name
	open OUTPUT,">$working/4-rename/$$table1[0]{'sample_name'}_assembly_$$table1[0]{'assembler_name'}_renamed.fasta" or die "The output file was not created";

	#set prefix name
	$prefix = $$table1[0]{'assembler_shortcut'};

	#reading file
	while ($line = <INPUT>){

		#extract the first character of each line
		$a = substr $line,0,1;

		#To do if $a is a >
		if ($a eq ">"){

			#print the new sequence name
			print OUTPUT ">$prefix"."_$read_number\n";

			#increase the read_number
			$read_number = $read_number + 1;
		}

		#To do if $a is not a >
		elsif ($a ne ">"){

			#Print directly the line in the output
			print OUTPUT $line;

		}	
	}

	#close input and output files
	close INPUT;
	close OUTPUT;

	#implement $c variable
	$c ++;
}

#reset $c variable to 0
$c = 0;

#create a folder to contain the transcripts clustering data and stats
########################################################################
system "mkdir $working/5-transcripts_clustering";

#intra_assembly transcripts clustering
######################################

#treat single assembler output one by one
while ($c < $b){

	#store the starting time
	my $clustering_start_time=time();

	#create table name
	$table1 = "assembler".$c;

	#clustering
	system "$cd_hit_est/cd-hit-est -i $working/4-rename/$$table1[0]{'sample_name'}_assembly_$$table1[0]{'assembler_name'}_renamed.fasta -o $working/5-transcripts_clustering/$$table1[0]{'sample_name'}_unigene_$$table1[0]{'assembler_name'}.fasta -c 1.00 -M 2549 -T $threads -G 0 -aS 1.00 -aL 0.005 >log_$$table1[0]{'sample_name'}_unigene_$$table1[0]{'assembler_name'}.txt";
	system "mv log_$$table1[0]{'sample_name'}_unigene_$$table1[0]{'assembler_name'}.txt $working/5-transcripts_clustering/";

	#print infos
	print LOG "\tIntra-assembly transcripts clustering done in ",time()-$clustering_start_time,"sec for $$table1[0]{'assembler_name'} assembler data\n";

	#implement $c variable
	$c ++;
}

#reset $c variable to 0
$c = 0;

#create a folder to contain the concatenated unigene data
#########################################################
system "mkdir $working/6-Merged_unigene";

#concatenate the data for the assemblers
########################################

#create a return file
open RETURN,">$working/return.fasta" or die "The return file was not created";
print RETURN "\n";
close RETURN;

#treat single assembler output one by one
while ($c < $b){

	#create table name
	$table1 = "assembler".$c;

	#create a list of element to concatenate
	$liste = $liste." $working/5-transcripts_clustering/$$table1[0]{'sample_name'}_unigene_$$table1[0]{'assembler_name'}.fasta $working/return.fasta";

	#create a list of assembler name
	$liste_assembler = $liste_assembler.$$table1[0]{'assembler_name'}."_";

	#implement $c variable
	$c ++;
}

#concatenate the file
system "cat $liste > $working/6-Merged_unigene/$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta";

#reset $c variable to 0
$c = 0;

#remove return file
system "rm $working/return.fasta";


#create a folder to contain the Transdecoder data and stats
###########################################################
system "mkdir $working/7-Transdecoder";

#Transdecoder ORF detection
###########################

#ORF detection

#store the starting time
my $Transdecoder_start_time=time();

#make the Transdecoder LongOrfs
system "$TransDecoder/TransDecoder.LongOrfs -t $working/6-Merged_unigene/$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta -m $size >log_$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene_Transdecoder_LongOrf.txt";

#make the Transdecoder Predict
system "$TransDecoder/TransDecoder.Predict -t $working/6-Merged_unigene/$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta >log_$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene_Transdecoder_Predict.txt";

#move the results files
system "mv $$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta.transdecoder.* $working/7-Transdecoder";
system "mv log_$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene_Transdecoder_LongOrf.txt $working/7-Transdecoder";
system "mv log_$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene_Transdecoder_Predict.txt $working/7-Transdecoder";
system "mv $$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta.transdecoder_dir $working/7-Transdecoder";

#create a folder to contain the transcripts clustering data and stats
########################################################################
system "mkdir $working/8-CDS_clustering";

#grouped CDS clustering
##############################

#store the starting time
my $clustering_start_time=time();

#clustering c0.98_aS1.00
system "$cd_hit_est/cd-hit-est -i $working/7-Transdecoder/$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta.transdecoder.cds -o $working/8-CDS_clustering/$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta.transdecoder_grouped_unicds_c$div"."_aS$aS.fasta -c $div -M 2549 -T $threads -G 0 -aS $aS -aL 0.005 -d 0 >log_$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta.transdecoder_grouped_unicds_c$div"."_aS$aS.txt";
system "mv log_$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta.transdecoder_grouped_unicds_c$div"."_aS$aS.txt $working/8-CDS_clustering/";

#create a folder to contain the files for annotation
########################################################################
system "mkdir $working/9-concatenated_assembly";

#store the starting time
my $creation_file_start_time=time();

#file creation

#open output file for transcript
open LOG_OUT,">$working/9-concatenated_assembly/Log_$$table1[0]{'sample_name'}_$liste_assembler"."for_annotation.txt" or die "The output log file was not created during creation of files for annotation";

#Trinity transcript table creation

#Starting time for global transcript table treatment
my $starting_time_transcript_storage = time();
#Variable
#########
my $h;
my $transcript_name;
my $transcript_seq;
my $transcript_size;

#treat single assembler output one by one
while ($c < $b){

	#create table name
	$table1 = "assembler".$c;

	#open renamed transcript file
	open TRANSCRIPTS,"<$working/4-rename/$$table1[0]{'sample_name'}_assembly_$$table1[0]{'assembler_name'}_renamed.fasta" or die "The $$table1[0]{'assembler_name'} renamed transcripts file was not found during creation of files for annotation";

	#create table name
	$table2 = $$table1[0]{'assembler_shortcut'};

	#reading file
	while ($line1 = <TRANSCRIPTS>){

		#extract the first character of each line
		$h = substr $line1,0,1;

		#if $e is a >
		if ($h eq ">"){

			#if the transcript name variable is not empty
			if ($transcript_name ne ""){

				#calculation of the transcript size
				$transcript_size = length($transcript_seq);

				#store the info into a table
				push @$table2, {	transcript_name =>$transcript_name,
							transcript_sequence => $transcript_seq,
							transcript_size => $transcript_size,
							transcript_freq => 0,
							transcript_picking => 0,
									};
			}

			#Clean the transcript name
			$line1 =~ s/\n//g;
			$line1 =~ s/\r//g;
			$line1 =~ s/\>//g;

			#Store the transcript name
			$transcript_name = $line1;
		}

		#if $e is not a >
		elsif ($h ne ">"){

			#Clean the transcript sequence
			$line1 =~ s/\n//g;
			$line1 =~ s/\r//g;
			$line1 =~ s/\ //g;

			#Store the transcript sequence
			$transcript_seq = $line1;
		}
	}

	#store the info into a table
	push @$table2, {	transcript_name =>$transcript_name,
				transcript_sequence => $transcript_seq,
				transcript_size => $transcript_size,
				transcript_freq => 0,
				transcript_picking => 0,
						};

	#reset $transcript_* variable
	$transcript_name = "";
	$transcript_seq = "";
	$transcript_size = "";

	#close the renamed transcript file
	close TRANSCRIPTS;

	#implement $c variable
	$c ++;
}

#reset $c variable to 0
$c = 0;


#print general computation time
print LOG_OUT "Transcript storage for the single assembler ended in ",time() - $starting_time_transcript_storage," sec.\n\n";

#print general computation time
print LOG_OUT "Transcript storage for the 3 assembler ended in ",time() - $starting_time_transcript_storage," sec.\n\n";

#Variable
#########
my $a;
my $size_ref = "";
my $name_ref = "";
my $b = 0;
my $c = 0;
my $nb_single_CDS_cluster = 0;
my $nb_transcript_to_pick = 0;
my $d = 0;
my $nb_cluster_multiple_CDS_without_multiple_ref_size = 0;
my $nb_cluster_multiple_CDS_with_multiple_ref_size = 0;
my $e = 0;
my $table_name = "multiple_CDS".$e;
my $global_transcript_table_Size;
my $nb_CDS_equal_size_to_ref = 0;
my $multiple_occured_transcript_Size;


#clustering information file treatment
#Starting time for clustering file treatment
my $starting_time_clustering_file_treatment = time();

#open clustering info file
open INFO,"<$working/8-CDS_clustering/$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta.transdecoder_grouped_unicds_c$div"."_aS$aS.fasta.clstr" or die "The clustering info file was not found during creation of files for annotation";

#reading file
while ($line = <INFO>){

	#extract the first characters of the line
	$a = substr $line,0,1;

	#to do if the first character is ">"
	if ($a eq ">" and $c ne "0"){

		#count the cluster number
		$nb_cluster ++;

		#to do if the CDS number per cluster is 1 ($c eq 1)
		if ($c eq "1"){

			#split the transcript name
			@decompose5 = split (/\_/,$name_ref);

			#Store the assembler and line info into variable
			$assembler = $decompose5[0];
			$line = $decompose5[1];

			#implement the transcript frequency
			$$assembler[$line]{'transcript_freq'} ++;

			#implement the number of cluster with single CDS
			$nb_single_CDS_cluster ++;

			#implement the number of transcript to pick
			$nb_transcript_to_pick ++;

			#store the transcript into the transcript list to pick
			push @table_transcript_to_pick, {	transcript_name =>$name_ref,
										};			
		}

		#to do if the CDS number per cluster is higher than 1 ($c > 1)
		if ($c > "1"){

			#read the cluster table
			foreach my $i (0..$#table){

				#if the CDS size is equal to the ref
				if ($table[$i]{'size'} eq $size_ref){

					#store the data in a table cluster specific
					push @$table_name, {	transcript_name =>$table[$i]{'transcript_name'},
									};

					#split the transcript name
					@decompose6 = split (/\_/,$table[$i]{'transcript_name'});

					#Store the assembler and line info into variable
					$assembler = $decompose6[0];
					$line = $decompose6[1];

					#implement the transcript frequency
					$$assembler[$line]{'transcript_freq'} ++;

					#implement the number of CDS with equal size to ref
					$nb_CDS_equal_size_to_ref ++;

					#implement $d variable
					$d ++;
				}
			}

			#to do if $d is equal to 0
			if ($d eq "0"){

				#implement the number of cluster without multiple CDS but not with length equal to the ref
				$nb_cluster_multiple_CDS_without_multiple_ref_size ++;

				#implement the number of transcript to pivk
				$nb_transcript_to_pick ++;

				#store the transcript into the transcript list to pick
				push @table_transcript_to_pick, {	transcript_name =>$name_ref,
											};

				#split the transcript name
				@decompose7 = split (/\_/,$name_ref);

				#Store the assembler and line info into variable
				$assembler = $decompose7[0];
				$line = $decompose7[1];

				#implement the transcript frequency
				$$assembler[$line]{'transcript_freq'} ++;
			}

			#to do if $d is superior to 0
			elsif ($d > "0"){

				#store the data in a table cluster specific
				push @$table_name, {	transcript_name =>$name_ref,
								};

				#split the transcript name
				@decompose8 = split (/\_/,$name_ref);

				#Store the assembler and line info into variable
				$assembler = $decompose8[0];
				$line = $decompose8[1];

				#implement the transcript frequency
				$$assembler[$line]{'transcript_freq'} ++;

				#implement the number of cluster with multiple CDS but not with length equal to the ref
				$nb_cluster_multiple_CDS_with_multiple_ref_size ++;

				#implement the number of CDS with equal size to ref
				$nb_CDS_equal_size_to_ref ++;

				#reset $d variable
				$d = 0;

				#implement $e variable
				$e ++;

				#create new table name
				$table_name = "multiple_CDS".$e;
			}			
		}

		#reset the $c value
		$c = 0;

		#clean the @table
		undef @table;
	}

	#to do if the first character is not a ">"
	elsif($a ne ">"){

		#count the CDS number
		$nb_CDS ++;

		#decompose line
		@decompose = split (/\t/,$line);

		#redecompose line
		@decompose1 = split (/\ /,$decompose[1]);
	
		#remove useless characters
		$decompose1[2] =~ s/\n//g;
		$decompose1[2] =~ s/\r//g;

		#to do if the character is a "*"
		if ($decompose1[2] eq "*"){

			#remove useless characters
			$decompose1[1] =~ s/\>//g;
			$decompose1[1] =~ s/\.//g;

			#decompose the line to get the transcript name
			@decompose2 = split (/\::/,$decompose1[1]);

			#Storing size_ref info
			$size_ref = substr $decompose1[0], 0, -3;
			$name_ref = $decompose2[1];

			#implement $c variable
			$c ++;
		}

		#to do if the character is not a "*"
		if ($decompose1[2] ne "*"){

			#remove useless characters
			$decompose1[1] =~ s/\>//g;
			$decompose1[1] =~ s/\.//g;

			#store the size into a variable
			$size =  substr $decompose1[0], 0, -3;

			#decompose the line to get the transcript name
			@decompose2 = split (/\::/,$decompose1[1]);

			#Storing cluster mapping contig info into a table
			push @table, {	transcript_name =>$decompose2[1],
						size =>$size,
								};

			#implement $c variable
			$c ++;
		}		
	}
}

#to store the info of the final cluster

#count the cluster number
$nb_cluster ++;

#to do if the CDS number per cluster is 1 ($c eq 1)
if ($c eq "1"){

	#implement the number of cluster with single CDS
	$nb_single_CDS_cluster ++;

	#implement the number of transcript to pivk
	$nb_transcript_to_pick ++;

	#store the transcript into the transcript list to pick
	push @table_transcript_to_pick, {	transcript_name =>$name_ref,
								};

	#split the transcript name
	@decompose9 = split (/\_/,$name_ref);

	#Store the assembler and line info into variable
	$assembler = $decompose9[0];
	$line = $decompose9[1];

	#implement the transcript frequency
	$$assembler[$line]{'transcript_freq'} ++;
}

#to do if the CDS number per cluster is higher than 1 ($c > 1)
if ($c > "1"){

	#store the cluster contig mapping name into a hash table
	foreach my $i (0..$#table){

		#if the CDS size is equal to the ref
		if ($table[$i]{'size'} eq $size_ref){

			#store the data in a table cluster specific
			push @$table_name, {	transcript_name =>$table[$i]{'transcript_name'},
							};

			#split the transcript name
			@decompose10 = split (/\_/,$table[$i]{'transcript_name'});

			#Store the assembler and line info into variable
			$assembler = $decompose10[0];
			$line = $decompose10[1];

			#implement the transcript frequency
			$$assembler[$line]{'transcript_freq'} ++;

			#implement the number of CDS with equal size to ref
			$nb_CDS_equal_size_to_ref ++;

			#implement $d variable
			$d ++;
		}
	}

	#to do if $d is equal to 0
	if ($d eq "0"){

		#implement the number of cluster without multiple CDS but not with length equal to the ref
		$nb_cluster_multiple_CDS_without_multiple_ref_size ++;

		#implement the number of transcript to pivk
		$nb_transcript_to_pick ++;

		#store the transcript into the transcript list to pick
		push @table_transcript_to_pick, {	transcript_name =>$name_ref,
									};

		#split the transcript name
		@decompose11 = split (/\_/,$name_ref);

		#Store the assembler and line info into variable
		$assembler = $decompose11[0];
		$line = $decompose11[1];

		#implement the transcript frequency
		$$assembler[$line]{'transcript_freq'} ++;
	}

	#to do if $d is superior to 0
	elsif ($d > "0"){

		#store the data in a table cluster specific
		push @$table_name, {	transcript_name =>$name_ref,
						};

		#split the transcript name
		@decompose12 = split (/\_/,$name_ref);

		#Store the assembler and line info into variable
		$assembler = $decompose12[0];
		$line = $decompose12[1];

		#implement the transcript frequency
		$$assembler[$line]{'transcript_freq'} ++;

		#implement the number of cluster with multiple CDS but not with length equal to the ref
		$nb_cluster_multiple_CDS_with_multiple_ref_size ++;

		#implement the number of CDS with equal size to ref
		$nb_CDS_equal_size_to_ref ++;

		#implement $e variable
		$e ++;

		#create new table name
		$table_name = "multiple_CDS".$e;
	}			
}

#close INFO file
close INFO;

#count the number of line in the transcript to pick table
$transcript_to_pick_table_size = $#table_transcript_to_pick + 1;

#print info lines
print LOG_OUT "There are $nb_cluster clusters and $nb_CDS CDS in the cluster info file\n\n";
print LOG_OUT "There are $nb_single_CDS_cluster clusters with single CDS in the cluster info file (category 1)\n";
print LOG_OUT "There are $nb_cluster_multiple_CDS_without_multiple_ref_size clusters with multiples CDS but no multiple CDS equal to the reference size (category 2)\n";
print LOG_OUT "There are $nb_transcript_to_pick transcripts to pick coming from these two categories (categories 1+2)\n\n";

print LOG_OUT "There are $nb_cluster_multiple_CDS_with_multiple_ref_size clusters with multiple CDS with size equal to the reference (category 3)\n";
print LOG_OUT "There are $nb_CDS_equal_size_to_ref CDS with a size equal to reference in the multiples CDS and multiple reference size clusters\n";

#print general computation time
print LOG_OUT "Cluster file treatment ended in ",time() - $starting_time_clustering_file_treatment," sec.\n\n";


#Storage of the CDS into a table per transcript
#Starting time for cds storage
my $starting_time_cds_storage = time();

#Variable
#########
my $nb_cds = 0;
my $cds_seq = "";
my $transcript_table_name = "";
my $cds_name = "";

#open the file containing the cds sequence
open CDS,"<$working/7-Transdecoder/$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta.transdecoder.cds" or die "cds file is not found during creation of files for annotation";

#Start to read the file
#######################

#file reading
while ($line4=<CDS>){

	#extract first character of each line
	$o = substr $line4,0,1;

	#to get the ORF name
	if ($o eq ">"){

		#implement sequence number variable
		$nb_cds ++;

		#if the name variable is not empty
		if ($cds_name ne ""){

			#split the transcript name
			@decompose13 = split (/\::/,$cds_name);

			#split the last part of the line
			@decompose16 = split (/\:/,$decompose13[5]);

			#replace the CDS name
			$cds_name = $decompose13[0]."::".$decompose13[1]."_c0_g1_i1::".$decompose13[2]."::".$decompose13[3]."::".$decompose13[4]."_c0_g1_i1::".$decompose16[0].":".$decompose16[1].":".$decompose16[2]."_c0_g1_i1:".$decompose16[3];

			#set the table name
			$transcript_table_name = $decompose13[1];

			#store the data in the table specific of the transcript
			push @$transcript_table_name, {	cds_name =>$cds_name,
							cds_seq =>$cds_seq,
								};

			#reset the sequence variable
			$cds_seq = "";		
		}

		#remove unusefull character and store the new contig name
		$line4 =~ s/\n//g;
		$line4 =~ s/\r//g;
		$line4 =~ s/\>//g;
		$cds_name = $line4;
	}

	#to get the sequence 
	elsif ($o ne ">"){

		#remove unusefull character
		$line4 =~ s/\n//g;
		$line4 =~ s/\r//g;
		$line4 =~ s/\ //g;

		#store the sequence into a variable
		$cds_seq = $cds_seq.$line4;
	}
}

#split the transcript name
@decompose13 = split (/\::/,$cds_name);

#split the last part of the line
@decompose16 = split (/\:/,$decompose13[5]);

#replace the CDS name
$cds_name = $decompose13[0]."::".$decompose13[1]."_c0_g1_i1::".$decompose13[2]."::".$decompose13[3]."::".$decompose13[4]."_c0_g1_i1::".$decompose16[0].":".$decompose16[1].":".$decompose16[2]."_c0_g1_i1:".$decompose16[3];

#set the table name
$transcript_table_name = $decompose13[1];

#store the data in the table specific of the transcript
push @$transcript_table_name, {	cds_name =>$cds_name,
				cds_seq =>$cds_seq,
					};

#reset the sequence variable
$cds_seq = "";

close CDS;

#print general computation time
print LOG_OUT "cds storage ended in ",time() - $starting_time_cds_storage," sec.\n\n";


#Storage of the peptide into a table per transcript
#Starting time for cds storage
my $starting_time_peptide_storage = time();

#Variable
#########
my $nb_peptide = 0;
my $peptide_seq = "";
my $transcript_table_name = "";
my $pep_name = "";

#open the file containing the cds sequence
open PEP,"<$working/7-Transdecoder/$$table1[0]{'sample_name'}_$liste_assembler"."merged_unigene.fasta.transdecoder.pep" or die "cds file is not found during creation of files for annotation";

#Start to read the file
#######################

#file reading
while ($line5=<PEP>){

	#extract first character of each line
	$p = substr $line5,0,1;

	#to get the ORF name
	if ($p eq ">"){

		#implement sequence number variable
		$nb_peptide ++;

		#if the name variable is not empty
		if ($pep_name ne ""){

			#split the transcript name
			@decompose14 = split (/\::/,$pep_name);

			#split the final segment of $decompose14 in order to also make an addition
			@decompose15 = split (/\:/,$decompose14[5]);

			#replace the CDS name
			$pep_name = $decompose14[0]."::".$decompose14[1]."_c0_g1_i1::".$decompose14[2]."::".$decompose14[3]."::".$decompose14[4]."_c0_g1_i1::".$decompose15[0].":".$decompose15[1].":".$decompose15[2]."_c0_g1_i1:".$decompose15[3];

			#set the table name
			$transcript_table_name = $decompose14[1]."pep";

			#store the data in the table specific of the transcript
			push @$transcript_table_name, {	pep_name =>$pep_name,
							pep_seq =>$pep_seq,
								};

			#reset the sequence variable
			$pep_seq = "";		
		}

		#remove unusefull character and store the new contig name
		$line5 =~ s/\n//g;
		$line5 =~ s/\r//g;
		$line5 =~ s/\>//g;
		$pep_name = $line5;
	}

	#to get the sequence 
	elsif ($p ne ">"){

		#remove unusefull character
		$line5 =~ s/\n//g;
		$line5 =~ s/\r//g;
		$line5 =~ s/\ //g;

		#store the sequence into a variable
		$pep_seq = $pep_seq.$line5;
	}
}

#split the transcript name
@decompose14 = split (/\::/,$pep_name);

#split the final segment of $decompose14 in order to also make an addition
@decompose15 = split (/\:/,$decompose14[5]);

#replace the CDS name
$pep_name = $decompose14[0]."::".$decompose14[1]."_c0_g1_i1::".$decompose14[2]."::".$decompose14[3]."::".$decompose14[4]."_c0_g1_i1::".$decompose15[0].":".$decompose15[1].":".$decompose15[2]."_c0_g1_i1:".$decompose15[3];

#set the table name
$transcript_table_name = $decompose14[1]."pep";

#store the data in the table specific of the transcript
push @$transcript_table_name, {	pep_name =>$pep_name,
				pep_seq =>$pep_seq,
					};

#reset the sequence variable
$pep_seq = "";	


close PEP;

#print general computation time
print LOG_OUT "peptide storage ended in ",time() - $starting_time_peptide_storage," sec.\n\n";

#Look into the list of each cluster in order to pick the longest transcript
######################################################################################################################################################################################################################################################################

#Starting time for cluster table treatment
my $starting_time_cluster_treatment = time();

#Variable
#########
my $k = 0;
my $assembler;
my $line;
my $cluster_ref_name = "";
my $cluster_ref_size = "";
my $cluster_ref_seq = "";
my $cluster_ref_freq = "";
my $nb_transcript_to_pick_from_cluster = 0;
my $transcript_to_pick_table_size = 0;


#while the $k variable is lower than $e
while ($k < $e){

	#creation of the table name
	$table_name = "multiple_CDS".$k;

	#read the table
	foreach my $l (0..$#$table_name){

		#to do if $l is equal to 0
		if ($l eq "0"){

			#split the transcript name
			@decompose3 = split (/\_/,$$table_name[$l]{'transcript_name'});

			#Store the assembler and line info into variable
			$assembler = $decompose3[0];
			$line = $decompose3[1];

			#store the info about the transcript into variable to compare with the other of the cluster
			$cluster_ref_name = $$assembler[$line]{'transcript_name'};
			$cluster_ref_size = $$assembler[$line]{'transcript_size'};
			$cluster_ref_freq = $$assembler[$line]{'transcript_freq'};
		}

		#to do if $l is superior to 0
		elsif ($l > "0"){

			#split the transcript name
			@decompose4 = split (/\_/,$$table_name[$l]{'transcript_name'});

			#Store the assembler and line info into variable
			$assembler = $decompose4[0];
			$line = $decompose4[1];

			#if the transcript freq is nul
			if ($$assembler[$line]{'transcript_freq'} eq "0"){

				#print an error message
				print "The frequency of the transcript $$assembler[$line]{'transcript_name'} is nul which is impossible\n";
			}

			#to do if the freq of the transcript is lower than the ref
			if ($$assembler[$line]{'transcript_freq'} < $cluster_ref_freq){

				#set ref variable with the new info
				$cluster_ref_name = $$assembler[$line]{'transcript_name'};
				$cluster_ref_size = $$assembler[$line]{'transcript_size'};
				$cluster_ref_freq = $$assembler[$line]{'transcript_freq'};
			}

			#to do if the freq of the transcript is equal than the ref
			elsif ($$assembler[$line]{'transcript_freq'} eq $cluster_ref_freq){

				#To do if the current transcript size is higher than the ref size
				if ($$assembler[$line]{'transcript_size'} > $cluster_ref_size){

					#set ref variable with the new info
					$cluster_ref_name = $$assembler[$line]{'transcript_name'};
					$cluster_ref_size = $$assembler[$line]{'transcript_size'};
					$cluster_ref_freq = $$assembler[$line]{'transcript_freq'};
				}

				#To do if the current transcript size is lower than the ref size
				elsif ($$assembler[$line]{'transcript_size'} < $cluster_ref_size){

				}
			}

			#to do if the freq of the transcript is higher than the ref
			elsif ($$assembler[$line]{'transcript_freq'} > $cluster_ref_freq){

			}
		}
	}

	#store the ref into the transcript list to pick
	push @table_transcript_to_pick, {	transcript_name =>$cluster_ref_name,
								};

	#implement the number of transcript to pick
	$nb_transcript_to_pick_from_cluster ++;

	#reset the ref variables
	$cluster_ref_name = "";
	$cluster_ref_size = "";
	$cluster_ref_freq = "";

	#implement $k variable
	$k ++;
}

#count the number of line in the transcript to pick table
$transcript_to_pick_table_size = $#table_transcript_to_pick + 1;

#print general computation time
print LOG_OUT "Treatment of the cluster table ended in ",time() - $starting_time_cluster_treatment," sec.\n";

print LOG_OUT "There are $nb_transcript_to_pick_from_cluster transcripts to pick coming from the cluster with multiple CDS with reference size treatment (category 3)\n\n";#$e\n\n$transcript_to_pick_table_size\n\n";

#Look into the list of transcript to pick in order to get the sequence
######################################################################################################################################################################################################################################################################

#variable
#########
my $cds_table_name = "";
my $pep_table_name = "";
my $nb_transcript_with_cds = 0;
my $nb_transcript_without_cds = 0;
my $nb_transcript_without_pep = 0;
my $nb_pick_CDS = 0;
my $nb_transcript_with_multiple_CDS = 0;

#open output file for transcript
open OUT,">$working/9-concatenated_assembly/$$table1[0]{'sample_name'}_$liste_assembler"."for_annotation.fasta" or die "The output file was not created";

#open output file for cds
open OUT_CDS,">$working/9-concatenated_assembly/$$table1[0]{'sample_name'}_$liste_assembler"."for_annotation"."_cds.fasta" or die "The output CDS file was not created";

#open output file for peptide
open OUT_PEP,">$working/9-concatenated_assembly/$$table1[0]{'sample_name'}_$liste_assembler"."for_annotation"."_pep.fasta" or die "The output peptide file was not created";

#Starting time for picking transcript
my $starting_time_transcript_picking = time();

#read the table
foreach my $m (0..$#table_transcript_to_pick){

	#split the transcript name
	@decompose4 = split (/\_/,$table_transcript_to_pick[$m]{'transcript_name'});

	#Store the assembler and line info into variable
	$assembler = $decompose4[0];
	$line = $decompose4[1];

	#to do if it is the first time that the transcript is meet
	if ($$assembler[$line]{'transcript_picking'} eq "0"){

		#implement the transcript_picking_frequency
		$$assembler[$line]{'transcript_picking'} ++;

		#set the cds_table_name
		$cds_table_name = $table_transcript_to_pick[$m]{'transcript_name'};

		#set the cds_table_name
		$pep_table_name = $table_transcript_to_pick[$m]{'transcript_name'}."pep";

		#to do if the CDS name of the transcript table is empty
		if ($$cds_table_name[0]{'cds_name'} eq ""){

			#implement the number of transcript without CDS
			$nb_transcript_without_cds ++;
		}

		#to do if the peptide name of the transcript table is empty
		if ($$pep_table_name[0]{'pep_name'} eq ""){

			#implement the number of transcript without CDS
			$nb_transcript_without_pep ++;
		}

		#to do if the first line of the transcript table is not empty
		elsif ($$pep_table_name[0]{'pep_name'} ne ""){

			#implement the number of transcript without CDS
			$nb_transcript_with_cds ++;

			#read the table
			foreach my $n (0..$#$pep_table_name){

				#count the number of CDS collected
				$nb_pick_CDS ++;

				#print the data in the CDS output file
				print OUT_CDS ">",$$cds_table_name[$n]{'cds_name'},"\n",$$cds_table_name[$n]{'cds_seq'},"\n";

				#print the data in the pep output file
				print OUT_PEP ">",$$pep_table_name[$n]{'pep_name'},"\n",$$pep_table_name[$n]{'pep_seq'},"\n";

				#if $n is equal to 1
				if ($n eq "1"){

					#implement the variable which count the number of transcript with multiple CDS
					$nb_transcript_with_multiple_CDS ++;
				}
			}
		}

		#print the name and sequence of the transcript of interest
		print OUT ">",$$assembler[$line]{'transcript_name'},"_c0_g1_i1\n",$$assembler[$line]{'transcript_sequence'},"\n";
	}
}

# if you want to clean the intermediate folders
if ($clean eq "1"){

	#remove all intermediate folder
	system "rm -R -f $working/2-* $working/3-* $working/4-* $working/5-* $working/6-* $working/7-* $working/8-*";
}


#close output file
close OUT;

#close output file
close OUT_PEP;
close OUT_CDS;

#control that there is not transcript without peptide
if ($nb_transcript_without_pep ne "0"){

	#if there is only one transcript without peptide
	if ($nb_transcript_without_pep eq "1"){

		#print a warning message if one transcript do not have peptide sequences
		print LOG_OUT "There is $nb_transcript_without_pep transcript without peptide sequence";
	}

	#if there is more than one transcript without peptide	
	elsif ($nb_transcript_without_pep > "1"){

		#print a warning message if one transcript do not have peptide sequences
		print LOG_OUT "There are $nb_transcript_without_pep transcripts without peptide sequence";
	}
}

#control that there is not transcript without peptide
if ($nb_transcript_without_cds ne "0"){

	#if there is only one transcript without peptide
	if ($nb_transcript_without_cds eq "1"){

		#print a warning message if one transcript do not have peptide sequences
		print LOG_OUT "There is $nb_transcript_without_cds transcript without cds sequence";
	}

	#if there is more than one transcript without peptide	
	elsif ($nb_transcript_without_cds > "1"){

		#print a warning message if one transcript do not have peptide sequences
		print LOG_OUT "There are $nb_transcript_without_cds transcripts without cds sequence";
	}
}

#if there is no transcript without peptide and CDS
if ($nb_transcript_without_pep eq "0" and $nb_transcript_without_cds eq "0"){

	#print general info
	print LOG_OUT "There are $nb_transcript_with_cds transcripts with at least one CDS in the concatenated assembly\n";
	print LOG_OUT "There are $nb_pick_CDS CDS in the concatenated assembly and $nb_transcript_with_multiple_CDS transcripts have multiple CDS\n\n";
	print LOG_OUT "Please note that some transcripts can have CDS in several categories (1, 2 and/or 3), so the number of transcripts in the concatenated assembly might be lower than the sum of the three categories\n";
}

#print general computation time
print LOG_OUT "Transcript picking ended in ",time() - $starting_time_transcript_picking," sec.\n";
