#!/usr/bin/env perl
################################################################################
#677968693279877869823280736982826932696779678132535748555051505432666947666976#
# Using a source database it creates a subset base on HGNC genenames.
#typical command line (see also getopt)
#perl '/home/user/CatchHGNC/CreateGenesDataset_v1i1.pl' --outputfolder '/home/user/CatchHGNC/CUGDtest01' --backupfolder '/media/user/Seagate Backup Plus Drive/CUGDtest01' --geneslist '/home/user/CHU_ULg_Genetic_Diseases/CUGD161231a.csv' --database '/home/user/CatchHGNC/HumGenDB'
#677968693279877869823280736982826932696779678132535748555051505432666947666976#
#some modules required (if not installed: 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';#to return the full path
my $codelocname =  abs_path($0);
my  $starttime = getdatetime();#for reporting sub routine fetching date and time
print "\n$starttime: $code ($codelocname) Process Started from folder $maindir.\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#overwrite
my $overwrite = $ARGV{'--overwrite'};
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#Output folders (main and backup)
my $mainoutputdir  = $ARGV{'--outputfolder'};
if (-e $mainoutputdir) {#exist
	if ($overwrite =~ m/^YES$/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;}}
	else {
		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
				$overwrite = "YES";
				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) {
	if ($overwrite =~ m/^YES$/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;}}
	else {
		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$starttime\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\nOVERWRITE\t$overwrite\n";
close LOG;
copymaintobackup($mainoutputdir,$backupoutputdir);
open (LOG,">>$runlog") ||
	die "Cannot open LOG file $runlog $!\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890	
#Checking source Database
my $database = $ARGV{'--database'};
if (!-e $database) {#does not exists
	die "Unable to find Database location (path) from $database\n";}
if (!-d $database) {#is not a folder
	die "Database location (path) is not a folder ($database)\n";}
#loading an array with all required files
my @tocompile;
#filenames
$tocompile[1][1] = "Genes.txt";
$tocompile[2][1] = "Codex.txt";
$tocompile[3][1] = "ExtendedCodex.txt";
$tocompile[4][1] = "Genes_hg38.bed";
$tocompile[5][1] = "Genes_hg19.bed";
$tocompile[6][1] = "Codex_hg38.bed";
$tocompile[7][1] = "Codex_hg19.bed";
$tocompile[8][1] = "ExtendedCodex_hg38.bed";
$tocompile[9][1] = "ExtendedCodex_hg19.bed";
$tocompile[10][1] = "Products.txt";
$tocompile[11][1] = "Synonymes.txt";
$tocompile[12][1] = "RejectedGenes.txt";
#info about header or not in files
$tocompile[1][2] = 1;
$tocompile[2][2] = 1;
$tocompile[3][2] = 1;
$tocompile[4][2] = 0;
$tocompile[5][2] = 0;
$tocompile[6][2] = 0;
$tocompile[7][2] = 0;
$tocompile[8][2] = 0;
$tocompile[9][2] = 0;
$tocompile[10][2] = 1;
$tocompile[11][2] = 1;
$tocompile[12][2] = 1;
#info about fields to use for ranking
$tocompile[1][3] = 4;
$tocompile[2][3] = 1;
$tocompile[3][3] = 1;
$tocompile[4][3] = 0;
$tocompile[5][3] = 0;
$tocompile[6][3] = 0;
$tocompile[7][3] = 0;
$tocompile[8][3] = 0;
$tocompile[9][3] = 0;
$tocompile[10][3] = 1;
$tocompile[11][3] = 1;
$tocompile[12][3] = 3;
#checking if all files are present
my $kcheck1 = 1;
while (defined($tocompile[$kcheck1][1])) {
	my $filename = "$database\/$tocompile[$kcheck1][1]";
	if (!-e $filename) {#does not exists
		die "Cannot identify file $tocompile[$kcheck1][1]\n";}
	if (!-f $filename) {
		die "File $tocompile[$kcheck1][1] is not a plain txt file\n";}
	++$kcheck1;}
print LOG "Source Database\t$database\n";
print "Source Database identified and validated ($database)\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#loading list of genes
my $geneslist =  $ARGV{'--geneslist'};
if (!-e $geneslist) {#does not exists
	die "Cannot find genes list file $geneslist\n";}
if (!-f $geneslist) {
		die "Genes List File $geneslist is not a plain txt file\n";}
open (GENESLIST,"<$geneslist") ||
	die "Cannot read from genes list file $geneslist\n";
my @tempgenesarr;
my $pointtempgenesarr = 0;
while (defined(my $line = <GENESLIST>)) {
	while($line =~ s/[\n\s\t\r\f]//) {}
	$tempgenesarr[$pointtempgenesarr] = uc $line;
	++$pointtempgenesarr;}
close GENESLIST;
print LOG "Genes List File\t$geneslist\n";
print LOG "Submitted Genes Names\t$pointtempgenesarr\n";
print "Genes List File\t$geneslist\n";
print "Submitted Genes Names\t$pointtempgenesarr\n";
#ranking by genes names (that is also the ranking in all tables from source database
my $field1 = 1;
my @tempranked = 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/]} @tempgenesarr;
undef @tempgenesarr;
my @genesarr;#array to store all submitted entries
#removing duplicates (if)
my $k1 = 0;
REMOVEDUPLI:
while (defined($tempranked[$k1])) {
	if (not defined($genesarr[0][0])) {#first entry
		$genesarr[++$genesarr[0][0]][1] = $tempranked[$k1];
		++$k1;
		next REMOVEDUPLI;}
	if ($genesarr[$genesarr[0][0]][1] ne $tempranked[$k1]) {
		$genesarr[++$genesarr[0][0]][1] = $tempranked[$k1];
		++$k1;
		next REMOVEDUPLI;}
	else {#duplicate
		++$k1;
		next REMOVEDUPLI;}}
undef @tempranked;
print LOG "APPROVED GENES\t$genesarr[0][0]\n";
print "APPROVED GENES\t$genesarr[0][0]\n";
#copying source genes list file
my $copyofgenelist = $geneslist;
if ($copyofgenelist =~ m/\/([^\/]+)$/) {
	$copyofgenelist = "$mainoutputdir\/$1";}
copy ($geneslist,$copyofgenelist) ||
	die "Genes list file: Copy failed! ($copyofgenelist)\n";
#subroutine that copies main to backup (used at different time points)
close LOG;
copymaintobackup($mainoutputdir,$backupoutputdir);
open (LOG,">>$runlog") ||
	die "Cannot open LOG file $runlog $!\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#starting generation of Dataset
#first genes managment to identify the one present or not in the dataset
#loading genes table / creating matched genes table
my $geneoutfile = "$mainoutputdir\/$tocompile[1][1]";
open(GENEOUT,">>$geneoutfile") ||
	die "Cannot create Genes Output Table $geneoutfile\n";
