#!/usr/bin/perl -w
#change the line above on Unix OS usually to: #!usr/bin/perl or wherever perl.exe lies
print "Content-type: text/html\n\n";

# ---------------------------- Script Summary ------------------------------------------------------------------
# This perl script receives the data sent by the input form in IncisiEnz.pl; the data includes:
# 	1) a sequence either from text boxes or as an uploaded file;
# 	2) A list of enzyme names; and
# 	3) Some option data for output view options , sequence type and sequence filtering of non-atcgi characters.

# The script will perform a restriction analysis with enzyme(s) selected by user on the submitted sequence and
# prints out the analysis in a summary and 3 views depending on the view options selected by user:
# 	1) Fragments Length View
# 	2) Sequence View
# 	3) Map View
# The views will be explained below.

# -------------------------------- 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.

# For a multiple-digest, I used some 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.

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


# -------------------------------------- Script dependancies ---------------------------------------------------------
# This script uses CGI module 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.

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


# ------------------------------------------------ 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;

# Get the raw data for analysis from the input form in IncisiEnz.pl
# Get sequence data
$seqName=$query->param("txtSeqName");
$seqContent=$query->param("txtSeqContent");
$uploadFileName=$query->param("fileSeq"); #gets the name of file uploaded
$uploadFileHandle = $query->upload("fileSeq"); #gets the file stream itself



# Get enzyme(s) selected
@enzymeNames=$query->param("cboEnzymes");

# Get PowerCleaner and analysis view options 
$powerCleaner=$query->param("ckCleaner");
$seqType=$query->param("ckCircular");

$viewLength=$query->param("ckViewLength");
$viewSeq=$query->param("ckViewSeq");
$viewMap=$query->param("ckViewMap");



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

<html><head>
<title>IncisiEnz Analysis Results</title>
<!-- A javascript function to use for making the map table in Map View no wrap. -->
<script language="javascript">
	function straighten(){
		var tds=document.getElementById("mapTable").cells;
		
		for(var i=0; i < tds.length ; i++){
				tds[i].noWrap="true";
		}
	}

</script>
<!-- End Script -->
</head>
<body 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_incisienz.gif" width="133" height="45">   </td>
  </tr>
  <tr>
    <td bgcolor="#0066FF"> <font color="#FFFFFF"><a href="../../rest_html/home.htm"><font color="#FFFFFF" size="2">Home</font> </a><font size="2"> | <a href="IncisiEnz.pl"><font color="#FFFFFF">New
       Restriction Analysis </font> </a> | <a href="../../rest_html/help_IncisiEnz.htm"><font color="#FFFFFF">Help</font></a> </font> </font></td>
  </tr>
</table>
 <!--  ==================================== End navigation bar ===============================  -->
 

MULTILNBLC


#-------------------------------- Input Data Validations ------------------------------------------------------
# Error messages will be thrown if a sequence does not have a description or content.
if ((length$seqContent<1) && (length$uploadFileName<1)){
	print "<p align='center'><BR><BR><b>Please submit a valid sequence. <BR><BR><a href='IncisiEnz.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='IncisiEnz.pl'>Back</a></b></p>";
	exit;

}


$uploadDir="/home/member/docreza/public_html/cgi-bin/restriction/usruploads";
#If both a sequence from seqContent textarea and a sequence file was uploaded by user 
# prefer the sequence from file and replace the values from text box.
if (length$uploadFileName > 0){
	# Save the uploaded file to disk with a UNIQUE name
	# This will prevent concurrent users' data interference.
	$fileTempName = time ."_". int(rand(1000)) .".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");

	# Get the first sequence (actually we need 1 sequence but even if user uploads a file containing
	# multiple sequences mistakenly, it seems wise to get the first sequence. 
	# If a file containing one sequence is uploaded, things are ok.	
	$mySeq=$SeqFile->next_seq;
	
	# Then a check is performed on the content and description of the sequence and 
	# errors will be thrown if either is empty. 	
	if (((length$mySeq->seq)<1) || ((length$mySeq->primary_id)<1)){
		print "<p align='center'><BR><BR><b>Please submit a valid sequence. The sequence may be empty or lacking a description.<BR><BR><a href='IncisiEnz.pl'>Back</a></b></p>";	
		exit;
	}

	
		
	
	# If anything goes ok, the sequence description (name) and content will be stored in 2 variables.
	else{
		$seqContent=$mySeq->seq;
		$seqName=$mySeq->primary_id;
		}
		

}


		
#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 - 300)){
 		 
