#!/usr/bin/perl -w
#-----------------------------------------------------------+
#                                                           | 
# Blast2Gb.pl                                               |
#                                                           |
#-----------------------------------------------------------+
# AUTHOR: Aurelien Chateigner                               |
# CONTACT: aurelien.chateigner_at_gmail.com                 |
# STARTED: 06/17/2011                                       |
# UPDATED:                                                  |
# LICENCE: CeCILL                                           |
# DESCRIPTION:                                              |
#  Converts Blast XML (5) output to Genebank format.        |
#                                                           |
#  Inspired by the James C. Estill's Blast2Gff.pl           |
#                                                           |
#-----------------------------------------------------------+

#-----------------------------+
# DATE 		                  |
#-----------------------------+
my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
$mday = sprintf("%02d", $mday);
my $dateUs =  ($mon+1)."/".$mday."/".(1900+$year);
my $MMM;
if ($mon+1 == 1){$MMM = JAN};
if ($mon+1 == 2){$MMM = FEB};
if ($mon+1 == 3){$MMM = MAR};
if ($mon+1 == 4){$MMM = APR};
if ($mon+1 == 5){$MMM = MAY};
if ($mon+1 == 6){$MMM = JUN};
if ($mon+1 == 7){$MMM = JUL};
if ($mon+1 == 8){$MMM = AUG};
if ($mon+1 == 9){$MMM = SEP};
if ($mon+1 == 10){$MMM = OCT};
if ($mon+1 == 11){$MMM = NOV};
if ($mon+1 == 12){$MMM = DEC};

#-----------------------------+
# INCLUDES                    |
#-----------------------------+
use strict;                   																		# Pragma to add restrictions to Perl rules
use Getopt::Long;              																		# To parse the command line
use POSIX;
use threads;																						# Mutithreading, Yeah !

#-----------------------------+
# VARIABLES                   |
#-----------------------------+

# Program variables
my $PROG_NAME = $0; 																				# $0 contains the name of the file containing the Perl script being executed  
my $VERSION   = "v3.1";  
my $VersionDate = "Decembre 16th 2011";
my $LastAdd = "works for blastn, tblastn, blastx, tblastx, in reverse and normal";

# User variables
my $InFile;                    																		# Full path to the blast-output input file
#my $OutFile;                   																	# Full path to the gb formatted output file
my $SourceFile;																						# Full path to the fasta input file
my $PrintHelp;                 																		# Boolean, print GBOUTthe Usage statement

# all infos needed in gb format
my $locus = "NOT_DEFINED";																			#ok
my $total_size;																						#ok
my $dna = "DNA";																					#ok
my $circular_or_linear = 'linear';																	#ok
my $vrl = "VRL";																					#ok	for Virus (see genbank format)
my $date = $mday."-".$MMM."-".(1900+$year);															#ok	format dd-MMM-yyyy
my $definition = "Undefined, please define in your blast fasta input";								#ok
my $accession = "AA000000";																			#ok
my $version = "AA000000.0";																			#ok
my $gi = "0000000";																					#ok
my $dblink = "NOT_DEFINED";																			#ok
my $keywords = ".";																					#ok
my $source = "Undefined";																			#ok
my $organism = "Undefined";																			#ok
my $reference = ".";																				#ok
my $authors = ".";																					#ok
my $title = ".";																					#ok
my $journal = ".";																					#ok
my $pubmed = ".";																					#ok
my $reference1 = ".";																				#ok
my $consrtm = ".";																					#ok
my $title2 = "."; 																					#ok
my $journal2 = ".";																					#ok 
my $reference2 = ".";																				#ok 
my $authors2 = ".";																					#ok 
my $title3 = ".";																					#ok 
my $journal3 = ".";																					#ok 
my $comment = ".";																					#ok
my $features = "Location/Qualifiers";																#ok
my $features_source;																				#ok
my $feature_source_organism;																		#ok
my $feature_source_mol_type = "genomic DNA";														#ok
my $feature_source_db_xref = ".";																	#ignoré
my $start;																							#ok
my $end;																							#ok
my $gene_position;																					#ok
my $gene_locus_tag = ".";																			#ok
my $gene_db_xref = ".";																				#ok
my $cds_position;																					#ok
my $cds_locus_tag = ".";																			#ok
my $cds_note = ".";																					#ok
my $cds_codon_start;																				#ok
my $cds_product;																					#ok
my $cds_protein_id;																					#ok
my $cds_protein_id_gi;																				#ok
my $cds_protein_id_geneid = ".";																	#ok
my $cds_translation;																				#ok
my $complete_sequence;																				#ok
my $hsp_num;																						#ok
my @p;																								#ok
my $j;																								#ok
my @count;																							#ok
my $z;																								#ok
my $count = 0;																						#ok
my $mol;																							#ok
my $test;																							#ok
my $hitdef;																							#ok
my $hit_format;																						#ok
my @hitdef;																							#ok
my $iteration_number = 0;																			#ok, number of fasta sequences blasted
my $fasta_count = 0;																				#ok, counts the number of fasta sequences, for the sequence print in the end
my $outf;																							#ok
my $count_fasta_61;																					#ok
my $query_def;																						#ok
my $target_start;
my $target_end;
my $color;
my $hit_count;
my ($first_def_line, $second_def_line);
my $gigbname;
my $query_gi;
my $query_accession;
my $query_name;
my $reverse;
my %scaffold_list = ();
my $hit_length;
my %seen = ();
my $print_all_hsp;
my $aftercoma;


# --------------------------------------------------+  
# SUB TO PRINT THE VERSION                          |
# --------------------------------------------------+  
sub Version {  
    print "
	$PROG_NAME version: $VERSION  
	Date  : $VersionDate  
	Author: Aurelien Chateigner
	Licence: CeCILL
	Last addition : $LastAdd\n\n";  
    exit;  
} 																									# end of Version  

