#!/usr/bin/env perl
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#This Perl code uses Getopt::Euclid module to describe its functionnalities
#For more details have from command line: perl FetchGeneData_v1i1.pl --man
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#Here are some modules required (if not installed use: cpan install <module>)
use Test::Simple; #used to check run (modules)
use Modern::Perl '2011'; # used as warnings/strict
use Smart::Comments '###'; # used to trace run
use Getopt::Euclid;#used to document perl code
use File::Basename; # used to capture file names and locations
use File::Copy 'copy'; # used to archive some files
use File::Compare; # used to compare files (archiving)
use Cwd; # used to catch current working directory
use DateTime;#to fetch date and time; remark: can be switched to alternate sub
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#initialization
my $code = $0;
#in case $0 return path+name
if ($code =~ m/\/([^\/]+)$/) {
	$code = $1;}
my $maindir = getcwd();#to catch current location
use Cwd 'abs_path';
my $codelocname =  abs_path($0);
my $rundatetime = getdatetime();#for reporting
print "\n$rundatetime: $code ($codelocname) Process Started from folder $maindir.\n";
my %codonhash = loadtranslation();#required to check for ATG, Stop and putative frameshift
my %tmhash = loadtmvalues();
my $showtime;
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#Output folders (main and backup)
my $mainoutputdir  = $ARGV{'--outputfolder'};
if (-e $mainoutputdir) {#exist
	ASKOVERWRITE1:
	while (1>0) {#for ever
		print "Selected main output folder $mainoutputdir already exists. Would you like to overwrite (Yes / No):";
		my $answer = <STDIN>;
		if ($answer =~ m/^YES\n$/i) {#so going to remove all files
			opendir (DIR,"$mainoutputdir") ||
				die "Cannot read from main output folder $mainoutputdir $!\n";
			my @catcharr = grep {$_ ne "." & $_ ne ".."} readdir DIR;
			closedir DIR;
			my $pointer = 0;
			while (defined($catcharr[$pointer])) {
				my $file = "$mainoutputdir\/$catcharr[$pointer]";
				if (-f $file) {#if plain txt file
					unlink $file ||
						die "Cannot erase file $file $!\n";}
				++$pointer;}
			last ASKOVERWRITE1;}
		elsif ($answer =~ m/^NO\n$/i) {
			die "OK. Select another main output folder. Bye! $!\n";}
		else {
			print "Bad answer. Try again (or CTRL-C to quit)\n";}}}
else {#creating
	mkdir $mainoutputdir ||
		die "Cannot create main output folder $mainoutputdir $!\n";}
my $backupoutputdir  = $ARGV{'--backupfolder'};
#NEW LAST MINUTE: should be different from $mainoutputdir
if ($backupoutputdir eq $mainoutputdir) {
	die "You provide a backup folder that is the main folder:\nMAIN:\n$mainoutputdir\nBACKUP:\n$backupoutputdir\nCorrect and try again!";}
if (-e $backupoutputdir) {
	ASKOVERWRITE2:
	while (1>0) {#for ever
		print "Selected backup output folder $backupoutputdir already exists. Would you like to overwrite (Yes / No):";
		my $answer = <STDIN>;
		if ($answer =~ m/^YES\n$/i) {#so going to remove all files
			opendir (DIR,"$backupoutputdir") ||
				die "Cannot read from backup output folder $backupoutputdir $!\n";
			my @catcharr = grep {$_ ne "." & $_ ne ".."} readdir DIR;
			closedir DIR;
			my $pointer = 0;
			while (defined($catcharr[$pointer])) {
				my $file = "$backupoutputdir\/$catcharr[$pointer]";
				if (-f $file) {#if plain txt file
					unlink $file ||
						die "Cannot erase file $file $!\n";}
				++$pointer;}
			last ASKOVERWRITE2;}
		elsif ($answer =~ m/^NO\n$/i) {
			die "OK. Select another backup output folder. Bye! $!\n";}
		else {
			print "Bad answer. Try again (or CTRL-C to quit)\n";}}}
else {#creating
	mkdir $backupoutputdir ||
		die "Cannot create backup output folder $backupoutputdir $!\n";}
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#creating a Run.log file (main output folder) and copying code to main output folder
my $runlog = "$mainoutputdir\/Run.Log";
open (LOG,">>$runlog") ||
	die "Cannot create LOG file $runlog $!\n";
print LOG "PERL_CODE\t$codelocname\nRUN_DATE_TIME\t$rundatetime\n";
print LOG "MAIN_OUTPUT_FOLDER\t$mainoutputdir\nBACKUP_OUTPUT_FOLDER\t$backupoutputdir\n";
#recording perl code
my $codecopy = "$mainoutputdir\/$code";
copy ($codelocname,$codecopy) ||
	die "Code Copy failed: $!\n";
#subroutine that copies main to backup (used at different time points)
print LOG "PERL_CODE_STATUS\tRecorded\n";
close LOG;
copymaintobackup($mainoutputdir,$backupoutputdir);
open (LOG,">>$runlog") ||
	die "Cannot open LOG file $runlog $!\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#loading genes for genes list file
my $geneslist = $ARGV{'--geneslist'};
if ((!-e $geneslist) || (!-f $geneslist)) {
	die "\nError: Proposed genes list file ($geneslist) does not exist or is not a plain text file. $!\n";}
#modification to have genes ranked by name and check if uniq
#also archiving input genes list file
my $geneslistname; 
if ($geneslist =~ m/\/([^\/]+)$/) {
	$geneslistname = $1;} 