my @notfoundarr;
my $geneinfile = "$database\/$tocompile[1][1]";
open(GENEIN,"<$geneinfile") ||
	die "Cannot read from Genes Input Table $geneinfile\n";
my $headerin = <GENEIN>;
print GENEOUT "$headerin";
my @expgenesinarr;
my @genesinarr;
my $pointgenesinarr = 0;
while (defined(my $line = <GENEIN>)) {
	chomp $line;
	$genesinarr[$pointgenesinarr] = $line;
	$line = $line . "\t";
	my $ctfield = 0;
	while($line =~ m/([^\t]*)\t/igc) {
		$expgenesinarr[$pointgenesinarr][++$ctfield] = uc $1;}
	++$pointgenesinarr;}
close GENEIN;
#so let's match one vs the other and collect the unmatching
my $km1 = 1;
my $km2 = 0;
my $countmatching = 0;
my @matchedarr;
my @unmatchedarr;
$unmatchedarr[0] = 0;
GETMATCHING1:
while(1>0) {
	#	die "if($expgenesinarr[$km2][$tocompile[1][3]] eq $genesarr[$km1][1])\n";
	if((not defined($genesinarr[$km2])) || (not defined($genesarr[$km1][1]))) {
		last GETMATCHING1;}
	if($expgenesinarr[$km2][$tocompile[1][3]] eq $genesarr[$km1][1]) {
		print GENEOUT "$genesinarr[$km2]\n";
		++$countmatching;
		$matchedarr[$countmatching] = $genesarr[$km1][1];
		++$km1;
		++$km2;
		next GETMATCHING1;}
	if ($expgenesinarr[$km2][$tocompile[1][3]] gt $genesarr[$km1][1]) {
		$unmatchedarr[++$unmatchedarr[0]] = $genesarr[$km1][1];
		++$km1;
		next GETMATCHING1;}
	if ($expgenesinarr[$km2][$tocompile[1][3]] lt $genesarr[$km1][1]) {
		++$km2;
		next GETMATCHING1;}}
close GENEOUT;
undef @expgenesinarr;
undef @genesinarr;
while(defined($genesarr[$km1][1])) {
	$unmatchedarr[++$unmatchedarr[0]] = $genesarr[$km1][1];
	++$km1;}
print "$countmatching and $unmatchedarr[0]\n";
#that's it if all are in the Database but
if ($unmatchedarr[0] == 0) {
	print LOG "Genes Identified from Database\t$countmatching\nGenes Not in Database\t0\nGenes Synonymes\tUseless\n";
	print "Genes Identified from Database\t$countmatching\nGenes Not in Database\t0\nGenes Synonymes\tUseless\n";
	print "All submitted Genes have been recovered from source Database.\n";}
else {#so problem with some genes going to stop but informative to correct
	#checking if these are not present in the database (from rejected table)
	my $rejectedgeneoutfile = "$mainoutputdir\/$tocompile[12][1]";
	open(REJGENEOUT,">>$rejectedgeneoutfile") ||
		die "Cannot create Rejected Genes Output Table $rejectedgeneoutfile\n";
	my @rejectednotfoundarr;
	my $rejectedgeneinfile = "$database\/$tocompile[12][1]";
	open(REJGENEIN,"<$rejectedgeneinfile") ||
		die "Cannot read from Rejected Genes Input Table $rejectedgeneinfile\n";
	$headerin = <REJGENEIN>;
	print REJGENEOUT "$headerin";
	my @exprejgenesinarr;
	my @rejgenesinarr;
	my $pointrejgenesinarr = 0;
	while (defined(my $line = <REJGENEIN>)) {
		chomp $line;
		$rejgenesinarr[$pointrejgenesinarr] = $line;
		$line = $line . "\t";
		my $ctfield = 0;
		while($line =~ m/([^\t]*)\t/igc) {
			$exprejgenesinarr[$pointrejgenesinarr][++$ctfield] = uc $1;}
		++$pointrejgenesinarr;}
	close REJGENEIN;
	#so let's match one vs the other and collect the unmatching
	$km1 = 1;
	$km2 = 0;
	my $rejcountmatching = 0;
	my @rejmatchedarr;
	my @rejunmatchedarr;
	$rejunmatchedarr[0] = 0;
	REJGETMATCHING1:
	while(1>0) {
		if((not defined($rejgenesinarr[$km2])) || (not defined($unmatchedarr[$km1]))) {
			last REJGETMATCHING1;}
		if($exprejgenesinarr[$km2][$tocompile[12][3]] eq $unmatchedarr[$km1]) {
			print REJGENEOUT "$rejgenesinarr[$km2]";
			++$rejcountmatching;
			$rejmatchedarr[$rejcountmatching] = $unmatchedarr[$km1];
			++$km1;
			++$km2;
			next REJGETMATCHING1;}
		if ($exprejgenesinarr[$km2][$tocompile[12][3]] gt $unmatchedarr[$km1]) {
			$rejunmatchedarr[++$rejunmatchedarr[0]] = $unmatchedarr[$km1];
			++$km1;
			next REJGETMATCHING1;}
		if ($exprejgenesinarr[$km2][$tocompile[12][3]] lt $unmatchedarr[$km1]) {
			++$km2;
			next REJGETMATCHING1;}}
	close REJGENEOUT;
	undef @exprejgenesinarr;
	undef @rejgenesinarr;
	while(defined($unmatchedarr[$km1])) {
		$rejunmatchedarr[++$rejunmatchedarr[0]] = $unmatchedarr[$km1];
		++$km1;}
	#that's it if the problem is not by using a not gene official symbol
	if ($rejunmatchedarr[0] == 0) {
		print LOG "Genes in Database\t$countmatching\nGenes not in Database\t$rejcountmatching\nGenes synonymes\tUseless\n";
		print "Genes in Database\t$countmatching\nGenes not in Database\t$rejcountmatching\nGenes synonymes\tUseless\n";
		die "Some submitted genes are not validated in the source Database.\nList of these is available from RejectedGenes.txt table.\nRemove these genes from your list and rerun please.\n";}
	else {#so here dealing with synonymes
		my $synfilein = "$database\/$tocompile[11][1]";
		open(SYNIN,"<$synfilein") ||
			die "Cannot read from  Genes Synonymes file $synfilein\n";
		my $synfileout = "$mainoutputdir\/BadHGNC_$tocompile[11][1]";
		open(SYNOUT,">>$synfileout") ||
			die "Cannot create Genes Synonymes file $synfileout\n";
		my @tempsynarr;
		my $pointtempsynarr =0;
		my $linein = <SYNIN>;
		print SYNOUT "$linein";
		while (defined($linein = <SYNIN>)) {
			$tempsynarr[$pointtempsynarr] = $linein;
			++$pointtempsynarr;}
		close SYNIN;
		#ranking (required)
		$field1 = 1;
		my @synarr = 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/]} @tempsynarr;
		undef @tempsynarr;
		#matching 
		my $ksyn1 = 0;
		my $ksyn2 = 1;
		GETSYN:
		while (1>0) {
			if ((not defined($synarr[$ksyn1])) || (not defined($rejunmatchedarr[$ksyn2]))) {
				last GETSYN;}
			my $synname;
			if ($synarr[$ksyn1] =~ m/^([^\t]+)\t/) {
				$synname = uc $1;}
			if ($synname eq $rejunmatchedarr[$ksyn2]) {
				print SYNOUT "$synarr[$ksyn1]";
				if (not defined($synarr[$ksyn1+1])) {
					++$ksyn2;
					last GETSYN;}
				if ($synarr[$ksyn1+1] =~ m/^$rejunmatchedarr[$ksyn2]\t/i) {
					++$ksyn1;
					next GETSYN;}
				else {
					++$ksyn1;
					++$ksyn2;
					next GETSYN;}}
			if ($synname gt $rejunmatchedarr[$ksyn2]) {
				print SYNOUT "$rejunmatchedarr[$ksyn2]\tUnknown Gene Symbol\n";
				++$ksyn2;
				next GETSYN;}
			if ($synname lt $rejunmatchedarr[$ksyn2]) {
				++$ksyn1;
				next GETSYN;}}
		undef @synarr;
		while(defined($rejunmatchedarr[$ksyn2])) {
			print SYNOUT "$rejunmatchedarr[$ksyn2]\tUnknown Gene Symbol\n";
			++$ksyn2;}
		close SYNOUT;
		if ($rejunmatchedarr[0] == 0) {#so only HGNC probs	
			print LOG "Genes in Database\t$countmatching\nGenes not in Database\t0\nGenes Synonymes\t$rejunmatchedarr[0]\n";
			print "Genes in Database\t$countmatching\nGenes not in Database\t0\nGenes Synonymes\t$rejunmatchedarr[0]\n";
			die "Some submitted genes are not using HGNC official gene symbol.\nReview content of the BadHGNC_$tocompile[11][1].\nRemove / Update and rerun.\n";}
		else {#so both not validated and HGNC probs	
			print LOG "Genes Identified from Database\t$countmatching\nGenes Not in Database\t$rejcountmatching\nGenes Synonymes review\t$rejunmatchedarr[0]\n";
			print "Genes Identified from Database\t$countmatching\nGenes Not in Database\t$rejcountmatching\nGenes Synonymes review\t$rejunmatchedarr[0]\n";
			die "Some submitted Genes are not validated in the source Database\nand some others do not use a valid HGNC gene symbol.\nReview both RejectedGenes.txt table and BadHGNC_$tocompile[11][1].\nRemove / Update and rerun.\n";}}}