sub PrintHelp {
	print "
	USAGE:
	perl Blast2Gb.pl -i InFile -s SourceFile [-c circular||linear] [-m 'molecular type'] [-r 0||1] [-a 0||1]
	
	-i:
	Your input has to be an xml blast output format (option 5).
	
	-r:
	It can be 0 for not reverse or 1 for reverse.
	If you want the reverse, the sourcefile is the db used for the blast.
	
	-a:
	It can be 0 for not all or 1 for all.
	It means that you want to print all the Hsp in a Hit (1), or just the best one (0).
	It is useful if you have exons (chose 1).
	
	-s:
	If your input has an odd name, please define your blast input fasta file this way :
	\">gi|0000000|gb|AA000000.0| Name of the organism, complete genome(or anything)\"
	And it works better if the sequence is on one line.
	
	EXAMPLE:
  	perl Blast2Gb.pl -i myinput.xml -s mydb.fasta -c circular -m 'genomic DNA' -r 0 -a 0
	\n";
	exit;
}


# Get the options
GetOptions('i=s' => \$InFile,
			's=s' => \$SourceFile,
#			'o=s' => \$OutFile,
			'c=s' => \$circular_or_linear,
			'h' => \&PrintHelp,
			'm=s' => \$feature_source_mol_type,
			'r=s' => \$reverse,
			'a=i' => \$print_all_hsp,
			"version" => \&Version);

# Open the different files
open (BLASTIN, "<", "$InFile") || die "Can not open BLAST input file.$InFile.\n";

open (SOURCE, "<", "$SourceFile") || die "Can not open Source fasta file.$SourceFile.\n";

#open (GBOUT, ">", "$OutFile") || die "Can not open GFF ouput file.$OutFile.\n";
if ($reverse eq "") {
	print "Please define -r\n";
	exit;
}
if ($print_all_hsp eq "") {
	print "Please define -a\n";
	exit;
}


#-------------------------------------------------+
# READ BLAST OUTPUT AND EXTRACT NEEDED FEATURES   |
#-------------------------------------------------+

(undef, $mol) = split (/ /, $feature_source_mol_type);
$feature_source_mol_type = "\/mol\_type\=\"$feature_source_mol_type\"";

