#!/usr/bin/env perl
################################################################################
#677968693279877869823280736982826932696779678132535748555051505432666947666976#
#some modules required (if not installed: cpan install <module>)
use Modern::Perl '2011'; # used as warnings/strict
use Smart::Comments '###'; # used to trace run
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
################################################################################
#677968693279877869823280736982826932696779678132535748555051505432666947666976#
#initialization
my $code = $0;
my @randomarr = loadrandom();
my %monthletterhash = monthletter();
my $maindir = getcwd();#to catch current location
my $rundatetime = getdatetimeold();#for reporting
print "\n$rundatetime: $code Process Started from folder $maindir.\n";
my %codonhash = loadtranslation();
my $showtime;
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
my $genelistfile = "/home/user/Exomes/NeuroGenesInfo.txt";
my @genelistarr;
if (-e $genelistfile) {
	open (LIST,"<$genelistfile") ||
		die "Cannot read from gene list file $genelistfile\n";
	my $header = uc <LIST>;
	if ($header ne "GENENAME\tGENEINFO\n") {
		die "Bad header ($header) in gene list file $genelistfile\n";}
	my @temparr;
	my $pointer = 0;
	while (defined(my $linein = <LIST>)) {
		if ($linein =~ m/^([^\t]+)\t([^\t]*)\n$/) {
			$temparr[$pointer] = uc $1 . "\t$2";
			++$pointer;}
		else {
			die "Bad format for entry in gene list: $linein\n";}}
	close LIST;
	$genelistarr[0][0][0] = $pointer;
	#now ranking
	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]}
	map {[$_, split /\t/]} @temparr;
	undef @temparr;
	my $k1 = 0;
	while (defined($rankarr[$k1])) {
		if ($rankarr[$k1] =~  m/^([^\t]+)\t([^\t]*)$/) {
			$genelistarr[$k1+1][0][1] = $1;
			$genelistarr[$k1+1][0][2] = $2;}
		else {
			die "Problem while review ranked gene list with $rankarr[$k1]\n";}
		++$k1;}
	undef @rankarr;
	print "Identified $k1 entries from gene list\n";}
else {
	die "A gene list file is required\n";}
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
my $genedir = "/home/user/Exomes/NeuroGenes";
if ((!-e $genedir) || (!-d $genedir)) {
	die "Cannot find gene folder or is not a folder $genedir";}
opendir (GENE,"$genedir")  ||
	die "Cannot read from gene folder $genedir";
my @catcharr = grep {$_ ne "." & $_ ne ".."} readdir GENE;
closedir GENE;
my $k1 = 0;
while (defined($catcharr[$k1])) {
	my $genefile = $genedir . "\/$catcharr[$k1]";
	if ((-e $genefile) && (-f $genefile)) {
		if ($catcharr[$k1] =~ m/^(.+)\.FASTA$/i) {
			my $genename = uc $1;
			my $k2 = 1;
			FOUNDFASTA:
			while (defined($genelistarr[$k2][0][1])) {
				if ($genename eq $genelistarr[$k2][0][1]) {
					$genelistarr[$k2][0][3] = $genefile;
					last FOUNDFASTA;}
				++$k2;}}
		elsif ($catcharr[$k1] =~ m/^(.+)\.GB$/i) {
			my $genename = uc $1;
			my $k2 = 1;
			FOUNDGB:
			while (defined($genelistarr[$k2][0][1])) {
				if ($genename eq $genelistarr[$k2][0][1]) {
					$genelistarr[$k2][0][4] = $genefile;
					last FOUNDGB;}
				++$k2;}}
		elsif ($catcharr[$k1] =~ m/^(.+)\.TXT$/i) {
			my $genename = uc $1;
			my $k2 = 1;
			FOUNDTXT:
			while (defined($genelistarr[$k2][0][1])) {
				if ($genename eq $genelistarr[$k2][0][1]) {
					$genelistarr[$k2][0][5] = $genefile;
					last FOUNDTXT;}
				++$k2;}}
		else {}}
	++$k1;}
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#Output Folder
my $outputdir = $maindir . "\/outputdir";
if ((-e $outputdir) && (-d $outputdir)) {}
elsif (-e $outputdir) {
	die "Output folder exists ($outputdir) but is not a folder\n";}
else {
	mkdir $outputdir ||
		die "Cannot create Output folder ($outputdir)\n";}
opendir (OUTDIR,"$outputdir")  ||
	die "Cannot read from Output folder $outputdir";
my @catcharr2 = grep {$_ ne "." & $_ ne ".."} readdir OUTDIR;
closedir OUTDIR;
$k1 = 0;
while (defined($catcharr2[$k1])) {
	my $outfile = $outputdir . "\/$catcharr2[$k1]";
	if (-f $outfile) {
		unlink $outfile ||
			die "Cannot erase output file $outfile";}
	++$k1;}
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#Output files	
#CODEX BED file OK
my $codexbed = $outputdir ."\/Codex.bed";
my $tempfile = $codexbed;
if (-e $tempfile) {
	unlink $tempfile ||
		die "Unable to erase exsiting temporary file $tempfile";}