################################################################################
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
#if here the list of genes submitted is OK for all genes; so let's create subset
#genes hg38 bed (remember bed files dont have header)
#structure: chr19	58346806	58353499	HUM_GEN_A1BG_5
my $bedinfile = "$database\/$tocompile[4][1]";
my @tempbedinarr;
my $pointtempbedinarr = 0;
open(BEDIN,"<$bedinfile") ||
	die "Cannot read from bed infile $bedinfile\n";
while (defined(my $line = <BEDIN>)) {
	if ($line =~ m/^[^\t]+\t[^\t]+\t[^\t]+\t([^\t]+)\n$/) {
		my $bedval = $1;
		if ($bedval =~ m/^[^_]+_[^_]+_([^_]+)_\d+$/) {
			$tempbedinarr[$pointtempbedinarr] = "$1\t$line";
			++$pointtempbedinarr;}
		else {
			die "What's this bed info $bedval\n";}}
	else {
		die "What's this bed entry $line\n";}}
close BEDIN;
print LOG "Genes HG38 Bed database entries\t$pointtempbedinarr\n";
print "Genes HG38 Bed database entries\t$pointtempbedinarr\n";
#ranking on genename
$field1 = 1;
my @bedinarr = 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/]} @tempbedinarr;
undef @tempbedinarr;
my $bedoutfile = "$mainoutputdir\/$tocompile[4][1]";
open(BEDOUT,">>$bedoutfile") ||
	die "Cannot create bed output file $bedoutfile\n";
#matching
my $kb1 = 1;
my $kb2 = 0;
my $bedctm = 0;
my $bedctu = 0;
BEDGENES38:
while (1>0) {
	if ((not defined($genesarr[$kb1][1])) || (not defined($bedinarr[$kb2]))) {
		last BEDGENES38;}
	my $valbedprint;
	my $valbedmatch;
	if ($bedinarr[$kb2] =~ m/^([^\t]+)\t(.+)$/) {
		$valbedprint = $2;
		$valbedmatch = $1;}
	if ($valbedmatch eq $genesarr[$kb1][1]) {
		++$bedctm;
		print BEDOUT "$valbedprint\n";
		++$kb2;
		++$kb1;
		next BEDGENES38;}
	if ($genesarr[$kb1][1] gt $valbedmatch) {
		++$kb2;
		next BEDGENES38;}
	if ($genesarr[$kb1][1] lt $valbedmatch) {
		++$bedctu;
		print "$genesarr[$kb1][1] not found\n";
		++$kb1;
		next BEDGENES38;}}
close BEDOUT;
undef @bedinarr;
while(defined($genesarr[$kb1][1])) {
	++$bedctu;
	print "$genesarr[$kb1][1] not found\n";
	++$kb1;}
print LOG "Genes bed Hg38 Matching\t$bedctm\nGenes bed Hg38 UnMatching\t$bedctu\n";
print "Genes bed Hg38 Matching\t$bedctm\nGenes bed Hg38 UnMatching\t$bedctu\n";
################################################################################
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
#genes hg19 bed (remember bed files dont have header)
#structure: chr19	58858172	58864865	HUM_GEN_A1BG_5chr19
$bedinfile = "$database\/$tocompile[5][1]";
$pointtempbedinarr = 0;
open(BEDIN,"<$bedinfile") ||
	die "Cannot read from bed infile $bedinfile\n";
while (defined(my $line = <BEDIN>)) {
	if ($line =~ m/^[^\t]+\t[^\t]+\t[^\t]+\t([^\t]+)\n$/) {
		my $bedval = $1;
		if ($bedval =~ m/^[^_]+_[^_]+_([^_]+)_\d+$/) {
			$tempbedinarr[$pointtempbedinarr] = "$1\t$line";
			++$pointtempbedinarr;}
		else {
			die "What's this bed info $bedval\n";}}
	else {
		die "What's this bed entry $line\n";}}
