#!/usr/bin/perl
use strict;
=Instructions about the tool
Author: Iliyas Rashid
Program name: HREFinder.pl
_______________________________________________________________________________
DESCRIPTION: Tool for the identification and localization of
	(I)  all canonical Hypoxia-Response Elements (HRE) 
	(II) all functional HRE (localization of Hypoxia-Anciliary Sequences (HAS) towards downstream from HRE on spaces of 7-15 nt)

PROGRAM EXECUTION:

	#perl HREFinder.pl <FASTA file>
	<FASTAfile>    Single file in FASTA format containing the sequence(s).
	
OUTPUT
	This program generates three output files  in the working directory with the prefix of the file name of the input sequence
	'sequence_AllMinedHRE.txt' contains all identified HRE in both strands
	'sequence_FunctionalHRE.txt' contains all predicted functional HRE in both strands
	'sequence_HREFrequencyTable.txt' contains all frequencies of HRE distribution in every 500 nt consecutive fragment set of supplied sequences
EXAMPLE
	The 'TestDataSet' directory contains input sequence files which were used in test studies
	The 'ResultsTestDataSet' directory contains output files of all example input sequences
 _______________________________________________________________________________

=cut
##########################################################################################################################################
my$infile=$ARGV[0] or die"Please provide input sequences in a FASTA format file in .text/.fasta file extension for example  #perl HREFinder.pl upstreamseqfile.fasta\n";
if($infile!~/(.+)\.(txt|fasta)/)
{
	die "\nProgram acceptes only input file with  '.fasta' or '.txt'  extension\n";
}
my$MinedHRE=$infile;
$MinedHRE=~s/^.+\/(\S+)$/$1/gi;
$MinedHRE=~s/(.+)\.(txt|fasta)/$1\_AllMinedHRE\.txt/g;
my$FunctionalHRE=$MinedHRE;
my$FrequencyTable=$MinedHRE;
$FunctionalHRE=~s/(.+)\_AllMinedHRE.+/$1\_FunctionalHRE\.txt/g;
$FrequencyTable=~s/(.+)\_AllMinedHRE.+/$1\_HREFrequencyTable\.txt/g;
open(HRE,">$MinedHRE") or die "'$MinedHRE' file cant open for writing: $1\n";
open(FUNCHRE,">$FunctionalHRE") or die "'$FunctionalHRE' file cant open for writing: $1\n";
open(FREQ,">$FrequencyTable") or die "'$FrequencyTable' file cant open for writing: $1\n";
my$found=0;
my$notfound=0;
open(IN,$infile) or die"$infile cant open for reading $!";
my@contents=<IN>;
close IN;
my$strings=join('',@contents);
$strings=~s/\>/>>>/g;
my@seqfiles=split('>>',$strings);
shift @seqfiles;
my%Consenscollect=();
my%RevConsenscollect=();
my$HREsubseqsize=30;
print "\nA total ".scalar @seqfiles." sequence(s) have been supplied in file '$infile' \n";
foreach my$seq(@seqfiles)
{$found++;
	my$header='';	
	if($seq=~/(\>.+)/g)
	{
		$header=$1;
		$header=~s/[\n\t\r\s]+/ /g;
	}
	$seq=~s/\n//g;
	$seq=~s/$header//g;
	$seq=~s/[\n\r\t\s]+//g;	
	my@consens1=motifpos($seq,'ACGTG');
	my@consens2=motifpos($seq,'GCGTG');
	foreach(@consens1)
	{	my$hreseq=substr($seq,$_,$HREsubseqsize);
		$Consenscollect{"$header\tACGTG\t$_\t$hreseq"}=$seq;
	}
	foreach(@consens2)
	{
		my$hreseq=substr($seq,$_,$HREsubseqsize);		
		$Consenscollect{"$header\tGCGTG\t$_\t$hreseq"}=$seq;
	}
	#mining HRE in reverse sequences
	my$reverseq=revcomplmnt($seq);
	my@revconsens1=motifpos($reverseq,'ACGTG'); 
	my@revconsens2=motifpos($reverseq,'GCGTG');
	foreach(@revconsens1)
	{	my$hreseq=substr($reverseq,$_,$HREsubseqsize);	
		$RevConsenscollect{"$header\tACGTG\t$_\t$hreseq"}=$reverseq;
	}
	foreach(@revconsens2)
	{
		my$hreseq=substr($reverseq,$_,$HREsubseqsize);
		$RevConsenscollect{"$header\tGCGTG\t$_\t$hreseq"}=$reverseq;
	}
}
my$hrecount=0;
my$revhrecount=0;
my$funchrecount=0;
my$revfunchrecount=0;
my@HREpos=();
my@RevHREpos=();
print HRE"S.No.\tSequence_definition\tHRE_location_on_sequence\tHRE sequence\tHRE (A|G)CGTG)\n";
print FUNCHRE"S.No.\tSequence_definition\tHRE_location_on_sequence\tHRE sequence\tHRE (A|G)CGTG)\tHAS (CA(C|G)(A|G)(T|G|C))\n";
foreach (sort keys %Consenscollect)
{$hrecount++;
	
	my($header,$hre,$hrestrpos,$hreseq)=split(/\t/,$_);
	my$hreend=$hrestrpos+$HREsubseqsize;
	push(@HREpos,$hrestrpos);
	print HRE"$hrecount\t$header\t$hrestrpos..$hreend\t$hreseq\t$hre\n";
	if($hreseq=~/((A|G)CGTG).{7,15}(CA(C|G)(A|G)(T|G|C))/g)
	{$funchrecount++;
		my$fhre=$1;
		my$has=$3;		
		print FUNCHRE"$funchrecount\t$header\t$hrestrpos..$hreend\t$hreseq\t$fhre\t$has\n";
	}
}
print HRE"\nHRE in reverse strand \n--------------------------------------------------\n";
print FUNCHRE"\nFunctional HRE in reverse strand \n--------------------------------------------------\n";
print HRE"S.No.\tSequence_definition\tHRE_location_on_sequence\tHRE sequence\tHRE (A|G)CGTG)\n";
print FUNCHRE"S.No.\tSequence_definition\tHRE_location_on_sequence\tHRE sequence\tHRE (A|G)CGTG)\tHAS (CA(C|G)(A|G)(T|G|C))\n";
foreach (sort keys %RevConsenscollect)
{$revhrecount++;
	
	my($header,$hre,$hrestrpos,$hreseq)=split(/\t/,$_);
	my$hreend=$hrestrpos+$HREsubseqsize;
	push(@RevHREpos,$hrestrpos);	
	print HRE"$revhrecount\t$header\t$hrestrpos..$hreend\t$hreseq\t$hre\n";
	if($hreseq=~/((A|G)CGTG).{7,15}(CA(C|G)(A|G)(T|G|C))/g)
	{$revfunchrecount++;
		my$fhre=$1;
		my$has=$3;		
		print FUNCHRE"$revfunchrecount\t$header\t$hrestrpos..$hreend\t$hreseq\t$fhre\t$has\n";
	}
}
print"\nHREs mined in original strand of input sequences: $hrecount\nHREs mined in reverse strand of input sequences: $revhrecount\n";
print"\nAll possible functional HREs in original strand: $funchrecount\nAll possible functional HREs in reverse strand: $revfunchrecount\n";
print"Three output files have been written in the working directory\n\t'$MinedHRE'\n\t'$FunctionalHRE' \n\t'$FrequencyTable' \n\n********************** End **********************\n";
my%HREPOS=HREFreq(@HREpos);
foreach (sort{$a <=> $b} keys %HREPOS)
{
	print FREQ"$_\t$HREPOS{$_}\n";

}
print FREQ"Reverse strand\n";
my%RevHREPOS=HREFreq(@RevHREpos);
foreach (sort{$a  <=> $b} keys %RevHREPOS)
{
	print FREQ"$_\t$RevHREPOS{$_}\n";

}
########################### methods ####################################
sub BreakSequence
{
	my$seq=shift;
	my $length = length $seq;
	my$breakseq='';
	for(my$i=0;$i<=$length; $i=$i+120)
	{
		my $substr=substr($seq,$i,120);
		$breakseq.="$substr\n";
	}
	return $breakseq;
}
sub HREFreq
{
	my@pos=@_;
	@pos=sort{ $b <=> $a }@pos;
	my$maxlength=$pos[0];
	@pos=sort{ $a <=> $b }@pos;	
	my%HREPOS=();
	my$start=1;
	my$range=500;
	my$i=0;
	for($start;$start<=$maxlength;$start+=500)
	{$i++;
		#print"\t\tLoop $i--- start:$start rang: $range\n";
		foreach my$strtpos(@pos)
		{	
			if($strtpos>=$start && $strtpos<=$range)
			{
				if(exists $HREPOS{"$start-$range"})
				{
					$HREPOS{"$start-$range"}++;
				}
				else{					
					$HREPOS{"$start-$range"}=1;
				}							
			}
		}
		$range+=500;
	}
	return (%HREPOS);
}

sub motifpos{
	my$string=shift;
	my$char = shift;
	my $offset = 0;
	my@poss=();
	my $result = index($string, $char, $offset);
	#print"$result";
	#Perl's index() will return -1 if there was no match
	while ($result != -1)
	{
		push(@poss,$result);
		#print "Found $char at $result\n";
		$offset = $result + 1;
		$result = index($string, $char, $offset);		
	}
	return @poss;
}
sub revcomplmnt
{
	my$seq=shift;
	$seq=~s/[\n\r\t\s]+//g;
	my$revscompstr=reverse($seq);
	$revscompstr=~tr/atcgATCG/tagcTAGC/;
	return ($revscompstr);
}