while (<BLASTIN>) {
	if (/<BlastOutput_program>blastx/ || /<BlastOutput_program>tblastx/){
		while (<BLASTIN>) {
			if (/<BlastOutput_query-def>/){
				(undef, $test) = split (/\>/);
				($test, undef) = split (/\</, $test);
# 				print "$test\n";
				if ($test ne "No definition line"){
					(undef, $gi, undef, $version, $definition) = split (/\|/);
					($definition, undef) = split (/\</, $definition);
					$definition = substr $definition, 1;
					($locus, undef) = split (/\./, $version);
					($source, undef) = split(/\,/, $definition);
					$feature_source_organism = "\/organism\=\"$source\"";
					$organism = "$source\n            Viruses; dsDNA viruses, no RNA stage; Poxviridae; Entomopoxvirinae;\n            Betaentomopoxvirus.";
				}
			}
			elsif (/<Iteration_iter-num>/){
				$hit_count = 0;
			}
			elsif (/<Iteration_query-len>/){
				(undef, $total_size) = split (/\>/);
				($total_size, undef) = split (/\</, $total_size);
				$features_source = "1..$total_size";
			}
			elsif (/<Iteration_query-def>/){
				chomp ($_);
				(undef, $query_def) = split (/\>/);
				($query_def, undef) = split (/\</, $query_def);
				if (/\,/){
					($gigbname, $query_def) = split (/\,/, $query_def);
					$query_def = substr $query_def, 1;
					(undef, $query_gi, undef, $query_accession, $query_name) = split (/\|/, $gigbname);
					$query_name = substr $query_name, 1;
					$version = $query_accession;
				}
				else {
					(undef, $query_gi, undef, $query_accession, $query_name) = split (/\|/, $query_def);
					$query_name = substr $query_name, 1;
					$version = $query_accession;
				}
				$first_def_line = $query_name;
				$second_def_line = "";
				if ($query_name =~ / /){
				($first_def_line, $second_def_line) = split (/ /, $query_name);
				}
				$second_def_line =~ s/ //;
				$query_def = "$first_def_line\n            $second_def_line";
			}
			elsif (/<Iteration_hits>/){
				if ($reverse != 1){
					$iteration_number += 1;
				}
			}
			elsif (/<Hit_num>/){
				$hit_count ++;
			}
			elsif (/<Hit_def>/){
				if (/<Hit_def>RecName/){
					(undef, $hitdef) = split (/\>/);
					($hitdef, undef) = split (/\</, $hitdef);
					$hit_format = 2;
				}
				elsif (/\&gt\;gi/){
					@hitdef = split (/\&gt\;/);
					$hitdef = $hitdef[@hitdef-1];
					chomp ($hitdef);
					(undef, $cds_protein_id_gi, undef, $cds_protein_id, $cds_product) = split (/\|/, $hitdef);
					($cds_product, undef) = split (/\</, $cds_product);
					$cds_product =~ s/^\s+//;
					$hit_format = 1;
				}
				elsif (/\>gi/){
					(undef, $cds_protein_id_gi, undef, $cds_protein_id, $cds_product) = split (/\|/);
					($cds_product, undef) = split (/\</, $cds_product);
					$cds_product = substr $cds_product, 1;
					($cds_product, $hitdef) = split (/\,/, $cds_product);
					$hitdef =~ s/ //;
					$hit_format = 1;
				}
				else {
					$hitdef = $_;
					(undef, $hitdef) = split (/\>/, $hitdef);
					($hitdef, undef) = split (/\</, $hitdef);
					$hit_format = 2;
				}
				if ($reverse == 1){					
					if ($scaffold_list{$hitdef}){
						my $iteration_file_number = $scaffold_list{$hitdef};
						$outf = "$InFile"."$iteration_file_number"."\.gb";
						print "$hitdef allreardy in the list, printing in file $outf\n";
					}
					else {
						$iteration_number += 1;
						$scaffold_list{$hitdef} = $iteration_number;
						$outf = "$InFile"."$iteration_number"."\.gb";
						print "$hitdef added to the list, new file created: $outf\n";
					}
				}
			}
			elsif (/<Hit_len>/){
				(undef, $hit_length) = split (/\>/);
				($hit_length, undef) = split (/\</, $hit_length);

				if ($reverse != 1){
					if (! $seen{$iteration_number}){
						$seen{$iteration_number} = 1;
						$outf = "$InFile"."$iteration_number"."\.gb";
						print "infos and CDS printed in $outf\n";
				
						open (GBOUT, ">", "$outf") || die "Can not open GB ouput file.$outf.\n";
				
						printf GBOUT "%-12s", "LOCUS";
						print GBOUT "$gi                $total_size bp    $mol     $circular_or_linear   VRL $date\n";
						printf GBOUT "%-12s", "DEFINITION";
		 				print GBOUT "$source, $query_def\.\n";
						printf GBOUT "%-12s", "ACCESSION";
						($accession, undef) = split (/\./, $version);
						print GBOUT "$accession\n";
						printf GBOUT "%-12s", "VERSION";
						print GBOUT "$version  GI:$gi\n";
						printf GBOUT "%-12s", "KEYWORDS";
						print GBOUT "$keywords\n";
						printf GBOUT "%-12s", "SOURCE";
						print GBOUT "$source\n";
						printf GBOUT "%-12s", "  ORGANISM";
						print GBOUT "$organism\n";
						printf GBOUT "%-12s", "REFERENCE";
						print GBOUT "$reference\n";
						printf GBOUT "%-12s", "  AUTHORS";
						print GBOUT "$authors\n";
						printf GBOUT "%-12s", "  TITLE";
						print GBOUT "$title\n";
						printf GBOUT "%-12s", "  JOURNAL";
						print GBOUT "$journal\n";
						printf GBOUT "%-12s", "  PUBMED";
						print GBOUT "$pubmed\n";
						printf GBOUT "%-12s", "COMMENT";
						print GBOUT "This is not a record yet, it has just been created by Blast2Gb.pl.\n";
						printf GBOUT "%-12s", "FEATURES";
						print GBOUT "         Location\/Qualifiers\n";
						printf GBOUT "%-12s", "     source";
						print GBOUT "         1\.\.$total_size\n";
						print GBOUT "                     \/organism=\"$source\"\n";
						print GBOUT "                     $feature_source_mol_type\n";
						print GBOUT "                     \/label=\"source $source\"\n";
						if ($feature_source_db_xref ne "."){
							print GBOUT "                    \/db_xref\=\"$feature_source_db_xref\"\n";
						}
					}
				}
				elsif ($reverse == 1){
					if ($scaffold_list{$hitdef} && ! $seen{$hitdef}){
						$seen{$hitdef} = 1;
						open (GBOUT, ">", "$outf") || die "Can not open GB ouput file.$outf.\n";
				
						printf GBOUT "%-12s", "LOCUS";
						print GBOUT "$cds_protein_id_gi                $hit_length bp    $mol     $circular_or_linear   VRL $date\n";
						printf GBOUT "%-12s", "DEFINITION";
						print GBOUT "$cds_product\.\n";
						printf GBOUT "%-12s", "ACCESSION";
						(my $cds_accession, undef) = split (/\./, $cds_protein_id);
						print GBOUT "$cds_accession\n";
						printf GBOUT "%-12s", "VERSION";
						print GBOUT "$cds_protein_id  GI:$cds_protein_id_gi\n";
						printf GBOUT "%-12s", "KEYWORDS";
						print GBOUT "$keywords\n";
						printf GBOUT "%-12s", "SOURCE";
						($source, undef) = split (/\,/, $cds_product);
						print GBOUT "$source\n";
						printf GBOUT "%-12s", "  ORGANISM";
						$organism = "$source\n            Viruses; dsDNA viruses, no RNA stage; Poxviridae; Entomopoxvirinae;\n            Betaentomopoxvirus.";
						print GBOUT "$organism\n";
						printf GBOUT "%-12s", "REFERENCE";
						print GBOUT "$reference\n";
						printf GBOUT "%-12s", "  AUTHORS";
						print GBOUT "$authors\n";
						printf GBOUT "%-12s", "  TITLE";
						print GBOUT "$title\n";
						printf GBOUT "%-12s", "  JOURNAL";
						print GBOUT "$journal\n";
						printf GBOUT "%-12s", "  PUBMED";
						print GBOUT "$pubmed\n";
						printf GBOUT "%-12s", "COMMENT";
						print GBOUT "This is not a record yet, it has just been created by Blast2Gb.pl.\n";
						printf GBOUT "%-12s", "FEATURES";
						print GBOUT "         Location\/Qualifiers\n";
						printf GBOUT "%-12s", "     source";
						print GBOUT "         1\.\.$hit_length\n";
						print GBOUT "                     \/organism=\"$query_name\"\n";
						print GBOUT "                     $feature_source_mol_type\n";
						print GBOUT "                     \/label=\"source $query_name\"\n";
						print "infos and CDS printed in $outf\n";
						if ($feature_source_db_xref ne "."){
							print GBOUT "                    \/db_xref\=\"$feature_source_db_xref\"\n";
						}
					}
				}
			}
			elsif (/<Hsp_num>/){
				(undef, $hsp_num) = split (/\>/);
				($hsp_num, undef) = split (/\</, $hsp_num);
			}

			elsif (/<Hsp_query-from>/){
				(undef, $start) = split (/\>/);
				($start, undef) = split (/\</, $start);
			}
			elsif (/<Hsp_query-to>/){
				(undef, $end) = split (/\>/);
				($end, undef) = split (/\</, $end);
			}
			elsif (/<Hsp_hit-from>/){
				(undef,$target_start) = split (/\>/);
				($target_start, undef) = split (/\</, $target_start);
			}
			elsif (/<Hsp_hit-to>/){
				(undef, $target_end) = split (/\>/);
				($target_end, undef) = split (/\</, $target_end);
				if ($reverse == 1){
					($start, $end, $target_start, $target_end) = ($target_start, $target_end, $start, $end);
				}	
				$gene_position = "$start\.\.$end";
				if ($start >= $end){
				$gene_position = "complement($end\.\.$start)";
				}
				$cds_position = $gene_position;
			}
			elsif (/<Hsp_query-frame>/){
				(undef, $cds_codon_start) = split (/\>/);
				($cds_codon_start, undef) = split (/\</, $cds_codon_start);
				if ($cds_codon_start < 0){
					$gene_position = "complement($gene_position)";
					$gene_position =~ s/complement\(complement/complement/;
					$gene_position =~ s/\)\)/\)/;
				}
			}
			elsif (/<Hsp_hseq>/){
				if ($reverse != 1){
					$outf = "$InFile$iteration_number\.gb";
					open (GBOUT, ">>"."$outf") || die "Can not open GB ouput file.$outf\n";
					if ($print_all_hsp == 0 && $hsp_num == 1) {
						(undef, $cds_translation) = split (/\>/);
						($cds_translation, undef) = split (/\</, $cds_translation);
						$cds_translation =~ s/\W//;
						$cds_translation = "\/translation\=\"$cds_translation\"";
						my @cds_translation = split (/(.{58})/, $cds_translation);
						printf GBOUT "%-12s", "     CDS";
						print GBOUT "         $gene_position\n";
						print GBOUT "                     \/locus_tag=\"$gene_locus_tag\"\n";
						print GBOUT "                     \/db_xref=\"$gene_db_xref\"\n";
						print GBOUT "                     \/locus_tag=\"$cds_locus_tag\"\n";
						print GBOUT "                     \/note=\"$cds_note\"\n";
						print GBOUT "                     \/codon_start=\"$cds_codon_start\"\n";
						if ($hit_format == 1){
							print GBOUT "                     \/product=\"$cds_product\_$hitdef\_$target_start\_$target_end\"\n";
							print GBOUT "                     \/protein_id=\"$cds_protein_id\"\n";
							print GBOUT "                     \/db_xref=\"GI:$cds_protein_id_gi\"\n";
							print GBOUT "                     \/label=\"hit$hit_count\_hsp$hsp_num\"\n";
						}
						if ($hit_format == 2){
							print GBOUT "                     \/product=\"$hitdef\"\n";
							print GBOUT "                     \/protein_id=\".\"\n";
							print GBOUT "                     \/db_xref=\"GI:.\"\n";
							print GBOUT "                     \/label=\"hit$hit_count\_hsp$hsp_num\"\n";
						}
						print GBOUT "                     \/db_xref=\"GeneID:$cds_protein_id_geneid\"\n";
						foreach my $cds (@cds_translation){
							if ($cds !~ m#^$#){
								chomp ($cds);
								print GBOUT "                     $cds\n";
							}
						}
						print "Hit printed in $outf\n";
					}
					elsif ($print_all_hsp == 1) {
						(undef, $cds_translation) = split (/\>/);
						($cds_translation, undef) = split (/\</, $cds_translation);
						$cds_translation =~ s/\W//;
						$cds_translation = "\/translation\=\"$cds_translation\"";
						my @cds_translation = split (/(.{58})/, $cds_translation);
						printf GBOUT "%-12s", "     CDS";
						print GBOUT "         $gene_position\n";
						print GBOUT "                     \/locus_tag=\"$gene_locus_tag\"\n";
						print GBOUT "                     \/db_xref=\"$gene_db_xref\"\n";
						print GBOUT "                     \/locus_tag=\"$cds_locus_tag\"\n";
						print GBOUT "                     \/note=\"$cds_note\"\n";
						print GBOUT "                     \/codon_start=\"$cds_codon_start\"\n";
						if ($hit_format == 1){
							print GBOUT "                     \/product=\"$cds_product\_$hitdef\_$target_start\_$target_end\"\n";
							print GBOUT "                     \/protein_id=\"$cds_protein_id\"\n";
							print GBOUT "                     \/db_xref=\"GI:$cds_protein_id_gi\"\n";
							print GBOUT "                     \/label=\"hit$hit_count\_hsp$hsp_num\"\n";
						}
						if ($hit_format == 2){
							print GBOUT "                     \/product=\"$hitdef\"\n";
							print GBOUT "                     \/protein_id=\".\"\n";
							print GBOUT "                     \/db_xref=\"GI:.\"\n";
							print GBOUT "                     \/label=\"hit$hit_count\_hsp$hsp_num\"\n";
						}
						print GBOUT "                     \/db_xref=\"GeneID:$cds_protein_id_geneid\"\n";
						foreach my $cds (@cds_translation){
							if ($cds !~ m#^$#){
								chomp ($cds);
								print GBOUT "                     $cds\n";
							}
						}
						print "Hit printed in $outf\n";
					}
					
				}
			}
			elsif (/<Hsp_qseq>/){
				if ($reverse == 1){
					open (GBOUT, ">>"."$outf") || die "Can not open GB ouput file.$outf\n";
					if ($print_all_hsp == 0 && $hsp_num == 1){
						(undef, $cds_translation) = split (/\>/);
						($cds_translation, undef) = split (/\</, $cds_translation);
						$cds_translation =~ s/\W//;
						$cds_translation = "\/translation\=\"$cds_translation\"";
						my @cds_translation = split (/(.{58})/, $cds_translation);
						printf GBOUT "%-12s", "     CDS";
						print GBOUT "         $gene_position\n";
						print GBOUT "                     \/locus_tag=\"$gene_locus_tag\"\n";
						print GBOUT "                     \/db_xref=\"$gene_db_xref\"\n";
						print GBOUT "                     \/locus_tag=\"$cds_locus_tag\"\n";
						print GBOUT "                     \/note=\"$cds_note\"\n";
						print GBOUT "                     \/codon_start=\"$cds_codon_start\"\n";
						if ($hit_format == 1){
							print GBOUT "                     \/product=\"$cds_product\_$hitdef\_$target_start\_$target_end\"\n";
							print GBOUT "                     \/protein_id=\"$cds_protein_id\"\n";
							print GBOUT "                     \/db_xref=\"GI:$cds_protein_id_gi\"\n";
							print GBOUT "                     \/label=\"hit$hit_count\_hsp$hsp_num\"\n";
						}
						if ($hit_format == 2){
							print GBOUT "                     \/product=\"$hitdef\"\n";
							print GBOUT "                     \/protein_id=\".\"\n";
							print GBOUT "                     \/db_xref=\"GI:.\"\n";
							print GBOUT "                     \/label=\"hit$hit_count\_hsp$hsp_num\"\n";
						}
						print GBOUT "                     \/db_xref=\"GeneID:$cds_protein_id_geneid\"\n";
						foreach my $cds (@cds_translation){
							if ($cds !~ m#^$#){
								chomp ($cds);
								print GBOUT "                     $cds\n";
							}
						}
						print "Hit printed in $outf\n";
					}
					elsif ($print_all_hsp == 1){
						(undef, $cds_translation) = split (/\>/);
						($cds_translation, undef) = split (/\</, $cds_translation);
						$cds_translation =~ s/\W//;
						$cds_translation = "\/translation\=\"$cds_translation\"";
						my @cds_translation = split (/(.{58})/, $cds_translation);
						printf GBOUT "%-12s", "     CDS";
						print GBOUT "         $gene_position\n";
						print GBOUT "                     \/locus_tag=\"$gene_locus_tag\"\n";
						print GBOUT "                     \/db_xref=\"$gene_db_xref\"\n";
						print GBOUT "                     \/locus_tag=\"$cds_locus_tag\"\n";
						print GBOUT "                     \/note=\"$cds_note\"\n";
						print GBOUT "                     \/codon_start=\"$cds_codon_start\"\n";
						if ($hit_format == 1){
							print GBOUT "                     \/product=\"$cds_product\_$hitdef\_$target_start\_$target_end\"\n";
							print GBOUT "                     \/protein_id=\"$cds_protein_id\"\n";
							print GBOUT "                     \/db_xref=\"GI:$cds_protein_id_gi\"\n";
							print GBOUT "                     \/label=\"hit$hit_count\_hsp$hsp_num\"\n";
						}
						if ($hit_format == 2){
							print GBOUT "                     \/product=\"$hitdef\"\n";
							print GBOUT "                     \/protein_id=\".\"\n";
							print GBOUT "                     \/db_xref=\"GI:.\"\n";
							print GBOUT "                     \/label=\"hit$hit_count\_hsp$hsp_num\"\n";
						}
						print GBOUT "                     \/db_xref=\"GeneID:$cds_protein_id_geneid\"\n";
						foreach my $cds (@cds_translation){
							if ($cds !~ m#^$#){
								chomp ($cds);
								print GBOUT "                     $cds\n";
							}
						}
						print "Hit printed in $outf\n";
					}
				}
			}
		}
	}

	# if it is a blastn output -------------------------------------------------------------------------------------------------------------------------
	elsif (/<BlastOutput_program>blastn/ || /<BlastOutput_program>tblastn/){
		while (<BLASTIN>) {
			if (/<BlastOutput_query-def>/){
				(undef, $test) = split (/\>/);
				($test, undef) = split (/\</, $test);
# 				print "$test\n";
				if ($test ne "No definition line"){
					(undef, $gi, undef, $version, $definition) = split (/\|/);
					($definition, undef) = split (/\</, $definition);
					$definition = substr $definition, 1;
					($locus, undef) = split (/\./, $version);
					($source, undef) = split(/\,/, $definition);
					$feature_source_organism = "\/organism\=\"$source\"";
					$organism = "$source\n            Viruses; dsDNA viruses, no RNA stage; Poxviridae; Entomopoxvirinae;\n            Betaentomopoxvirus.";
				}
			}
			elsif (/<Iteration_iter-num>/){
				$hit_count = 0;
			}
			elsif (/<Iteration_query-len>/){
				(undef, $total_size) = split (/\>/);
				($total_size, undef) = split (/\</, $total_size);
				$features_source = "1..$total_size";
			}
			elsif (/<Iteration_query-def>/){
				chomp ($_);
				(undef, $query_def) = split (/\>/);
				($query_def, undef) = split (/\</, $query_def);
				if (/\,/){
					($gigbname, $query_def) = split (/\,/, $query_def);
					$query_def = substr $query_def, 1;
					$aftercoma = $query_def;
					(undef, $query_gi, undef, $query_accession, $query_name) = split (/\|/, $gigbname);
					$query_name = substr $query_name, 1;
					$version = $query_accession;
				}
				else {
					(undef, $query_gi, undef, $query_accession, $query_name) = split (/\|/, $query_def);
					$query_name = substr $query_name, 1;
					$version = $query_accession;
				}
				$first_def_line = $query_name;
				$second_def_line = "";
				if ($query_name =~ / /){
				($first_def_line, $second_def_line) = split (/ /, $query_name);
				}
				$second_def_line =~ s/ //;
				$query_def = "$first_def_line\n            $second_def_line";
			}
			elsif (/<Iteration_hits>/){
				if ($reverse != 1){
					$iteration_number += 1;
				}
			}
			elsif (/<Hit_num>/){
				$hit_count ++;
			}
			elsif (/<Hit_def>/){
				if (/<Hit_def>RecName/){
					(undef, $hitdef) = split (/\>/);
					($hitdef, undef) = split (/\</, $hitdef);
					$hit_format = 2;
				}
				elsif (/\&gt\;gi/){
					@hitdef = split (/\&gt\;/);
					$hitdef = $hitdef[@hitdef-1];
					chomp ($hitdef);
					(undef, $cds_protein_id_gi, undef, $cds_protein_id, $cds_product) = split (/\|/, $hitdef);
					($cds_product, undef) = split (/\</, $cds_product);
					$cds_product =~ s/^\s+//;
					$hit_format = 1;
				}
				elsif (/\>gi/){
					(undef, $cds_protein_id_gi, undef, $cds_protein_id, $cds_product) = split (/\|/);
					($cds_product, undef) = split (/\</, $cds_product);
					$cds_product = substr $cds_product, 1;
					($cds_product, $hitdef) = split (/\,/, $cds_product);
					$hitdef =~ s/ //;
					$hit_format = 1;
				}
				else {
					$hitdef = $_;
					(undef, $hitdef) = split (/\>/, $hitdef);
					($hitdef, undef) = split (/\</, $hitdef);
					$hit_format = 2;
				}


				if ($reverse == 1){					
					if ($scaffold_list{$hitdef}){
						my $iteration_file_number = $scaffold_list{$hitdef};
						$outf = "$InFile"."$iteration_file_number"."\.gb";
						print "$hitdef allreardy in the list, printing in file $outf\n";
					}
					else {
						$iteration_number += 1;
						$scaffold_list{$hitdef} = $iteration_number;
						$outf = "$InFile"."$iteration_number"."\.gb";
						print "$hitdef added to the list, new file created: $outf\n";
					}
				}
			}
			elsif (/<Hit_len>/){
				(undef, $hit_length) = split (/\>/);
				($hit_length, undef) = split (/\</, $hit_length);

				if ($reverse != 1){
					if (! $seen{$iteration_number}){
						$seen{$iteration_number} = 1;
						$outf = "$InFile"."$iteration_number"."\.gb";
						print "infos and CDS printed in $outf\n";
				
						open (GBOUT, ">", "$outf") || die "Can not open GB ouput file.$outf.\n";
				
						printf GBOUT "%-12s", "LOCUS";
						print GBOUT "$query_gi                $total_size bp    $mol     $circular_or_linear   VRL $date\n";
						printf GBOUT "%-12s", "DEFINITION";
						print GBOUT "$query_name\.\n";
						printf GBOUT "%-12s", "ACCESSION";
						($accession, undef) = split (/\./, $version);
						print GBOUT "$accession\n";
						printf GBOUT "%-12s", "VERSION";
						print GBOUT "$version  GI:$query_gi\n";
						printf GBOUT "%-12s", "KEYWORDS";
						print GBOUT "$keywords\n";
						printf GBOUT "%-12s", "SOURCE";
						print GBOUT "$query_name\n";
						printf GBOUT "%-12s", "  ORGANISM";
						print GBOUT "$query_name\n";
						printf GBOUT "%-12s", "REFERENCE";
						print GBOUT "$reference\n";
						printf GBOUT "%-12s", "  AUTHORS";
						print GBOUT "$authors\n";
						printf GBOUT "%-12s", "  TITLE";
						print GBOUT "$title\n";
						printf GBOUT "%-12s", "  JOURNAL";
						print GBOUT "$journal\n";
						printf GBOUT "%-12s", "  PUBMED";
						print GBOUT "$pubmed\n";
						printf GBOUT "%-12s", "COMMENT";
						print GBOUT "This is not a record yet, it has just been created by Blast2Gb.pl.\n";
						printf GBOUT "%-12s", "FEATURES";
						print GBOUT "         Location\/Qualifiers\n";
						printf GBOUT "%-12s", "     source";
						print GBOUT "         1\.\.$total_size\n";
						print GBOUT "                     \/organism=\"$cds_product\"\n";
						print GBOUT "                     $feature_source_mol_type\n";
						print GBOUT "                     \/label=\"source $cds_product\"\n";
						if ($feature_source_db_xref ne "."){
							print GBOUT "                    \/db_xref\=\"$feature_source_db_xref\"\n";
						}
					}
				}
				elsif ($reverse == 1){
					if ($scaffold_list{$hitdef} && ! $seen{$hitdef}){
						$seen{$hitdef} = 1;
						open (GBOUT, ">", "$outf") || die "Can not open GB ouput file.$outf.\n";
				
						printf GBOUT "%-12s", "LOCUS";
						print GBOUT "$cds_protein_id_gi                $hit_length bp    $mol     $circular_or_linear   VRL $date\n";
						printf GBOUT "%-12s", "DEFINITION";
						print GBOUT "$cds_product\.\n";
						printf GBOUT "%-12s", "ACCESSION";
						(my $cds_accession, undef) = split (/\./, $cds_protein_id);
						print GBOUT "$cds_accession\n";
						printf GBOUT "%-12s", "VERSION";
						print GBOUT "$cds_protein_id  GI:$cds_protein_id_gi\n";
						printf GBOUT "%-12s", "KEYWORDS";
						print GBOUT "$keywords\n";
						printf GBOUT "%-12s", "SOURCE";
						($source, undef) = split (/\,/, $cds_product);
						print GBOUT "$source\n";
						printf GBOUT "%-12s", "  ORGANISM";
						$organism = "$source\n            Viruses; dsDNA viruses, no RNA stage; Poxviridae; Entomopoxvirinae;\n            Betaentomopoxvirus.";
						print GBOUT "$organism\n";
						printf GBOUT "%-12s", "REFERENCE";
						print GBOUT "$reference\n";
						printf GBOUT "%-12s", "  AUTHORS";
						print GBOUT "$authors\n";
						printf GBOUT "%-12s", "  TITLE";
						print GBOUT "$title\n";
						printf GBOUT "%-12s", "  JOURNAL";
						print GBOUT "$journal\n";
						printf GBOUT "%-12s", "  PUBMED";
						print GBOUT "$pubmed\n";
						printf GBOUT "%-12s", "COMMENT";
						print GBOUT "This is not a record yet, it has just been created by Blast2Gb.pl.\n";
						printf GBOUT "%-12s", "FEATURES";
						print GBOUT "         Location\/Qualifiers\n";
						printf GBOUT "%-12s", "     source";
						print GBOUT "         1\.\.$hit_length\n";
						print GBOUT "                     \/organism=\"$query_name\"\n";
						print GBOUT "                     $feature_source_mol_type\n";
						print GBOUT "                     \/label=\"source $query_name\"\n";
						print "infos and CDS printed in $outf\n";
						if ($feature_source_db_xref ne "."){
							print GBOUT "                    \/db_xref\=\"$feature_source_db_xref\"\n";
						}
					}
				}
			}
			elsif (/<Hsp_num>/){
				(undef, $hsp_num) = split (/\>/);
				($hsp_num, undef) = split (/\</, $hsp_num);
			}
			elsif (/<Hsp_query-from>/){
				(undef, $start) = split (/\>/);
				($start, undef) = split (/\</, $start);
			}
			elsif (/<Hsp_query-to>/){
				(undef, $end) = split (/\>/);
				($end, undef) = split (/\</, $end);
			}
			elsif (/<Hsp_hit-from>/){
				(undef,$target_start) = split (/\>/);
				($target_start, undef) = split (/\</, $target_start);
			}
			elsif (/<Hsp_hit-to>/){
				(undef, $target_end) = split (/\>/);
				($target_end, undef) = split (/\</, $target_end);
				if ($reverse == 1){
					($start, $end, $target_start, $target_end) = ($target_start, $target_end, $start, $end);
				}
				$gene_position = "$start\.\.$end";
				if ($start >= $end){
				$gene_position = "complement($end\.\.$start)";
				}
				$cds_position = $gene_position;
			}
			elsif (/<Hsp_query-frame>/){
				(undef, $cds_codon_start) = split (/\>/);
				($cds_codon_start, undef) = split (/\</, $cds_codon_start);
				if ($cds_codon_start < 0){
					$gene_position = "complement($gene_position)";
					$gene_position =~ s/complement\(complement/complement/;
					$gene_position =~ s/\)\)/\)/;
				}
			}
			elsif (/<Hsp_hseq>/){
				if ($reverse != 1){
					$outf = "$InFile$iteration_number\.gb";
					open (GBOUT, ">>"."$outf") || die "Can not open GB ouput file.$outf\n";
					if ($print_all_hsp == 0 && $hsp_num == 1){
						(undef, $cds_translation) = split (/\>/);
						($cds_translation, undef) = split (/\</, $cds_translation);
						$cds_translation =~ s/\W//;
						$cds_translation = "\/translation\=\"$cds_translation\"";
						my @cds_translation = split (/(.{58})/, $cds_translation);
						printf GBOUT "%-12s", "     BLASTN_HIT";
						print GBOUT "      $gene_position\n";
						$cds_note = "$hitdef\_$target_start\_$target_end";
						print GBOUT "                     \/note=\"$cds_note\"\n";
						print GBOUT "                     \/color=4\n";								# Here you change the color !!! ==============================================
						print GBOUT "                     \/label=\"hit$hit_count hsp$hsp_num $hitdef BLASTN HIT\"\n";
						print "Hit printed in $outf\n";
					}
					elsif ($print_all_hsp == 1){
						(undef, $cds_translation) = split (/\>/);
						($cds_translation, undef) = split (/\</, $cds_translation);
						$cds_translation =~ s/\W//;
						$cds_translation = "\/translation\=\"$cds_translation\"";
						my @cds_translation = split (/(.{58})/, $cds_translation);
						printf GBOUT "%-12s", "     BLASTN_HIT";
						print GBOUT "      $gene_position\n";
						$cds_note = "$hitdef\_$target_start\_$target_end";
						print GBOUT "                     \/note=\"$cds_note\"\n";
						print GBOUT "                     \/color=4\n";								# Here you change the color !!! ==============================================
						print GBOUT "                     \/label=\"hit$hit_count hsp$hsp_num $hitdef BLASTN HIT\"\n";
						print "Hit printed in $outf\n";
					}
				}
			}
			elsif (/<Hsp_qseq>/){
				if ($reverse == 1){
					open (GBOUT, ">>"."$outf") || die "Can not open GB ouput file.$outf\n";
					if ($print_all_hsp == 0 && $hsp_num == 1){
						(undef, $cds_translation) = split (/\>/);
						($cds_translation, undef) = split (/\</, $cds_translation);
						$cds_translation =~ s/\W//;
						$cds_translation = "\/translation\=\"$cds_translation\"";
						my @cds_translation = split (/(.{58})/, $cds_translation);
						printf GBOUT "%-12s", "     BLASTN_HIT";
						print GBOUT "      $gene_position\n";
	
						$cds_note = "$query_name\_$target_start\_$target_end";# 						}
						print GBOUT "                     \/note=\"$cds_note\"\n";
						print GBOUT "                     \/color=4\n";								# Here you change the color !!! ==============================================
						print GBOUT "                     \/label=\"hit$hit_count hsp$hsp_num $query_name" . ", " . "$aftercoma" . " BLASTN HIT\"\n";
						print "Hit printed in $outf\n";
					}
					elsif ($print_all_hsp == 1){
						(undef, $cds_translation) = split (/\>/);
						($cds_translation, undef) = split (/\</, $cds_translation);
						$cds_translation =~ s/\W//;
						$cds_translation = "\/translation\=\"$cds_translation\"";
						my @cds_translation = split (/(.{58})/, $cds_translation);
						printf GBOUT "%-12s", "     BLASTN_HIT";
						print GBOUT "      $gene_position\n";
	
						$cds_note = "$query_name\_$target_start\_$target_end";# 						}
						print GBOUT "                     \/note=\"$cds_note\"\n";
						print GBOUT "                     \/color=4\n";								# Here you change the color !!! ==============================================
						print GBOUT "                     \/label=\"hit$hit_count hsp$hsp_num $query_name" . ", " . "$aftercoma" . " BLASTN HIT\"\n";
						print "Hit printed in $outf\n";
					}
				}
			}
		}
	}
}