close BEDIN;
print LOG "Genes HG19 Bed database entries\t$pointtempbedinarr\n";
print "Genes HG19 Bed database entries\t$pointtempbedinarr\n";
#ranking on genename
$field1 = 1;
@bedinarr = 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/]} @tempbedinarr;
undef @tempbedinarr;
$bedoutfile = "$mainoutputdir\/$tocompile[5][1]";
open(BEDOUT,">>$bedoutfile") ||
	die "Cannot create bed output file $bedoutfile\n";
#matching
$kb1 = 1;
$kb2 = 0;
$bedctm = 0;
$bedctu = 0;
BEDGENES19:
while (1>0) {
	if ((not defined($genesarr[$kb1][1])) || (not defined($bedinarr[$kb2]))) {
		last BEDGENES19;}
	my $valbedprint;
	my $valbedmatch;
	if ($bedinarr[$kb2] =~ m/^([^\t]+)\t(.+)$/) {
		$valbedprint = $2;
		$valbedmatch = $1;}
	if ($valbedmatch eq $genesarr[$kb1][1]) {
		++$bedctm;
		print BEDOUT "$valbedprint\n";
		++$kb2;
		++$kb1;
		next BEDGENES19;}
	if ($genesarr[$kb1][1] gt $valbedmatch) {
		++$kb2;
		next BEDGENES19;}
	if ($genesarr[$kb1][1] lt $valbedmatch) {
		++$bedctu;
		print "$genesarr[$kb1][1] not found\n";
		++$kb1;
		next BEDGENES19;}}
close BEDOUT;
undef @bedinarr;
while(defined($genesarr[$kb1][1])) {
	++$bedctu;
	print "$genesarr[$kb1][1] not found\n";
	++$kb1;}
print LOG "Genes bed Hg19 Matching\t$bedctm\nGenes bed Hg19 UnMatching\t$bedctu\n";
print "Genes bed Hg19 Matching\t$bedctm\nGenes bed Hg19 UnMatching\t$bedctu\n";
################################################################################
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
#Codex hg38 bed (remember bed files dont have header)
#structure: chr19	58353404	58353437	HUM_CDX_A1BG_5_1_chr19_58353404_58353437
$bedinfile = "$database\/$tocompile[6][1]";
$pointtempbedinarr = 0;
open(BEDIN,"<$bedinfile") ||
	die "Cannot read from bed infile $bedinfile\n";
while (defined(my $line = <BEDIN>)) {
	if ($line =~ m/^[^\t]+\t[^\t]+\t[^\t]+\t([^\t]+)\n$/) {
		my $bedval = $1;
		if ($bedval =~ m/^[^_]+_[^_]+_([^_]+)_.+$/) {
			$tempbedinarr[$pointtempbedinarr] = "$1\t$line";
			++$pointtempbedinarr;}
		else {
			die "What's this bed info $bedval\n";}}
	else {
		die "What's this bed entry $line\n";}}
close BEDIN;
print LOG "Codex HG38 Bed database entries\t$pointtempbedinarr\n";
print "Codex HG38 Bed database entries\t$pointtempbedinarr\n";
#ranking on genename
$field1 = 1;
@bedinarr = 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/]} @tempbedinarr;
undef @tempbedinarr;
$bedoutfile = "$mainoutputdir\/$tocompile[6][1]";
open(BEDOUT,">>$bedoutfile") ||
	die "Cannot create bed output file $bedoutfile\n";
#matching
$kb1 = 1;
$kb2 = 0;
$bedctm = 0;
$bedctu = 0;
BEDCODEX38:
while (1>0) {
	if ((not defined($genesarr[$kb1][1])) || (not defined($bedinarr[$kb2]))) {
		last BEDCODEX38;}
	my $valbedprint;
	my $valbedmatch;
	if ($bedinarr[$kb2] =~ m/^([^\t]+)\t(.+)$/) {
		$valbedprint = $2;
		$valbedmatch = $1;}
	if ($valbedmatch eq $genesarr[$kb1][1]) {
		++$bedctm;
		print BEDOUT "$valbedprint\n";
		if (not defined($bedinarr[$kb2+1])) {
			last BEDCODEX38;}
		if ($bedinarr[$kb2+1] =~ m/^$genesarr[$kb1][1]\t/i) {#so next is the same gene
			++$kb2;
			next BEDCODEX38;}
		else {
			++$kb2;
			++$kb1;
			next BEDCODEX38;}}
	if ($genesarr[$kb1][1] gt $valbedmatch) {
		++$kb2;
		next BEDCODEX38;}
	if ($genesarr[$kb1][1] lt $valbedmatch) {
		++$bedctu;
		print "$genesarr[$kb1][1] not found\n";
		++$kb1;
		next BEDCODEX38;}}
close BEDOUT;
undef @bedinarr;
while(defined($genesarr[$kb1][1])) {
	++$bedctu;
	print "$genesarr[$kb1][1] not found\n";
	++$kb1;}
print LOG "Codex bed Hg38 Matching\t$bedctm\nCodex bed Hg38 UnMatching\t$bedctu\n";
print "Codex bed Hg38 Matching\t$bedctm\nCodex bed Hg38 UnMatching\t$bedctu\n";
################################################################################
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
#Codex hg19 bed (remember bed files dont have header)
#structure: chr19	58864770	58864803	HUM_CDX_A1BG_5_1_chr19_58353404_58353437
$bedinfile = "$database\/$tocompile[7][1]";
$pointtempbedinarr = 0;
open(BEDIN,"<$bedinfile") ||
	die "Cannot read from bed infile $bedinfile\n";
while (defined(my $line = <BEDIN>)) {
	if ($line =~ m/^[^\t]+\t[^\t]+\t[^\t]+\t([^\t]+)\n$/) {
		my $bedval = $1;
		if ($bedval =~ m/^[^_]+_[^_]+_([^_]+)_.+$/) {
			$tempbedinarr[$pointtempbedinarr] = "$1\t$line";
			++$pointtempbedinarr;}
		else {
			die "What's this bed info $bedval\n";}}
	else {
		die "What's this bed entry $line\n";}}
close BEDIN;
print LOG "Codex HG19 Bed database entries\t$pointtempbedinarr\n";
print "Codex HG19 Bed database entries\t$pointtempbedinarr\n";
#ranking on genename
$field1 = 1;
@bedinarr = 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/]} @tempbedinarr;
undef @tempbedinarr;
$bedoutfile = "$mainoutputdir\/$tocompile[7][1]";
open(BEDOUT,">>$bedoutfile") ||
	die "Cannot create bed output file $bedoutfile\n";