my @genesarr;
print LOG "GENES_LIST_FILE\t$geneslist\n";
if (open(GENES,"<$geneslist")) {
	my $genesheader = uc <GENES>;
	chomp $genesheader;
	if ($genesheader ne "GENENAME\tGENEINFO") {
		die "Error: Proposed genes list file does not have a correct header ($genesheader). Should be: GeneName GeneInfo. $!\n";}
	my $ctlines = 0;
	while (defined(my $line = <GENES>)) {
		chomp $line;
		if ($line =~ m/^([^\t]+)\t([^\t]*)$/) {
			$genesarr[$ctlines] = uc $1 . "\t$2";
			++$ctlines;}
		else {
			close GENES;
			die "Error: Line ($line) in submitted genes list file is not as expected (format 2 tab delimited fields) $!\n";}}
	close GENES;
	if ($ctlines == 0) {
		print LOG "GENES_LIST_FILE_STATUS\tFAILED\nGENES_COUNT\t\n";
		die "\nError: Cannot find any valid entry from genes list file $geneslist $!\n";}
	#now ranking by genes names
	my $field1 = 1;
	my @rankarr;
	@rankarr = map {$_->[0]} 
	sort {
		my @a_fields = @$a[1..$#$a];
		my @b_fields = @$b[1..$#$b];
		$a_fields[$field1-1] cmp $b_fields[$field1-1]}
		#	||
		#$a_fields[$field2-1] <=> $b_fields[$field2-1]
		#	||
		#$a_fields[$field3-1] <=> $b_fields[$field3-1]}
	map {[$_, split /\t/]} @genesarr;
	undef @genesarr;
	my $pointer = 0;
	while (defined($rankarr[$pointer])) {
		if ($rankarr[$pointer] =~ m/([^\t]+)\t([^\t]*)$/) {
			$genesarr[++$genesarr[0][0]][1] = $1;
			$genesarr[$genesarr[0][0]][2] = $2;}
		++$pointer;}
	undef @rankarr;}
else {
	print LOG "GENES_LIST_FILE_STATUS\tFAILED\nGENES_COUNT\t\n";
	die "\nError: Cannot read from genes list file $geneslist $!\n";}
#copying this file to main output folder
my $copyofgeneslist = "$mainoutputdir\/$geneslistname";
if (copy ($geneslist,$copyofgeneslist)) {}
else {
	print LOG "GENES_LIST_FILE_STATUS\tPASSED BUT ARCHIVING FAILURE\nGENES_COUNT\t$genesarr[0][0]\n";
	die "Unable to record genes list file in main output folder $!\n";}
print LOG "GENES_LIST_FILE_STATUS\tPASSED,ARCHIVED\nGENES_COUNT\t$genesarr[0][0]\n";
close LOG;
copymaintobackup($mainoutputdir,$backupoutputdir);
open (LOG,">>$runlog") ||
	die "Cannot open LOG file $runlog $!\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
##creating an error report file
my $errorfile = "$mainoutputdir\/ErrorReport.txt";
open (ERROR,">>$errorfile") ||
	die "Cannot create the error report file $errorfile $!\n";
print ERROR "ErrorLocation\tErrorInfo\n";
print LOG "ERROR_FILE\t$errorfile\n";
my $cterrors = 0;
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#preliminary check of required input files (Genbank and FASTA) from input folder
my $errorlocation = "Preliminary check of Genbank and FASTA input files";
my $inputdir  = $ARGV{'--inputfolder'};
if ((!-e $inputdir) || (!-d $inputdir)) {
	die "\nError: Proposed Input folder ($inputdir) does not exist or is not a folder. $!\n";}
print LOG "INPUT_FOLDER\t$inputdir\n";
my $kchk1 = 1;
while ($kchk1 <= $genesarr[0][0]) {
	my $gbfile = "$inputdir\/$genesarr[$kchk1][1]\.gb";
	if ((-e $gbfile) && (-f $gbfile)) {
		$genesarr[$kchk1][3] = $gbfile;}
	else {
		print ERROR "$errorlocation\tUnable to identify Genbank file for gene $genesarr[$kchk1][1]\n";
		++$cterrors;}
	my $fastafile = "$inputdir\/$genesarr[$kchk1][1]\.fasta";
	if ((-e $fastafile) && (-f $fastafile)) {
		$genesarr[$kchk1][4] = $fastafile;}
	else {
		print ERROR "$errorlocation\tUnable to identify FASTA file for gene $genesarr[$kchk1][1]\n";
		++$cterrors;}
	++$kchk1;}
close ERROR;
if ($cterrors != 0) {
	print LOG "STATUS_PRELIMINARY_CHECK_INPUT_FILES\tFAILED\n";
	close LOG;
	copymaintobackup($mainoutputdir,$backupoutputdir);
	die "Some errors have been identified for $errorlocation. Check error file. $!\n";}
else {
	print LOG "STATUS_PRELIMINARY_CHECK_INPUT_FILES\tPASSED\n";
	close LOG;
	copymaintobackup($mainoutputdir,$backupoutputdir);
	open (ERROR,">>$errorfile") ||
		die "Cannot open the error report file $errorfile $!\n";
	open (LOG,">>$runlog") ||
		die "Cannot open LOG file $runlog $!\n";}
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#fetching the blast extend value and checking for BlastDB and other thresholds
my $blastextend  = $ARGV{'--blastextend'};
if (not defined($blastextend)) {
	$blastextend = 21;}
my $blastdbgeneric = $ARGV{'--blastdb'};
my $blastfile = $blastdbgeneric . "\.nhr";
if (!-e $blastfile) {
	die "Cannot find Blast DB file $blastfile";}
$blastfile = $blastdbgeneric . "\.nin";
if (!-e $blastfile) {
	die "Cannot find Blast DB file $blastfile";}
$blastfile = $blastdbgeneric . "\.nog";
if (!-e $blastfile) {
	die "Cannot find Blast DB file $blastfile";}
$blastfile = $blastdbgeneric . "\.nsd";
if (!-e $blastfile) {
	die "Cannot find Blast DB file $blastfile";}
$blastfile = $blastdbgeneric . "\.nsi";
if (!-e $blastfile) {
	die "Cannot find Blast DB file $blastfile";}
$blastfile = $blastdbgeneric . "\.nsq";
if (!-e $blastfile) {
	die "Cannot find Blast DB file $blastfile";}
my $blastred  = $ARGV{'--blastredlimit'};
my $blastorange  = $ARGV{'--blastorangelimit'};
my $blastyellow  = $ARGV{'--blastyellowlimit'};
my $blastevalue  = $ARGV{'--blastevalue'};
my $blastpctidentity  = $ARGV{'--blastpctidentity'};
my $blastmaxwinsize = $ARGV{'--blastmaxwinsize'};
my $blastminwinsize = $ARGV{'--blastminwinsize'};
my $blastmaxhits = $ARGV{'--blastmaxhits'};
print LOG "BLASTN_FILES_STATUS\tPASSED\n";
print LOG "RED TM THRESHOLD\t$blastred\n";
print LOG "ORANGE TM THRESHOLD\t$blastorange\n";
print LOG "YELLOW TM THRESHOLD\t$blastyellow\n";
print LOG "BLAST EVALUE THRESHOLD\t$blastevalue\n";
print LOG "BLAST PCT IDENTITY THRESHOLD\t$blastpctidentity\n";
print LOG "MAX WINDOW SIZE\t$blastmaxwinsize\n";
print LOG "MIN WINDOW SIZE\t$blastminwinsize\n";
print LOG "BLAST MAX HITS THRESHOLD\t$blastmaxhits\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#creating all output files
my $genesbed = "$mainoutputdir\/genes.bed";
open(GENESBED,">>$genesbed") ||
	die "Cannot create genes bed file $genesbed $!\n";
my $codexbed = "$mainoutputdir\/codex.bed";
open(CODEXBED,">>$codexbed") ||
	die "Cannot create codex bed file $codexbed $!\n";
my $genestable = "$mainoutputdir\/genes.txt";
open(GENESTABLE,">>$genestable") ||
	die "Cannot create genes table file $genestable $!\n";
my $codextable = "$mainoutputdir\/codex.txt";
open(CODEXTABLE,">>$codextable") ||
	die "Cannot create codex table file $codextable $!\n";
print CODEXTABLE "CodexID\tGeneName\tCodexNum\tRefGeno\tRankingChr\tChr\tStrand";
print CODEXTABLE "\tInGenePosIN\tInGenePosOUT\tRefGenoPosIN\tRefGenoPosOUT\tSize\tSequence";
print CODEXTABLE "\tCodexUcscID\tCodexUcscLink";
print CODEXTABLE "\tNumProductList\tIDproductList\tGeneFullName\tGeneType\tGenePosIN";
print CODEXTABLE "\tGenePosOUT\tGeneSize\tGeneProductCount\tGeneCodexCount";
print CODEXTABLE "\tGeneHgnc\tGeneHgncLink\tGeneEntrez\tGeneEntrezLink\tGeneUcsc\tGeneUcscLink\tGeneMim\tGeneMimLink\n";
my $cdstable = "$mainoutputdir\/cds.txt";
open(CDSTABLE,">>$cdstable") ||
	die "Cannot create CDS table file $cdstable $!\n";	
my $genesfasta = "$mainoutputdir\/genes.fasta";
open(GENESFASTA,">>$genesfasta") ||
	die "Cannot create genes FASTA file $genesfasta $!\n";
my $codexfasta = "$mainoutputdir\/codex.fasta";
open(CODEXFASTA,">>$codexfasta") ||
	die "Cannot create codex FASTA file $codexfasta $!\n";
my $cdsfasta = "$mainoutputdir\/cds.fasta";
open(CDSFASTA,">>$cdsfasta") ||
	die "Cannot create CDS FASTA file $cdsfasta $!\n";
my $codexblast = "$mainoutputdir\/codex.blast";
open(CODEXBLAST,">>$codexblast") ||
	die "Cannot create codex BLAST file $codexblast $!\n";
my $synonymfile = "$mainoutputdir\/genesynonyms.txt";
open(SYN,">>$synonymfile") ||
	die "Cannot create Genes Synonyms file $synonymfile $!\n";
print LOG "OUTPUT_FILES_CREATION\tPASSED\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#loading FASTA data (and generating genes FASTA and BED file
#structure
#>gi|568815586:8822472-8887202 Homo sapiens chromosome 12, GRCh38.p7 Primary Assembly
#GTCCACCACCTAACCTGGTGGTTACAGAAGGAGGAACACAGTCCTCCTGTGTTCCAGGTGCCAGGAGGTG
$errorlocation = "Extracting FASTA genes sequences";
my $keeprefgeno;
my $kf1 = 1;
my @genesbedarr;
FETCHFASTA:
while ($kf1 <= $genesarr[0][0]) {
	my $fasta = $genesarr[$kf1][4];
	open (FASTA,"<$fasta") ||
		die "Unable to read from FASTA file $fasta for gene $genesarr[$kf1][1] $!\n";
	my $idline = <FASTA>;
	chomp $idline;
	if ($idline =~ m/\:[^\d]*(\d+)\-(\d+)[^\d].+CHROMOSOME\s+([^\,]+)[\,]\s+([^\s]+)\s/i) {
		my $posin = $1;
		my $posout = $2;
		my $chr = uc $3;
		my $strand;
		my $rankchr;
		my $refgeno = $4;
		if(not defined($keeprefgeno)) {
			$keeprefgeno = $refgeno;}
		if ($posin > $posout) {
			my $keepposin = $posin;
			$posin = $posout;
			$posout = $keepposin;
			$strand = "-";}
		else {
			$strand = "+";}
		if ($chr =~ m/^(\d)$/) {
			$rankchr = "Chr0$1";}
		else {
			$rankchr = "Chr$chr";}
		print GENESFASTA "\>$genesarr[$kf1][1]\_GENE chr$chr\:$posin\-$posout\_$strand\_$refgeno\n";	
		my $genesize = $posout - $posin + 1;
		$genesarr[$kf1][5] = $rankchr;
		$genesarr[$kf1][6] = $chr;
		$genesarr[$kf1][7] = $posin;
		$genesarr[$kf1][8] = $posout;
		$genesarr[$kf1][9] = $strand;
		$genesarr[$kf1][10] = $refgeno;
		$genesarr[$kf1][11] = $genesize;
		$genesbedarr[$kf1-1] = "$rankchr\tchr$chr\t$posin\t$posout\t$genesarr[$kf1][1]\_chr$chr\_$posin\_$posout\_$genesize\_$strand\_$refgeno";}
	else {
		++$cterrors;
		print ERROR "$errorlocation\t$genesarr[$kf1][1] ID line has not the expected format: $idline\n";}
	my $geneseq = "";
	my $ctnotiupac = 0;
	my $stringnotiupac = "";
	while (defined(my $line = <FASTA>)) {
		$line = uc $line;
		while ($line =~ s/[\s\t\r\f\n]//) {}
		while ($line =~ s/([^A^C^G^T^V^B^D^H^M^Y^K^R^S^W^N])//) {
			$stringnotiupac = $stringnotiupac . $1;
			++$ctnotiupac;}
		$geneseq = $geneseq . $line;}
	close FASTA;
	if ($ctnotiupac != 0) {
		++$cterrors;
		print ERROR "$errorlocation\t$genesarr[$kf1][1] sequence contains $ctnotiupac Not IUPAC chars ($stringnotiupac)\n";}
	else {
		my $checksize = length($geneseq);
		if ($genesarr[$kf1][11] == $checksize) {
			print GENESFASTA "$geneseq\n";
			$genesarr[$kf1][12] = $geneseq;}
		else {
			++$cterrors;
			print ERROR "$errorlocation\t$genesarr[$kf1][1] sequence has not expected length ($checksize vs. $genesarr[$kf1][11]\n";}}
	++$kf1;}
close ERROR;
close GENESFASTA;			
if ($cterrors != 0) {
	print LOG "STATUS_FASTA_FILES\tFAILED\n";
	close LOG;
	copymaintobackup($mainoutputdir,$backupoutputdir);
	die "Some errors have been identified for $errorlocation. Check error file. $!\n";}
else {
	#ranking the bed data
	my $field1 = 1;
	my $field2 = 3;
	my $field3 = 4;
	my @rankarr;
	@rankarr = map {$_->[0]} 
	sort {
		my @a_fields = @$a[1..$#$a];
		my @b_fields = @$b[1..$#$b];
		$a_fields[$field1-1] cmp $b_fields[$field1-1]
			||
		$a_fields[$field2-1] <=> $b_fields[$field2-1]
			||
		$a_fields[$field3-1] <=> $b_fields[$field3-1]}
	map {[$_, split /\t/]} @genesbedarr;
	undef @genesbedarr;
	my $klgb = 0;
	while (defined($rankarr[$klgb])) {
		if ($rankarr[$klgb] =~ m/^[^\t]+\t(.+)$/) {
			print GENESBED "$1\n";}
			++$klgb;}
	close GENESBED;
	print LOG "STATUS_FASTA_FILES\tPASSED\n";
	close LOG;
	copymaintobackup($mainoutputdir,$backupoutputdir);
	open (ERROR,">>$errorfile") ||
		die "Cannot open the error report file $errorfile $!\n";
	open (LOG,">>$runlog") ||
		die "Cannot open LOG file $runlog $!\n";}
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#Now dealing with genbank files
#LOCUS       NC_000012             144481 bp    DNA     linear   CON 06-JUN-2016
#DEFINITION  Homo sapiens chromosome 12, GRCh38.p7 Primary Assembly.
#ACCESSION   NC_000012 REGION: complement(21797389..21941869) GPC_000001304
$errorlocation = "ReviewingGenbankData";
my @codexbed;
my @genesynonym;
my $pointersyn = 0;
my @geneinfoarr;
my @keepinfoblast;
my @finalfinalcdx;
my $pointfinalfinalcdx = 0;
my $kgb1 = 1;
FETCHGB:
while ($kgb1 <= $genesarr[0][0]) {
	my $gb = $genesarr[$kgb1][3];
	open (GB,"<$gb") ||
		die "Unable to read from GB file $gb for gene $genesarr[$kgb1][1] $!\n";
	my $firstline = <GB>;
	if ($firstline =~ m/^LOCUS\s+[^\s]+\s+(\d+)\s+/i) {
		my $genexsize = $1;
		if ($genexsize != $genesarr[$kgb1][11]) {
			++$cterrors;
			print ERROR "$errorlocation\t$genesarr[$kgb1][1] GB firstline sequence has not expected length ($genexsize vs. $genesarr[$kgb1][11]\n";}}
	else {
		++$cterrors;
		print ERROR "$errorlocation\t$genesarr[$kgb1][1] GB firstline has not the expected format\n";}
	my $secondline = <GB>;
	if ($secondline =~ m/^.+CHROMOSOME\s+([^\,]+)[\,]\s+([^\s]+)\s/i) {
		my $chrx = $1;
		my $refgenox = $2;
		if (($chrx ne $genesarr[$kgb1][6]) || ($refgenox ne $genesarr[$kgb1][10])) {
			++$cterrors;
			 print ERROR "$errorlocation\t$genesarr[$kgb1][1] GB second line info not compatible with FASTA $chrx ne $genesarr[$kgb1][6] or $refgenox ne $genesarr[$kgb1][10]\n";}}
	else {
		++$cterrors;
		print ERROR "$errorlocation\t$genesarr[$kgb1][1] GB second line has not the expected format\n";}
	my $thirdline = <GB>;
	if ($thirdline =~ m/[^\d](\d+)[\.]+(\d+)[^\d]/i) {
		my $posinx = $1;
		my $posoutx = $2;
		if (($posinx != $genesarr[$kgb1][7]) || ($posoutx != $genesarr[$kgb1][8])) {
			++$cterrors;
			 print ERROR "$errorlocation\t$genesarr[$kgb1][1] GB third line info not compatible with FASTA $posinx != $genesarr[$kgb1][7] or $posoutx != $genesarr[$kgb1][8]\n";}}
	else {
		++$cterrors;
		print ERROR "$errorlocation\t$genesarr[$kgb1][1] GB third line has not the expected format\n";}
	#everything is cross checked vs FASTA so let's go now for features
	GOFEATURES:
	while (defined(my $line = <GB>)) {
		if ($line =~ m/^FEATURES\s+/) {
			last GOFEATURES;}}
	my @features;
	FETCHFEATURES:
	while (defined(my $line = <GB>)) {
		chomp $line;
		if ($line =~ m/^ORIGIN\s*/) {
			last FETCHFEATURES;}
		if ($line =~ m/^\s{5}([^\s]+)\s+([^\s].*)$/) {
			$features[++$features[0][0]][1] =$1;
			$features[$features[0][0]][2] = $2;}
		else {
			$features[$features[0][0]][2] = $features[$features[0][0]][2] . $line;}}
	close GB;
	#all features are recorded now let's review these
	#print "$genesarr[$kgb1][1] has $features[0][0] features\n";
	#first is to capture info from gene feature 
	#http://omim.org/entry/613726	
	#https://www.ncbi.nlm.nih.gov/gene/55129
	#http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=25519
	#     gene            1..144481
    #                 /gene="ABCC9"
    #                 /gene_synonym="ABC37; ATFB12; CANTU; CMD1O; SUR2"
    #                 /note="ATP binding cassette subfamily C member 9; Derived
    #                 by automated computational analysis using gene prediction
    #                 method: BestRefSeq,Gnomon."
    #                 /db_xref="GeneID:10060"
    #                 /db_xref="HGNC:HGNC:60"
    #                 /db_xref="MIM:601439"
	my $kft1 = 1;
	FOUNDGENEINFO:
	while ($kft1 <= $features[0][0]) {
		#here a removal of unvisible parameters in dupli
		while ($features[$kft1][2] =~ s/[\t\n\r\f]/ /) {}
		while ($features[$kft1][2] =~ s/  / /) {}
		if ($features[$kft1][1] =~ m/^GENE$/i) {
			if ($features[$kft1][2] =~ m/\/GENE\=\"([^\"]+)\"/i) {
				my $checkgene = uc $1;
				if ($checkgene eq $genesarr[$kgb1][1]) {
					$geneinfoarr[$kgb1][0] = $checkgene;
					if ($features[$kft1][2] =~ m/\/NOTE\=\"([^\;]+)\;/i) {
						$geneinfoarr[$kgb1][1] = $1;}
					if ($features[$kft1][2] =~ m/\/DB_XREF\=\"GENEID\:([^\"]+)\"/i) {
						$geneinfoarr[$kgb1][3] = $1;}
					if ($features[$kft1][2] =~ m/\/DB_XREF\=\"HGNC\:HGNC\:([^\"]+)\"/i) {
						$geneinfoarr[$kgb1][4] = $1;}
					if ($features[$kft1][2] =~ m/\/DB_XREF\=\"MIM\:([^\"]+)\"/i) {
						$geneinfoarr[$kgb1][5] = $1;}
					if ($features[$kft1][2] =~ m/\/GENE_SYNONYM\=\"([^\"]+)\"/i) {
						$geneinfoarr[$kgb1][2] = uc $1 . ";";
						while ($geneinfoarr[$kgb1][2] =~ s/[\s\t\n\r\f]//) {}}
					last FOUNDGENEINFO;}}}
		else {}
		++$kft1;}
	if (not defined($geneinfoarr[$kgb1][0])) {
		++$cterrors;
		print "Error $geneinfoarr[$kgb1][1]\n";
		print ERROR "$errorlocation\t$genesarr[$kgb1][1] GB: Unable to identify Gene Info\n";}
	else {
		#print "$geneinfoarr[$kgb1][0]\n";
		my $fillna = 1;
		while ($fillna <= 5) {
			if (not defined($geneinfoarr[$kgb1][$fillna])) {
				$geneinfoarr[$kgb1][$fillna] = "NA";}
			#print "\t$geneinfoarr[$kgb1][$fillna]\n";
			++$fillna;}
		#preparation of synomym list
		my $getsyn = $geneinfoarr[$kgb1][2];
		#https://www.ncbi.nlm.nih.gov/gene/55129
	#http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=25519
		$genesynonym[$pointersyn] = "$geneinfoarr[$kgb1][0]\t$geneinfoarr[$kgb1][0]\tHGNC:$geneinfoarr[$kgb1][4]\t=hyperlink\(\"http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=$geneinfoarr[$kgb1][4]\";\"HgncLink:$geneinfoarr[$kgb1][1]\"\)\tEntrez:$geneinfoarr[$kgb1][3]\t=hyperlink\(\"https://www.ncbi.nlm.nih.gov/gene/$geneinfoarr[$kgb1][3]\";\"EntrezLink:$geneinfoarr[$kgb1][1]\"\)";
		++$pointersyn;
		if ($getsyn ne "NA") {
			while ($getsyn =~ m/([^\;]+)\;+/igc) {
				$genesynonym[$pointersyn] = "$1\t$geneinfoarr[$kgb1][0]\tHGNC:$geneinfoarr[$kgb1][4]\t=hyperlink\(\"http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=$geneinfoarr[$kgb1][4]\";\"HgncLink:$geneinfoarr[$kgb1][1]\"\)\tEntrez:$geneinfoarr[$kgb1][3]\t=hyperlink\(\"https://www.ncbi.nlm.nih.gov/gene/$geneinfoarr[$kgb1][3]\";\"EntrezLink:$geneinfoarr[$kgb1][1]\"\)";
				++$pointersyn;}}}
	#now identification of CDS / ncRNA
#     CDS             join(5196..5337,7947..8088,15807..15928,24767..24933,
#                      135810..135872,140689..140826)
#                      /gene="ABCC9"
#                      /gene_synonym="ABC37; ATFB12; CANTU; CMD1O; SUR2"
#                      /note="Derived by automated computational analysis using
#                      gene prediction method: Gnomon."
#                      /codon_start=1
#                      /product="ATP-binding cassette sub-family C member 9
#                      isoform X5"
#                      /protein_id="XP_005253347.1"
#                      /db_xref="GI:530398974"
#                      /db_xref="GeneID:10060"
#                      /db_xref="HGNC:HGNC:60"
#                      /db_xref="MIM:601439"	
    my @cdsarr;
    my @ncrnaarr;
    $kft1 = 1;
	FOUNDGENEINFO:
	while ($kft1 <= $features[0][0]) {
		if ($features[$kft1][1] =~ m/^CDS$/i) {
			if ($features[$kft1][2] =~ m/^(.+)\/GENE\=\"([^\"]+)\"/i) {
				my $join = $1;
				my $checkgene = uc $2;
				if ($checkgene eq $genesarr[$kgb1][1]) {
					#processing CDS to validate it (not yet sequence but check if unique)
					if ($join !~ m/[\>\<]/) {#so not truncated CDS
						$join = $join .",";
						while ($join =~ s/[^\d^\.^\,]//) {}
						while ($join =~ s/\,\,/\,/) {}
						my $newjoin = "";
						while ($join =~ m/([^\,]+)\,/igc) {
							my $joinval = $1;
							if ($joinval =~ m/^(\d+)$/) {
								$newjoin = $newjoin . "$1\_$1,";}
							elsif($joinval =~ m/^(\d+)[\.]+(\d+)$/) {
								$newjoin = $newjoin . "$1\_$2,";}
							else {
								die "Oups: Not catched. $features[$kft1][2]";}}
						#we have the join valid and reformatted. Let's see if this is not already recorded (different mRNA could lead to the same protein
						my $protid = "NA";;
						if ($features[$kft1][2] =~ m/\/PROTEIN_ID\=\"([^\"]+)\"/i) {
							$protid = $1;}
						#next for unique based on join
						my $kc1 = 1;
						my $foundcds = 0;
						VIEWEDCDS:
						while (defined($cdsarr[$kc1][1])) {
							if ($cdsarr[$kc1][1] eq $newjoin) {
								$foundcds = 1;
								last VIEWEDCDS;}
							++$kc1;}
						if ($foundcds == 0) {
							++$cdsarr[0][0];
							$cdsarr[$kc1][1] = $newjoin;
							$cdsarr[$kc1][2] = $protid;}}}}}
			elsif ($features[$kft1][1] =~ m/^NCRNA$/i) {
				if ($features[$kft1][2] =~ m/^(.+)\/GENE\=\"([^\"]+)\"/i) {
					my $join = $1;
					my $checkgene = uc $2;
					if ($checkgene eq $genesarr[$kgb1][1]) {
						#processing CDS to validate it (not yet sequence but check if unique)
						if ($join !~ m/[\>\<]/) {#so not truncated CDS
							$join = $join .",";
							while ($join =~ s/[^\d^\.^\,]//) {}
							while ($join =~ s/\,\,/\,/) {}
							my $newjoin = "";
							while ($join =~ m/([^\,]+)\,/igc) {
								my $joinval = $1;
								if ($joinval =~ m/^(\d+)$/) {
									$newjoin = $newjoin . "$1\_$1,";}
								elsif($joinval =~ m/^(\d+)[\.]+(\d+)$/) {
									$newjoin = $newjoin . "$1\_$2,";}
								else {
									die "Oups: Not catched. $features[$kft1][2]";}}
							#we have the join valid and reformatted. Let's see if this is not already recorded
							my $transcriptid = "NA";;
							if ($features[$kft1][2] =~ m/\/TRANSCRIPT_ID\=\"([^\"]+)\"/i) {
								$transcriptid = $1;}
							#next for unique based on join
							my $kc1 = 1;
							my $foundncrna = 0;
							VIEWEDNCRNA:
							while (defined($ncrnaarr[$kc1][1])) {
								if ($ncrnaarr[$kc1][1] eq $newjoin) {
									$foundncrna = 1;
									last VIEWEDNCRNA;}
								++$kc1;}
							if ($foundncrna == 0) {
								++$ncrnaarr[0][0];
								$ncrnaarr[$kc1][1] = $newjoin;
								$ncrnaarr[$kc1][2] = $transcriptid;}}}}}
			else {}#other features not interesting
		++$kft1;}
	if ((not defined($cdsarr[0][0])) && (not defined($ncrnaarr[0][0]))) {#should be one or the other
		++$cterrors;
		print "Error $geneinfoarr[$kgb1][1]\n";
		print ERROR "$errorlocation\t$genesarr[$kgb1][1] GB: Unable to identify Gene type\n";}
	#reviewing features
	if (not defined($cdsarr[0][0])) {#here dealing with ncRNA gene and extracting codex
		#no verification for ncRNA
		print "$ncrnaarr[0][0] ncRNA for $genesarr[$kgb1][1]";
		my @temparr;
		my $kr1 = 1;
		while ($kr1 <= $ncrnaarr[0][0]) {
			my $xjoin = $ncrnaarr[$kr1][1];
			my $exonseq;
			my $exonblast;
			my $ncrnaseq = "";
			while ($xjoin =~ m/([^\,]+)\,/igc) {
				my $exon = $1;
				if ($exon =~ m/^(\d+)\_(\d+)$/) {
					my $start = $1;
					my $end = $2;
					my $ktmp = 1;
					my $foundexon = 0;
					FINDEXON:
					while (defined($temparr[$ktmp][1])) {#unique exon here without additionnal checks upfront on ncRNA
						if ($temparr[$ktmp][1] eq $exon) { #already recorded but adding info about transcript ID
							$foundexon = 1;
							$temparr[$ktmp][17] =  $temparr[$ktmp][17] . "$kr1;";
							$temparr[$ktmp][10] = $temparr[$ktmp][10] . "$ncrnaarr[$kr1][2];";#list of trancripts
							++$temparr[$ktmp][9];#count of transcripts
							last FINDEXON;}
						++$ktmp;}
					if ($foundexon == 0) {#new one
						++$temparr[0][0];#count of exons
						$temparr[$ktmp][1] = $exon;#identifier posin posout
						$temparr[$ktmp][10] = "$ncrnaarr[$kr1][2];";#transcript list
						++$temparr[$ktmp][9];#count of transcript
						$temparr[$ktmp][2] = $start;#in gene posin
						$temparr[$ktmp][3] = $end;#in gen posout
						$temparr[$ktmp][4] = $genesarr[$kgb1][5] ;#rank chr
						$temparr[$ktmp][5] = $genesarr[$kgb1][6] ;#chr
						$temparr[$ktmp][6] = $genesarr[$kgb1][7] ;#posin
						$temparr[$ktmp][7] = $genesarr[$kgb1][8] ;#posout
						$temparr[$ktmp][8] = $genesarr[$kgb1][9] ;#strand
						$temparr[$ktmp][17] = "$kr1;";
						#fetching codex sequence and genome coordinates and codex blastseq
						my $start = $1;
						my $end = $2;
						my $ctnuc = 0;
						my $exonseq = "";
						my $exonblast = "";
						my $startb = $start - $blastextend;
						my $endb = $end + $blastextend;
						if ($startb < 1) {
							$startb = 1;}
						my $geneseq = $genesarr[$kgb1][12];
						my $before = 0;
						my $after = 0;
						CATCHSEQ2:
						while ($geneseq =~ m/(.)/igc) {
							my $nuc = uc $1;
							++$ctnuc;
							if ($ctnuc < $startb) {
								next CATCHSEQ2;}
							if ($ctnuc > $endb) {
								last CATCHSEQ2;}
							$exonblast = $exonblast . $nuc;
							if ($ctnuc < $start) {
								++$before;}
							if ($ctnuc > $end) {
								++$after;}
							if (($ctnuc >= $start) && ($ctnuc <= $end)) {
								$exonseq = $exonseq . $nuc;
								$ncrnaseq = $ncrnaseq . $nuc;}}
						my $viewlen = length($exonseq);
						my $chklen = $end -$start + 1;
						if ($viewlen != $chklen) {
							die "bad perl code processing. $viewlen of $exonseq  $chklen = $end - $start+1; $exon Correct! $!\n";}
						$temparr[$ktmp][11] = $exonseq;
						$temparr[$ktmp][12] = $exonblast;		
						$temparr[$ktmp][13] = length($exonseq);
						$temparr[$ktmp][14] = length($exonblast);
						$temparr[$ktmp][18] = "$before\_$after";
						#now refgeno coordinate of exon
						#deducing refgeno positions
						if ($genesarr[$kgb1][9] eq "+") {
							$temparr[$ktmp][15] = $genesarr[$kgb1][7] + $start - 1;
							$temparr[$ktmp][16] = $genesarr[$kgb1][7] + $end - 1;}
						else {
							$temparr[$ktmp][15] = $genesarr[$kgb1][8] - $end + 1;
							$temparr[$ktmp][16] = $genesarr[$kgb1][8] - $start + 1;}
						if (($temparr[$ktmp][16] - $temparr[$ktmp][15] + 1) != length($exonseq)) {
							print "$exonseq\n";
							die "Bad fetching refgeno positions for $genesarr[$kgb1][1] in ncRNA $kr1 $ncrnaarr[$kr1][2] $!\n";}}}#strand	
				else {
					print "\tOups what's this exon $exon ($ncrnaarr[$kr1][1])\n";}}
			++$kr1;}
		#so this array needs to be ranked on in gene posin and posout n(to start from 5' end)
		my @toranktmp;
		my $finalcodexncrna = $temparr[0][0];
		$kr1 = 1;
		while ($kr1 <= $temparr[0][0]) {
			my $kr2 = 1;
			$toranktmp[$kr1-1] = "";
			while (defined($temparr[$kr1][$kr2])) {
				$toranktmp[$kr1-1] = $toranktmp[$kr1-1] . "$temparr[$kr1][$kr2]\t";
				++$kr2;}
			++$kr1;}
		#now ranking by in gene posin and posout (to give ID number starting from 5')
		my $field1z = 1;
		my $field2z = 3;
		my @rankarrz;
		@rankarrz = map {$_->[0]} 
		sort {
			my @a_fields = @$a[1..$#$a];
			my @b_fields = @$b[1..$#$b];
			$a_fields[$field1z-1] <=> $b_fields[$field1z-1]
				||
			$a_fields[$field2z-1] <=> $b_fields[$field2z-1]}
		map {[$_, split /\t/]} @toranktmp;
		#printing in various files
		undef @toranktmp;
		my @finalcds;
		my @finalcodex;
		$kr1 = 0;
		while (defined($rankarrz[$kr1])) {
			++$finalcodex[0][0];
			my $refcod = $kr1 +1;
			my @linearr;
			while ($rankarrz[$kr1] =~ m/([^\t]+)\t/igc) {
				$linearr[++$linearr[0]] = $1;
				$finalcodex[$finalcodex[0][0]][$linearr[0]] = $1;}
			if ($linearr[0] != 18) {
				die "Oups: bad number of fields $linearr[0] != 18 for $genesarr[$kgb1][1] in Codex $refcod = $rankarrz[$kr1]\n";}
			#Fasta Codex Unique ID and info as Chr:Posin-Posout_Strand_Refgeno
			print CODEXFASTA "\>$genesarr[$kgb1][1]\_CODEX_$refcod chr$linearr[5]\:$linearr[15]\-$linearr[16]\_$linearr[8]\_$genesarr[$kgb1][10]\n";
			print CODEXFASTA "$linearr[11]\n";
			#Blast Codex Unique ID ... But sequence is Blast one
			my $xstrand = $genesarr[$kgb1][9];
			if ($xstrand eq "+") {
				$xstrand = "P";}
			else {
				$xstrand = "M";}
			#print CODEXBLAST "\>$genesarr[$kgb1][1]\_CODEX_$refcod\_$linearr[13]\_$linearr[14]\_BLASTADD_$linearr[26]";
			#print CODEXBLAST " GENE_chr$genesarr[$kgb1][6]\_$genesarr[$kgb1][7]\_$genesarr[$kgb1][8]\_$xstrand\n";
			print CODEXBLAST "\>$genesarr[$kgb1][1]\_CODEX_$refcod\_$linearr[15]\_$linearr[16]\_BLASTADD_$linearr[18] GENE_chr$genesarr[$kgb1][6]\_$genesarr[$kgb1][7]\_$genesarr[$kgb1][8]\_$xstrand\n";
			$linearr[12] = uc $linearr[12];
			print CODEXBLAST "$linearr[12]\n";
			#Bed codex (still ranked by gene name)
			print CODEXBED "chr$linearr[5]\t$linearr[15]\t$linearr[16]\t$genesarr[$kgb1][1]\_CODEX_$refcod\_chr$linearr[5]\_$linearr[15]\_$linearr[16]\_$linearr[8]\_$genesarr[$kgb1][10]\n";
			#adding to ncRNA exons list (CODEXTABLE)
			#print CODEXTABLE "RankingCodexID\tCodexID\tRefGeno\tRankingChr\tChr\tStrand";
			#print CODEXTABLE \tInGenePosIN\tInGenePosOUT\tRefGenoPosIN\tRefGenoPosOUT\tSize\tSequence";
			#print CODEXTABLE \tCodexUcscID\tCodexUcscLink";
			#print CODEXTABLE \tNumProductList\tIDproductList\tGeneName\tGeneFullName\tGeneType\tGenePosIN";
			#print CODEXTABLE "\tGenePosOUT\tGeneSize\tGeneProductCount\tGeneCodexCount";
			#print CODEXTABLE "\tGeneHgnc\tGeneHgncLink\tGeneEntrez\tGeneEntrezLink\tGeneUcsc\tGeneUcscLink\tGeneMim\tGeneMimLink\n";
			my $rkcdx;
			if ($refcod < 10) {
				$rkcdx = "00$refcod";}
			elsif ($refcod < 100) {
				$rkcdx = "0$refcod";}
			else {
				$rkcdx = $refcod;}
			print CODEXTABLE "$genesarr[$kgb1][1]\_CODEX_$refcod\t$genesarr[$kgb1][1]\t$refcod\t$genesarr[$kgb1][10]\t$linearr[4]\tchr$linearr[5]\t$linearr[8]";
			print CODEXTABLE "\t$linearr[2]\t$linearr[3]\t$linearr[6]\t$linearr[7]\t$linearr[13]\t$linearr[11]";
			print CODEXTABLE "\tchr$linearr[5]\:$linearr[6]\-$linearr[7]\t=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr$linearr[5]\:$linearr[6]\-$linearr[7]\";\"UCSC_CODEX_LINK\"\)";
			print CODEXTABLE "\t$linearr[17]\t$linearr[10]\t$geneinfoarr[$kgb1][1]\tncRNA\t$genesarr[$kgb1][7]";
			print CODEXTABLE "\t$genesarr[$kgb1][8]\t$genesarr[$kgb1][11]\t$ncrnaarr[0][0]\t$finalcodexncrna";
			print CODEXTABLE "\tHGNC:$geneinfoarr[$kgb1][4]\t=hyperlink\(\"http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=$geneinfoarr[$kgb1][4]\";\"HGNC_GENE_LINK\"\)";
			print CODEXTABLE "\tEntrez:$geneinfoarr[$kgb1][3]\t=hyperlink\(\"https://www.ncbi.nlm.nih.gov/gene/$geneinfoarr[$kgb1][3]\";\"Entrez_GENE_LINK\"\)";
			print CODEXTABLE "\tUCSC:chr$linearr[5]\:$genesarr[$kgb1][7]\-$genesarr[$kgb1][8]\t=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr$linearr[5]\:$genesarr[$kgb1][7]\-$genesarr[$kgb1][8]\";\"UCSC_GENE_LINK\"\)";
			print CODEXTABLE "\tOMIM:$geneinfoarr[$kgb1][5]\t=hyperlink\(\"http://omim.org/entry/$geneinfoarr[$kgb1][5]\";\"OMIM_GENE_LINK\"\)\n";
			$finalfinalcdx[$pointfinalfinalcdx] = "$genesarr[$kgb1][1]\_CODEX_$refcod\t$genesarr[$kgb1][1]\t$refcod\t$genesarr[$kgb1][10]\t$linearr[4]\tchr$linearr[5]\t$linearr[8]";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\t$linearr[2]\t$linearr[3]\t$linearr[6]\t$linearr[7]\t$linearr[13]\t$linearr[11]";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tchr$linearr[5]\:$linearr[6]\-$linearr[7]\t=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr$linearr[5]\:$linearr[6]\-$linearr[7]\";\"UCSC_CODEX_LINK\"\)";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\t$linearr[17]\t$linearr[10]\t$geneinfoarr[$kgb1][1]\tncRNA\t$genesarr[$kgb1][7]";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\t$genesarr[$kgb1][8]\t$genesarr[$kgb1][11]\t$ncrnaarr[0][0]\t$finalcodexncrna";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tHGNC:$geneinfoarr[$kgb1][4]\t=hyperlink\(\"http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=$geneinfoarr[$kgb1][4]\";\"HGNC_GENE_LINK\"\)";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tEntrez:$geneinfoarr[$kgb1][3]\t=hyperlink\(\"https://www.ncbi.nlm.nih.gov/gene/$geneinfoarr[$kgb1][3]\";\"Entrez_GENE_LINK\"\)";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tUCSC:chr$linearr[5]\:$genesarr[$kgb1][7]\-$genesarr[$kgb1][8]\t=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr$linearr[5]\:$genesarr[$kgb1][7]\-$genesarr[$kgb1][8]\";\"UCSC_GENE_LINK\"\)";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tOMIM:$geneinfoarr[$kgb1][5]\t=hyperlink\(\"http://omim.org/entry/$geneinfoarr[$kgb1][5]\";\"OMIM_GENE_LINK\"\)";
			++$pointfinalfinalcdx;	
			++$kr1;}
		print "\tProcessed $temparr[0][0] ncRNA exons\n";}
	else {#here dealing with protein coding gene and extracting codex
		#step one: Building CDS, checking for ATG and STOP and Translation
		print "$cdsarr[0][0] CDS for $genesarr[$kgb1][1] ";
		my $kc1 = 1;
		my @finalcds;
		my @finalcodex;
		my $cdsok = 0;
		FETCHCDSSEQ:
		while ($kc1 <= $cdsarr[0][0]) {#for all unique CDS
			my $cdsjoin = $cdsarr[$kc1][1];
			my $cdsprotid = $cdsarr[$kc1][2];
			my $cdsseq = "";
			my @tempcodex;
			while ($cdsjoin =~ m/([^\,]+)\,/igc) {
				my $exonstring = $1;
				if ($exonstring =~ m/^(\d+)\_(\d+)$/) {
					my $start = $1;
					my $end = $2;
					my $ctnuc = 0;
					my $exonseq = "";
					my $exonblast = "";
					my $startb = $start - $blastextend;
					my $endb = $end + $blastextend;
					if ($startb < 1) {
						$startb = 1;}
					my $geneseq = $genesarr[$kgb1][12];
					my $before = 0;
					my $after = 0;
					CATCHSEQ1:
					while ($geneseq =~ m/(.)/igc) {
						my $nuc = uc $1;
						++$ctnuc;
						if ($ctnuc < $startb) {
							next CATCHSEQ1;}
						if ($ctnuc > $endb) {
							last CATCHSEQ1;}
						$exonblast = $exonblast . $nuc;
						if ($ctnuc < $start) {
							++$before;}
						if ($ctnuc > $end) {
							++$after;}
						if (($ctnuc >= $start) && ($ctnuc <= $end)) {
							$exonseq = $exonseq . $nuc;
							$cdsseq = $cdsseq . $nuc;}}
					my $viewlen = length($exonseq);
					my $chklen = $end -$start + 1;
					if ($viewlen != $chklen) {
						die "bad perl code processing. $viewlen of $exonseq  $chklen = $end - $start+1; $exonstring Correct! $!\n";}
					#recording temporarilly all info for codex (only validated if CDS validated)
					$tempcodex[++$tempcodex[0][0]][1] = $exonseq;
					$tempcodex[$tempcodex[0][0]][0] = $exonstring;#to identify all exons as unique
					$tempcodex[$tempcodex[0][0]][2] = $exonblast;
					$tempcodex[$tempcodex[0][0]][3] = $start;
					$tempcodex[$tempcodex[0][0]][4] = $end;
					$tempcodex[$tempcodex[0][0]][5] = length($exonseq);
					$tempcodex[$tempcodex[0][0]][15] = length($exonblast);
					$tempcodex[$tempcodex[0][0]][6] = $cdsprotid;
					$tempcodex[$tempcodex[0][0]][7] = $genesarr[$kgb1][5] ;#gene rank chr
					$tempcodex[$tempcodex[0][0]][8] = $genesarr[$kgb1][6] ;#gene chr
					$tempcodex[$tempcodex[0][0]][9] = $genesarr[$kgb1][7] ;#gene posin
					$tempcodex[$tempcodex[0][0]][10] = $genesarr[$kgb1][8] ;#gene posout
					$tempcodex[$tempcodex[0][0]][11] = $genesarr[$kgb1][9] ;#gene strand
					$tempcodex[$tempcodex[0][0]][14] = $genesarr[$kgb1][1] ;#gene name
					$tempcodex[$tempcodex[0][0]][16] = "$before\_$after"; #capture what was added for blast
					#deducing refgeno positions
					if ($genesarr[$kgb1][9] eq "+") {
						$tempcodex[$tempcodex[0][0]][12] = $genesarr[$kgb1][7] + $start - 1;
						$tempcodex[$tempcodex[0][0]][13] = $genesarr[$kgb1][7] + $end - 1;}
					else {
						$tempcodex[$tempcodex[0][0]][12] = $genesarr[$kgb1][8] - $end + 1;
						$tempcodex[$tempcodex[0][0]][13] = $genesarr[$kgb1][8] - $start + 1;}
					if (($tempcodex[$tempcodex[0][0]][13] - $tempcodex[$tempcodex[0][0]][12] + 1) != length($exonseq)) {
						die "Bad fetching refgeno positions for $genesarr[$kgb1][1] in CDS $kc1 $cdsarr[$kc1][2] $!\n";}}
				else {
					die "What's this exon string $exonstring. Bad processing for $genesarr[$kgb1][1] in CDS $kc1 $cdsarr[$kc1][2] $!\n";}}
			#here we have a cdsseq and all related codex info: remark this are not error. we just devalidate the CDS
			#verification length is a multiple of 3 (frameshift)
			my $lencds = length($cdsseq);
			my $rd3len = sprintf("%.0f",($lencds/3));
			if ($rd3len != ($lencds/3)) {#skipped
				print"Length not a multiple of 3 for $genesarr[$kgb1][1] in CDS $kc1 $cdsarr[$kc1][2]\n";
				++$kc1;
				next FETCHCDSSEQ;}
			if ($cdsseq !~ m/^ATG/) {#skipped
				print "First codon not ATG for $genesarr[$kgb1][1] in CDS $kc1 $cdsarr[$kc1][2]\n";
				++$kc1;
				next FETCHCDSSEQ;}
			#checkink if last is stop and fetching only aa coding seq
			if ($cdsseq =~ m/^(.+)(.{3})$/) {
				my $lastcodon = $2;
				$cdsseq = $1;
				if (not defined($codonhash{$lastcodon})) {#skipped
					print "First codon not ATG for $genesarr[$kgb1][1] in CDS $kc1 $cdsarr[$kc1][2]\n";
					++$kc1;
					next FETCHCDSSEQ;}
				if ($codonhash{$lastcodon} ne "B") {#skipped
					print "$lastcodon last codon is not a Stop codon for $genesarr[$kgb1][1] in CDS $kc1 $cdsarr[$kc1][2]\n";
					++$kc1;
					next FETCHCDSSEQ;}}
			#verification of translation without internal stop	
			my $checkstop = $cdsseq;
			my $cdsprot = "";
			my $ctaa = 0;
			my $cdsctrl = "";
			while ($checkstop =~ m/(.{3})/igc) {
				my $codon = $1;
				++$ctaa;
				if (not defined($codonhash{$codon})) {
					print "$codon codon at aa pos $ctaa or $genesarr[$kgb1][1] in CDS $kc1 $cdsarr[$kc1][2]: undefined\n$cdsctrl\n$cdsprot\n";
					++$kc1;
					next FETCHCDSSEQ;}
				if ($codonhash{$codon} eq "B") {
					print "$codon codon at aa pos $ctaa or $genesarr[$kgb1][1] in CDS $kc1 $cdsarr[$kc1][2]: stop codon\n$cdsctrl\n$cdsprot\n";
					++$kc1;
					next FETCHCDSSEQ;}
				$cdsctrl = $cdsctrl . $codon;
				$cdsprot = $cdsprot . $codonhash{$codon};}
			++$cdsok;
			print "$kc1 ";
			#print "\t$cdsok CDS ok\n";
			#here we have a CDS ok (and unique). But we cannot do anything else then recording info for this CDS and its Codex 
			#(only when all CDS and codex are known we can move on; example: Gene PTEN has a CDS not framed)
			$finalcds[++$finalcds[0][0]][1] = $cdsprotid;#product ID
			$finalcds[$finalcds[0][0]][2] = $cdsarr[$kc1][1];#join of exons
			$finalcds[$finalcds[0][0]][3] = $cdsseq;#nucleotide sequence (excluded stop)
			$finalcds[$finalcds[0][0]][4] = length($cdsseq);#length of nucleotide sequence
			$finalcds[$finalcds[0][0]][5] = $cdsprot;#amino acid sequence
			$finalcds[$finalcds[0][0]][6] = "$genesarr[$kgb1][1]\_CDS\_$finalcds[0][0]";#cds id
			$finalcds[$finalcds[0][0]][7] = $genesarr[$kgb1][1];#gene name 
			$finalcds[$finalcds[0][0]][8] = $genesarr[$kgb1][10];#ref geno
			$finalcds[$finalcds[0][0]][9] = $genesarr[$kgb1][5];#RankChr
			$finalcds[$finalcds[0][0]][10] = $genesarr[$kgb1][6];#Chr
			$finalcds[$finalcds[0][0]][11] = $genesarr[$kgb1][9];#strand
			$finalcds[$finalcds[0][0]][12] = $genesarr[$kgb1][7];#posin 
			$finalcds[$finalcds[0][0]][13] = $genesarr[$kgb1][8];#posout 
			$finalcds[$finalcds[0][0]][14] = $genesarr[$kgb1][11];#size
			$finalcds[$finalcds[0][0]][15] = "HGNC:$geneinfoarr[$kgb1][4]";
			$finalcds[$finalcds[0][0]][16] = "=hyperlink\(\"http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=$geneinfoarr[$kgb1][4]\";\"HgncLink:$geneinfoarr[$kgb1][1]\"\)";
			$finalcds[$finalcds[0][0]][17] = "Entrez:$geneinfoarr[$kgb1][3]";
			$finalcds[$finalcds[0][0]][18] = "=hyperlink\(\"https://www.ncbi.nlm.nih.gov/gene/$geneinfoarr[$kgb1][3]\";\"EntrezLink:$geneinfoarr[$kgb1][1]\"\)"; #last = missed previously
			$finalcds[$finalcds[0][0]][19] = $tempcodex[0][0];#number of codex
			$finalcds[$finalcds[0][0]][20] = "";#used to concatenate later all codex (once they are ranked);
			#now let us record codex (unique...)
			my $kcod1 = 1;
			while ($kcod1 <= $tempcodex[0][0]) {#for all recorded codex
				my $kcod2 = 1;
				my $alreadyincodex = 0;
				RECORDCODEX:
				while(defined($finalcodex[$kcod2][0])) {
					if ($tempcodex[$kcod1][0] eq $finalcodex[$kcod2][0]) {#same codex
						++$finalcodex[$kcod2][16]; #count of cds for this codex
						$finalcodex[$kcod2][17] = $finalcodex[$kcod2][17] ."$finalcds[0][0];"; #list  of num id of cds
						$finalcodex[$kcod2][18] = $finalcodex[$kcod2][18] . "$finalcds[$finalcds[0][0]][1];";# list of protein id of cds 
						$alreadyincodex = 1;
						last RECORDCODEX;}
					++$kcod2;}
				if ($alreadyincodex == 0) {
					++$finalcodex[0][0]; 
					my $kcod3 = 0;
					while ($kcod3 <= 15) {
						$finalcodex[$kcod2][$kcod3] = $tempcodex[$kcod1][$kcod3];
						++$kcod3;}
					++$finalcodex[$kcod2][16];
					$finalcodex[$kcod2][17] = "$finalcds[0][0];";
					$finalcodex[$kcod2][18] = "$finalcds[$finalcds[0][0]][1];";
					$finalcodex[$kcod2][19] = $genesarr[$kgb1][10];#refgeno 
					$finalcodex[$kcod2][20] = $genesarr[$kgb1][11];#gene size
					$finalcodex[$kcod2][21] = "HGNC:$geneinfoarr[$kgb1][4]";
					$finalcodex[$kcod2][22] = "=hyperlink\(\"http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=$geneinfoarr[$kgb1][4]\";\"HgncLink:$geneinfoarr[$kgb1][1]\"\)";
					$finalcodex[$kcod2][23] = "Entrez:$geneinfoarr[$kgb1][3]";
					$finalcodex[$kcod2][24] = "=hyperlink\(\"https://www.ncbi.nlm.nih.gov/gene/$geneinfoarr[$kgb1][3]\";\"EntrezLink:$geneinfoarr[$kgb1][1]\"\)";
					$finalcodex[$kcod2][25] = $tempcodex[$kcod1][16];} #last = missed previously
				++$kcod1;}
			++$kc1;}
		my $finalcodexcount = $finalcodex[0][0];
		#OK lets us deal with all unique codex
		#giving these codex their final ID (low number = ATG high number = Stop)
		my @rankcodex;
		my $kcod1 = 1;
		while ($kcod1 <= $finalcodex[0][0]) {
			my $kcod2 = 0;
			$rankcodex[$kcod1-1] = "";
			while ($kcod2 <= 25) {
				$rankcodex[$kcod1-1] = $rankcodex[$kcod1-1] . "$finalcodex[$kcod1][$kcod2]\t";
				++$kcod2;}
			++$kcod1;}
		undef @finalcodex;#going to be rebuilt
		#now ranking by in gene posin and posout (to give ID number starting from 5' of the gene)
		my $field1w = 4;
		my $field2w = 5;
		my @rankarrw;
		@rankarrw = map {$_->[0]} 
		sort {
			my @a_fields = @$a[1..$#$a];
			my @b_fields = @$b[1..$#$b];
			$a_fields[$field1w-1] <=> $b_fields[$field1w-1]
				||
			$a_fields[$field2w-1] <=> $b_fields[$field2w-1]}
		map {[$_, split /\t/]} @rankcodex;
		undef @rankcodex;
		#now we have all in hands to deal with codex and CDS
		$kcod1 = 0;
		while (defined($rankarrw[$kcod1])) {
			++$finalcodex[0][0];
			my $refcod = $kcod1 +1;
			my @linearr;
			while ($rankarrw[$kcod1] =~ m/([^\t]+)\t/igc) {
				$linearr[++$linearr[0]] = $1;
				$finalcodex[$finalcodex[0][0]][$linearr[0]] = $1;}
			if ($linearr[0] != 26) {
				die "Oups: bad number of fields $linearr[0] != 26 for $genesarr[$kgb1][1] in Codex $refcod = $rankarrw[$kcod1]\n";}
			#adding codex ID to respective CDS it is used in (using cds ref list in field 18 (in 7 of CDS final)
			my $listofcdsnum = $linearr[18];
			while ($listofcdsnum =~ m/(\d+)\;/igc) {
				$finalcds[$1][20] = $finalcds[$1][20] . "$genesarr[$kgb1][1]\_CODEX_$refcod;";}
			#printing in Codex table, Codex Blast, Codex Fasta, Codex Bed
			#Fasta Codex Unique ID and info as Chr:Posin-Posout_Strand_Refgeno
			print CODEXFASTA "\>$genesarr[$kgb1][1]\_CODEX_$refcod chr$linearr[9]\:$linearr[13]\-$linearr[14]\_$linearr[12]\_$linearr[20]\n";
			print CODEXFASTA "$linearr[2]\n";
			#Blast Codex Unique ID ... But sequence is Blast one
			my $xstrand = $genesarr[$kgb1][9];
			if ($xstrand eq "+") {
				$xstrand = "P";}
			else {
				$xstrand = "M";}
			print CODEXBLAST "\>$genesarr[$kgb1][1]\_CODEX_$refcod\_$linearr[13]\_$linearr[14]\_BLASTADD_$linearr[26] GENE_chr$genesarr[$kgb1][6]\_$genesarr[$kgb1][7]\_$genesarr[$kgb1][8]\_$xstrand\n";
			$linearr[3] = uc $linearr[3];
			print CODEXBLAST "$linearr[3]\n";
			#Bed codex (still ranked by gene name)
			print CODEXBED "chr$linearr[9]\t$linearr[13]\t$linearr[14]\t$genesarr[$kgb1][1]\_CODEX_$refcod\_chr$linearr[9]\_$linearr[13]\_$linearr[14]\_$linearr[12]\_$linearr[20]\n";
			#Codex Table (here below fields)
			#print CODEXTABLE "RankingCodexID\tCodexID\tRefGeno\tRankingChr\tChr\tStrand";
			#print CODEXTABLE \tInGenePosIN\tInGenePosOUT\tRefGenoPosIN\tRefGenoPosOUT\tSize\tSequence";
			#print CODEXTABLE \tCodexUcscID\tCodexUcscLink";
			#print CODEXTABLE \tNumProductList\tIDproductList\tGeneName\tGeneFullName\tGeneType\tGenePosIN";
			#print CODEXTABLE "\tGenePosOUT\tGeneSize\tGeneProductCount\tGeneCodexCount";
			#print CODEXTABLE "\tGeneHgnc\tGeneHgncLink\tGeneEntrez\tGeneEntrezLink\tGeneUcsc\tGeneUcscLink\tGeneMim\tGeneMimLink\n";
			my $rkcdx;
			if ($refcod < 10) {
				$rkcdx = "00$refcod";}
			elsif ($refcod < 100) {
				$rkcdx = "0$refcod";}
			else {
				$rkcdx = $refcod;}
			print CODEXTABLE "$genesarr[$kgb1][1]\_CODEX_$refcod\t$genesarr[$kgb1][1]\t$refcod\t$genesarr[$kgb1][10]\t$linearr[8]\tchr$linearr[9]\t$linearr[12]";
			print CODEXTABLE "\t$linearr[4]\t$linearr[5]\t$linearr[13]\t$linearr[14]\t$linearr[6]\t$linearr[2]";
			print CODEXTABLE "\tchr$linearr[9]\:$linearr[13]\-$linearr[14]\t=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr$linearr[9]\:$linearr[13]\-$linearr[14]\";\"UCSC_CODEX_LINK\"\)";
			print CODEXTABLE "\t$linearr[18]\t$linearr[19]\t$geneinfoarr[$kgb1][1]\tCDS\t$genesarr[$kgb1][7]";
			print CODEXTABLE "\t$genesarr[$kgb1][8]\t$genesarr[$kgb1][11]\t$finalcds[0][0]\t$finalcodexcount";
			print CODEXTABLE "\tHGNC:$geneinfoarr[$kgb1][4]\t=hyperlink\(\"http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=$geneinfoarr[$kgb1][4]\";\"HGNC_GENE_LINK\"\)";
			print CODEXTABLE "\tEntrez:$geneinfoarr[$kgb1][3]\t=hyperlink\(\"https://www.ncbi.nlm.nih.gov/gene/$geneinfoarr[$kgb1][3]\";\"Entrez_GENE_LINK\"\)";
			print CODEXTABLE "\tUCSC:chr$linearr[9]\:$genesarr[$kgb1][7]\-$genesarr[$kgb1][8]\t=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr$linearr[9]\:$genesarr[$kgb1][7]\-$genesarr[$kgb1][8]\";\"UCSC_GENE_LINK\"\)";
			print CODEXTABLE "\tOMIM:$geneinfoarr[$kgb1][5]\t=hyperlink\(\"http://omim.org/entry/$geneinfoarr[$kgb1][5]\";\"OMIM_GENE_LINK\"\)\n";
			$finalfinalcdx[$pointfinalfinalcdx] = "$genesarr[$kgb1][1]\_CODEX_$refcod\t$genesarr[$kgb1][1]\t$refcod\t$genesarr[$kgb1][10]\t$linearr[8]\tchr$linearr[9]\t$linearr[12]";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\t$linearr[4]\t$linearr[5]\t$linearr[13]\t$linearr[14]\t$linearr[6]\t$linearr[2]";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tchr$linearr[9]\:$linearr[13]\-$linearr[14]\t=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr$linearr[9]\:$linearr[13]\-$linearr[14]\";\"UCSC_CODEX_LINK\"\)";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\t$linearr[18]\t$linearr[19]\t$geneinfoarr[$kgb1][1]\tCDS\t$genesarr[$kgb1][7]";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\t$genesarr[$kgb1][8]\t$genesarr[$kgb1][11]\t$finalcds[0][0]\t$finalcodexcount";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tHGNC:$geneinfoarr[$kgb1][4]\t=hyperlink\(\"http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=$geneinfoarr[$kgb1][4]\";\"HGNC_GENE_LINK\"\)";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tEntrez:$geneinfoarr[$kgb1][3]\t=hyperlink\(\"https://www.ncbi.nlm.nih.gov/gene/$geneinfoarr[$kgb1][3]\";\"Entrez_GENE_LINK\"\)";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tUCSC:chr$linearr[9]\:$genesarr[$kgb1][7]\-$genesarr[$kgb1][8]\t=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr$linearr[9]\:$genesarr[$kgb1][7]\-$genesarr[$kgb1][8]\";\"UCSC_GENE_LINK\"\)";
			$finalfinalcdx[$pointfinalfinalcdx] = $finalfinalcdx[$pointfinalfinalcdx] . "\tOMIM:$geneinfoarr[$kgb1][5]\t=hyperlink\(\"http://omim.org/entry/$geneinfoarr[$kgb1][5]\";\"OMIM_GENE_LINK\"\)";
			++$pointfinalfinalcdx;
			++$kcod1;}
		#Now CDS TABLE and FASTA
		my $kcds1 = 1;
		while ($kcds1 <= $finalcds[0][0]) {
			print CDSFASTA "\>$genesarr[$kgb1][1]_CDS_$kcds1 $finalcds[$kcds1][1]\_$finalcds[$kcds1][4]\n$finalcds[$kcds1][3]\n";
			++$kcds1;}
		print " CDS OK = $finalcds[0][0] and Codex = $finalcodex[0][0]\n";}
	++$kgb1;}
close CDSFASTA;
close CODEXTABLE;
close CODEXFASTA;
close CODEXBLAST;
close CODEXBED;
close ERROR;
print LOG "FINAL CODEX RECORDED\t$pointfinalfinalcdx\n";
if ($cterrors != 0) {
	print LOG "STATUS_GB_FILES\tFAILED\n";
	close LOG;
	copymaintobackup($mainoutputdir,$backupoutputdir);
	die "Some errors ($cterrors) have been identified for $errorlocation. Check error file. $!\n";}
else {
	print LOG "STATUS_GB_FILES\tPASSED\n";
	close LOG;
	copymaintobackup($mainoutputdir,$backupoutputdir);
	open (ERROR,">>$errorfile") ||
		die "Cannot open the error report file $errorfile $!\n";
	open (LOG,">>$runlog") ||
		die "Cannot open LOG file $runlog $!\n";}
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
################################################################################
#big problem with Blast
#Query= HRAS_CODEX_1_534212_534322_BLASTADD_30_30 GENE_chr11_532242_535567_M = too short not 80 chars+
        #677968693279877869823280736982826932766967796781326669765357485550515055495154#
#Query= HPRT1_CODEX_9_134500030_134500077_BLASTADD_30_30 GENE_chrX_134460145_134500668_P = OK 80+
#to solve check length of ID and add chars up to 80 or have a cross check about next line empty ;-)
#now feeding genes synonym
my $field1s = 1;
my @rankarrs;
@rankarrs = map {$_->[0]} 
sort {
	my @a_fields = @$a[1..$#$a];
	my @b_fields = @$b[1..$#$b];
	$a_fields[$field1s-1] cmp $b_fields[$field1s-1]}
map {[$_, split /\t/]} @genesynonym;
undef @genesynonym;
print SYN "GeneName\tGeneHgncName\tHgncRef\tEntrezRef\n";
my $ksyn = 0;
while (defined($rankarrs[$ksyn])) {
	print SYN "$rankarrs[$ksyn]\n";
	++$ksyn;}
close SYN;
print LOG "GENES_SYNONYMS_STATUS\tPASSED\n"; 
#last but not least: the blast of codex to identify sequence similarities AND summarization by gene
#Leonor commandline for blast
#/src/ncbi-blast-2.5.0+/bin/blastn -query CodexSeq_fromPierre.fasta -outfmt "6" -db hg19 -max_target_seqs 5 -num_threads 4 > CodexSeq_fromPierre_hg19.blast 
my $outputblast  = "$mainoutputdir\/output.blast";
system "blastn -query $codexblast -outfmt 0 -db $blastdbgeneric -perc_identity $blastpctidentity -num_alignments $blastmaxhits -max_hsps $blastmaxhits -num_threads 4 -evalue $blastevalue > $outputblast";
my @blastarr;
my $summaryblast  = "$mainoutputdir\/summary.blast";
open (SUMB,">>$summaryblast") ||
	die "Unable to create Summary Blast Output file $summaryblast";
open (OUTB,"<$outputblast") ||
	die "Unable to read from Blast Output file $outputblast";
my @xcodexarr;
my $hitdb;
#is to capture matches on same chromosome
my $firstsubject = 0;
my $keepquery;
my $keepinfo;
my $keepsizequery;
my $keepsubject;
my $usethissubject = 1;
FETCHBLAST:
while (defined(my $linein = <OUTB>)) {
	chomp $linein;
	if ($linein =~ m/^Database\:\s+([^\s]+)/) {
		$hitdb = $1;}
		#Database: hg38.fa
	if ($linein =~ m/^Query\=\s*([^\s]+)\s*([^\s]*)$/) {
		#modified for short ID + Info and also modified for more then pone match per chr
		#Query= A2ML1_CODEX_1_chr12_8822652_8822713
        #BlastAt21_chr12:8822652-8822713_+_GRCh38.p7
		#$blastarr[++$blastarr[0][0][0]][0][1] = $1;#ID of query
		$keepquery = $1;
		my $putativeinfo = $2;
		if ($putativeinfo =~ m/[A-Za-z0-9]/) {
			$keepinfo = $putativeinfo;}
			#$blastarr[$blastarr[0][0][0]][0][2] = $putativeinfo;}
		else {
			$linein = <OUTB>;
			chomp $linein;
			$keepinfo = $linein;}}
			#$blastarr[$blastarr[0][0][0]][0][2] = $linein;}}#info of query
			#Length=104
	elsif ($linein =~ m/^Length\=(\d+)/i) {
		$keepsizequery = $1;}
#		$blastarr[$blastarr[0][0][0]][0][3] = $1;}#size of query
	else {
		#relates to subject hits
		if ($linein =~ m/^\>\s*([^\s]+)/) {
			#>chr12
			my $valsubject = $1;
			if ($valsubject !~ m/^chr[0-9A-Za-z]{1,2}$/i) {
				$usethissubject = 0;}
			else {
				$keepsubject = $valsubject;
				$usethissubject = 1;}}
		if (($linein =~ m/^\s+Score\s+\=\s+([^\s]+)\s+bits\s+\((\d+)\)\s*\,\s+Expect\s+\=\s+([^\s]+)/i) && ($usethissubject == 1)) {
			# Score = 193 bits (104),  Expect = 2e-47
			$blastarr[++$blastarr[0][0][0]][0][1] = $keepquery;#ID of query
			$blastarr[$blastarr[0][0][0]][0][2] = $keepinfo;
			$blastarr[$blastarr[0][0][0]][0][3] = $keepsizequery;
			++$blastarr[$blastarr[0][0][0]][0][4];#
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][1] = $keepsubject;
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][2] = $1;#bit score
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][3] = $2;#length of alignement
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][4] = $3;#Evalue
			#next to initialization of query & subject sequence and coordinates
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][13] = "";#query sequence
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][14] = "";#query join coordinates low
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][15] = "";#query join coordinates high
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][16] = "";#subject sequence
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][17] = "";#subject join coordinates low
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][18] = "";}#subject join coordinates high}#ratio gap
		if (($linein =~ m/^\s+Identities\s+\=\s+(\d+)\/(\d+)\s+.+\,\s+Gaps\s+\=\s+(\d+)\/(\d+)\s/i) && ($usethissubject == 1)){
			# Identities = 104/104 (100%), Gaps = 0/104 (0%)
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][5] = $1;#identities
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][6] = $2;#over...
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][7] = $1/$2;#ratio identities
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][8] = $3;#gaps
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][9] = $4;#over...
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][10] = $3/$4;}
		if (($linein =~ m/^\s+Strand\=([A-Za-z]+)\/([A-Za-z]+)/i) && ($usethissubject == 1)){
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][11] = $1;#Query Strand
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][12] = $2;}#Subject Strand
		if (($linein =~ m/^Query\s+(\d+)\s+([^\s]+)\s+(\d+)/i) && ($usethissubject == 1)) {
			#Query  1        AAATCTGTCTCACCCACAAAGATGTGGGCTCAGCTCCTTCTAGGAATGTTGGCCCTATCA  60
			#    			 ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
			#Sbjct  8822631  AAATCTGTCTCACCCACAAAGATGTGGGCTCAGCTCCTTCTAGGAATGTTGGCCCTATCA  8822690
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][13] = $blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][13] . uc $2;#query sequence
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][14] = $blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][14] . "$1;";#query join coordinates low
			$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][15] = $blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][15] . "$3;";#query join coordinates high
			$linein = <OUTB>;
			$linein = <OUTB>;
			chomp $linein;
			if ($linein =~ m/^Sbjct\s+(\d+)\s+([^\s]+)\s+(\d+)/i) {
				$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][16] = $blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][16] . uc $2;#subject sequence
				$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][17] = $blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][17] . "$1;";#subject join coordinates low
				$blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][18] = $blastarr[$blastarr[0][0][0]][$blastarr[$blastarr[0][0][0]][0][4]][18] . "$3;";}}}}#subject join coordinates high
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#here all relevant info in blast output is captured
#lets review and print in a blast summary file
print SUMB "EntryNum\tCodexUniqueID\tGeneName\tCodexNum\tCodexPosIN\tCodexPosOUT\tBlastBefore\tBlastAfter";
#                0         1             2         3         4          5            6             7
print SUMB "\tchr\tGenePosIN\tGenePosOUT\tGeneStrand";
#              8     9          10           11
print SUMB "\tHitID\tHitChr\tTypeHitID";
#                12    13         14
print SUMB "\tBitScore\tSizeAlign\tEvalue\tNbrIdentities\tIdentitiesOver\tPctIdentities\tNbrGaps\tGapsOver\tPctGaps";
#                15          16        17        18         19              20            21           22     23
print SUMB "\tHitPosIN\tHitPosOUT\tHitStrand";
#                24        25         26
print SUMB "\tRefGeno\tHitDB\tMaxTm\tFlag\tColor\tRemark\tUcscGeneID\tUcscGeneLink\tUcscCodexID\tUcscCodexLink\tUcscHitID\tUcscHitLink\n";
#               27       28     29     30   31      32       33         34             35         36              37        38
my @genesumblast;
my $kob1 = 1;
while ($kob1 <= $blastarr[0][0][0]) {
	my $kob2 = 1;
	while ($kob2 <= $blastarr[$kob1][0][4]) {
		my $kob3 = 1;
		#print SUMB "$blastarr[$kob1][0][1]\t$blastarr[$kob1][0][2]\t$kob2";
		my $seqq = $blastarr[$kob1][$kob2][13];
		my $seqs = $blastarr[$kob1][$kob2][16];
		if (length($seqq) != length($seqs)) {
			die "Perl code badly written to recover sequence from blast output $seqq vs $seqs $!\n";}
		my $longuestmatch = 0;
		my $longueststring = "";
		my $ctlonguest = 0;
		my @calcultm;
		$calcultm[0] = 0;
		my $maxtmvalue = 0;
		while ($seqq =~ m/(.)/igc) {
			my $nucq = $1;
			my $nucs;
			if ($seqs =~ /^(.)(.*)$/) {
				$nucs = $1;
				$seqs = $2;}
			if ($nucs ne $nucq) {
				if ($calcultm[0] >= $blastminwinsize) {#minimum size to start evaluation (windowed up to 26 values)
					if ($calcultm[0] <= $blastmaxwinsize) {#so using everything
						my $kcalc = 1;
						my $tmvalue = 0;
						while (defined($calcultm[$kcalc])) {
							$tmvalue = $tmvalue + $calcultm[$kcalc];
							++$kcalc;}
						if ($tmvalue > $maxtmvalue) {
							$maxtmvalue = $tmvalue;}}
					else {#so longer then $blastmaxwinsize so windowed approach at 26
						my $kcalc1 = $blastmaxwinsize;
						while (defined($calcultm[$kcalc1])) {
							my $kcalc2 = $kcalc1 - $blastmaxwinsize + 1;
							my $tmvalue = 0;
							while ($kcalc2 <= $kcalc1) {
								$tmvalue = $tmvalue + $calcultm[$kcalc2];
								++$kcalc2;}
							if ($tmvalue > $maxtmvalue) {
								$maxtmvalue = $tmvalue;}
							++$kcalc1;}}}
				undef @calcultm;
				$calcultm[0] = 0;}
			else {
				$calcultm[++$calcultm[0]] = $tmhash{$nucs};}}
		if ($calcultm[0] >= $blastminwinsize) {#minimum size to start evaluation (windowed up to 26 values)
			if ($calcultm[0] <= $blastmaxwinsize) {#so using everything
				my $kcalc = 1;
				my $tmvalue = 0;
				while (defined($calcultm[$kcalc])) {
					$tmvalue = $tmvalue + $calcultm[$kcalc];
					++$kcalc;}
				if ($tmvalue > $maxtmvalue) {
					$maxtmvalue = $tmvalue;}}
			else {#so longer then 26 so windowed approach at 26
				my $kcalc1 = $blastmaxwinsize;
				while (defined($calcultm[$kcalc1])) {
					my $kcalc2 = $kcalc1 - $blastminwinsize +1;
					my $tmvalue = 0;
					while ($kcalc2 <= $kcalc1) {
						$tmvalue = $tmvalue + $calcultm[$kcalc2];
						++$kcalc2;}
					if ($tmvalue > $maxtmvalue) {
						$maxtmvalue = $tmvalue;}
						++$kcalc1;}}}
			undef @calcultm;
			$calcultm[0] = 0;
			$blastarr[$kob1][$kob2][20] = $maxtmvalue;
			#new = review, print and store in arrays for adding to codex / gene tables
			my $idvalue = $blastarr[$kob1][0][1];
			if ($idvalue =~ m/^(([^\_]+)_CODEX_(\d+))_(\d+)_(\d+)_BLASTADD_(\d+)_(\d+)$/i) {
				#A2ML1_CODEX_1_8822652_8822713_BLASTADD_30_30
				#print SUMB "CodexUniqueID\tGeneName\tCodexNum\tCodexPosIN\tCodexPosOUT\tBlastBefore\tBlastAfter"; 
				$xcodexarr[++$xcodexarr[0][0]][1] = $1;#codex unique ID
				$xcodexarr[$xcodexarr[0][0]][2] = $2;#codex genename
				$xcodexarr[$xcodexarr[0][0]][3] = $3;#codex num
				$xcodexarr[$xcodexarr[0][0]][4] = $4;#codex posin
				$xcodexarr[$xcodexarr[0][0]][5] = $5;#codex posout
				$xcodexarr[$xcodexarr[0][0]][6] = $6;#codex blast nuc added before
				$xcodexarr[$xcodexarr[0][0]][7] = $7;}#codex blast nuc added after
			else {
				die "Bad code processing for id value $idvalue (Blast output 1)";}
			my $infovalue = $blastarr[$kob1][0][2];
			if ($infovalue =~ m/^GENE_(chr[^\_]+)_(\d+)_(\d+)_([PM])$/i) {
				#GENE_chr12_8822472_8887202_P
				#print SUMB "\tchr\tGenePosIN\tGenePosOUT\tGeneStrand";
				$xcodexarr[$xcodexarr[0][0]][8] = $1;#gene chr
				$xcodexarr[$xcodexarr[0][0]][9] = $2;#gene posin
				$xcodexarr[$xcodexarr[0][0]][10] = $3;#gene posout
				$xcodexarr[$xcodexarr[0][0]][11] = uc $4;#gene strand
				if ($xcodexarr[$xcodexarr[0][0]][11] eq "P") {
					$xcodexarr[$xcodexarr[0][0]][11] = "+";}
				else {
					$xcodexarr[$xcodexarr[0][0]][11] = "-";}}
			else {
				die "Bad code processing for info value $infovalue (Blast output 1)";}
			$xcodexarr[$xcodexarr[0][0]][12] =  $blastarr[$kob1][$kob2][1]; #Hit ID
			if ($xcodexarr[$xcodexarr[0][0]][12] =~ m/^(chr[^\_]+)\_/i) {
				 $xcodexarr[$xcodexarr[0][0]][13] = $1;
				 $xcodexarr[$xcodexarr[0][0]][14] = 0;}
			else {
				 $xcodexarr[$xcodexarr[0][0]][13] =  $xcodexarr[$xcodexarr[0][0]][12];
				 $xcodexarr[$xcodexarr[0][0]][14] = 1;}
			#print SUMB "\tHitID\tHitChr\tTypeHitID";
			#                  12    13      14
			#print SUMB "\tBitScore\tSizeAlign\tEvalue\tNbrIdentities\tIdentitiesOver\tPctIdentities\tNbrGaps\tGapsOver\tPctGaps";
			#                 15        16        17       18            19             20             21           22    23
			my $kload = 2;
			while ($kload <= 10) {
				$xcodexarr[$xcodexarr[0][0]][$kload+13] = $blastarr[$kob1][$kob2][$kload];
				++$kload;} 
			#hit coordinates Pointer start here at 24
			#print SUMB "\tHitPosIN\tHitPosOUT\tHitStrand";
			#                  24        25       26
			my $joinstart = $blastarr[$kob1][$kob2][17];
			my $joinend = $blastarr[$kob1][$kob2][18];
			my $subjectstrand = $blastarr[$kob1][$kob2][12];
		if ($joinstart =~ m/^(\d+)\;/) {
			$joinstart = $1;}
		if ($joinend =~ m/(\d+)\;$/) {
			$joinend = $1;}
		if (uc $subjectstrand eq "PLUS") {
			$xcodexarr[$xcodexarr[0][0]][26] = "+";
			$xcodexarr[$xcodexarr[0][0]][24] = $joinstart;
			$xcodexarr[$xcodexarr[0][0]][25] = $joinend;}
		else {
			$xcodexarr[$xcodexarr[0][0]][26] = "-";
			$xcodexarr[$xcodexarr[0][0]][25] = $joinstart;
			$xcodexarr[$xcodexarr[0][0]][24] = $joinend;}
		#print SUMB "\tRefGeno\tHitDB\tHitFlag\tHitMaxTm\tHitRemark\n";
		#                 27      28   29      30           31
		$xcodexarr[$xcodexarr[0][0]][27] = $keeprefgeno;
		$xcodexarr[$xcodexarr[0][0]][28] = $hitdb;
		$xcodexarr[$xcodexarr[0][0]][29] = $blastarr[$kob1][$kob2][20];
		if ($xcodexarr[$xcodexarr[0][0]][29] > 72) {
			$xcodexarr[$xcodexarr[0][0]][29] = 72;}
		if ($xcodexarr[$xcodexarr[0][0]][29] < 48) {
			$xcodexarr[$xcodexarr[0][0]][29] = 48;}
		#now evaluation of hit to check if same gene / exon or something else
		if (($xcodexarr[$xcodexarr[0][0]][13] eq $xcodexarr[$xcodexarr[0][0]][8])  && ($xcodexarr[$xcodexarr[0][0]][11] eq $xcodexarr[$xcodexarr[0][0]][26])) {#corresponding chr and strand)
			my $realcdxstart = $xcodexarr[$xcodexarr[0][0]][24] + $xcodexarr[$xcodexarr[0][0]][6];
			my $realcdxend =  $xcodexarr[$xcodexarr[0][0]][25] - $xcodexarr[$xcodexarr[0][0]][7];
			if ($xcodexarr[$xcodexarr[0][0]][26] eq "-") {
				$realcdxstart = $xcodexarr[$xcodexarr[0][0]][24] + $xcodexarr[$xcodexarr[0][0]][7];
				$realcdxend =  $xcodexarr[$xcodexarr[0][0]][25] - $xcodexarr[$xcodexarr[0][0]][6];}
			if (($realcdxstart == $xcodexarr[$xcodexarr[0][0]][4]) && ($realcdxend == $xcodexarr[$xcodexarr[0][0]][5])) {#exactly expected coordinates
				if ($xcodexarr[$xcodexarr[0][0]][20] == 1) {#so perfect
					$xcodexarr[$xcodexarr[0][0]][30] = 0;
					$xcodexarr[$xcodexarr[0][0]][31] = "Violet";
					$xcodexarr[$xcodexarr[0][0]][32] = "Perfect!";}
				else {#bizarre ??? Error GB vs Blast???
					$xcodexarr[$xcodexarr[0][0]][30] = 1;
					$xcodexarr[$xcodexarr[0][0]][31] = "Blue";
					$xcodexarr[$xcodexarr[0][0]][32] = "Unexpected mismatch(es)!";}}
			else {#so same chromosome and strand but not perfectly the same coordinates
				if (($realcdxstart >= $xcodexarr[$xcodexarr[0][0]][4]) && ($realcdxend <= $xcodexarr[$xcodexarr[0][0]][5])) {
					#probable split of codex during the blast but inner part 
					$xcodexarr[$xcodexarr[0][0]][30] = 1;
					$xcodexarr[$xcodexarr[0][0]][31] = "Blue";
					$xcodexarr[$xcodexarr[0][0]][32] = "Split of Codex!";}
				else {#so another location and using tm to categorise
					if ($xcodexarr[$xcodexarr[0][0]][29] >= $blastred) {
						$xcodexarr[$xcodexarr[0][0]][30] = 5;
						$xcodexarr[$xcodexarr[0][0]][31] = "Red";
						$xcodexarr[$xcodexarr[0][0]][32] = "To reject!";}
					elsif ($xcodexarr[$xcodexarr[0][0]][29] >= $blastorange) {
						$xcodexarr[$xcodexarr[0][0]][30] = 4;
						$xcodexarr[$xcodexarr[0][0]][31] = "Orange";
						$xcodexarr[$xcodexarr[0][0]][32] = "To review carefully!";}
					elsif ($xcodexarr[$xcodexarr[0][0]][29] >= $blastyellow) {
						$xcodexarr[$xcodexarr[0][0]][30] = 3;
						$xcodexarr[$xcodexarr[0][0]][31] = "Yellow";
						$xcodexarr[$xcodexarr[0][0]][32] = "May require attention!";}
					else {
						$xcodexarr[$xcodexarr[0][0]][30] = 2;
						$xcodexarr[$xcodexarr[0][0]][31] = "Green";
						$xcodexarr[$xcodexarr[0][0]][32] = "No problem!";}}}}
		 else {#so definitively another location so another gene and use of Tm to categorize
			 if ($xcodexarr[$xcodexarr[0][0]][29] >= $blastred) {
				$xcodexarr[$xcodexarr[0][0]][30] = 5;
				$xcodexarr[$xcodexarr[0][0]][31] = "Red";
				$xcodexarr[$xcodexarr[0][0]][32] = "To reject!";}
			elsif ($xcodexarr[$xcodexarr[0][0]][29] >= $blastorange) {
				$xcodexarr[$xcodexarr[0][0]][30] = 4;
				$xcodexarr[$xcodexarr[0][0]][31] = "Orange";
				$xcodexarr[$xcodexarr[0][0]][32] = "To review carefully!";}
			elsif ($xcodexarr[$xcodexarr[0][0]][29] >= $blastyellow) {
				$xcodexarr[$xcodexarr[0][0]][30] = 3;
				$xcodexarr[$xcodexarr[0][0]][31] = "Yellow";
				$xcodexarr[$xcodexarr[0][0]][32] = "May require attention!";}
			else {
				$xcodexarr[$xcodexarr[0][0]][30] = 2;
				$xcodexarr[$xcodexarr[0][0]][31] = "Green";
				$xcodexarr[$xcodexarr[0][0]][32] = "No problem!";}}
		#UCSC http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr1:2228125-2282759
		#Works also http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr11_KI270832v1_alt:206431-206738
		$xcodexarr[$xcodexarr[0][0]][33] = "$xcodexarr[$xcodexarr[0][0]][8]\:$xcodexarr[$xcodexarr[0][0]][9]\-$xcodexarr[$xcodexarr[0][0]][10]";
		$xcodexarr[$xcodexarr[0][0]][34] = "=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=$xcodexarr[$xcodexarr[0][0]][33]\";\"UCSC_GENE_LINK\"\)";
		$xcodexarr[$xcodexarr[0][0]][35] = "$xcodexarr[$xcodexarr[0][0]][8]\:$xcodexarr[$xcodexarr[0][0]][4]\-$xcodexarr[$xcodexarr[0][0]][5]";
		$xcodexarr[$xcodexarr[0][0]][36] = "=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=$xcodexarr[$xcodexarr[0][0]][35]\";\"UCSC_CODEX_LINK\"\)"; 
		$xcodexarr[$xcodexarr[0][0]][37] = "$xcodexarr[$xcodexarr[0][0]][12]\:$xcodexarr[$xcodexarr[0][0]][24]\-$xcodexarr[$xcodexarr[0][0]][25]";
		$xcodexarr[$xcodexarr[0][0]][38] = "=hyperlink\(\"http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=$xcodexarr[$xcodexarr[0][0]][37]\";\"UCSC_HIT_LINK\"\)"; 
		print SUMB "$xcodexarr[0][0]";
		my $kload3 = 1;		
		while (defined($xcodexarr[$xcodexarr[0][0]][$kload3] )){
			print SUMB "\t$xcodexarr[$xcodexarr[0][0]][$kload3]";
			++$kload3;}
		print SUMB "\n";
		++$kob2;}
	++$kob1;}
