#!/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 t_DistinctiEnz.pl; the data includes:
# 	1) An uploaded file containing a set of FASTA sequences;
# 	2) A list of enzyme names; and
# 	3) Whether a multiple digest with graphical output must be performed or 'test each enzyme mode' is selected. 
# 	4) Some option data for T-RFLP analysis, sequence filtering of non-atcg characters, and graphical output adjustment.
# The script will perform a T-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 2 tagged fragments and their sizes from T-RFLP of each sequence.
# 2) Each fragment size actually will appear as peak (vertical line) in T-RFLP sequencer plot.
# 3) 3'- and 5'-tagged peaks will be of distinct colors in sequencer plot
# 4) If the fragment pairs of each sequence are unique (interms of size of similarly-tagged fragments), resolution
# 	 wiil be 100%. (example: F20-T30 from seq1 and F45-T240 from seq2 are unique
# 						but F20-T30 from seq1 and F20-T30 from seq2 is NOT unique
# 						even F20-T30 from seq1 and F20-T40 from seq2 is unique
# 						also F30-T20 from seq1 and F20-T30 from seq2 is unique
# 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.
# the $uploadDir="/home/member/docreza/public_html/cgi-bin/restriction/usruploads"; must be set as appropriate on your server.
# the $rootDir="/home/member/docreza/public_html/cgi-bin/restriction"; must be set as appropriate on your server.


# ---------------------------------------------------------------------------------------------------------------
# ----------------------------------------- Acknowledgements -------------------------------------------------------------
# Bioperl modules from http://www.bioperl.org was used to:
# 1) make a reference list of restriction enzymes;
# 2) perform restriction enzyme analyses; and
# 3) make a sequence object model from file or directly from strings.
# I also wish to thank valuable comments from my classmates in Biocomputing Module
# in University of Manchester. 
# My special thanks goes to Caroline Reiff my classmate in Biocomputing Module who was the first one
# to suggest the idea of making such an appication.
# For a multiple-digest, I used some code examples taught by Robert Edwards [redwards@pub.open-bio.org] in 
# [ http://bioperl.org/pipermail/bioperl-guts-l/2003-December/015775.html ] in addition to documentation 
# for Bio::Restriction::Analysis.
# --------------------------------------------------------------------------------------------------
# ------------------------------------------------ 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 type of T-RFLP, and PowerCleaner options
$analysis_type=$query->param("cboAction");
$powerCleaner=$query->param("ckCleaner");
# $seqType=$query->param("ckCircular");
#Define whether a multiple digest with graphical output must be performed or 'test each enzyme mode' is selected. 
$doBatchAnalysis=$query->param("ckBatch");
if($doBatchAnalysis || (scalar@enzymeNames == 1)){
	$digest_mode='single';	
}
# use print<< command for out-printing multiline strings in Perl. We produce HTML headers and a navigation bar
print <<MULTILNBL;
<html><head><title>T-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_v');
			var label=document.getElementById('x_offset');
			document.getElementById('x_offset').value = (offset_x - 50) + ' bp';
			if(offset_x > 50){ 
				marker.style.visibility='visible';			
				label.style.visibility='visible';						
			}
			else{
				marker.style.visibility='hidden';			
				label.style.visibility='hidden';				
			}
			document.getElementById('marker_v').style.height=pic.height;
			document.getElementById('marker_v').style.width=1;
			document.getElementById('marker_v').style.left= offset_x;
			document.getElementById('x_offset').style.top= offset_y;
			document.getElementById('x_offset').style.left= offset_x;
		}
		return true; 
	}	
	function hideReport(){
		if(document.all.reporting){	document.all.reporting.style.visibility='hidden';}
	}
	function hideMarkers(){
		var marker=document.getElementById('marker_v');
		var label=document.getElementById('x_offset');
		marker.style.visibility='hidden';			
		label.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_t_distinctienz.gif" width="133" height="45">  </td>
  </tr>
  <tr>
    <td bgcolor="#CCCCFF"> <font color="#000099"><a href="home.htm"><font color="#000099" size="2">
            Home</font> </a><font size="2"> | <a href="t_DistinctiEnz.pl"><font color="#000099">New
        T-RFLP Analysis</font> </a> | <a href="help_t_DistinctiEnz.htm"><font color="#000099">Help</font></a> </font> </font></td>
  </tr>
</table>
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='t_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='t_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";
eval { 
	# 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;
	# ===============
};
if($@){
	print "An error occured in saving uploaded file.<BR><a href='t_DistinctiEnz.pl'>Back</a>" ;
}
# 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='t_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='t_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 abundance in [0-9{*}] pattern
%seq_abundances=();
foreach $seq(@sequences){
	my $seq_name=$seq->primary_id;
	# A default 
	$abundance=100;
	#Look for [\d+] pattern
	if ($seq_name =~ /\|(\d*)\|/ ){
		$abundance = $1;
	}
	$seq_abundances{$seq->primary_id} = $abundance;
}
#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();
}
#-----------------------------------------------------------------
# @abkeys=keys%seq_abundances;
# foreach $k(@abkeys){
# print "<br>-$k :" . $seq_abundances{$k}	;
# }
#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>". ($sequence->length())." bp </i>" ;
}
print "</OL><BR><BR>";
# Define a 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, $t_rflp_type, $digest_mode)=@_;
	#%AnalysisResult will contain a hash with "sequence_primary_id" => @fragments from its analysis
	my %AnalysisResult=();
	# 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;	
	foreach $seq(@$seqArray){
		$analyzer=Bio::Restriction::Analysis->new(-seq => $seq);
		# 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=();
		if($digest_mode ne 'single'){		
			#Cut the sequence in multiple mode
			$analyzer->cut('multiple',$customEnzymeCollection);
		}
		if($digest_mode eq 'single'){
			my $enz=$enzymeNames[0];
			@frags=$analyzer->fragments($enz);
		}
		else{
			@frags=$analyzer->fragments('multiple_digest');
		}
		if($t_rflp_type eq '5-prime'){
			push(@sizes, length($frags[0]));
		}
		if($t_rflp_type eq '3-prime'){
			push(@sizes, length($frags[-1]));
		}
		if($t_rflp_type eq 'dual'){
			push(@sizes, 'F' . length($frags[0]));
			push(@sizes, 'T' . length($frags[-1]));
		}
		# For 5' or 3' t-RFLP, @sizes only contains one value the length of 3'-tagged or 5'-tagged fragment
		#For dual t-RFLP @sizes contains 2 values sizes of 3'-tagged and 5'-tagged fragment.
		#Pack the @sizes array in a referenced variable and store it in the Hash
		my $sizesRef = \@sizes;
		$AnalysisResult{$seq->primary_id}=$sizesRef ;
	}
	#Deliver the output of this subroutine
	#The results will be a Hash like this: (sq_1_Sequence_a => (T12, F15), sq_2_Sequence_b => (T34, F50),...)
	# or a Hash like this: (sq_1_Sequence_a => (12), sq_2_Sequence_b => (34),...)
	return %AnalysisResult;
}
# ------------------------------------------------------------------------------------------
# ------------------------------------ CalculateResolution Sub--------------------------------------------------------
#This gets a Hash of T-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,$t_rflp_type)=@_;
	#Dereference the input
	my %myRFLPAnalysis=%$HashRFLPAnalysis;
	my $resolutionResult=0;