#matching
$kb1 = 1;
$kb2 = 0;
$bedctm = 0;
$bedctu = 0;
BEDCODEX19:
while (1>0) {
	if ((not defined($genesarr[$kb1][1])) || (not defined($bedinarr[$kb2]))) {
		last BEDCODEX19;}
	my $valbedprint;
	my $valbedmatch;
	if ($bedinarr[$kb2] =~ m/^([^\t]+)\t(.+)$/) {
		$valbedprint = $2;
		$valbedmatch = $1;}
	if ($valbedmatch eq $genesarr[$kb1][1]) {
		++$bedctm;
		print BEDOUT "$valbedprint\n";
		if (not defined($bedinarr[$kb2+1])) {
			last BEDCODEX19;}
		if ($bedinarr[$kb2+1] =~ m/^$genesarr[$kb1][1]\t/i) {#so next is the same gene
			++$kb2;
			next BEDCODEX19;}
		else {
			++$kb2;
			++$kb1;
			next BEDCODEX19;}}
	if ($genesarr[$kb1][1] gt $valbedmatch) {
		++$kb2;
		next BEDCODEX19;}
	if ($genesarr[$kb1][1] lt $valbedmatch) {
		++$bedctu;
		print "$genesarr[$kb1][1] not found\n";
		++$kb1;
		next BEDCODEX19;}}
close BEDOUT;
undef @bedinarr;
while(defined($genesarr[$kb1][1])) {
	++$bedctu;
	print "$genesarr[$kb1][1] not found\n";
	++$kb1;}
print LOG "Codex bed Hg19 Matching\t$bedctm\nCodex bed Hg19 UnMatching\t$bedctu\n";
print "Codex bed Hg19 Matching\t$bedctm\nCodex bed Hg19 UnMatching\t$bedctu\n";
################################################################################
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
#Extended Codex hg38 bed (remember bed files dont have header)
#structure: chr19	58353394	58353447	EXT10_HUM_CDX_A1BG_5_1_chr19_58353394_58353447
$bedinfile = "$database\/$tocompile[8][1]";
$pointtempbedinarr = 0;
open(BEDIN,"<$bedinfile") ||
	die "Cannot read from bed infile $bedinfile\n";
while (defined(my $line = <BEDIN>)) {
	if ($line =~ m/^[^\t]+\t[^\t]+\t[^\t]+\t([^\t]+)\n$/) {
		my $bedval = $1;
		if ($bedval =~ m/^[^_]+_[^_]+_[^_]+_([^_]+)_.+$/) {
			$tempbedinarr[$pointtempbedinarr] = "$1\t$line";
			++$pointtempbedinarr;}
		else {
			die "What's this bed info $bedval\n";}}
	else {
		die "What's this bed entry $line\n";}}
close BEDIN;
print LOG "Extended Codex HG38 Bed database entries\t$pointtempbedinarr\n";
print "Extended Codex HG38 Bed database entries\t$pointtempbedinarr\n";
#ranking on genename
$field1 = 1;
@bedinarr = 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/]} @tempbedinarr;
undef @tempbedinarr;
$bedoutfile = "$mainoutputdir\/$tocompile[8][1]";
open(BEDOUT,">>$bedoutfile") ||
	die "Cannot create bed output file $bedoutfile\n";
#matching
$kb1 = 1;
$kb2 = 0;
$bedctm = 0;
$bedctu = 0;
BEDEXTCODEX38:
while (1>0) {
	if ((not defined($genesarr[$kb1][1])) || (not defined($bedinarr[$kb2]))) {
		last BEDEXTCODEX38;}
	my $valbedprint;
	my $valbedmatch;
	if ($bedinarr[$kb2] =~ m/^([^\t]+)\t(.+)$/) {
		$valbedprint = $2;
		$valbedmatch = $1;}
	if ($valbedmatch eq $genesarr[$kb1][1]) {
		++$bedctm;
		print BEDOUT "$valbedprint\n";
		if (not defined($bedinarr[$kb2+1])) {
			last BEDEXTCODEX38;}
		if ($bedinarr[$kb2+1] =~ m/^$genesarr[$kb1][1]\t/i) {#so next is the same gene
			++$kb2;
			next BEDEXTCODEX38;}
		else {
			++$kb2;
			++$kb1;
			next BEDEXTCODEX38;}}
	if ($genesarr[$kb1][1] gt $valbedmatch) {
		++$kb2;
		next BEDEXTCODEX38;}
	if ($genesarr[$kb1][1] lt $valbedmatch) {
		++$bedctu;
		print "$genesarr[$kb1][1] not found\n";
		++$kb1;
		next BEDEXTCODEX38;}}
close BEDOUT;
undef @bedinarr;
while(defined($genesarr[$kb1][1])) {
	++$bedctu;
	print "$genesarr[$kb1][1] not found\n";
	++$kb1;}
print LOG "Extended Codex bed Hg38 Matching\t$bedctm\nExtended Codex bed Hg38 UnMatching\t$bedctu\n";
print "Extended Codex bed Hg38 Matching\t$bedctm\nExtended Codex bed Hg38 UnMatching\t$bedctu\n";
################################################################################
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
#Extended Codex hg19 bed (remember bed files dont have header)
#structure: chr19	58864760	58864813	EXT10_HUM_CDX_A1BG_5_1_chr19_58353394_58353447
$bedinfile = "$database\/$tocompile[9][1]";
$pointtempbedinarr = 0;
open(BEDIN,"<$bedinfile") ||
	die "Cannot read from bed infile $bedinfile\n";
while (defined(my $line = <BEDIN>)) {
	if ($line =~ m/^[^\t]+\t[^\t]+\t[^\t]+\t([^\t]+)\n$/) {
		my $bedval = $1;
		if ($bedval =~ m/^[^_]+_[^_]+_[^_]+_([^_]+)_.+$/) {
			$tempbedinarr[$pointtempbedinarr] = "$1\t$line";
			++$pointtempbedinarr;}
		else {
			die "What's this bed info $bedval\n";}}
	else {
		die "What's this bed entry $line\n";}}
close BEDIN;
print LOG "Extended Codex HG19 Bed database entries\t$pointtempbedinarr\n";
print "Extended Codex HG19 Bed database entries\t$pointtempbedinarr\n";
#ranking on genename
$field1 = 1;
@bedinarr = 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/]} @tempbedinarr;
undef @tempbedinarr;
$bedoutfile = "$mainoutputdir\/$tocompile[9][1]";
open(BEDOUT,">>$bedoutfile") ||
	die "Cannot create bed output file $bedoutfile\n";
