#!/usr/bin/env perl
################################################################################
#677968693279877869823280736982826932696779678132535748555051505432666947666976#
# This code compares nuc by nuc 2 bed files and delivers coverages one by the other
################################################################################
#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
#Output Folder
my $outputdir = $maindir . "\/BTB_CUGDC_CSID";
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;
my $k1 = 0;
while (defined($catcharr2[$k1])) {
	my $outfile = $outputdir . "\/$catcharr2[$k1]";
	if (-f $outfile) {
		unlink $outfile ||
			die "Cannot erase output file $outfile";}
	++$k1;}
my $outfile = $outputdir . "\/Coverage.txt";
open(INFO,">>$outfile") ||
	die "Cannot create ouput file $outfile";
print INFO "RUN_DATE_TIME\t$rundatetime\nPERL_CODE\t$code\nMAIN_DIR\t$maindir\nOUTPUT_DIR\t$outputdir\n";
#print INFO "Chromosome\tStart\tEnd\tCheckOK\tLength\tEntryInfo\tMatchCount\tPctMatch\tMatchType\tMaxUncover\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#loading covered target file
#Input files used:
#CSID Improved Regions "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/CSID_Improved_Regions.bed";
#CSID Agilent Regions "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/CSID_Regions.bed";
#Fake1 "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/Fake1.bed";
#Fake2 "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/Fake2.bed";
#CUGDG HG19 "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/CUGDG_HG19.bed";
#CUDGC HG19 "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/CUGDC_HG19.bed";
#CUGDC HG38 "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/CUGDC_HG38.bed";
#V5 Agilent Regions  "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/V5_Regions.bed";
#V6r2 Agilent Regions  "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/V6r2_Regions.bed";
#Agilent SureDesign 798 "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/AgilentSureDesign_798_Regions.bed";
my $csidfile = "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/CUGDC_HG19.bed";
my @tempcsidarr;
my $pointer = 0;
my $ctskippedcsid = 0;
if (-e $csidfile) {
	open (CSID,"<$csidfile") ||
		die "Cannot read from CSID file $csidfile";
		while (defined(my $linein = <CSID>)) {
			#$linein = uc $linein;
			chomp $linein;
			if ($linein =~ m/CHR([^\t]+)\t(\d+)\t(\d+)\t([^\t]{1,128})/i) {
				my $infogene = $4;
				my $chr = uc $1;
				my $posin = $2;
				my $posout = $3;
				while ($infogene =~ s/[\t\n\f\r]//) {}
				my $strand = "+";
				if ($posout < $posin) {
					my $posinv = $posout;
					$posout = $posin;
					$strand = "-";
					$posin = $posinv;}
				if ($posout < $posin) {
					die "What's this line: $linein";}
				my $len = $posout - $posin +1;
				my $rankchr;
				if ($chr =~ m/^\d$/i) {
					$rankchr = "CHR0$chr";}
				else {
					$rankchr =  "CHR$chr";}
				$rankchr = uc $rankchr;
				if (($rankchr !~ m/^CHR\d{2}$/) && ($rankchr !~ m/^CHR[XY]$/)) {
					die "What's this Rank Chr $rankchr in $linein";}
				$tempcsidarr[$pointer] = "$rankchr\t$posin\t$posout\t$strand\t$len\tchr$chr:$posin\-$posout,$infogene\t";
				++$pointer;}
			else {#skipped
				++$ctskippedcsid;
				}}
		close CSID;}
else {
	die "Cannot find CSID file $csidfile";}
#ranking the array
my $field1 = 1;
my $field2 = 2;
my $field3 = 3;
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/]} @tempcsidarr;
undef @tempcsidarr;
@tempcsidarr = @rankarr;
undef @rankarr;
#loading data
my @csidarr;
$k1 = 0;
while (defined($tempcsidarr[$k1])) {
	if ($tempcsidarr[$k1] =~ m/^([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t$/) {
		$csidarr[++$csidarr[0][0]][1] = $1;
		$csidarr[$csidarr[0][0]][2] = $2;
		$csidarr[$csidarr[0][0]][3] = $3;
		$csidarr[$csidarr[0][0]][4] = $4;
		$csidarr[$csidarr[0][0]][5] = $5;
		$csidarr[$csidarr[0][0]][6] = $6;}
	else {
		die "What's this raned entry: $tempcsidarr[$k1]";}
	++$k1;}
undef @tempcsidarr;
print "Found and ranked $csidarr[0][0] (from not ranked $pointer) entries from CSID file $csidfile\n";
print INFO "COVERED_BED\t$csidfile\nENTRIES\t$csidarr[0][0]\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#loading the other Bed file data (covering)
my $bedfile =  "/media/user/Seagate Backup Plus Drive/ValidatedBedToBed/InputFiles/CSID_Regions.bed";
my @tempbedarr;
my $ctskippedbed = 0;
$pointer = 0;
if (-e $bedfile) {
	open (BED,"<$bedfile") ||
		die "Cannot read from BED file $bedfile";
		while (defined(my $linein = <BED>)) {
			#$linein = uc $linein;
			chomp $linein;
			if ($linein =~ m/CHR([^\t]+)\t(\d+)\t(\d+)\t([^\t]{1,128})/i) {
				my $infogene = $4;
				my $chr = uc $1;
				my $posin = $2;
				my $posout = $3;
				while ($infogene =~ s/[\t\n\f\r]//) {}
				my $strand = "+";
				if ($posout < $posin) {
					my $posinv = $posout;
					$posout = $posin;
					$strand = "-";
					$posin = $posinv;}
				if ($posout < $posin) {
					die "What's this line: $linein";}
				my $len = $posout - $posin +1;
				my $rankchr;
				if ($chr =~ m/^\d$/i) {
					$rankchr = "CHR0$chr";}
				else {
					$rankchr =  "CHR$chr";}
				$rankchr = uc $rankchr;
				if (($rankchr !~ m/^CHR\d{2}$/) && ($rankchr !~ m/^CHR[XY]$/)) {
					die "What's this Rank Chr $rankchr in $linein";}
				$tempbedarr[$pointer] = "$rankchr\t$posin\t$posout\t$strand\t$len\tchr$chr:$posin\-$posout,$infogene\t";
				++$pointer;}
			else {
				++ $ctskippedbed;
				#just skipped
				}}
		close BED;}
else {
	die "Cannot find BED file $bedfile";}
undef @rankarr;
#ranking the array
$field1 = 1;
$field2 = 2;
$field3 = 3;
@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/]} @tempbedarr;
undef @tempbedarr;
@tempbedarr = @rankarr;
undef @rankarr;
#loading data
my @bedarr;
$k1 = 0;
while (defined($tempbedarr[$k1])) {
	if ($tempbedarr[$k1] =~ m/^([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t([^\t]+)\t$/) {
		$bedarr[++$bedarr[0][0]][1] = $1;
		$bedarr[$bedarr[0][0]][2] = $2;
		$bedarr[$bedarr[0][0]][3] = $3;
		$bedarr[$bedarr[0][0]][4] = $4;
		$bedarr[$bedarr[0][0]][5] = $5;
		$bedarr[$bedarr[0][0]][6] = $6;
		#die "$bedarr[$bedarr[0][0]][1] $bedarr[$bedarr[0][0]][2] $bedarr[$bedarr[0][0]][3] $bedarr[$bedarr[0][0]][4] $bedarr[$bedarr[0][0]][5] $bedarr[$bedarr[0][0]][6]\n";
		}
	else {
		die "What's this ranked entry: $tempbedarr[$k1]";}
	++$k1;}
undef @tempbedarr;
print "Found and ranked $bedarr[0][0] (from not ranked $pointer) entries from BED file $bedfile\n";
print INFO "COVERING_BED\t$bedfile\nENTRIES\t$bedarr[0][0]\n";
################################################################################
#2345678901234567890123456789012345678901234567890123456789012345678901234567890
#starting comparison of bed entries for coverage
print INFO "Chromosome\tStart\tEnd\tCheckOK\tLength\tEntryInfo\tMatchCount\tPctMatch\tMatchType\tMaxUncover\n";
$showtime = getdatetimeold();#for reporting
print "\n$showtime: Starting evaluation.\n";
$k1 = 1;
while ($k1 <= $csidarr[0][0]) {
	my $chrc  = $csidarr[$k1][1];
	my $posinc = $csidarr[$k1][2];
	my $posoutc = $csidarr[$k1][3];
	my $size = $csidarr[$k1][5];
	my @collectarr;
	my $k2 = 1;
	while ($k2 <= $size) {#initialization at 0
		$collectarr[$k2] = 0;
		++$k2;}
	$k2 = 1;
	FETCHMATCH:
	while ($k2 <= $bedarr[0][0]) {
		if ($bedarr[$k2][1] lt $chrc) {
			++$k2;
			next FETCHMATCH;}
		elsif ($bedarr[$k2][1] gt $chrc) {
			last FETCHMATCH;}
		else {
			if ($bedarr[$k2][3] < $posinc) {
				++$k2;
				next FETCHMATCH;}
			elsif ($bedarr[$k2][2] > $posoutc) {
				last FETCHMATCH;}
			else {
				if (($bedarr[$k2][2] <= $posinc) && ($bedarr[$k2][3] >= $posoutc)) {#fully encompassing so done in one shoot
						my $k4 = 1;
						while ($k4 <= $size) {
							 ++$collectarr[$k4];
							 ++$k4;}
						last FETCHMATCH;}
				#to fetch counter in
				my $counter = 0;
				my $startpos = $bedarr[$k2][2];
				if ($posinc <= $startpos) {#inner start
					$counter = $startpos - $posinc;}
				else {
					while ($startpos >= $posinc) {
						++$startpos;}}
				#print "$startpos and $bedarr[$k2][3] vs  $posinc and $posoutc\n"; 
				MATCHPOS:
				while ($startpos <= $bedarr[$k2][3]) {
					if ($startpos < $posinc) {
						++$startpos;}
					elsif ($startpos > $posoutc) {
						last MATCHPOS;}
					else {
						++$collectarr[++$counter];
						++$startpos;}}
					++$k2;
					next FETCHMATCH;}}}
	#reviewing collectarr;
	my $matched = 0;
	my $uncover = 0;
	my $maxuncover = 0;
	$k2 = 1;
	while ($k2 <= $size) {
		if ($collectarr[$k2] > 0) {
			++$matched;
			if ($uncover > $maxuncover) {
				$maxuncover = $uncover;}
			$uncover = 0;}
		else {
			++$uncover;}
		++$k2;}
	if ($uncover > $maxuncover) {
		$maxuncover = $uncover;}
	my $dif = $maxuncover;
	my $matchtype;
	if ($matched == 0) {
		$matchtype = "Missing";}
	else {
		if ($dif == 0) {
			$matchtype = "Covered";}
		elsif ($dif <= 25) {
			$matchtype = "Padded25";}
		elsif ($dif <= 50) {
			$matchtype = "Padded50";}
		elsif ($dif <= 75) {
			$matchtype = "Padded75";}
		elsif ($dif <= 100) {
			$matchtype = "Padded100";}
		else {
			$matchtype = "LowCoverage";}}
	my $ratio = $matched/$size;
	$ratio = sprintf("%.2f",$ratio*100);
	print INFO "$csidarr[$k1][1]\t$csidarr[$k1][2]\t$csidarr[$k1][3]\t$csidarr[$k1][4]\t$csidarr[$k1][5]\t$csidarr[$k1][6]";
	print INFO "\t$matched\t$ratio\t$matchtype\t$maxuncover\n";
	++$k1;}
$showtime = getdatetimeold();#for reporting
print "\n$showtime: Evaluation completed\n";
print INFO "SKIPPED_COVERED\t$ctskippedcsid\nSKIPPED_COVERING\t$ctskippedbed\nRUN_COMPLETED_BY\t$showtime\n";
close INFO;
################################################################################
#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#
################################################################################