#+++++++++++++++++++++++++++++++++++++++++	Preliminary calculation of Resolution power +++++++++++++++++++++++++++++
	#Retrieve the size-arrays from the input hash and store them in another array
	# the array will include the data like: ((12), (34), (500)).
	#or in dual t-RFLP ((F12, T50), (F34,T45), (F50,T50))
	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;			
	my $RFLPUniqueSetOfFragSizes=0;
	#Get frag sizes from each one-member array in non-dual t-rflp
	if (($t_rflp_type eq '5-prime') || ($t_rflp_type eq '3-prime')){
		my @sizes=();
		foreach $v(@hashValues){
			my @arr_v=@$v;
			#get the first array member (the only array member)
			push(@sizes, $arr_v[0] );
		}
		#@sizes=(8,9,9);
		#Make unique; this will be compared with total in the end of this sub
		$RFLPUniqueSetOfFragSizes=scalar(&MakeUnique(\@sizes));
 	}
	# dual t-RFLP
	if ($t_rflp_type eq 'dual'){
		#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 7 members; after MakeUnique this will be:
		# (12, 15, 129, 34, 50, 279) with 6 members. Since there are one size difference 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(join(',','@arr1') eq join(',','@arr2')){	
					#If both of F and T fragments in 2 size-arrays are equal, then delete it.
					splice(@hashValues, $j,1);
				}
			}
		}
		}
		#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 t-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></div></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,$analysis_type);
		my $eachRflp=\%RFLPResults_1;
		my $eachResolution=&CalculateResolution($eachRflp,$analysis_type);
		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>";	