#matching
$kb1 = 1;
$kb2 = 0;
$bedctm = 0;
$bedctu = 0;
BEDEXTCODEX38:
while (1>0) {
	if ((not defined($genesarr[$kb1][1])) || (not defined($bedinarr[$kb2]))) {
		last BEDEXTCODEX38;}
	my $valbedprint;
	my $valbedmatch;
	if ($bedinarr[$kb2] =~ m/^([^\t]+)\t(.+)$/) {
		$valbedprint = $2;
		$valbedmatch = $1;}
	if ($valbedmatch eq $genesarr[$kb1][1]) {
		++$bedctm;
		print BEDOUT "$valbedprint\n";
		if (not defined($bedinarr[$kb2+1])) {
			last BEDEXTCODEX38;}
		if ($bedinarr[$kb2+1] =~ m/^$genesarr[$kb1][1]\t/i) {#so next is the same gene
			++$kb2;
			next BEDEXTCODEX38;}
		else {
			++$kb2;
			++$kb1;
			next BEDEXTCODEX38;}}
	if ($genesarr[$kb1][1] gt $valbedmatch) {
		++$kb2;
		next BEDEXTCODEX38;}
	if ($genesarr[$kb1][1] lt $valbedmatch) {
		++$bedctu;
		print "$genesarr[$kb1][1] not found\n";
		++$kb1;
		next BEDEXTCODEX38;}}
close BEDOUT;
undef @bedinarr;
while(defined($genesarr[$kb1][1])) {
	++$bedctu;
	print "$genesarr[$kb1][1] not found\n";
	++$kb1;}
print LOG "Extended Codex bed Hg19 Matching\t$bedctm\nExtended Codex bed Hg19 UnMatching\t$bedctu\n";
print "Extended Codex bed Hg19 Matching\t$bedctm\nExtended Codex bed Hg19 UnMatching\t$bedctu\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#Done for BED files now Products table
#structure: HUM_CDS_A1BG_5_NP_570602.2	NP_570602.2	GRCh38.p7	20161228_092056	HUM_GEN_A1BG_5	A1BG	PRT	chr19	19Q13.43	chr19	-	8	1;2;3;4;5;6;7;8;	63_96,173_208,303_572,945_1217,1813_2109,2849_3130,5860_6147,6471_6478,	=HYPERLINK("http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=HGNC:5";"HGNC:5_GENE_link")	=HYPERLINK("https://www.ncbi.nlm.nih.gov/nuccore/NP_570602.2";"NP_570602.2_REFSEQ_link")	=HYPERLINK("https://www.ncbi.nlm.nih.gov/projects/sviewer/?id=NP_570602.2";"NP_570602.2_NCBI_link")	1485	ATGTCCATGCTCGTGGTCTTTCTCTTGCTGTGGGGTGTCACCTGGGGCCCAGTGACAGAAGCAGCCATATTTTATGAGACGCAGCCCAGCCTGTGGGCAGAGTCCGAATCACTGCTGAAACCCTTGGCCAATGTGACGCTGACGTGCCAGGCCCACCTGGAGACTCCAGACTTCCAGCTGTTCAAGAATGGGGTGGCCCAGGAGCCTGTGCACCTTGACTCACCTGCCATCAAGCACCAGTTCCTGCTGACGGGTGACACCCAGGGCCGCTACCGCTGCCGCTCGGGCTTGTCCACAGGATGGACCCAGCTGAGCAAGCTCCTGGAGCTGACAGGGCCAAAGTCCTTGCCTGCTCCCTGGCTCTCGATGGCGCCAGTGTCCTGGATCACCCCCGGCCTGAAAACAACAGCAGTGTGCCGAGGTGTGCTGCGGGGTGTGACTTTTCTGCTGAGGCGGGAGGGCGACCATGAGTTTCTGGAGGTGCCTGAGGCCCAGGAGGATGTGGAGGCCACCTTTCCAGTCCATCAGCCTGGCAACTACAGCTGCAGCTACCGGACCGATGGGGAAGGCGCCCTCTCTGAGCCCAGCGCTACTGTGACCATTGAGGAGCTCGCTGCACCACCACCGCCTGTGCTGATGCACCATGGAGAGTCCTCCCAGGTCCTGCACCCTGGCAACAAGGTGACCCTCACCTGCGTGGCTCCCCTGAGTGGAGTGGACTTCCAGCTACGGCGCGGGGAGAAAGAGCTGCTGGTACCCAGGAGCAGCACCAGCCCAGATCGCATCTTCTTTCACCTGAACGCGGTGGCCCTGGGGGATGGAGGTCACTACACCTGCCGCTACCGGCTGCATGACAACCAAAACGGCTGGTCCGGGGACAGCGCGCCGGTCGAGCTGATTCTGAGCGATGAGACGCTGCCCGCGCCGGAGTTCTCCCCGGAGCCGGAGTCCGGCAGGGCCTTGCGGCTGCGGTGCCTGGCGCCCCTGGAGGGCGCGCGCTTCGCCCTGGTGCGCGAGGACAGGGGCGGGCGCCGCGTGCACCGTTTCCAGAGCCCCGCTGGGACCGAGGCGCTCTTCGAGCTGCACAACATTTCCGTGGCTGACTCCGCCAACTACAGCTGCGTCTACGTGGACCTGAAGCCGCCTTTCGGGGGCTCCGCGCCCAGCGAGCGCTTGGAGCTGCACGTGGACGGACCCCCTCCCAGGCCTCAGCTCCGGGCGACGTGGAGTGGGGCGGTCCTGGCGGGCCGAGATGCCGTCCTGCGCTGCGAGGGACCCATCCCCGACGTCACCTTCGAGCTGCTGCGCGAGGGCGAGACGAAGGCCGTGAAGACGGTCCGCACCCCCGGGGCCGCGGCGAACCTCGAGCTGATCTTCGTGGGGCCCCAGCACGCCGGCAACTACAGGTGCCGCTACCGCTCCTGGGTGCCCCACACCTTCGAATCGGAGCTCAGCGACCCTGTGGAGCTCCTGGTGGCAGAAAGC	MSMLVVFLLLWGVTWGPVTEAAIFYETQPSLWAESESLLKPLANVTLTCQAHLETPDFQLFKNGVAQEPVHLDSPAIKHQFLLTGDTQGRYRCRSGLSTGWTQLSKLLELTGPKSLPAPWLSMAPVSWITPGLKTTAVCRGVLRGVTFLLRREGDHEFLEVPEAQEDVEATFPVHQPGNYSCSYRTDGEGALSEPSATVTIEELAAPPPPVLMHHGESSQVLHPGNKVTLTCVAPLSGVDFQLRRGEKELLVPRSSTSPDRIFFHLNAVALGDGGHYTCRYRLHDNQNGWSGDSAPVELILSDETLPAPEFSPEPESGRALRLRCLAPLEGARFALVREDRGGRRVHRFQSPAGTEALFELHNISVADSANYSCVYVDLKPPFGGSAPSERLELHVDGPPPRPQLRATWSGAVLAGRDAVLRCEGPIPDVTFELLREGETKAVKTVRTPGAAANLELIFVGPQHAGNYRCRYRSWVPHTFESELSDPVELLVAES
$bedinfile = "$database\/$tocompile[10][1]";
$pointtempbedinarr = 0;
open(BEDIN,"<$bedinfile") ||
	die "Cannot read from bed infile $bedinfile\n";
