#!/usr/bin/perl
#change the line above on Unix OS usually to: !usr/bin/perl or wherever perl.exe lies
$| = 1;

print "Content-type: text/html\n\n";
# ----------------------------------------- Script Summary ------------------------------------------------------------------
# This perl script receives the data sent by the input form in DistinctiEnz.pl; the data includes:
# 	1) An uploaded file containing a set of FASTA sequences;
# 	2) A list of enzyme names; and
# 	3) Some option data for RFLP Gel dilation factor and sequence filtering of non-atcg characters.

# The script will perform an RFLP analysis with enzyme(s) selected by user on the submitted set of sequences and
# prints out the analysis as a Resolution Profile and a virtual RFLP Gel view.
# The main logic is that:
# 1) We get some fragments and fragment sizes from digest of each sequence.
# 2) Each fragment size actually will appear as aband on RFLP Gel.
# 3) Each sequence will have a column in RFLP Gel.
# 4) If each column has a unique band feature then the resolution will be 100%.
# 5) A unique band feature depends on 2 factors. One is the unique number of bands and the other is 
#    the unique set of fragment sizes when the number of fragments prouced by some sequences are equal so that
#    the bands produced will appear in different distances and therefore producing distinguishable patterns.

# Three subroutines were defined to perform repetitive tasks and make the program more readable.

# -------------------------------------- Script dependencies ---------------------------------------------------------
# This script uses CGI and GD modules of Perl and the following Bioperl modules:
# 1) Bio::Restriction::Analysis
# 2) Bio::Restriction::EnzymeCollection
# 3) Bio::Restriction::Enzyme
# 4) Bio::PrimarySeq
# 5) Bio::SeqIO
# However some basic Bioperl modules may also be needed; therefore, it is recommended that Bioperl
# be installed on your server.
# It is strongly recommended to study Bioperl documentation on the above-mentioned modules in http://www.bioperl.org.


# ---------------------------------------------------------------------------------------------------------------





# ------------------------------------------------ Main Script ----------------------------------------------
# Use CGI module of Perl; this module is the standard and reliable way of handling data from forms
#  and producing html output. First a new instance of CGI object is created as $query.
use CGI;
$query=new CGI;
#In Apache 2.0 on Windows I learned it is better to get a non-file-object parameter from the input form
#This will also ensue that no CGItemp file will remain in the folder. Strange but it happened!!!
#It is why I get enzymeNames before fileSeq object data !!!!
# Get enzyme(s) selected
@enzymeNames=$query->param("cboEnzymes");

# Get the raw data for analysis from the input form in IncisiEnz.pl
# Get sequence data

$uploadFileHandle = $query->upload("fileSeq");
$uploadFileName=$query->param("fileSeq");




# Get PowerCleaner, sequence type (circular), and dilation factor options 
$dilationFactor=$query->param("cboDilate");
$powerCleaner=$query->param("ckCleaner");
$seqType=$query->param("ckCircular");

$doBatchAnalysis=$query->param("ckBatch");
$ladderSet=$query->param("cboLadder");


# use print<< command for out-printing multiline strings in Perl. We produce HTML headers and a navigation bar
print <<MULTILNBL;

<html><head><title>DistinctiEnz Analysis Results</title>
<script language='javascript'>
	var IE = document.all? true : false; 
	if (!IE) document.captureEvents(Event.MOUSEMOVE); 
	document.onmousemove = getMousexy; 

	function getMousexy(e) { 
		if(!e) { 
			if( window.event ) { 
				//DOM 
				e = window.event; 
			} 
			else { 
				//TOTAL FAILURE, WE HAVE NO WAY of REFERENCING THE EVENT 
				return; 
			} 
		} 
		var tempX = 0; 
		var tempY = 0; 
		
		if(typeof(e.offsetX) == 'number'){ 
			var offset_x = e.offsetX; 
			var offset_y = e.offsetY; 
		} else if( typeof(e.clientx) == 'number') { 
			var offset_x = e.clientx; 
			var offset_y = e.clientY; 
		} 
		if(e.srcElement.id=='pic'){
			var pic=document.getElementById('pic');
			var marker=document.getElementById('marker_h');
			var label=document.getElementById('y_offset');
			var y = offset_y - 27;
			var bp=Math.pow(10, (4 - (y/$dilationFactor)));
			
			document.getElementById('y_offset').value = Math.round(bp) + ' bp';
			if(offset_y > 27){ 
				marker.style.visibility='visible';			
				label.style.visibility='visible';						
			}
			else{
				marker.style.visibility='hidden';			
				label.style.visibility='hidden';				
			}
			document.getElementById('marker_h').style.width=pic.width;
			document.getElementById('marker_h').style.top= offset_y;
			document.getElementById('y_offset').style.top= offset_y;
			document.getElementById('y_offset').style.left= offset_x;
		}
		return true; 
	} 
	
	function hideReport(){
		if(document.all.reporting){	document.all.reporting.style.visibility='hidden';}
	}