open(CODEXBED,">>$tempfile") ||
	die "Unable to create temporary file $tempfile";
#Genes BED file OK
my $genesbed = $outputdir ."\/Genes.bed";
$tempfile = $genesbed;
if (-e $tempfile) {
	unlink $tempfile ||
		die "Unable to erase exsiting temporary file $tempfile";}
open(GENESBED,">>$tempfile") ||
	die "Unable to create temporary file $tempfile";
#CODEX Table OK
my $codextable = $outputdir ."\/Codex.txt";
$tempfile = $codextable;
if (-e $tempfile) {
	unlink $tempfile ||
		die "Unable to erase exsiting temporary file $tempfile";}
open(CODEXTABLE,">>$tempfile") ||
	die "Unable to create temporary file $tempfile";
print CODEXTABLE "Codex_ID\tHGNC_GeneSymbol\tHGNC_GeneFullname\tNCBI_GeneID\tHGNC_GeneID\tGene_User_Info\tNCBI_Summary\tCodex_Num";
print CODEXTABLE "\tAlignementSourceReference\tASR_Chromosome\tASR_StartPosition\tASR_EndPosition\tASR_Strand\tCodingStrandStartPosition\tCodingStrandEndPosition";
print CODEXTABLE "\tLinked_CDS_Count\tLinked_CDS_NumList\tGene_CDS_Count\tCodexSequence\n";
#Genes Table OK
my $genestable = $outputdir ."\/Genes.txt";
$tempfile = $genestable;
if (-e $tempfile) {
	unlink $tempfile ||
		die "Unable to erase exsiting temporary file $tempfile";}
open(GENESTABLE,">>$tempfile") ||
	die "Unable to create temporary file $tempfile";
print GENESTABLE "HGNC_GeneSymbol\tHGNC_GeneFullname\tHGNC_GeneID\tNCBI_GeneID\tGene_User_Info\tNCBI_Summary\tGene_Type\tGene_RelatedLinks\tGeneOtherNames";
print GENESTABLE "\tAlignementSourceReference\tASR_GeneLocus\tASR_Chromosome\tASR_StartPosition\tASR_EndPosition\tASR_Strand\tASR_CDS_Count\tASR_Codex_Count\tASR_GeneSequence\n";
#CDS Table OK
my $cdstable = $outputdir ."\/Cds.txt";
$tempfile = $cdstable;
if (-e $tempfile) {
	unlink $tempfile ||
		die "Unable to erase exsiting temporary file $tempfile";}
open(CDSTABLE,">>$tempfile") ||
	die "Unable to create temporary file $tempfile";
print CDSTABLE "CDS_ID\tGENE_ID\tCDS_NBR\tGENE_TYPE\tCODEX_LIST\tCDS_SEQ\tCDS_PRT\n";
my $codexfasta = $outputdir ."\/Codex.fasta";
$tempfile = $codexfasta;
if (-e $tempfile) {
	unlink $tempfile ||
		die "Unable to erase exsiting temporary file $tempfile";}
open(CODEXFASTA,">>$tempfile") ||
	die "Unable to create temporary file $tempfile";
#Genes FASTA OK
my $genesfasta = $outputdir ."\/Genes.fasta";
$tempfile = $genesfasta;
if (-e $tempfile) {
	unlink $tempfile ||
		die "Unable to erase exsiting temporary file $tempfile";}
open(GENESFASTA,">>$tempfile") ||
	die "Unable to create temporary file $tempfile";
#CDS FASTA OK
my $cdsfasta = $outputdir ."\/Cds.fasta";
$tempfile = $cdsfasta;
if (-e $tempfile) {
	unlink $tempfile ||
		die "Unable to erase exsiting temporary file $tempfile";}
open(CDSFASTA,">>$tempfile") ||
	die "Unable to create temporary file $tempfile";
#Also known as table
my $genesnames = $outputdir ."\/GenesNames.txt";
$tempfile = $genesnames;
if (-e $tempfile) {
	unlink $tempfile ||
		die "Unable to erase exsiting temporary file $tempfile";}
open(GENESNAMES,">>$tempfile") ||
	die "Unable to create temporary file $tempfile";

################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890	
my $listofmissing = "";
my $countofmissing  = 0;
$k1 = 1;
while (defined($genelistarr[$k1][0][1])) {
	if ((not defined($genelistarr[$k1][0][3])) || (not defined($genelistarr[$k1][0][4])) || (not defined($genelistarr[$k1][0][5]))) {
		$listofmissing =$listofmissing . "$genelistarr[$k1][0][1]\n";
		++$countofmissing;}
	++$k1;}
if ($countofmissing > 0) {
	die "Missing files for $countofmissing genes (list here below):\n$listofmissing";}