my $headerprod = <BEDIN>;
while (defined(my $line = <BEDIN>)) {
	if ($line =~ m/^([^\t]+)\t.+\n$/) {
		my $bedval = $1;
		if ($bedval =~ m/^[^_]+_[^_]+_([^_]+)_.+$/) {
			$tempbedinarr[$pointtempbedinarr] = "$1\t$line";
			++$pointtempbedinarr;}
		else {
			die "What's this bed info $bedval\n";}}
	else {
		die "What's this bed entry $line\n";}}
close BEDIN;
print LOG "Products database entries\t$pointtempbedinarr\n";
print "Products database entries\t$pointtempbedinarr\n";
#ranking on genename
$field1 = 1;
@bedinarr = 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/]} @tempbedinarr;
undef @tempbedinarr;
$bedoutfile = "$mainoutputdir\/$tocompile[10][1]";
open(BEDOUT,">>$bedoutfile") ||
	die "Cannot create bed output file $bedoutfile\n";
print BEDOUT "$headerprod";
#matching
$kb1 = 1;
$kb2 = 0;
$bedctm = 0;
$bedctu = 0;
PRODUCTS:
while (1>0) {
	if ((not defined($genesarr[$kb1][1])) || (not defined($bedinarr[$kb2]))) {
		last PRODUCTS;}
	my $valbedprint;
	my $valbedmatch;
	if ($bedinarr[$kb2] =~ m/^([^\t]+)\t(.+)$/) {
		$valbedprint = $2;
		$valbedmatch = $1;}
	if ($valbedmatch eq $genesarr[$kb1][1]) {
		++$bedctm;
		print BEDOUT "$valbedprint\n";
		if (not defined($bedinarr[$kb2+1])) {
			last PRODUCTS;}
		if ($bedinarr[$kb2+1] =~ m/^$genesarr[$kb1][1]\t/i) {#so next is the same gene
			++$kb2;
			next PRODUCTS;}
		else {
			++$kb2;
			++$kb1;
			next PRODUCTS;}}
	if ($genesarr[$kb1][1] gt $valbedmatch) {
		++$kb2;
		next PRODUCTS;}
	if ($genesarr[$kb1][1] lt $valbedmatch) {
		++$bedctu;
		print "$genesarr[$kb1][1] not found\n";
		++$kb1;
		next PRODUCTS;}}
close BEDOUT;
undef @bedinarr;
while(defined($genesarr[$kb1][1])) {
	++$bedctu;
	print "$genesarr[$kb1][1] not found\n";
	++$kb1;}
print LOG "Products Matching\t$bedctm\nproducts UnMatching\t$bedctu\n";
print "Products Matching\t$bedctm\nproducts UnMatching\t$bedctu\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#Now Codex table
#structure: HUM_CDX_A1BG_5_1	GRCh38.p7	20161228_092056	HUM_GEN_A1BG_5	A1BG	PRT	chr19	19Q13.43	chr19	-	0	1	58353404	58353437	63	96	1	NP_570602.2;	=HYPERLINK("http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=HGNC:5";"HGNC:5_GEN_link")	=HYPERLINK("http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr19:58353404-58353437";"UCSC_CDX_link")	ATGTCCATGCTCGTGGTCTTTCTCTTGCTGTGGG	GRCH37hg19	1	chr19	58864770	58864803	34	1	1	0	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
$bedinfile = "$database\/$tocompile[2][1]";
$pointtempbedinarr = 0;
open(BEDIN,"<$bedinfile") ||
	die "Cannot read from bed infile $bedinfile\n";
my $headerprod = <BEDIN>;
while (defined(my $line = <BEDIN>)) {
	if ($line =~ m/^([^\t]+)\t.+\n$/) {
		my $bedval = $1;
		if ($bedval =~ m/^[^_]+_[^_]+_([^_]+)_.+$/) {
			$tempbedinarr[$pointtempbedinarr] = "$1\t$line";
			++$pointtempbedinarr;}
		else {
			die "What's this bed info $bedval\n";}}
	else {
		die "What's this bed entry $line\n";}}
close BEDIN;
print LOG "Codex database entries\t$pointtempbedinarr\n";
print "Codex database entries\t$pointtempbedinarr\n";
#ranking on genename
$field1 = 1;
@bedinarr = 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/]} @tempbedinarr;
undef @tempbedinarr;
$bedoutfile = "$mainoutputdir\/$tocompile[2][1]";
open(BEDOUT,">>$bedoutfile") ||
	die "Cannot create bed output file $bedoutfile\n";
print BEDOUT "$headerprod";
#matching
$kb1 = 1;
$kb2 = 0;
$bedctm = 0;
$bedctu = 0;
PRODUCTS:
while (1>0) {
	if ((not defined($genesarr[$kb1][1])) || (not defined($bedinarr[$kb2]))) {
		last PRODUCTS;}
	my $valbedprint;
	my $valbedmatch;
	if ($bedinarr[$kb2] =~ m/^([^\t]+)\t(.+)$/) {
		$valbedprint = $2;
		$valbedmatch = $1;}
	if ($valbedmatch eq $genesarr[$kb1][1]) {
		++$bedctm;
		print BEDOUT "$valbedprint\n";
		if (not defined($bedinarr[$kb2+1])) {
			last PRODUCTS;}
		if ($bedinarr[$kb2+1] =~ m/^$genesarr[$kb1][1]\t/i) {#so next is the same gene
			++$kb2;
			next PRODUCTS;}
		else {
			++$kb2;
			++$kb1;
			next PRODUCTS;}}
	if ($genesarr[$kb1][1] gt $valbedmatch) {
		++$kb2;
		next PRODUCTS;}
	if ($genesarr[$kb1][1] lt $valbedmatch) {
		++$bedctu;
		print "$genesarr[$kb1][1] not found\n";
		++$kb1;
		next PRODUCTS;}}
close BEDOUT;
undef @bedinarr;
while(defined($genesarr[$kb1][1])) {
	++$bedctu;
	print "$genesarr[$kb1][1] not found\n";
	++$kb1;}