</script>

</head>
<body onLoad="javascript:hideReport();" style="font-family:Arial; font-size:13px" bgcolor="#FFFFFF" link="#0033CC" vlink="#990099" alink="#CC3300" leftmargin="3" topmargin="0" marginwidth="4" marginheight="0">

<!-- ================= A simple navigation bar to browse around the web application ===================== -->
<table width="100%"  style="font-family:Tahoma, Arial; font-size:14px; font-weight:bold" border="0" cellpadding="2" cellspacing="0" bordercolor="#CC9999" bgcolor="#FFFFFF">
  <tr>
    <td valign='middle' bgcolor="#FFFFFF"><img src="../../rest_html/images/logo_distinctienz.gif" width="133" height="45">  </td>
  </tr>
  <tr>
    <td bgcolor="#000099"> <font color="#FFFFFF"><a href="../../rest_html/home.htm"><font color="#FFFFFF" size="2">
            Home</font> </a><font size="2"> | <a href="DistinctiEnz.pl"><font color="#FFFFFF">New
        RFLP Analysis</font> </a> | <a href="../../rest_html/help_DistinctiEnz.htm"><font color="#FFFFFF">Help</font></a> </font> </font></td>
  </tr>
</table>
 <!--  ==================================== End navigation bar ===============================  -->

MULTILNBL


#-------------------------------- Input Data Validations ------------------------------------------------------
# An error message will be thrown when no file is uploaded.
if (length$uploadFileName<1){
	print "<p align='center'><BR><BR><b>Please submit a valid sequence. <BR><BR><a href='DistinctiEnz.pl'>Back</a></b></p>";	
	exit;
}

# An error message will be thrown if no enzyme is selected.
if (scalar@enzymeNames==0){
	print "<p align='center'><BR><BR><b>Please select an enzyme. <BR><BR><a href='DistinctiEnz.pl'>Back</a></b></p>";
	exit;

}


#Create an array container for all sequences uploaded in the file
@sequences=();
$uploadDir="/home/member/docreza/public_html/cgi-bin/restriction/usruploads";
	# Save the uploaded file to disk with a UNIQUE name
	# This will prevent concurrent users' data interference.
	$fileTempName = time ."_". int(rand(1000)) ."d.fasta";
	
	open UPLOADFILE, ">$uploadDir/$fileTempName";
	binmode UPLOADFILE;
	while(<$uploadFileHandle>){
		print UPLOADFILE;
	}
	close UPLOADFILE;

	# ===============


# Using Bio::SeqIO module of Bioperl make a sequence directly out of the filehandle (above).
# No files are saved to the hard disk of the server.
# This Bioperl module gets all sequences from a file in the given format (FASTA) and stores them in hash.
use Bio::SeqIO;
$SeqFile=Bio::SeqIO->new(-file=>"<$uploadDir/$fileTempName",-format=>"Fasta");

#This counter will keep the No of seq and is used in creating a sq_No format prefix for each sequence
my $seqCounter=0; 