# 	$t0=time();
# 	print $t0 . '<br>';	
	#Perform RFLP and save results in a Hash
	%RFLPResults=&RFLPAnalyzer($sequencesArray,$selEnzCollection,$analysis_type,$digest_mode);
# 	#$delta=time() - $t0;
# 	#print "After RFLP:" . time() . "  -Delta: $delta	seconds<br>";
# 	$t1=  time();
# 	print $t1 . '<br>';
	print "<b>Type of t-RFLP analysis: $analysis_type </b><br><br>";
	#Print out the resolution profile
	print "<TABLE width='400px' border='1'><tr bgcolor='pink'><th>";
	print "<center>Resolution Profile <a target='_blank' href='help_DistinctiEnz.htm#whatis'>[?]</a></th></tr><td>";
	print &CalculateResolution(\%RFLPResults,$analysis_type);
	print "</b></td></tr></Table><hr>";
# 	$delta=time() - $t1;
# 	print "After Resolution:" . time() . "  -Delta: $delta	seconds<br>";
	#Ouput a list of sequences and the UNIQUE set of fragment sizes each produce in t-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 in order of appearing in the uploaded file and their tagged fragments :</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'>Tagged fragments by size</th><th valign='top' bgcolor='pink'>Alerts</th></tr>";
	#an array to store all sizes for next sorting for getting longest frag for pic width
	@all_sizes=();
	foreach $key(@seqKeysSorted){
		$value=$RFLPResults{$key};
		print "<tr><td valign='top'><b>$key </b></td><td valign='top'>";
		print join " - ", @$value ;
		foreach $sz(@$value){
			$sz_num=$sz;
			$sz_num =~ s/[A-Z]//g;
			push(@all_sizes, $sz_num); 
		}
		@all_sizes_sorted=sort{$a <=> $b}@all_sizes;
		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($sz_num == $seq_sizes{$key}){print " No enzyme cut this sequence.";}else{print "--";}
		print "</center></td></tr>";
	}
	print "</Table><br><hr>";
	#print $all_sizes_sorted[-1];
	#================================================= 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=1;
	my $max_horiz_frag_size=$all_sizes_sorted[-1];
	my $startPoint=50; #Sets from which horizontal distance to the left border of the picture RFLP bands Start to appear.
	my $picWidth= 800;
	if($max_horiz_frag_size > 800){$picWidth=$startPoint + ($dilation * ($max_horiz_frag_size)) + 50;}
	my $picHeight=$query->param("cboPicHeight");
	my $max_vertical_abundance=$query->param("cboMaxVertValue");
	my $vertical_compression=sprintf("%.3f",($picHeight/($max_vertical_abundance * 2)));
	my $baseline=($picHeight/2);#baseline is in the middle of picture vertically and bands will be drawn upwards or dawnwards off that line
	#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);
	my $pink = $im->colorAllocate(230, 150, 230);
	$im->transparent($white);
	$im->interlaced('true');
	#------------------------
	#Draw a blue border
	$im->rectangle(0, 0, $picWidth - 2, $picHeight - 2, $blue); 
	#Middle line 
	$im->line($startPoint, $baseline , ($picWidth - 2), $baseline, $black);	
	#Vertical line 
	$im->line($startPoint,$picHeight, $startPoint, 0, $black);	
	#Arrow for 3-prime	
	$poly3 = new GD::Polygon;
	#Add point (x,y) to the polygon. 
    $poly3->addPt(2,$baseline);
    $poly3->addPt(2,$baseline/2);
    $poly3->addPt(22,$baseline/2 - 20);
    $poly3->addPt(22,$baseline);
    $im->filledPolygon($poly3,$red);
    #Write 3-prime vertically
    $im->stringUp(gdSmallFont,($startPoint - 45), ($baseline -10), '3-prime', $black);
	#Arrow for 5-prime	
	$poly5 = new GD::Polygon;
	#Add point (x,y) to the polygon. 
    $poly5->addPt(2,$baseline + 3);
    $poly5->addPt(2,(3/2) * $baseline);
    $poly5->addPt(22,((3/2) * $baseline) + 20);
    $poly5->addPt(22,$baseline + 3);
    $im->filledPolygon($poly5,$green);
    #Write 3-prime vertically
    $im->stringUp(gdSmallFont,($startPoint - 45), ($baseline + 50), '5-prime', $black);
		#---------------- Draw markers (ladders,...) -------------
	#Draw horizontal markers 
	my $i=0;
	while($i <= $picWidth){
		my $markDist = ($dilation * $i);
		$im->line(($startPoint + $markDist) , $baseline + 1, ($startPoint + $markDist) , $baseline - 1, $black);
		$im->string(gdTinyFont,($startPoint + $markDist), $baseline + 1, $i, $black);
		$i+=25;
	}
	#Draw vertical markers 
	$i=0;
	while($i <= $max_vertical_abundance){
		my $markDist = sprintf("%u",($vertical_compression * $i));
		$im->line(($startPoint - 1) , ($baseline + $markDist), ($startPoint + 1) , ($baseline + $markDist), $black);
		$im->line(($startPoint - 1) , ($baseline - $markDist), ($startPoint + 1) , ($baseline - $markDist), $black);
		$im->string(gdTinyFont,($startPoint - 25), ($baseline + $markDist), $i, $black);
		$im->string(gdTinyFont,($startPoint - 25), ($baseline - $markDist), $i, $black);
		$i += 50;
	}
	#-------------	
	#Now start drawing the bands for each fragment set of each sequence.
	# Use @seqKeysSorted created in analysis above.
	#Create a hash to store each size and its relative cumulative abundance
	my %abundance4sizes=();
	foreach $key(@seqKeysSorted){
		$value=$RFLPResults{$key};
		foreach my $size(@$value){
			$consider_abundances=$query->param("ckConsiderAbundances");		
			if($consider_abundances eq 'yes'){
				#Default height of line is the abundance of each frag size (its corresponding seq. abundance).
				$abundance_height=$seq_abundances{$key};
				# test: print 'abundance_height: ' . $seq_abundances{$key} . '--';
				#If this size is being repeated (there has already been a size for this) add its abundance to previous ones
				$abundance_height= ($seq_abundances{$key}) + ($abundance4sizes{$size});
				#add this size and the new TOTAL abundance (up to now) for this size to hash			
				$abundance4sizes{$size}=$abundance_height;
				#test: print 'abundance4sizes: ' . $abundance4sizes{$size}	. '<br>';				
			}
			else{
				$abundance_height=$max_vertical_abundance/2 ;
			}
			#remove F or T characters
			$size_num = $size;
			$size_num =~ s/[A-Z]//g;		
			my $distance=$size_num * $dilation;
			#Here, draw the band for this fragment.
			#apply vertical compression to abundance_height
			$abundance_height=sprintf("%u", ($abundance_height * $vertical_compression));
			#Depending on analysis type set y coordinate of the vert. line above or below baseline			
			if($analysis_type eq '3-prime'){
				$height_end=$baseline - $abundance_height;
				$line_color=$red;
			}
			if($analysis_type eq '5-prime'){
				$height_end=$baseline + $abundance_height;
				$line_color=$green;
			}
			if($analysis_type eq 'dual'){
				if($size =~ /T\d+/ ){
					$height_end=$baseline - $abundance_height;
					$line_color=$red;
				}
				if($size =~ /F\d+/ ){
					$height_end=$baseline + $abundance_height;
					$line_color=$green;
				}
			}	
			#If abundances cumulate to exceed height of picture in either direction, cut them short.
			#In addition, put a mark at the end of the line to show its abundance was cut.
			if($height_end < 0){
				$abund= sprintf("%u", $height_end);
				$height_end = 0;
				$im->stringUp(gdTinyFont, $startPoint +  $distance - 4, ($baseline - $height_end)/2 , $abund . ' >|', $blue);
			}
			if($height_end > $picHeight){
				$height_end = $picHeight;
				$mark_ceiling='T';
				$im->stringUp(gdTinyFont, $startPoint +  $distance - 4, ($baseline + $height_end)/2 , '|< ' . $height_end , $blue);
			}
			#Now draw the vertical line for this band
			$im->line($startPoint +  $distance, $baseline , $startPoint +  $distance, $height_end, $line_color);
			#$im->stringUp(gdTinyFont, $startPoint +  $distance, 10, $height_end/2, '$height_end >|', $black);
	}
		#Also write the name of each sequence at the top of its band column
		#$im->string(gdTinyFont	, $startPoint , 2, substr($key,0,12), $black);
	}
	#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>t-RFLP graphical output<center><a target='_blank' href='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' onMouseOut='javascript:hideMarkers();'/>
    <img src='../../rest_html/images/vert_marker.gif' id='marker_v' style='position:absolute; top:0px; left:10px; visibility:hidden' color='#990000' />
    <hr size='1'  id='marker_h' style='position:absolute; top:0px; left:0px; visibility:hidden' color='#990000' />
	<input name='x_offset' type='text' id='x_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='t_DistinctiEnz.pl'> << New T-RFLP</a></p><br></body></html>";