unlink $fileName;
 	 }
}

# If user has NOT turned on PowerCleaner option, check if 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);
	$seqContent =~ s/[ \t15\n\r]//g ;
	if($seqContent =~ m/[^ATCGatcg]/){
		print "<p align='center'><BR><BR><b>The sequence may contain non-atcg characters.<BR>You can turn on <i>PowerCleaner</i> and try again.<BR><BR><a href='IncisiEnz.pl'>Back</a></p>";	
		exit;	
	}
}

# If user has turned on PowerCleaner option, remove all non-atcg characters. 
if ((length$seqContent > 0) && ($powerCleaner eq 'doClean')){
	$seqContent =uc($seqContent);
	$seqContent =~ s/[^ATCGatcg]//g ;
	# Perform another check after PowerCleaning because all sequence characters may be non-atcg and so
	# an empty sequence may have been produced after replacement.
	if (length$seqContent < 1){
		print "<p align='center'><BR><BR><b>All sequence characters were non-atcg! Please submit a valid sequence.<BR><BR><a href='IncisiEnz.pl'>Back</a></b></p>";	
		exit;
	}

}

# ------------------------------------- End Validations -----------------------------------------------------------

# Now use the cleaned and validated sequence description and content stored in $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 whether we are getting a seq from
# file or from text boxes in the input form. From here on, we have one valid sequence to work on.
use Bio::PrimarySeq;
$seq1=Bio::PrimarySeq->new(-seq => $seqContent,
				-primary_id=> $seqName,
				-molecule=> 'dna'
				);
#If sequence is circular define it.
if ($seqType eq 'isCircular'){
	$seq1->is_circular(1);
}

# IncisiEnz performs a multiple-digest , that is the sequence can be simultaneously cut by more than one enzyme 
# though the issues such as buffer compatibility and others should be considered by user and not the application.

# To perform a multiple digest, a sequence and a cusytomly defined collection of selected enzymes are needed as input.

# Make available necessary Bioperl modules
use Bio::Restriction::EnzymeCollection;
use Bio::Restriction::Enzyme;
use Bio::Restriction::Analysis;

# 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
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.
my @selectedEnzymes=();
foreach $enz(@enzymeNames){
 push(@selectedEnzymes, ($enzCollection->get_enzyme($enz)));
}

# 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);

# Using Bio::Restriction::Analysis of Bioperl, make an analyzer object to perform cuts
# using the sequence and the custom enzyme collection ($selEnzCollection) created above.
my $analyzer=Bio::Restriction::Analysis->new(-seq => $seq1);
#Cut the sequence in multiple mode
$analyzer->cut('multiple',$selEnzCollection);

# Get the product of the multiple-digest analysis and store the fragments, cut positions and fragment sizes.
@resultFragments=$analyzer->fragments('multiple_digest');
my @positions=$analyzer->positions('multiple_digest');
my @sizes=$analyzer->sizes('multiple_digest');

# Sort fragment sizes and fragments themselves by length.
@sizesSorted=sort {$a<=>$b} @sizes;
@fragmentsSorted=sort {length($a)<=>length($b)} @resultFragments;