close SUMB;
print LOG "BLAST FOR HOMOLOGS\tSTAGE 1 COMPLETED\n";
print LOG "BLAST SUMMARY FILE\t$summaryblast\n";
print LOG "BLAST HITS\t$xcodexarr[0][0]\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#now reviewing and summarizing per exon and per gene.
my @exonblast;
my @tempexonb;
my @tempgeneb;
my $novalidblast = 0;
my $validblast = 0;
my $krx1 = 1;
FORSKIP:
while ($krx1 <= $xcodexarr[0][0]) {
	if ((not defined($tempexonb[0][0])) || ($xcodexarr[$krx1][1] eq $tempexonb[1][1])) {#first entry or same exon and not last entry of array
		++$tempexonb[0][0];
		my $krx2 = 1;
		while (defined($xcodexarr[$krx1][$krx2])) {
			$tempexonb[$tempexonb[0][0]][$krx2] = $xcodexarr[$krx1][$krx2];
			++$krx2;}}
	else {
		#processing the exon data
		#Loading values ready for ranking
		my @rankarrk;
		my $kk1 = 1;
		while ($kk1 <= $tempexonb[0][0]) {
			my $point = $kk1 -1;
			my $kk2 = 2;
			$rankarrk[$point] = "$tempexonb[$kk1][1]\t";
			while (defined($tempexonb[$kk1][$kk2])) {
				$rankarrk[$point] = $rankarrk[$point] . "$tempexonb[$kk1][$kk2]\t";
			++$kk2;}
		++$kk1;}
		#cases with only one entry
		if (not defined($rankarrk[1])) {
			if (($rankarrk[0] !~ m/\tVIOLET\t/i) && ($rankarrk[0] !~ m/\tBLUE\t/i)) {#this single entry doesnot match within expected coordinates
				#better to skip this one as no blue, no violet
				++$novalidblast;
				++$krx1;
				next FORSKIP;}
			else {
				#exploding entry to recover only interesting values.
				++$validblast;
				my @explode;
				undef @explode;
				my $best = $rankarrk[0];
				while($best =~ m/([^\t]+)\t/igc) {
					$explode[++$explode[0]] = $1;}
				$exonblast[++$exonblast[0][0]][1] = $explode[1];#Exon Unique ID
				$exonblast[$exonblast[0][0]][2] = $explode[2];#genename
				if ($explode[1] =~ m/CODEX\_(\d+)$/) {
					$exonblast[$exonblast[0][0]][3] = $1;}#exon number
				$exonblast[$exonblast[0][0]][4] = 1;#overall exon Flag
				$exonblast[$exonblast[0][0]][5] = "NA";#flag of worst hit
				$exonblast[$exonblast[0][0]][6] = "NA";#Evalue of worst hit
				$exonblast[$exonblast[0][0]][7] = "NA";#Pct Identity of worst hit
				$exonblast[$exonblast[0][0]][8] = "NA";#UCSC ID of worst hit
				$exonblast[$exonblast[0][0]][9] = "NA";#UCSC Link of worst hit
				$exonblast[$exonblast[0][0]][10] = "NA";#flag of second worst hit 			
				$exonblast[$exonblast[0][0]][11] = "NA";#Evalue of second worst hit
				$exonblast[$exonblast[0][0]][12] = "NA";#Pct Identity second of worst hit
				$exonblast[$exonblast[0][0]][13] = "NA";#UCSC ID of second worst hit
				$exonblast[$exonblast[0][0]][14] = "NA";}}#UCSC Link of second worst hit
		else {
			#need first to search for a Violet / Blue
			my $best;
			my $bestref;
			my $kb1 = 0;
			my $foundbest = 0;
			SEARCHBESTV:#looking for perfect match
			while (defined($rankarrk[$kb1])) {
				if ($rankarrk[$kb1] =~ m/\tVIOLET\t/i) {
					$best = $rankarrk[$kb1];
					$bestref = $kb1;
					$foundbest = 1;
					last SEARCHBESTV;}
				++$kb1;}
			if ($foundbest == 0) {#if no perfect looking for blue one
				my $kb2 = 0;
				SEARCHBESTB:
				while (defined($rankarrk[$kb2])) {
					if ($rankarrk[$kb2] =~ m/\tBLUE\t/i) {
						$best = $rankarrk[$kb2];
						$bestref = $kb2;
						$foundbest = 1;
						last SEARCHBESTB;}
					++$kb2;}}
			if ($foundbest == 0) {#better to skip this one as no blue, no violet
				++$novalidblast;
				++$krx1;
				next FORSKIP;}
			++$validblast;
			my @explode;
			undef @explode;
			$best = $rankarrk[$bestref];
			while($best =~ m/([^\t]+)\t/igc) {
				$explode[++$explode[0]] = $1;}
			$exonblast[++$exonblast[0][0]][1] = $explode[1];#Exon Unique ID
			$exonblast[$exonblast[0][0]][2] = $explode[2];#genename
			if ($explode[1] =~ m/CODEX\_(\d+)$/) {
				$exonblast[$exonblast[0][0]][3] = $1;}#exon number
			$exonblast[$exonblast[0][0]][4] = 1;#overall exon Flag means passed perfect and nothign else (for right now)
			#let us remove this entry from array
			my @rankarrk2;
			my $krm1 = 0;
			my $pointb = 0;
			while (defined($rankarrk[$krm1])) {
				if ($krm1 == $bestref) {}
				else {
					$rankarrk2[$pointb] = $rankarrk[$krm1];
					++$pointb;}
				++$krm1;}
			undef @rankarrk;
			#let us rank dicreasing on flag then increasing on evalue
			my $field1k = 30;
			my $field2k = 17;
			my @rankarrk = map {$_->[0]} 
			sort {
				my @a_fields = @$a[1..$#$a];
				my @b_fields = @$b[1..$#$b];
				$b_fields[$field1k-1] <=> $a_fields[$field1k-1]
					||
				$b_fields[$field2k-1] <=> $a_fields[$field2k-1]}
			map {[$_, split /\t/]} @rankarrk2;
			undef @rankarrk2;
			my $wcs1 = $rankarrk[0];
			undef @explode;
			while($wcs1 =~ m/([^\t]+)\t/igc) {
				$explode[++$explode[0]] = $1;}
			my $valfin = $explode[30];
			if ($valfin <= 2) {
				$exonblast[$exonblast[0][0]][4] = 1;}
			else {
				$exonblast[$exonblast[0][0]][4] = $explode[30]-1;}
			$exonblast[$exonblast[0][0]][5] = $valfin;#flag of worst hit
			$exonblast[$exonblast[0][0]][6] = $explode[17];#Evalue of worst hit
			$exonblast[$exonblast[0][0]][7] = $explode[20];#Pct Identity of worst hit
			$exonblast[$exonblast[0][0]][8] = $explode[37];#UCSC ID of worst hit
			$exonblast[$exonblast[0][0]][9] = $explode[38];#UCSC Link of worst hit
			if (not defined($rankarrk[1])) {	
				$exonblast[$exonblast[0][0]][10] = "NA";#flag of second worst hit 			
				$exonblast[$exonblast[0][0]][11] = "NA";#Evalue of second worst hit
				$exonblast[$exonblast[0][0]][12] = "NA";#Pct Identity second of worst hit
				$exonblast[$exonblast[0][0]][13] = "NA";#UCSC ID of second worst hit
				$exonblast[$exonblast[0][0]][14] = "NA";}#UCSC Link of second worst hit
			else {
				my $wcs2 = $rankarrk[1];
				my @explode2;
				undef @explode2;
				while ($wcs2 =~ m/([^\t]+)\t/igc) {
					$explode2[++$explode2[0]] = $1;}
				if ($explode2[30] <= 2) {
					$exonblast[$exonblast[0][0]][10] = 1;#flag of second worst hit
					$exonblast[$exonblast[0][0]][11] = $explode2[17];#Evalue of second worst hit
					$exonblast[$exonblast[0][0]][12] = $explode2[20];#Pct Identity second of worst hit
					$exonblast[$exonblast[0][0]][13] = $explode2[37];#UCS of second worst hit
					$exonblast[$exonblast[0][0]][14] = $explode2[38];}#UCSC Link of second worst hit
				else {
					$exonblast[$exonblast[0][0]][10] = $explode2[30]-1;#flag of second worst hit		
					$exonblast[$exonblast[0][0]][11] = $explode2[17];#Evalue of second worst hit
					$exonblast[$exonblast[0][0]][12] = $explode2[20];#Pct Identity second of worst hit
					$exonblast[$exonblast[0][0]][13] = $explode2[37];#UCS of second worst hit
					$exonblast[$exonblast[0][0]][14] = $explode2[38];}}}#UCSC Link of second worst hit
		#what do we keep?
		#print SUMB "EntryNum\tCodexUniqueID\tGeneName\tCodexNum\tCodexPosIN\tCodexPosOUT\tBlastBefore\tBlastAfter";
		#                0         1             2         3         4          5            6             7
		#print SUMB "\tchr\tGenePosIN\tGenePosOUT\tGeneStrand";
		#              8     9          10           11
		#print SUMB "\tHitID\tHitChr\tTypeHitID";
		#                12    13         14
		#print SUMB "\tBitScore\tSizeAlign\tEvalue\tNbrIdentities\tIdentitiesOver\tPctIdentities\tNbrGaps\tGapsOver\tPctGaps";
		#                15          16        17        18         19              20            21           22     23
		#print SUMB "\tHitPosIN\tHitPosOUT\tHitStrand";
		#                24        25         26
		#print SUMB "\tRefGeno\tHitDB\tMaxTm\tFlag\tColor\tRemark\tUcscGeneID\tUcscGeneLink\tUcscCodexID\tUcscCodexLink\tUcscHitID\tUcscHitLink\n";
		#               27       28     29     30   31      32       33         34             35         36              37        38
		undef @rankarrk;
		undef @tempexonb;
		++$tempexonb[0][0];
		my $krx2 = 1;
		while (defined($xcodexarr[$krx1][$krx2])) {
			$tempexonb[$tempexonb[0][0]][$krx2] = $xcodexarr[$krx1][$krx2];
			++$krx2;}}		
	++$krx1;}
#next are for last codex
if (defined($tempexonb[0][0])) {
	##processing the exon data
		#Loading values ready for ranking
		my @rankarrk;
		my $kk1 = 1;
		while ($kk1 <= $tempexonb[0][0]) {
			my $point = $kk1 -1;
			my $kk2 = 2;
			$rankarrk[$point] = "$tempexonb[$kk1][1]\t";
			while (defined($tempexonb[$kk1][$kk2])) {
				$rankarrk[$point] = $rankarrk[$point] . "$tempexonb[$kk1][$kk2]\t";
			++$kk2;}
		++$kk1;}
		#cases with only one entry
		if (not defined($rankarrk[1])) {
			if (($rankarrk[0] !~ m/\tVIOLET\t/i) && ($rankarrk[0] !~ m/\tBLUE\t/i)) {#this single entry doesnot match within expected coordinates
				#better to skip this one as no blue, no violet
				++$novalidblast;
				++$krx1;
				next FORSKIP;}
			else {
				#exploding entry to recover only interesting values.
				++$validblast;
				my @explode;
				undef @explode;
				my $best = $rankarrk[0];
				while($best =~ m/([^\t]+)\t/igc) {
					$explode[++$explode[0]] = $1;}
				$exonblast[++$exonblast[0][0]][1] = $explode[1];#Exon Unique ID
				$exonblast[$exonblast[0][0]][2] = $explode[2];#genename
				if ($explode[1] =~ m/CODEX\_(\d+)$/) {
					$exonblast[$exonblast[0][0]][3] = $1;}#exon number
				$exonblast[$exonblast[0][0]][4] = 1;#overall exon Flag
				$exonblast[$exonblast[0][0]][5] = "NA";#flag of worst hit
				$exonblast[$exonblast[0][0]][6] = "NA";#Evalue of worst hit
				$exonblast[$exonblast[0][0]][7] = "NA";#Pct Identity of worst hit
				$exonblast[$exonblast[0][0]][8] = "NA";#UCSC ID of worst hit
				$exonblast[$exonblast[0][0]][9] = "NA";#UCSC Link of worst hit
				$exonblast[$exonblast[0][0]][10] = "NA";#flag of second worst hit 			
				$exonblast[$exonblast[0][0]][11] = "NA";#Evalue of second worst hit
				$exonblast[$exonblast[0][0]][12] = "NA";#Pct Identity second of worst hit
				$exonblast[$exonblast[0][0]][13] = "NA";#UCSC ID of second worst hit
				$exonblast[$exonblast[0][0]][14] = "NA";}}#UCSC Link of second worst hit
		else {
			#need first to search for a Violet / Blue
			my $best;
			my $bestref;
			my $kb1 = 0;
			my $foundbest = 0;
			SEARCHBESTV:#looking for perfect match
			while (defined($rankarrk[$kb1])) {
				if ($rankarrk[$kb1] =~ m/\tVIOLET\t/i) {
					$best = $rankarrk[$kb1];
					$bestref = $kb1;
					$foundbest = 1;
					last SEARCHBESTV;}
				++$kb1;}
			if ($foundbest == 0) {#if no perfect looking for blue one
				my $kb2 = 0;
				SEARCHBESTB:
				while (defined($rankarrk[$kb2])) {
					if ($rankarrk[$kb2] =~ m/\tBLUE\t/i) {
						$best = $rankarrk[$kb2];
						$bestref = $kb2;
						$foundbest = 1;
						last SEARCHBESTB;}
					++$kb2;}}
			if ($foundbest == 0) {#better to skip this one as no blue, no violet
				++$novalidblast;
				++$krx1;
				next FORSKIP;}
			++$validblast;
			my @explode;
			undef @explode;
			$best = $rankarrk[$bestref];
			while($best =~ m/([^\t]+)\t/igc) {
				$explode[++$explode[0]] = $1;}
			$exonblast[++$exonblast[0][0]][1] = $explode[1];#Exon Unique ID
			$exonblast[$exonblast[0][0]][2] = $explode[2];#genename
			if ($explode[1] =~ m/CODEX\_(\d+)$/) {
				$exonblast[$exonblast[0][0]][3] = $1;}#exon number
			$exonblast[$exonblast[0][0]][4] = 1;#overall exon Flag means passed perfect and nothign else (for right now)
			#let us remove this entry from array
			my @rankarrk2;
			my $krm1 = 0;
			my $pointb = 0;
			while (defined($rankarrk[$krm1])) {
				if ($krm1 == $bestref) {}
				else {
					$rankarrk2[$pointb] = $rankarrk[$krm1];
					++$pointb;}
				++$krm1;}
			undef @rankarrk;
			#let us rank dicreasing on flag then increasing on evalue
			my $field1k = 30;
			my $field2k = 17;
			my @rankarrk = map {$_->[0]} 
			sort {
				my @a_fields = @$a[1..$#$a];
				my @b_fields = @$b[1..$#$b];
				$b_fields[$field1k-1] <=> $a_fields[$field1k-1]
					||
				$b_fields[$field2k-1] <=> $a_fields[$field2k-1]}
			map {[$_, split /\t/]} @rankarrk2;
			undef @rankarrk2;
			my $wcs1 = $rankarrk[0];
			undef @explode;
			while($wcs1 =~ m/([^\t]+)\t/igc) {
				$explode[++$explode[0]] = $1;}
			my $valfin = $explode[30];
			if ($valfin <= 2) {
				$exonblast[$exonblast[0][0]][4] = 1;}
			else {
				$exonblast[$exonblast[0][0]][4] = $explode[30]-1;}
			$exonblast[$exonblast[0][0]][5] = $valfin;#flag of worst hit
			$exonblast[$exonblast[0][0]][6] = $explode[17];#Evalue of worst hit
			$exonblast[$exonblast[0][0]][7] = $explode[20];#Pct Identity of worst hit
			$exonblast[$exonblast[0][0]][8] = $explode[37];#UCSC ID of worst hit
			$exonblast[$exonblast[0][0]][9] = $explode[38];#UCSC Link of worst hit
			if (not defined($rankarrk[1])) {	
				$exonblast[$exonblast[0][0]][10] = "NA";#flag of second worst hit 			
				$exonblast[$exonblast[0][0]][11] = "NA";#Evalue of second worst hit
				$exonblast[$exonblast[0][0]][12] = "NA";#Pct Identity second of worst hit
				$exonblast[$exonblast[0][0]][13] = "NA";#UCSC ID of second worst hit
				$exonblast[$exonblast[0][0]][14] = "NA";}#UCSC Link of second worst hit
			else {
				my $wcs2 = $rankarrk[1];
				my @explode2;
				undef @explode2;
				while ($wcs2 =~ m/([^\t]+)\t/igc) {
					$explode2[++$explode2[0]] = $1;}
				if ($explode2[30] <= 2) {
					$exonblast[$exonblast[0][0]][10] = 1;#flag of second worst hit
					$exonblast[$exonblast[0][0]][11] = $explode2[17];#Evalue of second worst hit
					$exonblast[$exonblast[0][0]][12] = $explode2[20];#Pct Identity second of worst hit
					$exonblast[$exonblast[0][0]][13] = $explode2[37];#UCS of second worst hit
					$exonblast[$exonblast[0][0]][14] = $explode2[38];}#UCSC Link of second worst hit
				else {
					$exonblast[$exonblast[0][0]][10] = $explode2[30]-1;#flag of second worst hit		
					$exonblast[$exonblast[0][0]][11] = $explode2[17];#Evalue of second worst hit
					$exonblast[$exonblast[0][0]][12] = $explode2[20];#Pct Identity second of worst hit
					$exonblast[$exonblast[0][0]][13] = $explode2[37];#UCS of second worst hit
					$exonblast[$exonblast[0][0]][14] = $explode2[38];}}}#UCSC Link of second worst hit
		#what do we keep?
		#print SUMB "EntryNum\tCodexUniqueID\tGeneName\tCodexNum\tCodexPosIN\tCodexPosOUT\tBlastBefore\tBlastAfter";
		#                0         1             2         3         4          5            6             7
		#print SUMB "\tchr\tGenePosIN\tGenePosOUT\tGeneStrand";
		#              8     9          10           11
		#print SUMB "\tHitID\tHitChr\tTypeHitID";
		#                12    13         14
		#print SUMB "\tBitScore\tSizeAlign\tEvalue\tNbrIdentities\tIdentitiesOver\tPctIdentities\tNbrGaps\tGapsOver\tPctGaps";
		#                15          16        17        18         19              20            21           22     23
		#print SUMB "\tHitPosIN\tHitPosOUT\tHitStrand";
		#                24        25         26
		#print SUMB "\tRefGeno\tHitDB\tMaxTm\tFlag\tColor\tRemark\tUcscGeneID\tUcscGeneLink\tUcscCodexID\tUcscCodexLink\tUcscHitID\tUcscHitLink\n";
		#               27       28     29     30   31      32       33         34             35         36              37        38
		undef @rankarrk;
		undef @tempexonb;}
print LOG "BLAST PROCESSED EXONS\t$exonblast[0][0]\n";
#now matching the source exon table data
#for matching best is to use the ranking codex ID (added position 15 with this value)
#First file to record exon data final
print "VALID BLAST CODEX: $validblast NotOK: $novalidblast\n";
my $finalcodextable  = "$mainoutputdir\/finalcodex.txt";
open (FINALCODEX,">>$finalcodextable") ||
	die "Unable to create Final Codex table $finalcodextable";
print FINALCODEX "CodexID\tCodexGene\tCodexNum\tRefGeno\tRankingChr\tChr\tStrand";
#                    1       2           3         4        5         6     7
print FINALCODEX "\tInGenePosIN\tInGenePosOUT\tRefGenoPosIN\tRefGenoPosOUT\tSize\tSequence";
#                    8              9             10            11             12     13
print FINALCODEX "\tCodexUcscID\tCodexUcscLink";
#                     14              15
print FINALCODEX "\tNumProductList\tIDproductList\tGeneFullName\tGeneType\tGenePosIN";
#                      16              17              18            19        20
print FINALCODEX "\tGenePosOUT\tGeneSize\tGeneProductCount\tGeneCodexCount";
#                      21         22              23              24
print FINALCODEX "\tGeneHgnc\tGeneHgncLink\tGeneEntrez\tGeneEntrezLink\tGeneUcsc\tGeneUcscLink\tGeneMim\tGeneMimLink";
#                      25          26          27            28           29         30         31           32
print FINALCODEX "\tExonBlast\tExonFlag\tWCS1_Flag\tWCS1_Evalue\tWCS1_PctIdentities\tWCS1_UCSC_ID\tWCS1_UCSC_LINK\tWCS2_Flag\tWCS2_Evalue\tWCS2_PctIdentities\tWCS2_UCSC_ID\tWCS2_UCSC_LINK\n";
#                       33       34            35          36                  37           38             39           40          41                    42          43           44
#now matching values
#first ranking the 2 datasets on genename and codexnum
#finalcdsarr
#let us rank dicreasing on flag then increasing on evalue
my $field1f = 2;
my $field2f = 3;
my @rankedfinalcdx = map {$_->[0]} 
sort {
	my @a_fields = @$a[1..$#$a];
	my @b_fields = @$b[1..$#$b];
	$a_fields[$field1f-1] cmp $b_fields[$field1f-1]
		||
	$a_fields[$field2f-1] <=> $b_fields[$field2f-1]}
map {[$_, split /\t/]} @finalfinalcdx;
undef @finalfinalcdx;
my @exfinalfinalcdx;
my $veriffieldex;
#now exploding ranked values
my $kex1 = 0;
while (defined($rankedfinalcdx[$kex1])) {
	my $rkex = "$rankedfinalcdx[$kex1]\t";
	$finalfinalcdx[$kex1+1] = $rankedfinalcdx[$kex1];
	my $ctrkex = 0;
	while ($rkex =~ m/([^\t]+)\t+/igc) {
		$exfinalfinalcdx[$kex1+1][++$ctrkex] = $1;}
	if (not defined($veriffieldex)) {
		$veriffieldex = $ctrkex;}
	else {
		if ($veriffieldex != $ctrkex) {
			die "Big code bug ct fields in ranked final final Codex $veriffieldex != $ctrkex for entry $rankedfinalcdx[$kex1]\n";}}
	++$kex1;}
undef @rankedfinalcdx;
#same with exons blast output
my @rankedblast;
my $rk1 = 1;
while ($rk1 <= $exonblast[0][0]) {
	my $point = $rk1-1;
	$rankedblast[$point] = $exonblast[$rk1][1];
	my $rk2 = 2;
	while(defined($exonblast[$rk1][$rk2])) {
		$rankedblast[$point] = $rankedblast[$point] . "\t$exonblast[$rk1][$rk2]";
		++$rk2;}
	++$rk1;}
my $field1b = 2;
my $field2b = 3;
my @rankedfinalblast = map {$_->[0]} 
sort {
	my @a_fields = @$a[1..$#$a];
	my @b_fields = @$b[1..$#$b];
	$a_fields[$field1b-1] cmp $b_fields[$field1b-1]
		||
	$a_fields[$field2b-1] <=> $b_fields[$field2b-1]}
map {[$_, split /\t/]} @rankedblast;
undef @rankedblast;
#now exploding ranked values
my @exrankedblast;
my $veriffieldex2;
$kex1 = 0;
while (defined($rankedfinalblast[$kex1])) {
	my $rkex = "$rankedfinalblast[$kex1]\t";
	$rankedblast[$kex1+1] = $rankedfinalblast[$kex1];
	my $ctrkex = 0;
	while ($rkex =~ m/([^\t]+)\t+/igc) {
		$exrankedblast[$kex1+1][++$ctrkex] = $1;}
	if (not defined($veriffieldex2)) {
		$veriffieldex2 = $ctrkex;}
	else {
		if ($veriffieldex2 != $ctrkex) {
			die "Big code bug ct fields in ranked final final Blast $veriffieldex2 != $ctrkex for entry $rankedfinalblast[$kex1]\n";}}
	++$kex1;}
undef @rankedfinalblast;
#now lets match both table
my @globalexon;
my $ctfoundfinal=0;
my $ctmissedfinal=0;
my $kfdc = 1;
my $kfdb = 1;
FINDFINAL:
while (1>0) {#until I go out of this looop using another way
	if (not defined($exfinalfinalcdx[$kfdc][1])) {#no more exons
		last FINDFINAL;}
	my $cdxgene = $exfinalfinalcdx[$kfdc][2];
	my $cdxnum =  $exfinalfinalcdx[$kfdc][3];
	if (not defined($exrankedblast[$kfdb][1])) {#no more blast
		++$ctmissedfinal;
		print FINALCODEX "$finalfinalcdx[$kfdc]\t0";
		my $stringexon = "$finalfinalcdx[$kfdc]\t0";
		my $klb = 4;
		while ($klb <= $veriffieldex2) {
			print  FINALCODEX "\tNA";
			$stringexon = $stringexon . "\tNA";
			++$klb;}
		$stringexon = $stringexon . "\t";
		++$globalexon[0][0];
		my $ctglob = 0;
		while ($stringexon =~ m/([^\t]+)\t+/igc) {
			$globalexon[$globalexon[0][0]][++$ctglob] = $1;}
		print FINALCODEX "\n";
		++$kfdc;
		next FINDFINAL;}
	my $blastgene = $exrankedblast[$kfdb][2];
	my $blastnum = $exrankedblast[$kfdb][3];
	if (($cdxgene eq $blastgene) && ($cdxnum == $blastnum)) {
		++$ctfoundfinal;
		print FINALCODEX "$finalfinalcdx[$kfdc]\t1";
		my $stringexon = "$finalfinalcdx[$kfdc]\t1";
		my $klb = 4;
		while ($klb <= $veriffieldex2) {
			print  FINALCODEX "\t$exrankedblast[$kfdb][$klb]";
			$stringexon = $stringexon . "\t$exrankedblast[$kfdb][$klb]";
			++$klb;}
		$stringexon = $stringexon . "\t";
		++$globalexon[0][0];
		my $ctglob = 0;
		while ($stringexon =~ m/([^\t]+)\t+/igc) {
			$globalexon[$globalexon[0][0]][++$ctglob] = $1;}
		print FINALCODEX "\n";
		++$kfdc;
		++$kfdb;
		next FINDFINAL;}
	if (($cdxgene lt $blastgene) || (($cdxgene eq $blastgene) && ($cdxnum < $blastnum))) {#no blast for this codex
		++$ctmissedfinal;
		print FINALCODEX "$finalfinalcdx[$kfdc]\t0";
		my $stringexon = "$finalfinalcdx[$kfdc]\t0";
		my $klb = 4;
		while ($klb <= $veriffieldex2) {
			$stringexon = $stringexon . "\tNA";
			print  FINALCODEX "\tNA";
			++$klb;}
		$stringexon = $stringexon . "\t";
		++$globalexon[0][0];
		my $ctglob = 0;
		while ($stringexon =~ m/([^\t]+)\t+/igc) {
			$globalexon[$globalexon[0][0]][++$ctglob] = $1;}
		print FINALCODEX "\n";
		++$kfdc;
		next FINDFINAL;}
	else {
		++$kfdb;
		next FINDFINAL;}}
close FINALCODEX;
print LOG "FINAL CODEX STATUS\tPassed\nBLAST MATCHED CODEX\t$ctfoundfinal\nBLAST MISSED CODEX\t$ctmissedfinal\n";
#now generating the gene table
my $kgn = 1;
my @globalgene;
while ($kgn <= $globalexon[0][0]) {
	if ($kgn == 1) {
		++$globalgene[0][0];}
	elsif ($globalexon[$kgn][2] ne $globalexon[$kgn-1][2]) {
		++$globalgene[0][0];}
	else {}
	++$globalgene[$globalgene[0][0]][0];#counting exons (should match at the end)
	if (not defined($globalgene[$globalgene[0][0]][1])) {#we record all interesting info
		$globalgene[$globalgene[0][0]][1] = $globalexon[$kgn][2];#genename
		$globalgene[$globalgene[0][0]][2] = $globalexon[$kgn][4];#refgeno
		$globalgene[$globalgene[0][0]][3] = $globalexon[$kgn][6];#chr
		$globalgene[$globalgene[0][0]][4] = $globalexon[$kgn][7];#Strand
		$globalgene[$globalgene[0][0]][5] = $globalexon[$kgn][20];#PosIN
		$globalgene[$globalgene[0][0]][6] = $globalexon[$kgn][21];#PosOUT
		$globalgene[$globalgene[0][0]][7] = $globalexon[$kgn][22];#Size
		$globalgene[$globalgene[0][0]][8] = $globalexon[$kgn][19];#Type
		$globalgene[$globalgene[0][0]][9] = $globalexon[$kgn][23];#Product Count
		$globalgene[$globalgene[0][0]][10] = $globalexon[$kgn][24];#Exons Count
		$globalgene[$globalgene[0][0]][11] = $globalexon[$kgn][25];#hgnc
		$globalgene[$globalgene[0][0]][12] = $globalexon[$kgn][26];#hgnc link
		$globalgene[$globalgene[0][0]][13] = $globalexon[$kgn][27];#entrez
		$globalgene[$globalgene[0][0]][14] = $globalexon[$kgn][28];#entrez link
		$globalgene[$globalgene[0][0]][15] = $globalexon[$kgn][29];#ucsc
		$globalgene[$globalgene[0][0]][16] = $globalexon[$kgn][30];#ucsc link
		$globalgene[$globalgene[0][0]][17] = $globalexon[$kgn][31];#OMIM
		$globalgene[$globalgene[0][0]][18] = $globalexon[$kgn][32];#OMIM link
		$globalgene[$globalgene[0][0]][32] = $globalexon[$kgn][18];#gene full name
		$globalgene[$globalgene[0][0]][19] = "";#init list of products (unique)
		$globalgene[$globalgene[0][0]][20] = 0;#init count blasted
		$globalgene[$globalgene[0][0]][29] = 0;#init count not blasted
		$globalgene[$globalgene[0][0]][21] = 1;#init global flag (looks perfet)
		$globalgene[$globalgene[0][0]][22] = 0;#init ct perfect
		$globalgene[$globalgene[0][0]][23] = 0;#init ct low
		$globalgene[$globalgene[0][0]][24] = 0;#init ct mid
		$globalgene[$globalgene[0][0]][25] = 0;#init ct high
		$globalgene[$globalgene[0][0]][28] = 0;#init ct NA
		$globalgene[$globalgene[0][0]][30] = 0;#init count of WCS products
		$globalgene[$globalgene[0][0]][31] = ""; #init codex WCS list
		$globalgene[$globalgene[0][0]][27] = 0;#count of products (should match at end)
		$globalgene[$globalgene[0][0]][26] = "";}#init list of WCS products
	my $prodlist = $globalexon[$kgn][17];
	while ($prodlist =~ m/([^\;]+)\;+/igc) {
		my $product = $1;
		if ($globalgene[$globalgene[0][0]][19] !~ m/$product\;/) {
			++$globalgene[$globalgene[0][0]][27];
			$globalgene[$globalgene[0][0]][19] = $globalgene[$globalgene[0][0]][19] . "$product;";}}
	if ($globalexon[$kgn][33] == 1) {#blast data available
		++$globalgene[$globalgene[0][0]][20];
		if ($globalexon[$kgn][34] == 1) {
			++$globalgene[$globalgene[0][0]][22];}
		elsif ($globalexon[$kgn][34] == 2) {
			 ++$globalgene[$globalgene[0][0]][23];}
		elsif ($globalexon[$kgn][34] == 3) {
			 ++$globalgene[$globalgene[0][0]][24];}
		elsif ($globalexon[$kgn][34] == 4) {
			 ++$globalgene[$globalgene[0][0]][25];}
		else  {
			 ++$globalgene[$globalgene[0][0]][28];}
		if ($globalexon[$kgn][34] > $globalgene[$globalgene[0][0]][21]) {
			$globalgene[$globalgene[0][0]][30] = 0;#this is a reset of values
			$globalgene[$globalgene[0][0]][31] = ""; #init codex WCS list
			$globalgene[$globalgene[0][0]][26] = "";}
		if ($globalexon[$kgn][34] >= $globalgene[$globalgene[0][0]][21]) {
			$globalgene[$globalgene[0][0]][31] = $globalgene[$globalgene[0][0]][31] . "$globalexon[$kgn][3];";
			$globalgene[$globalgene[0][0]][21] = $globalexon[$kgn][34];
			my $prodlist = $globalexon[$kgn][17];
			while ($prodlist =~ m/([^\;]+)\;+/igc) {
				my $product = $1;
				if ($globalgene[$globalgene[0][0]][26] !~ m/$product\;/) {
					++$globalgene[$globalgene[0][0]][30];
					$globalgene[$globalgene[0][0]][26] = $globalgene[$globalgene[0][0]][26] . "$product;";}}}}
	else {#no blast available
		++$globalgene[$globalgene[0][0]][29];}
	++$kgn;}
#ready to generate the gene table
print "GENES FROM EXONS TABLE\t$globalgene[0][0]\n";
print LOG "GENES FROM EXONS TABLE\t$globalgene[0][0]\n";
print GENESTABLE "GeneName\tRefGeno\tGeneChromosome\tGeneStrand\tGenePosIN\tGenePosOUT\tGeneSize";	
print GENESTABLE "\tGeneType\tGeneFullName\tGeneHgncID\tGeneHgncLink\tGeneEntrezID\tGeneEntrezLink\tGeneUcscID\tGeneUcscLink\tGeneOmimID\tGeneOmimLink";
print GENESTABLE "\tGeneCodexCount\tGeneProductCount\tGeneProductList\tCodexWithBlast\tCodexWithoutBlast";
print GENESTABLE "\tGlobalBlastFlag\tCountCodexNoHomology\tCountCodexLowHomology\tCountCodexMidHomology\tCountCodexHighHomology\tCountCodexOtherHomology";
print GENESTABLE "\tWcsProductList\tWcsProductCount\tWcsCodexList\n";
my $kgg = 1;
while ($kgg <= 	$globalgene[0][0]) {
	print GENESTABLE "$globalgene[$kgg][1]\t$globalgene[$kgg][2]\t$globalgene[$kgg][3]\t$globalgene[$kgg][4]\t$globalgene[$kgg][5]\t$globalgene[$kgg][6]\t$globalgene[$kgg][7]\t";
	print GENESTABLE "$globalgene[$kgg][8]\t$globalgene[$kgg][32]\t$globalgene[$kgg][11]\t$globalgene[$kgg][12]\t$globalgene[$kgg][13]\t$globalgene[$kgg][14]\t$globalgene[$kgg][15]\t";	
	print GENESTABLE "$globalgene[$kgg][16]\t$globalgene[$kgg][17]\t$globalgene[$kgg][18]\t$globalgene[$kgg][0]\t$globalgene[$kgg][27]\t$globalgene[$kgg][19]\t$globalgene[$kgg][20]\t";	
	print GENESTABLE "$globalgene[$kgg][29]\t$globalgene[$kgg][21]\t$globalgene[$kgg][22]\t$globalgene[$kgg][23]\t$globalgene[$kgg][24]\t$globalgene[$kgg][25]\t$globalgene[$kgg][28]\t";		
	print GENESTABLE "$globalgene[$kgg][26]\t$globalgene[$kgg][30]\t$globalgene[$kgg][31]\n";
	++$kgg;}
close GENESTABLE;
#now generating the product table
		




################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#last copy of everything
$showtime = getdatetime();
print LOG "PROCESS_STATUS\tCompleted by $showtime";
close LOG;
copymaintobackup($mainoutputdir,$backupoutputdir);
die "Congratulations. Process went to completion without any problem.\nProcess started by $rundatetime and ended by $showtime\n";
	#http://omim.org/entry/613726	
	#https://www.ncbi.nlm.nih.gov/gene/55129
	#http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=25519
#####################	
sub loadtmvalues {
	my %xtmhash = ( "A" => 2, "C" => 4, "G" => 4, "T" => 2,);
	return (%xtmhash);}
#####################	
sub loadtranslation {
	my %xcodonhash = ('AAA'=>'K','K'=>'AAA','AAC'=>'N','N'=>'AAC','AAG'=>'K','K'=>'AAG','AAT'=>'N','N'=>'AAT','ACA'=>'T','T'=>'ACA','ACC'=>'T','T'=>'ACC','ACG'=>'T','T'=>'ACG','ACT'=>'T','T'=>'ACT',
	'AGA'=>'R','R'=>'AGA','AGC'=>'S','S'=>'AGC','AGG'=>'R','R'=>'AGG','AGT'=>'S','S'=>'AGT','ATA'=>'I','I'=>'ATA','ATC'=>'I','I'=>'ATC','ATG'=>'M','M'=>'ATG','ATT'=>'I','I'=>'ATT',
	'CAA'=>'Q','Q'=>'CAA','CAC'=>'H','H'=>'CAC','CAG'=>'Q','Q'=>'CAG','CAT'=>'H','H'=>'CAT','CCA'=>'P','P'=>'CCA','CCC'=>'P','P'=>'CCC','CCG'=>'P','P'=>'CCG','CCT'=>'P','P'=>'CCT',
	'CGA'=>'R','R'=>'CGA','CGC'=>'R','R'=>'CGC','CGG'=>'R','R'=>'CGG','CGT'=>'R','R'=>'CGT','CTA'=>'L','L'=>'CTA','CTC'=>'L','L'=>'CTC','CTG'=>'L','L'=>'CTG','CTT'=>'L','L'=>'CTT',
	'GAA'=>'E','E'=>'GAA','GAC'=>'D','D'=>'GAC','GAG'=>'E','E'=>'GAG','GAT'=>'D','D'=>'GAT','GCA'=>'A','A'=>'GCA','GCC'=>'A','A'=>'GCC','GCG'=>'A','A'=>'GCG','GCT'=>'A','A'=>'GCT',
	'GGA'=>'G','G'=>'GGA','GGC'=>'G','G'=>'GGC','GGG'=>'G','G'=>'GGG','GGT'=>'G','G'=>'GGT','GTA'=>'V','V'=>'GTA','GTC'=>'V','V'=>'GTC','GTG'=>'V','V'=>'GTG','GTT'=>'V','V'=>'GTT',
	'TAA'=>'B','B'=>'TAA','TAG'=>'B','B'=>'TAG','TAC'=>'Y','Y'=>'TAC','TAT'=>'Y','Y'=>'TAT','TCA'=>'S','S'=>'TCA','TCC'=>'S','S'=>'TCC','TCG'=>'S','S'=>'TCG','TCT'=>'S','S'=>'TCT',
	'TGA'=>'B','B'=>'TGA','TGC'=>'C','C'=>'TGC','TGG'=>'W','W'=>'TGG','TGT'=>'C','C'=>'TGT','TTA'=>'L','L'=>'TTA','TTC'=>'F','F'=>'TTC','TTG'=>'L','L'=>'TTG','TTT'=>'F','F'=>'TTT',
	'OOO'=>'O','O'=>'OOO',);
	return (%xcodonhash);}	
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#subroutine that copies content of one folder to another one
sub copymaintobackup {
	my $main = shift;
	my $backup = shift;
	opendir (DIR,"$main") ||
		die "Cannot read from main output folder $main\n";
	my @catchmain = grep {$_ ne "." & $_ ne ".."} readdir DIR;
	closedir DIR;
	opendir (DIR,"$backup") ||
		die "Cannot read from backup output folder $backup\n";
	my @catchbackup = grep {$_ ne "." & $_ ne ".."} readdir DIR;
	closedir DIR;
	my $kmain = 0;
	while (defined($catchmain[$kmain])) {
		my $mainfile = "$main\/$catchmain[$kmain]";
		my $found = 0;
		my $kbackup = 0;
		FOUNDCOPY:
		while (defined($catchbackup[$kbackup])) {
			if ($catchmain[$kmain] eq $catchbackup[$kbackup]) {
				my $backupfile = "$backup\/$catchbackup[$kbackup]";
				if (compare("$mainfile","$backupfile") == 0) {
					$found = 1;
					last FOUNDCOPY;}
				else {
					unlink $backupfile ||
						die "Cannot erase backup file $backupfile";
					copy("$mainfile","$backupfile") or die "Copy failed: $!";
					$found = 1;
					last FOUNDCOPY;}}
			++$kbackup;}
		if ($found == 0) {
			my $backupfile = "$backup\/$catchmain[$kmain]";
			copy("$mainfile","$backupfile") or die "Copy failed: $!";}
		++$kmain;}
	print "\nBackup performed for $kmain files\n";}
#################
sub getdatetime {
#for fetching current date and time
	my $today = DateTime->now();
	return $today->ymd('') . '_' . $today->hms('');}
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
################################################################################

=head1 NAME

FetchGeneData_v1i1.pl (version 1, issue 1)

=head1 HELP

This code has been developped on a Linux (Ubuntu 14.04) 64bit platform.

This codes requires Perl to be installed (and some Perl Modules).

This code is part of the CHU ULg Genetic Diseases "pipeline".
It uses a tab delimited txt file (with header; 2 fields = 
GeneName & GeneInfo) to capture for each gene from 3 related txt files (= 
FASTA sequence file, Genbank Full Info file and NCBI Gene Result file) 
any relevant information and to create a primary database (tables).
Remark: All of these subsets of 3 files should be located in the same folder.

Important: This code uses blastn. So check that this is available from the platform
	
For additionnal info. contact the code Owner/Author (or from the WEB: www.ngyx.eu)

=head1 VERSION

This documentation refers to FetchGeneData_v1i1.pl (version 1, issue 1)

=head1  USAGE

	perl FetchGeneData_v1i1.pl [options]  
	
Example (no optional argument):
	
perl FetchGeneData_v1i1.pl \ 
--outputfolder '/home/user/TestOutputBlastx1' \
--backupfolder '/media/user/Seagate Backup Plus Drive/TestBackupBlastx1' \
--geneslist '/home/user/CHU_ULg_Genetic_Diseases/CUGD_v1i1_UserGeneListTest.txt' \
--inputfolder '/home/user/CHU_ULg_Genetic_Diseases/Infiles/CUGD20161031' \
--blastdb '/home/user/CHU_ULg_Genetic_Diseases/BlastDBs/HG38BlastDB/hg38'

=head1 REQUIRED ARGUMENTS

=over

=item --outputfolder [=] <outputdir>

Path of folder to record outputs

=for Euclid:
	outputdir.type: string

=item --backupfolder [=] <backupdir>

Path of folder to backup outputs

=for Euclid:
	backupdir.type: string
	
=item --geneslist [=] <genesfile>

Path+Name of file (Tab delimited Txt with header; 2 fields: GeneName & GeneInfo)

=for Euclid:
	genesfile.type: string

=item --inputfolder [=] <inputdir>

Path to input folder where NCBI, Genbank and FASTA files should be located

=for Euclid:
	inputdir.type: string

=item --blastdb [=] <blastdb>

Path + generic name (formatdb) of the genomic db used to fetch gene homologies

=for Euclid:
	blastdb.type: string
	

=back

=head1 OPTIONS

=over

=item --blastmaxhits [=] <blastmaxhits>

Integer positive <= 1000 (default = 3) used to restrict the number of hits returned

=for Euclid:
	blastmaxhits.type: int > 0
	blastmaxhits.default: 3
	
=item --blastextend [=] <blastextend>

Integer positive >= 12 (default = 36) used to extend codex both sides by adding nuc for Blast

=for Euclid:
	blastextend.type: int >= 12
	blastextend.default: 30

=item --blastredlimit [=] <blastred>

Integer positive >= 48 (default = 66) used as T° anneal threshold to flag for sure gene with pseudogene or...

=for Euclid:
	blastred.type: int >= 48
	blastred.default: 66

=item --blastorangelimit [=] <blastorange>

Integer positive >= 48 (default = 60)  used as T° anneal threshold to flag probable gene with pseudogene or...

=for Euclid:
	blastorange.type: int >= 48
	blastorange.default: 60
	
=item --blastyellowlimit [=] <blastyellow>

Integer positive >= 48 (default = ) used as T° anneal threshold to flag putative gene with pseudogene or...

=for Euclid:
	blastyellow.type: int >= 48
	blastyellow.default: 54
	
=item --blastevalue [=] <blastevalue>

Real positive for evalue <= 1.0E-06 (default) used to select blast hits to review

=for Euclid:
	blastevalue.type: +number <= 1.0E-06
 	blastevalue.default: 1.0E-06

=item --blastpctidentity [=] <blastpctidentity>

Real positive for pct identity <= 100 (default 90) used to select blast hits to review

=for Euclid:
	blastpctidentity.type: +number <= 100 
 	blastpctidentity.default: 90

=item --blastminwinsize [=] <blastminwinsize>

Integer positive for min of windowed approach calculating Tm (default 18) used to derive max Tm

=for Euclid:
	blastminwinsize.type: +integer <= 24 
 	blastminwinsize.default: 18
 	
=item --blastmaxwinsize [=] <blastmaxwinsize>

Integer positive for max of windowed approach calculating Tm (default 24) used to derive max Tm

=for Euclid:
	blastmaxwinsize.type: +number <= 30 
 	blastmaxwinsize.default: 24
	
=back
	
=item --version

=item --usage

=item --help

=item --man

=back

=head1 AUTHOR

LECOCQ Pierre 59072327136 (Pierre.Lecocq@NGYX.eu)

=head1 BUGS

In case a bug would be identified please contact code Owner/Author.

=head1 COPYRIGHT

Copyright (c) 2016JUL23, code Owner: LECOCQ Pierre 59072327136.
All Rights Reserved.

This program is not free software. It may not be used , redistributed
and/or modified without code Owner/Author written approval.

677968693279877869823280736982826932766967796781326669765357485550515055495154