#Loop through the set of sequences received from the filehandle
# First a prefix (e.g. sq_1_..., sq_2_..., etc) will be added to the description of each sequence to have a good order 
# of sequences as appearing in the uploaded file and creating an alphabetic order for future use.
while($mySeqFromFile=$SeqFile->next_seq){
		$seqCounter += 1;
  		$mySeq=$mySeqFromFile;
		
  		# Add prefix sq_No	
		my $seqContent=$mySeq->seq;	
 		my $seqName="sq_". $seqCounter ."_". ($mySeq->primary_id);
 		
 		# If user has NOT turned on PowerCleaner option, check if each sequence contains non-atcg characters. 
		# However, first, space, tab, and newline characters are filtered.
		# If sequence contains non-atcg characters, an error message is thrown and the user is recommended 
		# to use PowerCleaner.
		if ((length$seqContent > 0) && ($powerCleaner ne 'doClean')){
			$seqContent =uc($seqContent);
			
				eval{$seqContent =~ s/[ \t15\n\r]//g ;};
				if($@){print "An error occured in processing sequence content.<BR>Try again with checking PowerCleaner in input.<BR><a href='t_DistinctiEnz.pl'>Back</a>" ;}

			if($seqContent =~ m/[^ATCGatcg]/){
				print "<p align='center'><BR><BR><b>One of the sequences in your uploaded file may contain non-atcg characters.<BR>You can turn on <i>PowerCleaner</i> and try again.<BR><BR><a href='DistinctiEnz.pl'>Back</a></p>";	
				exit;	
			}
		}
		
		# If user has turned on PowerCleaner option, remove all non-atcg characters from each sequence. 
 		if ((length$seqContent > 0) && ($powerCleaner eq 'doClean')){
		$seqContent =uc($seqContent);
		$seqContent =~ s/[^ATCGatcg]//g ;
		}
		
		# Perform a check after PowerCleaning because all sequence characters may be non-atcg and so
		# an empty sequence may have been produced after replacement.
		# Alternatively, one of the sequence may have been empty or lacking a description from the beginning.
		if (((length$seqContent)<1) || ((length$seqName)<1)){
		print "<p align='center'><BR><BR><b>Please submit valid sequences. A sequence in your uploaded file may be empty, containing only non-atcg characters, or lacking a description. <BR><BR><a href='DistinctiEnz.pl'>Back</a></p>";	
		exit;
		}
		# ------------------------------------- End Validations -----------------------------------------------------------
		
		# Now use EACH cleaned and validated sequence description and content stored in a $seqName and $seqContent
		# variables to make a new sequence object with Bio::PrimarySeq module of Bioperl.
		# The above operations makes it possible to have a robust sequence validation for each of the sequences uploaded
		# in the file.
		use Bio::PrimarySeq;
		my $seqTemp=Bio::PrimarySeq->new(-seq => $seqContent,
				-primary_id=> $seqName,
				-molecule=> 'dna'
				);
		#If sequence is circular define it.
		if ($seqType eq 'isCircular'){
			$seqTemp->is_circular(1);
			}

		# Add the valid sequence to sequence container
		push(@sequences, $seqTemp);
			
#End while
}

#Make a Hash to store seq name and seq size 
%seq_sizes=();
foreach $seq(@sequences){
	my $seq_name=$seq->primary_id;
	$seq_sizes{$seq->primary_id} = $seq->length();
}

#Delete all Fasta files created before to save disk
my @fileList=<$uploadDir/*.fasta>;

#Check each file name and take out the time label part of file name (e.g see above)
foreach $fileName(@fileList){
	 my $tmpFN=$fileName;
	 $tmpFN =~ s/_.*//g;
 	 $tmpFN =~ s/\/.*\///g;
	#Delete file if time label says that it was created 5 seconds before.	 	
	 if ($tmpFN < (time - 5)){
 		 unlink $fileName;
 	 }
}


# Define an enzyme collection using Bio::Restriction::EnzymeCollection module of Bioperl. 
# This collection usually contains around 500 type II enzymes.
# Instead, you can define your enzyme collection using Bio::Restriction::IO module of Bioperl and 
# an enzyme file from REBASE at http://rebase.neb.com/rebase
use Bio::Restriction::EnzymeCollection;
use Bio::Restriction::Enzyme;
my $enzCollection=Bio::Restriction::EnzymeCollection->new();

# Make an array of enzyme OBJECTS using the array containing enzyme names (@enzymeNames) received from input form and
# the enzyme collection defined above. In fact, for each enzyme name, the corresponding full enzyme OBJECT is 
# retrieved from the enzyme collection defined above.
# Simultaneously print out the list of enzymes selected by user
print "<H5>RFLP resolution results for the effect of the following enzyme(s):</H5>";
my @selectedEnzymes=();

my $enz_list_to_show='';
my $j=1;
foreach $enz(@enzymeNames){
 push(@selectedEnzymes, ($enzCollection->get_enzyme($enz)));
 $enz_list_to_show .= "$j) $enz, ";
 $j++;
}
if(scalar@enzymeNames > 1){$enz_list_to_show = substr($enz_list_to_show,0,length($enz_list_to_show)-2);}
print "$enz_list_to_show</br><br>";