print LOG "Codex Matching\t$bedctm\ncodex UnMatching\t$bedctu\n";
print "Codex Matching\t$bedctm\ncodex UnMatching\t$bedctu\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#Now Extended Codex table
#structure: HUM_CDX_A1BG_5_1	GRCh38.p7	20161228_092056	HUM_GEN_A1BG_5	A1BG	PRT	chr19	19Q13.43	chr19	-	0	1	58353404	58353437	63	96	1	NP_570602.2;	=HYPERLINK("http://www.genenames.org/cgi-bin/gene_symbol_report?hgnc_id=HGNC:5";"HGNC:5_GEN_link")	=HYPERLINK("http://genome.cse.ucsc.edu/cgi-bin/hgTracks?db=hg38&position=chr19:58353404-58353437";"UCSC_CDX_link")	ATGTCCATGCTCGTGGTCTTTCTCTTGCTGTGGG	GRCH37hg19	1	chr19	58864770	58864803	34	1	1	0	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA	NA
$bedinfile = "$database\/$tocompile[3][1]";
$pointtempbedinarr = 0;
open(BEDIN,"<$bedinfile") ||
	die "Cannot read from bed infile $bedinfile\n";
my $headerprod = <BEDIN>;
while (defined(my $line = <BEDIN>)) {
	if ($line =~ m/^([^\t]+)\t.+\n$/) {
		my $bedval = $1;
		if ($bedval =~ m/^[^_]+_[^_]+_([^_]+)_.+$/) {
			$tempbedinarr[$pointtempbedinarr] = "$1\t$line";
			++$pointtempbedinarr;}
		else {
			die "What's this bed info $bedval\n";}}
	else {
		die "What's this bed entry $line\n";}}
close BEDIN;
print LOG "Extended Codex database entries\t$pointtempbedinarr\n";
print "Extended Codex database entries\t$pointtempbedinarr\n";
#ranking on genename
$field1 = 1;
@bedinarr = 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/]} @tempbedinarr;
undef @tempbedinarr;
$bedoutfile = "$mainoutputdir\/$tocompile[3][1]";
open(BEDOUT,">>$bedoutfile") ||
	die "Cannot create bed output file $bedoutfile\n";
print BEDOUT "$headerprod";
#matching
$kb1 = 1;
$kb2 = 0;
$bedctm = 0;
$bedctu = 0;
PRODUCTS:
while (1>0) {
	if ((not defined($genesarr[$kb1][1])) || (not defined($bedinarr[$kb2]))) {
		last PRODUCTS;}
	my $valbedprint;
	my $valbedmatch;
	if ($bedinarr[$kb2] =~ m/^([^\t]+)\t(.+)$/) {
		$valbedprint = $2;
		$valbedmatch = $1;}
	if ($valbedmatch eq $genesarr[$kb1][1]) {
		++$bedctm;
		print BEDOUT "$valbedprint\n";
		if (not defined($bedinarr[$kb2+1])) {
			last PRODUCTS;}
		if ($bedinarr[$kb2+1] =~ m/^$genesarr[$kb1][1]\t/i) {#so next is the same gene
			++$kb2;
			next PRODUCTS;}
		else {
			++$kb2;
			++$kb1;
			next PRODUCTS;}}
	if ($genesarr[$kb1][1] gt $valbedmatch) {
		++$kb2;
		next PRODUCTS;}
	if ($genesarr[$kb1][1] lt $valbedmatch) {
		++$bedctu;
		print "$genesarr[$kb1][1] not found\n";
		++$kb1;
		next PRODUCTS;}}
close BEDOUT;
undef @bedinarr;
while(defined($genesarr[$kb1][1])) {
	++$bedctu;
	print "$genesarr[$kb1][1] not found\n";
	++$kb1;}
print LOG "Extended Codex Matching\t$bedctm\nExtended codex UnMatching\t$bedctu\n";
print "Extended Codex Matching\t$bedctm\nExtended codex UnMatching\t$bedctu\n";
my  $endtime = getdatetime();
print LOG "COMPILATION COMPLETED BY\t$endtime\n";
close LOG;
copymaintobackup($mainoutputdir,$backupoutputdir);
print "Congratulations. Process started by $starttime and completed by $endtime\n";
die;
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#subroutine that copies files 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;
	NOTADIR:
	while (defined($catchmain[$kmain])) {
		my $mainfile = "$main\/$catchmain[$kmain]";
		if (-d $mainfile) {
			++$kmain;
			next NOTADIR;} 
		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";}
################################################################################
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
sub getdatetime {
#for fetching current date and time
	my $today = DateTime->now();
	return $today->ymd('') . '_' . $today->hms('');}
################################################################################
#677968693279877869823280736982826932766967796781326669765357485550515055495154#
=head1 NAME

CreateGenesDataset_v1i1.pl (version 1, issue 1)

=head1 HELP

This code has been developped on a Linux (Ubuntu 14.04) 64bit PC platform.
(RAM 16 Gb, 4 cores / 8 CPUs)

This codes requires Perl to be installed and some Perl Modules:

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

This code create a Dataset from a source Database (generated using perl code CatchHGNC_v1i1.pl)

It requires (provided as required arguments):
	a Tab delimited csv file  (argument --geneslist) containing (one per line) all genes HGNC official symbol that should be in the Dataset
	a main output folder (argument --outputfolder) and a backup output folder (argument --backupfolder) to record output files for the Dataset
	a Database folder (argument --database) with all expected files listed below:
		Genes.txt
		Codex.txt
		ExtendedCodex.txt
		Genes_hg38.bed
		Genes_hg19.bed
		Codex_hg38.bed
		Codex_hg19.bed
		ExtendedCodex_hg38.bed
		ExtendedCodex_hg19.bed
		Products.txt
		Synonymes.txt
		RejectedGenes.txt

For additionnal info. contact the code Owner/Author (or from WEB: www.ngyx.eu)

=head1 VERSION

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

=head1  USAGE

	perl CreateGenesDataset_v1i1.pl [options]  
	
Example (no optional argument):
	
perl '/home/user/CatchHGNC/CreateGenesDataset_v1i1.pl' \
--outputfolder '/home/user/CatchHGNC/CUGDtest01' \
--backupfolder '/media/user/Seagate Backup Plus Drive/CUGDtest01' \
--geneslist '/home/user/CHU_ULg_Genetic_Diseases/CUGD161231a.csv'
--database '/home/user/CatchHGNC/HumGenDB'

=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 [=] <geneslist>

Path of genes list file

=for Euclid:
	geneslist.type: string

=item --database [=] <database>

Path to folder with all required Database input files

=for Euclid:
	database.type: string
	
=back

=head1 OPTIONS

=over

=item --overwrite [=] <overwrite>

Yes or No; Default is No.

=for Euclid:
	overwrite.type: string 
 	overwrite.default: "No"
 	
=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