#If no cut happened exit and warn user. 
if ((scalar@resultFragments) < 2){
	print "<p align='center'><BR><BR><b>Enzymes you selected produced no cuts. <BR><BR><a href='IncisiEnz.pl'>Back</a></b></p>";
	exit;	
} 
#------------------------Print a summary of analysis----------------------
# This uses Bio::Restriction::Enzyme to get overhang, cut site and No of cuts for each enzyme OBJECT for the sequence.
print "<p>The digest results for <b>". $seq1->primary_id ."</b> with the following enzyme(s):<br><font size='1'>The enzymes with no cuts will have a blue background.</font><p>";
print "<Table style='font-family:Tahoma; font-size:12px' border='1' cellspacing='1'><tr bgcolor='pink'><th>Enzyme</th><th>Overhang</th><th>Overhang sequence</th><th>Cut site</th><th>No of cuts</th></tr>";
foreach my $enz(@selectedEnzymes){
	my $overhangSeq="--";
	if (length($enz->overhang_seq) > 0){
		$overhangSeq = $enz->overhang_seq;
	}
	my $bgCol="";
	#If an enzyme does not cut the sequence give its row a skyblue background.
	if ($analyzer->cuts_by_enzyme($enz->name) == 0){
		$bgCol="bgColor='#CCCCFF'";
	}
	print "<tr $bgCol><td><a target='_blank' href='http://rebase.neb.com/rebase/enz/". $enz->name .".html'>". $enz->name ."</a></td><td><center>". $enz->overhang ."</center></td><td><center> $overhangSeq </center></td><td>". $enz->site ."</td><td><center>". $analyzer->cuts_by_enzyme($enz->name) ."</center></td></tr>";
}
print "<tr><td colspan='5'><center><b>Number of fragments produced: ". scalar @resultFragments ."</b></center></td></Table ></p>";


#-------------------------------- show Length view -----------------------------------
# Length view gives a quick virtual image of the lengths of the fragments produced in the digest sorted by length.
# This is particularly handy when the fragments are too large and getting a sequence view alone can not produce
# a good virtual image effect.
# Lengths are produced so that the width of a length view will always be 400 pixels so that the page will be 
# more PRINTER friendly. A $sizeCoefficient is calculated using the longest fragment length which is retrieved from
# the $sizesSorted array as its last member having [-1] index.
# In short, smaller-than-400-bp fragments may expand while the bigger-than-400 bp ones shrink!
# The No 400 was chosen because when I tested it in IE 6 with higher figures I got wrong print previews though 
# Netscape 4.0 did not have this bug.

if ($viewLength eq 'showLength'){
	print "<HR><BR><CENTER><H4>Fragments Length View <a target='_blank' href='../../rest_html/help_IncisiEnz.htm#output'>[?]</a></H4></CENTER><br>The lengths are adjusted by sequence sizes to fit a width of 400 pixels on page. Therefore, lengths may vary for each set of fragments from different digests depending on the size of the longest fragment produced.<br><br>";
	
	print "<Table style='font-family:Tahoma; font-size:12px' ><tr><td>No</td><td align='center'>Size (bps)</td><td></td></tr>";
	my $sizeCoefficient=400/$sizesSorted[-1];
	
	for ($i=0; $i<scalar(@fragmentsSorted);$i++){
	print "<tr><td width='50'>". ($i+1), "-> </td><td><b>". length($fragmentsSorted[$i]) ." </b></td><td><HR align='left' size='5' noshade color='blue' width='". ($sizeCoefficient) * ($sizesSorted[$i]) ."'></td></tr> " ;
	}
	print "</Table>";
}
#-------------------------------- End Length view -----------------------------------

#-------------------------------- show Sequence view -----------------------------------
# This is useful to enable user to see the content of sequences and copy them by double-click to select and then copy.
if ($viewSeq eq 'showSeq'){
	print "<HR><BR><CENTER><H4>Sequence View <a target='_blank' href='../../rest_html/help_IncisiEnz.htm#output'>[?]</a></H4></CENTER>";
	print "<Table align='center' border ='1' style='font-family:Courier New; font-size:12px' ><tr bgcolor='pink' style='font-family:Tahoma; font-size:12px'><td>No</td><td>Length</td><td>Fragment sequence (To copy a sequence, double click to select.)</td></tr>";
	for ($i=0; $i<scalar(@fragmentsSorted);$i++){
	print "<tr><td>".($i+1)."</td><td><b>". length($fragmentsSorted[$i]) ."</b></td><td> ". ($fragmentsSorted[$i])."</td></tr>" ;
	}
	print "</Table>";
}
#-------------------------------- End Sequence view -----------------------------------