# Also print out the list of sequences in order of appearing in the file (alphabetical ordered by prefixes sq_1,...)

print "<H5>on the following sequence(s) in order of appearing in the uploaded file:</H5>";
#Print the name of sequences
print "<OL>";
foreach $sequence(@sequences){
	print "<LI>". $sequence->primary_id ." : <i>". length($sequence->seq)." bp </i>" ;
}
print "</OL><BR><BR>";


# Define an custom enzyme collection. This is performed by:

# 1) Creating an EMPTY enzyme collection and 
my $selEnzCollection=Bio::Restriction::EnzymeCollection->new(-empty=>1);
# 2) Filling it with the enzyme objects array (@selectedEnzymes) created above.
$selEnzCollection->enzymes(@selectedEnzymes);




#======================================= Subroutines used in this script ================================================


# -------------------- MakeUnique Sub--------------------------------------------------------
#This gets a reference to an array of strings or numbers and returns an array with unique (nonrepetitious) set of members.
#None of your array members must be "D~E~L"; this is a reserved string. You can change it with a very rare string!!
sub MakeUnique{
	my ($anyArray)=@_;
	#dereference the input
	my @myArray=@$anyArray;
	#Find repetition members and replace them with "D~E~L" a rare string.
	for ($i=0; $i<scalar@myArray; $i++){
		for($j=0; $j<scalar@myArray; $j++){
			if (($i!=$j) && ($myArray[$i] eq $myArray[$j])){$myArray[$j]="D~E~L";}
		}
	}
	#Remove "D~E~L" occurances;
	my @uniqeArray=();
	foreach my $m(@myArray){
		if ($m ne "D~E~L"){push(@uniqeArray, $m);}
	}
	return @uniqeArray; 
}

# --------------------------------------------------------------------------------------------


# -------------------- RFLPAnalyzer Sub--------------------------------------------------------
# Gets a "reference" to an array of sequences (Bio::Seq) and a collection of Enzymes (Bio::Restriction::EnzymeCollection)
# as input and returns a Hash including Sequence descriptions (primry_id) as keys and a reference to the array of the 
# set of UNIQUE fragment sizes (e.g: (12, 15, 129)) from each sequence as values.

sub RFLPAnalyzer{
	my($seqArray,$customEnzymeCollection)=@_;
	
	#%AnalysisResult will contain a hash with "sequence_primary_id" => @fragments from its analysis
	my %AnalysisResult=();
	
	foreach $seq(@$seqArray){
		
		# Using Bio::Restriction::Analysis of Bioperl, make an analyzer object to perform cuts
		# using the sequence and the custom enzyme collection ($selEnzCollection) created above.	
		use Bio::Restriction::Analysis;
		my $analyzer=Bio::Restriction::Analysis->new(-seq => $seq);
		
		#Cut the sequence in multiple mode
		$analyzer->cut('multiple',$customEnzymeCollection);
		
		# Get the product of the multiple-digest analysis and store the fragment sizes.
		#I used fragments array and then calculated the sizes from length of each fragment
		# I could have used #my @sizes=$analyzer->sizes('multiple_digest'); BUT I encountered
		#some bugs when the number of enzymes goes up and very short fragments and overlapping cuts occur.
		my @sizes=();
		my @frags=$analyzer->fragments('multiple_digest');
		foreach my $frag(@frags){
			push(@sizes, length($frag));
		}
		
		
		#Use MakeUnique subroutine to remove all but one of similar-length fragments and produce
		# a set of UNIQUE fragment lengths for each sequence (we put \@resultFragments to reference @resultFragments)
		my @uniqueSizes=&MakeUnique(\@sizes) ;
		
		#Sort the unique-fragment sizes
		my @uniqueSizesSorted=sort{$a<=>$b}@uniqueSizes;
		
		#Pack the @uniqueSizesSorted array in a referenced variable and store it in the Hash
		my $sizesRef = \@uniqueSizesSorted;
		$AnalysisResult{$seq->primary_id}=$sizesRef ;
	}
	#Deliver the output of this subroutine
	#The results will be a Hash like this: (sq_1_Sequence_a => (12, 15, 129), sq_2_Sequence_b => (34, 50, 50, 279),...)
	
	
	return %AnalysisResult;
}
# ------------------------------------------------------------------------------------------