# --------------------------+
# Do what Gbrenamer.pl does |
# --------------------------+

my %versions;
for (my $i = 1 ; $i <= $iteration_number ; $i++) {
# print "$InFile\n";
	my $version_to_mv;
	my $name;
	my $name_to_mv = "$InFile"."$i"."\.gb";

	open (BLASTIN, "<", "$name_to_mv") || die "Can not open file to rename file $name_to_mv\.\n";

	while (<BLASTIN>) {
	    if (/VERSION/) {
	        (undef, $version_to_mv) = split (/\s+/);
	        $versions{$version_to_mv} = $i;
	    }
	}
	$name = $InFile;
	$name =~ s/\.xml//;
	$name =~ s/\d\.gb//;
	my $file_mved = "$name\_$version_to_mv\.gb";
	system ("bash -c 'mv $name_to_mv $file_mved'");

	print "file $i renamed\n";
}


#----------------------------------------------------+
# READ SOURCE FILE, EXTRACT SEQUENCE AND FORMAT IT   |
#----------------------------------------------------+
my @list;
my $current_version;
$InFile =~ s/\.xml//;

if ($reverse != 1) {
	while (<SOURCE>) {
		chomp ($_);
		if (/^>/){
			(undef, undef, undef, $current_version, undef) = split (/\|/);
			push @list, $current_version;
			$count_fasta_61 = 0;
			$count = 0;
			$fasta_count ++;
			my $outf = "$InFile"."_"."$current_version"."\.gb";
			open (GBOUT, ">>", "$outf") || die "Can not open GFF ouput file.$outf\n";
			print GBOUT "ORIGIN      \n";
			print "sequence printed in $outf\n";
		}
		elsif (/^[A-Za-z]/) {
			my $lineseq = $_;
			$lineseq =~ s/^\s+//;
			$complete_sequence = $lineseq;
			$complete_sequence =~ tr/A-Z/a-z/;
			chomp($complete_sequence);
			my $length = length($complete_sequence);
			if ($length <= 61){
				$complete_sequence =~ s/(.{10})/$1 /g;
				if (length ($complete_sequence) == 67){
					chop ($complete_sequence);
				}
				my $number_of_nt = ($count_fasta_61 * 60) + 1;
				printf GBOUT "%9s", "$number_of_nt";
				chop ($complete_sequence);
				printf GBOUT " $complete_sequence\n";	
				$count_fasta_61 ++;
			}
			else {
				my $lines = $length / 60;
				$lines = ceil($lines);
				@p = 0..$lines-1;
				foreach (@p){
					$_ *= 60;
				}
				foreach (@p){
					$_ += 1;
				}
				$complete_sequence =~ s/(.{10})/$1 /g;
				chomp ($complete_sequence);
				my @sequence_lines = split (/(.{66})/, $complete_sequence);
				foreach my $line (@sequence_lines){
					if ($line !~ m/^$/){
						my $number_of_nt = ($count * 60) + 1;
						printf GBOUT "%9s", "$number_of_nt";
						print GBOUT " $line\n";
						$count ++;
					}
				}
			}
		}
	}
	foreach my $i (@list) {
		my $outf = "$InFile"."_"."$i"."\.gb";
		open (GBOUT, ">>", "$outf") || die "Can not open GFF ouput file.$outf\n";
		print GBOUT"\/\/\n";
	}
}
elsif ($reverse == 1) {
	while (<SOURCE>) {
		chomp ($_);
		if (/^>/){
			(undef, undef, undef, $current_version, undef) = split (/\|/);
			if ($versions{$current_version}){
				$count_fasta_61 = 0;
				$count = 0;
				$InFile =~ s/.xml//;
				my $outf = "$InFile"."_"."$current_version"."\.gb";
				open (GBOUT, ">>", "$outf") || die "Can not open GFF ouput file.$outf\n";
				print GBOUT "ORIGIN      \n";
				print "sequence printed in $outf\n";
			}
		}
		elsif (/^[A-Za-z]/) {
			if ($versions{$current_version}){
				my $lineseq = $_;
				$lineseq =~ s/^\s+//;
				$complete_sequence = $lineseq;
				$complete_sequence =~ tr/A-Z/a-z/;
				chomp($complete_sequence);
				my $length = length($complete_sequence);
				if ($length <= 61){
					$complete_sequence =~ s/(.{10})/$1 /g;
					if (length ($complete_sequence) == 67){
						chop ($complete_sequence);
					}
					my $number_of_nt = ($count_fasta_61 * 60) + 1;
					printf GBOUT "%9s", "$number_of_nt";
					chop ($complete_sequence);
					printf GBOUT " $complete_sequence\n";	
					$count_fasta_61 ++;
				}
				else {
					my $lines = $length / 60;
					$lines = ceil($lines);
					@p = 0..$lines-1;
					foreach (@p){
						$_ *= 60;
					}
					foreach (@p){
						$_ += 1;
					}
					$complete_sequence =~ s/(.{10})/$1 /g;
					chomp ($complete_sequence);
					my @sequence_lines = split (/(.{66})/, $complete_sequence);
					foreach my $line (@sequence_lines){
						if ($line !~ m/^$/){
							my $number_of_nt = ($count * 60) + 1;
							printf GBOUT "%9s", "$number_of_nt";
							print GBOUT " $line\n";
							$count ++;
						}
					}
				}
			}
		}
	}
	@list = keys (%versions);
	foreach my $i (@list) {
		my $outf = "$InFile"."_"."$i"."\.gb";
		open (GBOUT, ">>", "$outf") || die "Can not open GFF ouput file.$outf\n";
		print GBOUT"\/\/\n";
	}
}