$showtime = getdatetimeold();
print "$showtime: Identified $genelistarr[0][0][0] genes and FASTA, GB and NCBI putative files\n";
#loading NCBI data
$k1 = 1;
while (defined($genelistarr[$k1][0][1])) {
	open (NCBI,"<$genelistarr[$k1][0][5]") ||
		die "Cannot read from NCBI file of $genelistarr[$k1][0][1] ($genelistarr[$k1][0][5])\n";
	FETCHNCBI:
	while(defined(my $linein = <NCBI>)) {
		if ($linein =~ m/^Official Symbol\:\s+([^\(]+)\s+\(/i) {
			$genelistarr[$k1][1][1] = uc $1;
			if ($genelistarr[$k1][1][1] ne $genelistarr[$k1][0][1]) {
				die "Problem NCBI Official name ($genelistarr[$k1][1][1] ne $genelistarr[$k1][0][1])\n";}}
		if ($linein =~ m/^Gene ID\:\s+([^\,]+)\,/i) {
			$genelistarr[$k1][1][9] = $1;}
		if ($linein =~ m/^Official full name\:\s+(.+)\s+\(\s*PROVIDED/i) {
			$genelistarr[$k1][1][2] = $1;}
		if ($linein =~ m/^Primary source\:\s+(.+)\n$/i) {
			$genelistarr[$k1][1][3] = $1;}
		if ($linein =~ m/^See related\:\s+(.+)\n$/i) {
			$genelistarr[$k1][1][4] = $1;}
		if ($linein =~ m/^Gene type\:\s+(.+)\n$/i) {
			$genelistarr[$k1][1][5] = $1;
			if ($genelistarr[$k1][1][5] =~ m/protein coding/i) {
				$genelistarr[$k1][1][5] = "PRT";}
			else {
				$genelistarr[$k1][1][5] = "RNA";}}
		if ($linein =~ m/^Also known as\:\s+(.+)\n$/i) {
			$genelistarr[$k1][1][6] = $1;
			FETCHRELATED:
			while(defined($linein = <NCBI>)) {
				if ($linein =~ m/^\s*\n$/) {
					$genelistarr[$k1][1][6] = $genelistarr[$k1][1][6] . ';';
					last FETCHRELATED;}
				else {
					chomp $linein;
					$genelistarr[$k1][1][6] = $genelistarr[$k1][1][6] . " $linein;";}}}
		if ($linein =~ m/^Summary\:\s+(.+)\n$/i) {
			$genelistarr[$k1][1][7] = $1;
			FETCHSUMMARY:
			while(defined($linein = <NCBI>)) {
				if ($linein =~ m/^\s*\n$/) {
					last FETCHSUMMARY;}
				else {
					chomp $linein;
					$genelistarr[$k1][1][7] = $genelistarr[$k1][1][7] . " $linein";}}}
		if ($linein =~ m/^Gene type\:\s+(.+)\n$/i) {
			$genelistarr[$k1][1][8] = $1;}
		if ($linein =~ m/^GENOMIC CONTEXT\n$/i) {
			last FETCHNCBI;}}
	close NCBI;
	if (not defined($genelistarr[$k1][1][4])) {
		$genelistarr[$k1][1][4] = "No Additionnal References Available";}
	else {
		while ($genelistarr[$k1][1][4] =~ s/[\t\n\r\f\s]//) {}}
	if (not defined($genelistarr[$k1][1][7])) {
		$genelistarr[$k1][1][7] = "No NCBI Summary Available";}
	else {
		while ($genelistarr[$k1][1][7] =~ s/[\t\n\r\f]/ /) {}
		while ($genelistarr[$k1][1][7] =~ s/\s\s/ /) {}
		if ($genelistarr[$k1][1][7] =~ m/^\s+$/) {
			$genelistarr[$k1][1][7] = "No NCBI Summary Available";}}
	if (defined($genelistarr[$k1][1][3])) {
		while ($genelistarr[$k1][1][3] =~ s/[\t\n\r\f]//) {}}
	++$k1;}
$showtime = getdatetimeold();
print "$showtime: NCBI data loaded\n";
#loading FASTA data
#>gi|568815591:69598475-70793068 Homo sapiens chromosome 7, GRCh38.p7 Primary Assembly
#AGTGCCCGGGCTCCCGCTGCAGCCGCCCGGGGACACGCCGTGCACCCTCCGGCTCGGGGCTTTCTCGGCG
$k1 = 1;
while (defined($genelistarr[$k1][0][1])) {
	open (FASTA,"<$genelistarr[$k1][0][3]") ||
		die "Cannot read from FASTA file of $genelistarr[$k1][0][1] ($genelistarr[$k1][0][3])\n";
	my $firstline = <FASTA>;
	#>gi|568815591:69598475-70793068 Homo sapiens chromosome 7, GRCh38.p7 Primary Assembly
	if ($firstline =~ m/^\>[^\:]+\:(C?)(\d+)\-(\d+).+chromosome\s+([^\,]+)\,\s+([^\s]+)\s+/i) {
		my $ori = $1;
		$genelistarr[$k1][3][1] = "chr$4";
		my $num1 = $2;
		my $num2 = $3;
		$genelistarr[$k1][3][5] = $5;
		if ($num1 > $num2) {
			$genelistarr[$k1][3][2] = $num2;
			$genelistarr[$k1][3][3] = $num1;}
		else {
			$genelistarr[$k1][3][2] = $num1;
			$genelistarr[$k1][3][3] = $num2;}
		if ((defined($ori)) && (uc $ori)) {
			$genelistarr[$k1][3][4] = '-';}
		else {
			$genelistarr[$k1][3][4] = '+';}}
	else {
		die "Problem with FASTA file of $genelistarr[$k1][0][1] first lien $firstline";}
	$genelistarr[$k1][3][6] = "";
	while (my $linein = <FASTA>) {
		chomp $linein;
		$genelistarr[$k1][3][6] = $genelistarr[$k1][3][6] . uc ($linein);}
	close FASTA;
	if ($genelistarr[$k1][3][6] !~ m/^[ACGTVBDHYRSWKMN]+$/) {
		die "Problem with FASTA sequence for $genelistarr[$k1][0][1]";}
	$genelistarr[$k1][3][7] = length($genelistarr[$k1][3][6]);
	my $reflen = ($genelistarr[$k1][3][3] - $genelistarr[$k1][3][2] + 1);
	if ($reflen != $genelistarr[$k1][3][7]) {
		die "Problem with gene ref. $reflen != $genelistarr[$k1][3][7] for $genelistarr[$k1][0][1] $k1 ($genelistarr[$k1][3][3] - $genelistarr[$k1][3][2] $firstline)";}
	++$k1;}
$showtime = getdatetimeold();
print "$showtime: FASTA data loaded\n";
print GENESNAMES "UsedGeneName\tHGNC_Flag\tHGNC_GeneSymbol\tHGNC_GeneID\tHGNC_GeneFullname\tNCBI_GeneID\tNCBI_Summary\n";
my @usednamearr;
my $pointername = 0;
$k1 = 1;
while (defined($genelistarr[$k1][0][1])) {
	my $stringname = "$genelistarr[$k1][0][1]\t$genelistarr[$k1][1][3]\t$genelistarr[$k1][1][2]\t$genelistarr[$k1][1][9]\t$genelistarr[$k1][1][7]";
	$usednamearr[$pointername] = "$genelistarr[$k1][0][1]\t1\t$stringname";
	++$pointername;
	if (defined($genelistarr[$k1][1][6])) {
		my $stringofused = uc $genelistarr[$k1][1][6];
		while ($stringofused =~ s/[\n\f\r\s\t]//) {}
		while ($stringofused =~ s/\;\;/\;/) {}
		$genelistarr[$k1][1][6] = $stringofused;
		while ($stringofused =~ m/([^\;]+)\;+/igc) {
			$usednamearr[$pointername] = "$1\t0\t$stringname";
			++$pointername;}}
	else {
		$genelistarr[$k1][1][6] = 'None';}
	++$k1;}
#now ranking
my $field1x = 1;
my $field2x = 2;
my @newarr;
@newarr = map {$_->[0]} 
sort {
	my @a_fields = @$a[1..$#$a];
	my @b_fields = @$b[1..$#$b];
	$a_fields[$field1x-1] cmp $b_fields[$field1x-1]
	||
	$b_fields[$field2x-1] <=> $a_fields[$field2x-1]}
map {[$_, split /\t/]} @usednamearr;
undef @usednamearr;
$k1 = 0;
while (defined($newarr[$k1])) {
	print GENESNAMES "$newarr[$k1]\n";
	++$k1;}
undef @newarr;
close GENESNAMES;
#loading gb data
$k1 = 1;
while (defined($genelistarr[$k1][0][1])) {
	open (GB,"<$genelistarr[$k1][0][4]") ||
		die "Cannot read from GB file of $genelistarr[$k1][0][1] ($genelistarr[$k1][0][4])\n";
	my $type;
	if ($genelistarr[$k1][1][5] eq 'PRT') {
		$type = 'CDS';}
	else {
		$type = 'ncRNA';}
	my @cdsarr;
	undef @cdsarr;
	FETCHGB:
	while (defined(my $linein = <GB>)) {
		#DEFINITION  Homo sapiens chromosome 12, GRCh38.p7 Primary Assembly.
		#ACCESSION   NC_000012 REGION: 8822472..8887202 GPC_000001304
		#ACCESSION   NC_000012 REGION: complement(21797389..21941869) GPC_000001304
		#VERSION     NC_000012.12  GI:568815586
		#LOCUS       NC_000012              64731 bp
		if ($linein =~ m/^LOCUS\s+[^\s]+\s+(\d+)\s+/i) {
			$genelistarr[$k1][2][10] = $1;}
		if ($linein =~ m/^DEFINITION\s+(.+)\s+chromosome\s+([0-9A-X]+)\,\s+([^\s]+)\s*/i) {
			$genelistarr[$k1][2][1] = $1;
			$genelistarr[$k1][2][2] = "chr$2";
			$genelistarr[$k1][2][3] = $3;}
		if ($linein =~ m/^ACCESSION\s+.+\s+REGION\:\s+([^\s]+)\s*/i) {
			my $location = $1;
			if ($location =~ m/complement/i) {
				$genelistarr[$k1][2][4] = "-";}
			else {
				$genelistarr[$k1][2][4] = "+";}
			if ($location =~ m/(\d+)\.+(\d+)/i) {
				$genelistarr[$k1][2][5] = $1;
				$genelistarr[$k1][2][6] = $2;}}
		if ($linein =~ m/^VERSION\s+(.+)\s+GI\:(.+)\n/i) {
				$genelistarr[$k1][2][7] = $1;
				$genelistarr[$k1][2][8] = $2;}
		if (($genelistarr[$k1][1][5] eq "PRT") && ($linein =~ m/^\s+$type\s+(.+)\n/i)) {
			$cdsarr[++$cdsarr[0][0]][1] = "$1";
			FETCHCDS:
			while (defined($linein = <GB>)) {
				if ($linein =~ m/^\s+\/gene\=\"(.+)\"/) {
					my $genecheck = uc $1;
					if ($genecheck ne $genelistarr[$k1][0][1]) {
						$cdsarr[$cdsarr[0][0]][3] =  0;}
					else {
						$cdsarr[$cdsarr[0][0]][3] =  1;}
					last FETCHCDS;}
				chomp $linein;
				$cdsarr[$cdsarr[0][0]][1] = $cdsarr[$cdsarr[0][0]][1] . "$linein";}}
		if (($genelistarr[$k1][1][5] eq "RNA") && ($linein =~ m/^\s+$type\s+(.+)\n/i)) {
			$cdsarr[++$cdsarr[0][0]][1] = "$1";
			FETCHCDS2:
			while (defined($linein = <GB>)) {
				if ($linein =~ m/^\s+\/gene\=\"(.+)\"/) {
					my $genecheck = uc $1;
					if ($genecheck ne $genelistarr[$k1][0][1]) {
						$cdsarr[$cdsarr[0][0]][3] =  0;}
					else {
						$cdsarr[$cdsarr[0][0]][3] =  1;}
					last FETCHCDS2;}
				chomp $linein;
				$cdsarr[$cdsarr[0][0]][1] = $cdsarr[$cdsarr[0][0]][1] . "$linein";}}
		if ($linein =~ m/^ORIGIN/i) {
			last FETCHGB;}}
	close GB;
	if (($genelistarr[$k1][2][6] != $genelistarr[$k1][3][3]) || ($genelistarr[$k1][2][5]  != $genelistarr[$k1][3][2])) {
		die "Compatibility prob GB / FASTA for $genelistarr[$k1][0][1] $genelistarr[$k1][2][6] != $genelistarr[$k1][3][3] || $genelistarr[$k1][2][5]  != $genelistarr[$k1][3][2]";}
	if ($genelistarr[$k1][2][4] ne $genelistarr[$k1][3][4]) { 
		die "Compatibility prob GB / FASTA (strand) for  $genelistarr[$k1][0][1] ($genelistarr[$k1][2][4] != $genelistarr[$k1][3][4])";}
	#codex cleaning and filtering
	my @uniqcds;
	my $k2 = 1;
	REVIEWCDS:
	while ($k2 <= $cdsarr[0][0]) {#may be empty??
		if (not defined($cdsarr[$k2][3])) {
			die "Unable to identify gene for $genelistarr[$k1][0][1]";}
		if ($cdsarr[$k2][3] != 1) {
			++$k2;
			next REVIEWCDS;}
		my $cds = "$cdsarr[$k2][1],";
		if ($cds =~ m/[\>\<]/) {#incomplete; skipped
			++$k2;
			next REVIEWCDS;}
		while($cds =~ s/[^\.^\d^\,]//) {}
		while ($cds =~ s/,,/,/) {}
		if ($cds =~ m/^(\d+)(,.+)$/) {
			$cds = "$1\.\.$1" . $2;}
		if ($cds =~ m/^(.+,)(\d+),$/) {
			$cds = $1 . "$2\.\.$2,";}
		if ($cds =~ m/^(.+,)(\d+)(,.+)$/) {
			$cds = $1 . "$2\.\.$2" . $3;}
		my $kc = 1;
		my $foundcds = 0;
		while (defined($uniqcds[$kc])) {
			if ($cds eq $uniqcds[$kc]){
				$foundcds = 1;}
			++$kc;}
		if ($foundcds == 0) {
			$uniqcds[++$uniqcds[0][0]][1] = $cds;
			if (defined($cdsarr[$k2][2])) {
				$uniqcds[$uniqcds[0][0]][2] = $cdsarr[$k2][2];}}
		++$k2;}
	undef @cdsarr;
	@cdsarr = @uniqcds;
	undef @uniqcds;
	$k2 = 1;
	my @tmpcds;
	my $point = 0;
	my $ctallcds = $cdsarr[0][0];
	if (not defined($cdsarr[0][0])) {
		$cdsarr[0][0] = 0;
		print "Problem no valid cds for $genelistarr[$k1][0][1]\n";
		goto SKIPTHISENTRY;}
		
	#print "codex for $genelistarr[$k1][0][1] ($ctallcds CDS)\n";
	while ($k2 <= $cdsarr[0][0]) {
		my $cds = $cdsarr[$k2][1];
		while ($cds =~ m/(\d+)\.+(\d+)\,+/igc) {
			my $string = "$1\t$2";
			my $k3 = 1;
			my $foundit = 0;
			FOUNDIT:
			while (defined($tmpcds[$k3][1])) {
				if ($tmpcds[$k3][1] eq $string) {
					$foundit = 1;
					$tmpcds[$k3][2] = $tmpcds[$k3][2] . "$k2,";
					++$tmpcds[$k3][3];
					last FOUNDIT;}
				++$k3;}
			if ($foundit == 0) {
				++$tmpcds[0][0];
				$tmpcds[$k3][1] = $string;
				$tmpcds[$k3][2] = "$k2,";
				++$tmpcds[$k3][3];}}
		++$k2;}
	#now ranking
	#print " records $tmpcds[0][0]\n";
	my $codexnum = $tmpcds[0][0];
	$k2 = 1;
	my @torank;
	my $pointertorank = 0;
	while ($k2 <= $tmpcds[0][0]) {
		$torank[$pointertorank] = $tmpcds[$k2][1] . "\t$tmpcds[$k2][2]\t$tmpcds[$k2][3]";
		#print "$pointertorank $torank[$pointertorank]\n";
		++$pointertorank;
		++$k2;}
	undef @tmpcds;
	#now ranking
	my $field1 = 1;
	my $field2 = 2;
	my @rankarr;
	undef @rankarr;
	@rankarr = map {$_->[0]} 
	sort {
		my @a_fields = @$a[1..$#$a];
		my @b_fields = @$b[1..$#$b];
		$a_fields[$field1-1] <=> $b_fields[$field1-1]
			||
		$a_fields[$field2-1] <=> $b_fields[$field2-1]}
	map {[$_, split /\t/]} @torank;
	undef @torank;
	#Genes BED file 
	print GENESBED "$genelistarr[$k1][3][1]\t$genelistarr[$k1][3][2]\t$genelistarr[$k1][3][3]\t";
	print GENESBED "$genelistarr[$k1][0][1]\_GENE\@$genelistarr[$k1][2][3]\_$genelistarr[$k1][3][1]\_$genelistarr[$k1][3][2]\_$genelistarr[$k1][3][3]\_$genelistarr[$k1][2][4]\n";
	#Genes FASTA file
	print GENESFASTA "\>$genelistarr[$k1][0][1]\_GENE\t$genelistarr[$k1][2][3]\_$genelistarr[$k1][3][1]\_$genelistarr[$k1][3][2]\_$genelistarr[$k1][3][3]\_$genelistarr[$k1][2][4]\n";
	print GENESFASTA "$genelistarr[$k1][3][6]\n";
	#Codex BED file and FASTA file
	$k2 = 0;
	my $genein = $genelistarr[$k1][3][2];
	my $geneout = $genelistarr[$k1][3][3];
	my $genefullseq = $genelistarr[$k1][3][6];
	my @joincodex;
	while (defined($rankarr[$k2])) {
		if ($rankarr[$k2] =~ m/^(\d+)\t(\d+)\t([^\t]+)\t(.+)$/) {
			my $cdsin = $1;
			my $cdsout = $2;
			my $ctcds = $4;
			my $cdslist = $3;
			my $refcodex = $k2+1;
			my $refinx;
			my $refoutx;
			if ($genelistarr[$k1][2][4] eq '+') {
				$refinx = $cdsin + $genein - 1;
				$refoutx = $cdsout + $genein - 1;}
			else {
				$refinx = $geneout - $cdsout + 1;
				$refoutx = $geneout - $cdsin + 1;}
			if ($refinx > $refoutx) {
				die "Bad conversion of absolute references  ($refinx > $refoutx) gene $genelistarr[$k1][0][1] ";}
			#Codex BED file
			print CODEXBED "$genelistarr[$k1][3][1]\t$refinx\t$refoutx\t";
			print CODEXBED "$genelistarr[$k1][0][1]\_CODEX_$refcodex\@$genelistarr[$k1][2][3]\_$genelistarr[$k1][3][1]\_$refinx\_$refoutx\_$genelistarr[$k1][2][4]\n";
			#Codex FASTA file
			my $codexseq = $genefullseq;
			my $fetchseq = "";
			my $ctcodexnuc = 0;
			FETCHCODEXSEQ:
			while ($codexseq =~ m/(.)/igc) {
				++$ctcodexnuc;
				my $codexnuc = $1;
				if (($ctcodexnuc >= $cdsin) && ($ctcodexnuc <= $cdsout)) {
					$fetchseq = $fetchseq . $codexnuc;
					if ($ctcodexnuc == $cdsout) {
						last FETCHCODEXSEQ;}}}
			if (length($fetchseq) != ($cdsout - $cdsin + 1)) {
				die "Problem fetchiung codex sequence for $genelistarr[$k1][0][1]\_CODEX_$refcodex";}
			print CODEXFASTA "\>$genelistarr[$k1][0][1]\_CODEX_$refcodex\t$genelistarr[$k1][2][3]\_$genelistarr[$k1][3][1]\_$refinx\_$refoutx\_$genelistarr[$k1][2][4]\_$cdsin\_$cdsout\n";
			print CODEXFASTA "$fetchseq\n";
			#building the CDS sequence and info
			while ($cdslist =~ m/(\d+)\,/igc) {
				my $refcds = $1;
				if (not defined($joincodex[$refcds][1])) {
					$joincodex[$refcds][1] = $fetchseq;
					$joincodex[$refcds][2] = "$refcodex,";}
				else {
					$joincodex[$refcds][1] = $joincodex[$refcds][1] . $fetchseq;
					$joincodex[$refcds][2] = $joincodex[$refcds][2] . "$refcodex,";}}
			
			#here Codex table
			#codex frame if protein coding ERROR: Ref in Gene
			#my $frame = 0;
			#if ($genelistarr[$k1][1][5] eq "PRT") {
			#	$frame = sprintf("%.0f",((($cdsin/3 - int($cdsin/3))*3)+1));}		THIS CAN ONLY BE ADDED FROM CDS and DEPENDING on CDS framing so not solvable at codex level so array and then add???
			print CODEXTABLE "$genelistarr[$k1][0][1]\_CODEX_$refcodex\t$genelistarr[$k1][0][1]\t$genelistarr[$k1][1][2]\t$genelistarr[$k1][1][9]\t$genelistarr[$k1][1][3]\t$genelistarr[$k1][0][2]\t$genelistarr[$k1][1][7]";
			print CODEXTABLE "\t$refcodex\t$genelistarr[$k1][2][3]\t$genelistarr[$k1][3][1]\t$refinx\t$refoutx\t$genelistarr[$k1][2][4]";
			print CODEXTABLE "\t$cdsin\t$cdsout\t$ctcds\t$cdslist\t$ctallcds\t$fetchseq\n";
			
			}
		else {
			die "what's this codex entry $rankarr[$k2] for gene $genelistarr[$k1][0][1]\n";}
		++$k2;}
	#CDS FASTA
	$k2 = 1;
	while ($k2 <= $cdsarr[0][0]) {
		if ($genelistarr[$k1][1][5] eq "PRT") {
			if ($joincodex[$k2][1] !~ m/^ATG/i) {
				print "ATG not found \>$genelistarr[$k1][0][1]\_CDS_$k2\n$joincodex[$k2][2]\n";
				$joincodex[$k2][3] = 0; }
			else {
				$joincodex[$k2][3] = 1; }
				
			if ($joincodex[$k2][1] =~ m/(.{3})$/) {
				my $endcds = $1;
				if (($endcds ne 'TGA') && ($endcds ne 'TAG') && ($endcds ne 'TAA')) {
					print "STOP not found \>$genelistarr[$k1][0][1]\_CDS_$k2\n$joincodex[$k2][2]\n";
					$joincodex[$k2][4] = 0; }
				else {
					$joincodex[$k2][4] = 1; }}}
		print CDSFASTA "\>$genelistarr[$k1][0][1]\_CDS_$k2\t$joincodex[$k2][2]\n";
		print CDSFASTA "$joincodex[$k2][1]\n";
		#here Translation and CDS table
		if ($genelistarr[$k1][1][5] eq "PRT") {
			my $seqtotrans = $joincodex[$k2][1];
			my $lenseqtotrans = length($seqtotrans);
			my $pep = "";
			my $oktrans = 1;
			if (($lenseqtotrans/3) != int(($lenseqtotrans/3))) {
				$oktrans = 0;
				print "Problem with coding sequence (length/3) : $genelistarr[$k1][0][1]\_CDS_$k2\n";}
			FETCHAA:
			while ($seqtotrans =~ m/(.{3})/igc) {
				my $codon = $1;
				if (not defined($codonhash{$codon})) {
					print "Problem with coding sequence (unknown codon $codon) : $genelistarr[$k1][0][1]\_CDS_$k2\n";
					$oktrans = 0;
					last FETCHAA;}
				$pep = $pep . $codonhash{$codon};}
			if ($pep =~ m/B.+$/) {
				print "Problem with coding sequence (internal stop codon) : $genelistarr[$k1][0][1]\_CDS_$k2\n";
				$oktrans = 0;}
			print CDSTABLE "$genelistarr[$k1][0][1]\_CDS_$k2\t$genelistarr[$k1][0][1]\t$k2\tPRT\t$joincodex[$k2][2]\t$joincodex[$k2][1]";
			print CDSTABLE "\t$pep\n";}
		else {
			print CDSTABLE "$genelistarr[$k1][0][1]\_CDS_$k2\t$genelistarr[$k1][0][1]\t$k2\tRNA\t$joincodex[$k2][2]\t$joincodex[$k2][1]\t\n";}
		++$k2;}
	#CDS TABLE
	SKIPTHISENTRY:
	print GENESTABLE "$genelistarr[$k1][0][1]\t$genelistarr[$k1][1][2]\t$genelistarr[$k1][1][3]\t$genelistarr[$k1][1][9]\t$genelistarr[$k1][0][2]\t$genelistarr[$k1][1][7]";
	print GENESTABLE "\t$genelistarr[$k1][1][5]\t$genelistarr[$k1][1][4]\t$genelistarr[$k1][1][6]\t$genelistarr[$k1][2][3]\t$genelistarr[$k1][2][7]\t$genelistarr[$k1][3][1]";
	print GENESTABLE "\t$genelistarr[$k1][3][2]\t$genelistarr[$k1][3][3]\t$genelistarr[$k1][3][4]\t$cdsarr[0][0]\t$codexnum\t$genelistarr[$k1][3][6]\n";
	undef @rankarr;			
	undef @cdsarr;
	#if ($k1 == 12) {
	#	die "check stop\n";}
	#print "$k1 ";
	++$k1;}


$showtime = getdatetimeold();
close CODEXBED;
close GENESBED;
close CODEXTABLE;
close GENESTABLE;
close CDSTABLE;
close CODEXFASTA;
close GENESFASTA;
close CDSFASTA;
print "$showtime: GB data loaded\n";

################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
sub monthletter {
	my %xmonthletterhash = ('JAN'=>'01','FEB'=>'02','MAR'=>'03','APR'=>'04','MAY'=>'05','JUN'=>'06','JUL'=>'07','AUG'=>'08','SEP'=>'09','OCT'=>'10','NOV'=>'11','DEC'=>'12',);
	return (%xmonthletterhash);}
sub revcomp {}
sub getdatetime {
    #sub to capture GMT date and time; use it if module DateTime cannot be installed
    my ($sec,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = gmtime;#catch date starting process
    $year = $year + 1900;#$year was number of years since 1900
    $mon = $mon + 1;#first day is zero!!!
    if ($mon < 10) {
	$mon = "0" . $mon;}
    if ($mday < 10) {
	$mday = "0" . $mday;}
    if ($hour < 10) {
	$hour = "0" . $hour;}
    if ($min < 10) {
	$min = "0" . $min;}
    if ($sec < 10) {
	$sec = "0" . $sec;}
    my $date = "$year$mon$mday";
    my $time = "$hour$min$sec";
    my $xdatetime = $date ."_$time";
    return $xdatetime;}
####################
sub getdatetimeold {
#for fetching current date and time
	my $today = DateTime->now();
	return $today->ymd('') . '_' . $today->hms('');}
###############
sub loadcells {
#for loading in msexcel columns
	my %xcellhash = (1 => 'A',2 => 'B',3 => 'C',4 => 'D',5 => 'E',6 => 'F',7 => 'G',8 => 'H',9 => 'I',10 => 'J',11 => 'K',12 => 'L',13 => 'M',14 => 'N',15 => 'O',16 => 'P',17 => 'Q',18 => 'R',19 => 'S',20 => 'T',21 => 'U',22 => 'V',23 => 'W',24 => 'X',25 => 'Y',26 => 'Z',);
	return (%xcellhash);}
#####################	
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);}
################	
sub loadrandom {
my $xpi	= '31415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337867831652712019091456485669234603486104543266482133936072602491412737245870066063155881748815209209628292540917153643678925903600113305305488204665213841469519415116094330572703657595919530921861173819326117931051185480744623799627495673518857527248912279381830119491298336733624406566430860213949463952247371907021798609437027705392171762931767523846748184676694051320005681271452635608277857713427577896091736371787214684409012249534301465495853710507922796892589235420199561121290219608640344181598136297747713099605187072113499999983729780499510597317328160963185950244594553469083026425223082533446850352619311881710100031378387528865875332083814206171776691473035982534904287554687311595628638823';
my @xpiarr;
while ($xpi =~ m/(\d{6})/igc) {
	$xpiarr[++$xpiarr[0]] = $1;}
return (@xpiarr);}
#677968693279877869823280736982826932696779678132535748555051505432666947666976#
################################################################################