# ------------------------------------ CalculateResolution Sub--------------------------------------------------------
#This gets a Hash of RFLP analysis  results from from RFLPAnalyzer subroutine as input and returns the Resolution Profile
# to make distinction among a set of analyzed sequences by RFLP.
sub CalculateResolution{
	my ($HashRFLPAnalysis)=@_;
		
	#Dereference the input
	my %myRFLPAnalysis=%$HashRFLPAnalysis;
	#Ignore n, this was just a marker!
	my $n=0;
	my $resolutionResult=0;
#+++++++++++++++++++++++++++++++++++++++++	Preliminary calculation of Resolution power +++++++++++++++++++++++++++++
	# Keep a count of all "size arrays" in this array showing the number of input "size arrays".
	# Each contains the set of sizes for fragments from a sequence in an array.
	# First see if the size arrays produced are of varying length. In other words, check if 
	# each sequence has produced a unique NUMBER of unique fragment sizes. If so, resolution is 100 percent
	
	#Retrieve the size-arrays from the input hash and store them in another array
	# the array will include the data like: ((12,15,129), (34,50, 50, 279), (500, 580, 1200)).
	my @hashValues=values%myRFLPAnalysis;
	
	# Get the number of sequences in the digest, this will make the denominator of the Resolution Power ratio
	my $RFLPSetOfFragSizes=scalar@hashValues;			
	
	# Store the fragment size data of fragments (including nonunique ones) from each sequence digest in an array
	# the array will include the data like: (3, 4, 3) corresponding to the example array : ((12,15,129), (34,50, 50, 279), (500, 580, 1200))
	my @FragCount=();
	foreach $v(@hashValues){
		push(@FragCount,scalar@$v);
	}
	
	# Make a preliminary calculation to see if each sequence has produced a unique NUMBER of unique fragment sizes. 
	# If so, resolution is 100 percent. MakeUnique subroutine is used to remove similar fragment sizes 
	# from each size-arra. For example it will get (3, 4, 3) and return (3, 4).
	# This preliminary calculation will make the restof the Sub unnecessary if a 100% resolution is reached and
	# so will save time and program resources.
	my $resolutionFirstTest=scalar(&MakeUnique(\@FragCount))/scalar(@FragCount);
	if ($resolutionFirstTest==1){
		$resolutionResult="100%";
		$RFLPUniqueSetOfFragSizes=$RFLPSetOfFragSizes;
 	}
#+++++++++++++++++ Secondary calculation of Resolution power when a 100% resolution was not reached above ++++++++++

	# If a 100% resolution is NOT reached , we need to compare each 2 size-arrays which contain equal number of 
	# fragments to see if the sizes are different in thaose 2 size-arrays .
	else{
		
		#Use a dual loop to make a comparison of each size-array with the other size-arrays in $hashValues
		#This is performed by first merging each 2 nonself size arrays and then using MakeUnique sub to get a 
		#unique-sized version of each merged set. Comparing these 2 sets (by getting the length of each merged set
		# before and after MakeUnique process) will test whether there are at least one different fragment size 
		# in 2 size-arrays before merging. If all fragment sizes are similar then a 0.5 ratio will ensue. This shows that
		# not only the NUMBER of fragments produced from the 2 compared sequences were equal but ALL the fragment SIZES
		# were equal too making the 2 sequences producing completely SIMILAR and INDISTINGUISHABLE bands. Such a score
		# will DELETE one of the totally similar size arrays from the collection ($hashValues). And this will continue
		# untill one of the 2 totally similar size-arrays have been deleted from the collection.  
		# For example merging 2 size-arrays: (12,15,129) and (34,50, 50, 279) will yield:
		# (12, 15, 129, 34, 50, 50, 279) having a 7 members; after MakeUnique this will be:
		# (12, 15, 129, 34, 50, 279) with 6 members. Since there are many sizes different in 2 sets so will not get
		# a ratio=6/7 and not 0.5. 

		#Start the dual loop	
		for (my $i=0; $i<scalar@hashValues; $i++){
		for(my $j=0; $j<scalar@hashValues; $j++){
			#Ignore self-self comparison of arrays
			if($i!=$j){ 
				my $ref1=$hashValues[$i];
				my $ref2=$hashValues[$j];
				#Seperately get and dereference each size-array
				my @arr1=@$ref1;
				my @arr2=@$ref2;
				#Merge 2 size-arrays
				foreach $it(@arr1){
					push(@arr2,$it);
				}
				@merge2Arrays=@arr2;
				
				# Compare these 2 sets (by getting the length of each merged set before and after MakeUnique process)
				# and test whether there are at least one different fragment size in 2 size-arrays 
				# before merging. If all fragment sizes are similar then a 0.5 ratio will ensue. 
				my $compare4UniqueSets=scalar(&MakeUnique(\@merge2Arrays))/scalar(@merge2Arrays);
				if($compare4UniqueSets==0.5){
					#If all of the fragments in 2 size-arrays are equal, then delete it.
					splice(@hashValues, $j,1);
					#Ignore n, this was a test marker only.
					$n++;
				}
			}
		}
		}
		#After the totally similar size-arrays have been deleted count the array members which are now unique
		$RFLPUniqueSetOfFragSizes=scalar@hashValues;
		
		#if all sequences produce only one unique-fragment-size pattern then the resolution power is 0
		if($RFLPUniqueSetOfFragSizes == 1){
			$resolutionResult=0;
		}
		else{
			#Compare unique sets and the input set (the total number of sequences)
			$resolutionResult=$RFLPUniqueSetOfFragSizes/$RFLPSetOfFragSizes;
		}
		#Format the Resolution power in % with 3 digits
		$resolutionResult=(100 * sprintf("%.3f",$resolutionResult))."%";
	}
return "<font color='blue'><UL><li>Resolution: $resolutionResult</li><li>Unique-fragment RFLP sets: $RFLPUniqueSetOfFragSizes </li><li>Total No of sequences: $RFLPSetOfFragSizes</UL></font>";
}
# ------------------------------------------------------------------------------------------