#-------------------------------- show Map view -----------------------------------

# Map view is especially useful to get a virtual image of what exactly happened to the sequence by showing cut sites with
# enzyme name. It will also enable the user to get a printer friendly view of the digest-map. This map view was completely
# written out of scratch by myself and I did not use any sample codes.
# The map wraps at every 50 base pair length of the sequence.
if ($viewMap eq 'showMap'){	
	print "<HR><BR><CENTER><H4>Map View <a target='_blank' href='../../rest_html/help_IncisiEnz.htm#output'>[?]</a></H4></CENTER>";
	
	# First get an array containing all cut positions by ALL enzymes selected.
	my @AllPositions=$analyzer->positions('multiple_digest');

	# Make an empty hash and fill it with Keys: "cut positions" and values: "enzyme name" for 
	# the corresponding cut position.
	%PositionEnzymeNames=();
	foreach $enz(@enzymeNames){
		foreach $pos($analyzer->positions($enz)){
			$PositionEnzymeNames{$pos}=$enz;	
		}
	}
	# Store cut positions in an array to be used in map production
	@enzKeys=keys %PositionEnzymeNames;
	
	print "<BR>";
	
	# Loop through the sequence and detect the positions of wrapping (multiples of 50) and add a table row there to
	# wrap in HTML Table output
	$strMap="<Table id='mapTable' border='0' cellspacing='1' cellpadding='2' align='center' style='font-family:Courier New, Courier; font-size: 12px'><tr><td>start</td><td></td><td>end</td></tr><tr><td valign='top'>0</td><td>";
	#Start looping through sequence and synthesis of the map string
	for ($i=0; $i<length($seqContent); $i++){
		#Detect 50 bp milestones and add a table row where each was found.
		if((($i+1) % 50 == 0)){
			$strMap .= "<td valign='bottom'>".($i)."</td></tr><tr><td valign='top'>".($i+1)."</td><td>";
		}
		$strMap .= (substr $seqContent,$i,1) ;	
		
		#Whenever a cut position is reached, add the name of the corresponding cutter enzyme from tha hash %PositionEnzymeNames.
		if($PositionEnzymeNames{$i + 1}){
		$strMap .= "<span id='$PositionEnzymeNames{$i + 1}'><b><SUP><font style='font-size:10px'>".($i+1)."^</font></SUP><font id='enz' style='background-color:pink; font-size:11px'>[<a target='_blank' href='http://rebase.neb.com/rebase/enz/". $PositionEnzymeNames{$i + 1} .".html'>$PositionEnzymeNames{$i + 1}</a>]</font></b></span>";
		}
		
	}
	
	#Output the map to HTML
	print "$strMap</td><td valign='bottom'>". length($seqContent). "</td></tr></Table>";
}

#If sequence is circular remind the user.
if ($seqType eq 'isCircular'){
	print "<p align='center'><font color='blue'>Notice: the sequence is CIRCULAR.</font></p><BR>";
}
print "<p align='center'>If the map rows has wrapped click here: <a href='javascript:straighten()'>Make map No Wrap</a></p>";
#-------------------------------- End Map-----------------------------------
# print "<select name='cboEnzList' id='cboEnzList'>";
# foreach $enz(@enzymeNames){
# 	print "<option value='$enz'>$enz</option>";	
# }
# print "</select><input type='button' onClick='javascript:sh(document.all.cboEnzList.value);' name='btnShow'>";		
# Print a small footer to allow the user to go back to the top of page or to start a new analysis.
print "<p align='center'><BR><BR><a href='#'>^Top</a> | <a href='IncisiEnz.pl'>New Incision</a></p><br>";	




print "</body></html>";