#========================================== End of Subroutines =============================================================



#------------- use the subroutines to perform main analysis -----------------------------------------------------------

#Reference the sequences container (array) to be used as the input to RFLPAnalyzer subroutine.
$sequencesArray=\@sequences;
if ($doBatchAnalysis eq "batch"){

	print "<b>Resolution power for each enzyme:<div id='reporting'><br><font color='red'><b> Calculating resolution power for each enzyme....</b></font><br><br><br></div><br></b>";	
	print "<OL>";
	foreach my $enz(@selectedEnzymes){
		
		# Define a custom enzyme collection. This is performed by:

		# 1) Creating an EMPTY enzyme collection and 
		my $selEnzCol=Bio::Restriction::EnzymeCollection->new(-empty=>1);
		# 2) Filling it with the enzyme objects array (@selectedEnzymes) created above.
		my @enz=($enz);
		$selEnzCol->enzymes(@enz);

		my %RFLPResults_1=&RFLPAnalyzer($sequencesArray,$selEnzCol);
		my $eachRflp=\%RFLPResults_1;
		my $eachResolution=&CalculateResolution($eachRflp);
		print "<li><b>" . $enz->name .": $eachResolution</b></li>";
		#push(@allResolutions, $enz->name .": $eachResolution");
	}
	print "</OL>";
	#print "@allResolutions";
}
else{
	print "<div id='reporting'><br><font color='red'><b> Calculating resolution power ....<br></b></font></div>";	
		
	#Perform RFLP and save results in a Hash
	%RFLPResults=&RFLPAnalyzer($sequencesArray,$selEnzCollection);
	

	#Print out the resolution profile
	print "<TABLE width='400px' border='1'><tr bgcolor='pink'><th>";
	print "<center>Resolution Profile <a target='_blank' href='../../rest_html/help_DistinctiEnz.htm#whatis'>[?]</a></th></tr><td>";
	print &CalculateResolution(\%RFLPResults);
	print "</b></td></tr></Table><hr>";
	
	#Ouput a list of sequences and the UNIQUE set of fragment sizes each produce in RFLP
	@seqKeys=keys%RFLPResults;
	
	#Sort them by No of sq_No_SeqName before outprinting.
	@seqKeysSorted=sort{(substr($a,3,(index($a,'_',3) - 3)))  <=> (substr($b,3,(index($b,'_',3)) - 3))}(@seqKeys);
	print "<H5>Sequences and unique set of fragments in order of appearing in the uploaded file:</H5><BR>";
	print "<Table  style='font-size:11px' cellspacing='1'  border='1'><tr><th valign='top' bgcolor='pink'>Sequence</th><th valign='top' bgcolor='pink'>Unique fragments(No.)</th><th valign='top' bgcolor='pink'>Fragments by size</th><th valign='top' bgcolor='pink'>Alerts</th></tr>";
	my @v=@$value;
	foreach $key(@seqKeysSorted){
		$value=$RFLPResults{$key};
		print "<tr><td valign='top'><b>$key </b></td><td valign='top' ><center>". (scalar@$value) ."</center></td><td valign='top'>";
		print join " - ", @$value ;
		print "</td><td><center>";
		#Determine and warn user if no cut occured for this enzyme. (size of either cut would be equal to seq length.		
		if((scalar@$value == 1) && ($v[0] = $seq_sizes{$key})){print " No enzyme cut this sequence.";}else{print "--";}
		print "</center></td></tr>";

	}
	print "</Table><br><hr>";

	#================================================= RFLP virtual Gel Graphic Output =====================================================
	
	# The basis of this graphical ouput is Bio::Tools::Gel module of Bioperl; Copyright Allen Day <allenday@ucla.edu>.
	# However, I did not directly use this module. Instead, I used the formula used there for calculation of 
	#  migration distances using:   
	#            distance = dilation * (4 - log10(length(dna))
	# Source: Molecular Cloning, a Laboratory Manual. Sambrook, Fritsch, Maniatis. 
	# CSHL Press, 1989.
	# As documented in Bio::Tools::Gel, "the formula calculates migration distances based solely on
	# the length of the nucleotide sequence.  Secondary or tertiary structure, 
	# curvature, and other biophysical attributes of a sequence are currently 
	# not considered." 
	
	# Since the formula above generates negative values when log10(length(dna) is more than 4, I also made a correction
	# and for sequences meeting this condition no migration will occur. 
	
	# ------------------------------------------------------------------------------------------------
	
	
	#Make available GD Graphical Design module of Perl
	use GD;
	
	#Set some default properties of the image 
	my $dilation=$dilationFactor;
	my $startPoint=20; #Sets from which horizontal distance to the left border of the picture RFLP bands Start to appear.
	my $ladderWidth=100;#Sets width of ladder
	my $lineLengh=50; #Sets the width of each RFLP band on Gel. 
	my $spaceBetweenLines=20; #Sets the space between band columns on RFLP Gel. 
	my $upperSpace=27; #Sets the beginning of DNA band movement off the upper border of image.
	my $picWidth=(2 * $startPoint) + $ladderWidth + (scalar@seqKeysSorted * ($lineLengh + $spaceBetweenLines)); #Dynamically sets the width of picture
	my $picHeight=($upperSpace + ($dilation * 4) + 50); #Dynamically sets the height of picture
	
	
	#The image output files are named by time they are created with a suffix of a random number so that overwrites for simultaneous
	# users will not occur. The image files will be created in usr_images directory which MUST have WRITE permission.
	
	# An example of image file names produced below: 1118388978_127.png.
	# The part before underscore : 1118388978 is the time lable and the part after is the random number produced by rand(1000).
	# produced by time() function of Perl which "returns the number of non-leap seconds since 
	# whatever time the system considers to be the 
	# epoch (that's 00:00:00, January 1, 1904 for MacOS, and 00:00:00 UTC, January 1, 1970 for most other systems)".
	 
	$rootDir="/home/member/docreza/public_html/rest_html/";
	# To save web space of the application, all image files created 300 seconds (5 minutes before) by any previous users
	# are deleted. Five minutes are allowed so that if a user is connected through a slow connection can still load the
	# image in his browser.
	
	#Get a list of png image files in usr_images directory and store them in an array.
	my @fileList=<$rootDir/usr_images/*.png>;
	
	#Check each file name and take out the time label part of file name (e.g see above)
	foreach $fileName(@fileList){
		 my $tmpFN=$fileName;
		 $tmpFN =~ s/\/.*\///g;
		 $tmpFN =~ s/_.*//g;
		#Delete file if time label says that it was created 5 minutes (300 secs) before.	 	
		 if ($tmpFN < (time - 300)){
	 		 unlink $fileName;
	 	 }
	}
	
	#Start creating a random image file name with time(0 and rand() functions (see above)
	$nowString = time ."_". int(rand(1000));
	$imgFileName="$nowString.png";
	
	#Make a file handle  (stream for image file).
	open PICFILE, ">$rootDir/usr_images/$imgFileName" or die "Couldn't open image file: $!<BR>";
	
	#Define image size, some colors and the preliminary defaults of this image.
	my $im = new GD::Image($picWidth, $picHeight);
	
	my $white = $im->colorAllocate(255,255,255);
	my $black = $im->colorAllocate(0, 0, 0);
	my $red = $im->colorAllocate(255, 0, 0);
	my $blue = $im->colorAllocate(0, 0, 255);
	my $green = $im->colorAllocate(50, 200, 0);
	
	$im->transparent($white);
	$im->interlaced('true');
	
	#Draw a blue border
	$im->rectangle(0, 0, $picWidth - 2, $picHeight - 2, $blue); 
	
	#Draw two stripes, one at the top delineating Minus area and one at the bottom defining Plus area of RFLP Gel.
	
	#Upper stripe (lower border)
	$im->line(0, ($upperSpace -1) , ($picWidth - 2), ($upperSpace - 1) , $blue);
	#Minus string
	$im->string(gdTinyFont	,2 , ($upperSpace - 9), 'Minus (-)', $black);
	#Upper stripe (upper border)
	$im->line(0, ($upperSpace - 10) , ($picWidth - 2), ($upperSpace - 10), $blue);
	
	#Lower stripe (upper border)
	$im->line(0, ($picHeight - 10) , ($picWidth - 1), ($picHeight - 10), $blue);
	#Plus string
	$im->string(gdTinyFont, 2,($picHeight - 9), 'Plus (+)', $black);

	
	
	#----------------
	@ladder=split(/-/,$ladderSet);
	for($i=0; $i < scalar@ladder ; $i++){
		my $ladDist=0;
		#Calculate Log 10 of $i
		my $distLogVal=log($ladder[$i])/log(10);
		$ladDist= $upperSpace + ($dilation * (4 - $distLogVal));
		$im->line(($startPoint + 40) , $ladDist , $startPoint + 60 , $ladDist, $black);
		#For odd 
		$startLaddString=0;
		if(($i%2)== 0){$startLaddString=$startPoint + 15;} else{$startLaddString=$startPoint + 65};
		$im->string(gdTinyFont	, $startLaddString, $ladDist - 3, $ladder[$i], $black);

	}
	$startPoint += $ladderWidth;
	#-------------	
	#Now start drawing the bands for each fragment set of each sequence. See the explanations on the formula used 
	# for migration distances of dna fragments above.   distance = dilation * (4 - log10(length(dna));
	# Use @seqKeysSorted created in analysis above.
	foreach $key(@seqKeysSorted){
		$value=$RFLPResults{$key};
		foreach my $size(@$value){
			my $distance=0;
			
			#Calculate Log 10 of $size
			my $logVal=log($size)/log(10);
			$distance = $upperSpace + ($dilation * (4 - $logVal));
			
			# When log10(length(dna) is more than 4, migration distance wil be negative, I set it 
			# so that no migration occurs.
	 		if ($distance < $upperSpace){
				$distance = $upperSpace ;
			}
			#Here, draw the band for this fragment.
			$im->line($startPoint , $distance , $startPoint + $lineLengh, $distance, $red);
		}
		#Also write the name of each sequence at the top of its band column
		$im->string(gdTinyFont	, $startPoint , 2, substr($key,0,12), $black);
		
		#Make the start point go to right to be ready for drawing bands of the next column.
		$startPoint += ($lineLengh + $spaceBetweenLines); 
	}
	
	#Output the image in png format.  
	binmode PICFILE;
	print PICFILE $im->png;
	close PICFILE;      
	
	#Outprint the HTML to show the image in browser.
	print "<p><H3><center>See RFLP on GEL<center><a target='_blank' href='../../rest_html/help_DistinctiEnz.htm#output'>[?]</a></H3>
	
		<br>
	<div style='position:relative; top:auto' >
   <img src='../../rest_html/usr_images/$imgFileName' style='cursor:crosshair' id='pic'/>
    <hr size='1'  id='marker_h' style='position:absolute; top:0px; left:0px; visibility:hidden' color='#990000' />
	<input name='y_offset' type='text' id='y_offset'  style='position:absolute; visibility:hidden; width:50px; font-size:9px; background-color:#FFCC66'/>
	</div>";    
} #End of else 
#Print out a footer for navigation.
print "<p align='center'><a href='#'>^Top</a> | <a href='DistinctiEnz.pl'> << New RFLP</a></p><br></body></html>";



