#!/usr/bin/perl
package RestrictionMap;

#######################################################################
##### Author :	    Shuai Weng
##### Date   :      Jan. 2002
##### Description : This package contains all necessary methods for 
#####               performing the restriction analysis for a given 
#####               sequence name or a raw DNA sequence. 
#####              
#######################################################################
use strict;
use DBI;
use GD;
use CGI qw/:all/;
use Bio::SeqIO;
use Bio::Seq;
use CGI::Carp qw(fatalsToBrowser);

use lib "/usr/local/dicty/www_dictybase/db/lib/common";
use Login qw (ConnectToDatabase);
use TextUtil qw (DeleteUnwantedChar);
use lib "/usr/local/dicty/www_dictybase/db/lib/dictyBase";
use dictyBaseCentralMod qw(:formatPage);
use lib "/usr/local/dicty/www_dictybase/db/lib/dictyBase/Objects";
use SeqParamTranslator;
use ConfigURLdictyBase;
use ConfigPathdictyBase;
use Feature;
use Sequence;
use CreateGD;
use GCG;

#######################################################################
#################### global variables #################################
#######################################################################

my $dbh;
my $dblink; 
my $configUrl;
my $configPath;
my $scan4matches;
my $width = 620;
my ($y1, $x1, $x2, $diff);

#######################################################################
sub new {      ############ constructor ###############################
#######################################################################
       
	my ($self, %args) = @_;

	$self = {};

	bless $self;

      	$self->{'_database'} = $args{'database'};
	$self->{'_help'} = $args{'help'};
	$self->{'_title'} = $args{'title'};
	$self->{'_query'} = $args{'query'};
	$self->{'_sequence'} = $args{'sequence'};

    	return $self;
}

sub help { $_[0]->{_help} }
sub database { $_[0]->{_database} }
sub title { $_[0]->{_title} }


######################################################################
sub DESTROY {   ############ destructor ##############################
######################################################################
    	if (defined $dbh) {
		$dbh->disconnect;
    	}
}

######################################################################
sub start{
######################################################################

    	my ($self) = @_;

	$configUrl = ConfigURLdictyBase->new;
	$dblink = $configUrl->dblink($self->database);
	$configPath = ConfigPathdictyBase->new;

	if (param('id') && param('enzyme')) {
	    $self->displayFragments;
	}
	elsif (param('id') && param('seqlength')) {
	    $self->displayResult; 
	}
	elsif (!$self->{'_query'} && !$self->{'_sequence'}) {
	    $self->printEntryForm;
	}
	else {
	    $self->displayResult;
	}
	exit;
}

#######################################################################
sub printEntryForm {
#######################################################################
        my ($self) = shift;

	&printStartPage($self->database, $self->title, $self->help);

	if (param) {
	    $self->getSequence;
	}

	print blockquote("This form allows you to perform a restriction analysis by entering a sequence name or arbitrary DNA sequence."), p;
	print start_form, 
	      table(Tr(td(table({-width=>'100%',
				 -border=>'0',
				 -cellpadding=>'10'},
				Tr(td({-bgcolor=>'#b7d8e4',
				       -width=>'40%',
				       -valign=>'top'}, 
				      b(big("Enter a Sequence Name")).br.
				      "or the first few characters followed by *:".br.
				      textfield(-name=>'seqname',
						-value=>$self->{'_seqname'},
						-size=>'12').br.
				      b("Examples").":".br.
				      "Gene name - act1, Act*".br.
				      "ORF name - YHR023W, yal007*".br.
				      "GenBank name - SC8339, YSCACT*".br.
				      "Clone name - 70353").
				   td({-valign=>'top',
				       -align=>'center'},
				      b(font({-size=>'5',
					      -color=>'red'}, "OR"))).
				   td({-bgcolor=>'#b7d8e4',
				       -width=>'40%',
				       -valign=>'top'},
				      b(big("Type or Paste a DNA Sequence:")).br.
				      textarea(-name=>'sequence',
					       -value=>$self->{'_sequence'},
					       -cols=>'25',
					       -rows=>'7').br.
				      "The sequence can be provided in raw, fasta or GCG format without comments (numbers are okay)."))))).
		    Tr(td({-align=>'left'},
			  b(big("Choose Restriction Enzymes:")).br.
			  popup_menu(-name=>'cut',
				     -values=>["all", 
					      "3'overhang", 
					      "5'overhang",
					      "blunt end",
					      "cut once",
					      "cut twice",
					      "Six-base cutters"]).br.
			  font({-color=>'red'},
			       b("Note")).": To find enzymes that do not cut, choose 'all' and see the resulting list at bottom.".hr.
			  submit(-name=>'Display Map')." or ".
			  reset(-name=>'Reset Form'))));

	&printEndPage;
}

########################################################################
sub displayResult {
########################################################################
    my ($self) = @_;
    
    ###############################################################

    $dbh = &ConnectToDatabase($self->database); 

    $self->setVariables;

    ###############################################################

    $self->{'_title'} = "Genome Restriction Map";
    if (!$self->{'_showNm'}) {
	$self->{'_showNm'} = $self->{'_query'};
    }
    if ($self->{'_showNm'} && $self->{'_showNm'} !~ /^unknown/i) {
	$self->{'_title'} .= " of ".uc($self->{'_showNm'});
    }

    &printStartPage($self->database, $self->title, $self->help);

    ################################################

    ####### create tmp sequence file

    my $seqfile = $self->createSeqFile;

    ####### do the search
    my $outfile = $self->doSearch($seqfile);
    
    ####### process output 
    my ($dataHashRef, $notCutEnzymeArrayRef, 
	$offsetHashRef, $overhangHashRef) 
	 = $self->processData($outfile);
	
    ###############################################################

    ####### display graph
    $self->createGraph($dataHashRef, $offsetHashRef, $overhangHashRef);

    ####### display enzymes that do not cut
    $self->printEnzymes($notCutEnzymeArrayRef);

    ####### display legend graph
    $self->createLegendGraph;

    $dbh->disconnect;

    &printEndPage;
}

########################################################################
sub createSeqFile {
########################################################################
    my ($self) = @_;

    if (param('id')) { return; }
	
    my $seqfile = $configPath->tmpDir."restriction.tmp.$$"; 

    open(OUT, ">$seqfile") || 
	die "Can't open '$seqfile' for writing:$!\n";

    if ($self->{'_query'}) {

	my $seq;
	if ($self->{'_type'} eq 'gb') {
	    
	    my $gcg = GCG->new(map=>$self->{'_map'},
			       seqname=>$self->{'_query'},
			       rev=>$self->{'_rev'});
	    my $tmpseqfile = $gcg->getSequence;

	    open(IN, $tmpseqfile) ||
		die "Can't open '$tmpseqfile' for reading:$!\n";
	    while(<IN>) {
		chomp;
		$seq .= $_;
	    }
	    close(IN);
	
	}
	elsif ($self->{'_type'} =~ /^clone/i) { ## it is clone
	    my $seqObj = Sequence->new(dbh=>$dbh,
				       chrnum=>$self->{'_chr'},
				       beg=>$self->{'_beg'},
				       end=>$self->{'_end'});
	    $seq = $seqObj->getSequence;
	}
	else {

	    if (!$self->{'_feature_name'}) {
		print "No sequence found for ".$self->{'_query'}, p;
		exit;
	    }
	    
	    my $seqObj = Sequence->new(dbh=>$dbh,
				       feature_name=>$self->{'_feature_name'});
	    $seq = $seqObj->getSequence;
	}
	$self->{'_seqlength'} = length($seq);
	print OUT ">", $self->{'_query'}, "\n";
	print OUT $seq, "\n";
    }
    else {
	$self->{'_sequence'} =~ s/[^A-Za-z]//g;
	$self->{'_seqlength'} = length($self->{'_sequence'});
	print OUT ">unknown\n";
	print OUT $self->{'_sequence'}, "\n";	
    }
    close(OUT);

    return $seqfile;

}

########################################################################
sub doSearch {
########################################################################

    my ($self, $seqfile) = @_;

    if (param('id')) { 
	return $configPath->tmpDir."restriction.out.".param('id').".tmp";
    }

    my $patfile = $self->setPatFile;

    my $outfile = $configPath->tmpDir."restriction.out.$$.tmp";

    open(PAT, "$patfile") || die "Can't open '$patfile' for reading:\n";

    while(<PAT>) {

	chomp;

	my ($enzyme, $offset, $pat, $overhang) = split(/ /);

	my $patfile = $configPath->tmpDir."tmp.pat.$$";
	open(TMP,">$patfile");
	print TMP $pat;
	close(TMP);

	open(OUT, ">>$outfile") || 
	    die "RestrictionMapper: Can't open '$outfile' for appending:$!\n"; 
	print OUT "\>\>$enzyme: $offset $overhang $pat\n";
	close (OUT);

	system "/usr/bin/nice -10 $scan4matches -c $patfile < $seqfile >> $outfile";    

    }
    close(PAT);

    return $outfile;
}

########################################################################
sub processData {
########################################################################
    my ($self, $datafile) = @_;

    my (%dataHash, %offset, %overhang);
    my @notCutEnzyme;
 
    open(IN, $datafile) || die "Can't open '$datafile' for reading:$!\n";

    my $preLine;
    my $enzyme;
    while(<IN>) {
	chomp;
	if (/^>>([^\:]+)\: ([^ ]+) ([^ ]+) ([^ ]+)$/) {
	    $enzyme = $1;
	    $offset{$enzyme} = $2;
	    $overhang{$enzyme} = $3;
	    if ($preLine =~ /^>>([^\:]+)\: ([^ ]+) ([^ ]+) ([^ ]+)$/) {
		push(@notCutEnzyme, $1);
	    }
	}
	elsif (/\>.+\[([0-9]+\,[0-9]+)\]$/) {
	    if ($dataHash{$enzyme}) { $dataHash{$enzyme} .= ":"; }
	    $dataHash{$enzyme} .= $1;
	}
	$preLine = $_;
    }
    close(IN);
    if ($preLine =~ /^\>\>([^\:]+)\: ([^ ]+) ([^ ]+) ([^ ]+)$/) {
	push(@notCutEnzyme, $1);
    }

    if (param('cut') =~ /^(cut once|cut twice)/i) {

	my $cutLimit;
	if (param('cut') =~ /^cut once/i) { $cutLimit = 1; }
	else { $cutLimit = 2; } 
	    
	my %newDataHash;
	foreach my $key (keys %dataHash) {
	    my $coords = $dataHash{$key};
	    my @coords = split(/:/, $coords);
	    my $wCut;
	    my $cCut;
	    foreach my $coordPair ( @coords) {
		my ($beg, $end) = split(/\,/, $coordPair);
		if ($beg < $end) { $wCut++; }
		else { $cCut++; }
	    }
	    if (($cCut == $cutLimit && $wCut <= $cutLimit) ||
		($wCut == $cutLimit && $cCut <= $cutLimit)) {
	        $newDataHash{$key} = $coords;  
	    }
	}
	return (\%newDataHash, \@notCutEnzyme, \%offset, \%overhang);
    }
    else {
	return (\%dataHash, \@notCutEnzyme, \%offset, \%overhang);
    }
}

########################################################################
sub createGraph {
########################################################################
    my ($self, $dataHashRef, $offsetHashRef, $overhangHashRef) = @_;

    my $hitNum = keys(%$dataHashRef);

    my $beg = param('beg');
    if (!$beg) { $beg = 1; }

    my $displayNum;
    if ($hitNum <= 100) { 
	$displayNum = $hitNum; 
    }
    elsif ($hitNum <= 199) { 
	if ($beg == 1) { $displayNum = 100; }
	else { $displayNum = $hitNum - $beg + 1; }
    }
    else {
	if ($beg <= 100) {
	    $displayNum = 100;
	}
	else { $displayNum = $hitNum - $beg + 1; }
    }

    my $height = 2*$y1 + $displayNum*$diff + 5;

    my $gd = CreateGD->new(width=>$width,
			   height=>$height);

    print "<map name='RestGif'>";

    ### Put a blue frame around the graph 
    $gd->im->rectangle(0, 0, $width-1, $height-1, $gd->blue);
     
    ### Draw a dictyBase label and date on the bottom of the map
    my $date = `date`;
    chomp $date;
    $date =~ s/^[A-Za-z]+ ([A-Za-z]+) +([0-9]+) .+ ([0-9]+)$/$1 $2, $3/;
    $gd->im->string(gdSmallFont, 5, $height-15, "dictyBase", $gd->red);
    $gd->im->string(gdSmallFont, $width-77, $height-15, $date, $gd->red);

    ### Draw tick marks and units on the top ruler   
    $gd->im->string(gdMediumBoldFont, $width/2-15, 4, "bp", $gd->black);
    $gd->im->line($x1, $y1, $x2, $y1, $gd->black);

    my $tickSize;
    if ($self->{'_seqlength'} > 10000) {
	$tickSize = 5000;
    }
    elsif ($self->{'_seqlength'} > 5000) {
	$tickSize = 2000;
    }
    elsif ($self->{'_seqlength'} > 2000) {
	$tickSize = 1000;
    }
    elsif ($self->{'_seqlength'} > 500) {
	$tickSize = 500;
    }
    else { $tickSize = 100; }
    
    my $tickNum = $self->{'_seqlength'}/$tickSize;
    my $rulerWidth = $x2 - $x1;
    my $pixelNum = $rulerWidth*$tickSize/$self->{'_seqlength'};
    my $i;
    for ($i = 0; $i <= $tickNum; $i++) {
	$gd->im->line($x1+$i*$pixelNum+1, $y1, 
		      $x1+$i*$pixelNum+1, $y1+5,
		      $gd->black);
	$gd->im->string(gdSmallFont, $x1+$i*$pixelNum, 
			$y1-14, $i*$tickSize, $gd->black);
	if ($i > 0) {
	    $gd->im->dashedLine($x1+$i*$pixelNum+1, $y1+10, 
				$x1+$i*$pixelNum+1, $height-$y1,
				$gd->magenta);
	}
    }
    if ($self->{'_seqlength'} - ($i-1)*$tickSize > 30) {
	$gd->im->line($x1+$rulerWidth, $y1, 
		      $x1+$rulerWidth, $y1+5,
		      $gd->black);
	$gd->im->string(gdSmallFont, $x1+$rulerWidth-1, 
			$y1-17, $self->{'_seqlength'},
			$gd->black);	
    }

    ### Draw the lines 
    $y1 += 5;
    for (my $i = 1; $i <= $displayNum; $i++) {
	$gd->im->line($x1, $y1+$diff*$i, $x2, $y1+$diff*$i, 
		      $gd->black);
    }

    my $id;
    if (param('id')) { $id = param('id');}
    else { $id = $$; }

    ### Draw cut tick and enzyme name
    my ($count, $i);
    my %enzymeBYindex;
    foreach my $enzyme (sort (keys %$dataHashRef) ) {

	$count++;

	$enzymeBYindex{$count} = $enzyme;

	if ($count < $beg || $count >= $beg + 100) { 
	    next; 
	}
	
	$i++;

	my $offset = $$offsetHashRef{$enzyme};
	my $overhang = $$overhangHashRef{$enzyme};

	my $nameColor;
	if ($overhang == 0) {
	    $nameColor = $gd->orange;
	}
	elsif ($overhang < 0 ) {
	    $nameColor = $gd->darkGreen;
	}
	else {
	    $nameColor = $gd->magenta;
	}
	    
	$gd->im->string(gdSmallFont, $x2+5, 
			$y1+$diff*$i-5, $enzyme, $nameColor);
	
	my $NmX1 = $x2 + 5;  
	my $NmY1 = $y1+$diff*$i-5;
        my $NmX2 = $x2 + 5 + length($enzyme)*6;
        my $NmY2 = $y1+$diff*$i+3;
	my $showNm = $self->{'_showNm'} || $self->{'_query'};
	my $url = url."?id=$id&enzyme=$enzyme&showNm=$showNm&seqlength=".$self->{'_seqlength'}."&cut=".param('cut');

	##### for handling the genBank seqname
        ##### for some reason, there is an extra '(gb)' in the url... 
	$url =~ s/\(gb\).*(http.+)$/$1/;

	print "<area shape=\"rect\" coords=\"$NmX1,$NmY1,$NmX2,$NmY2\" href=\"$url\" target='infowin'>";

	my $coords = $$dataHashRef{$enzyme};
	my @coords = split(/:/, $coords);

	foreach my $coords (@coords) {

	    my ($beg, $end) = split(/\,/, $coords);
	    if ($beg < $end) { ### watson
		$gd->im->filledRectangle(
			 $x1+($x2-$x1)*($beg+$offset)/$self->{'_seqlength'}, 
			 $y1+$diff*$i-5, 
			 $x1+($x2-$x1)*($beg+$offset)/$self->{'_seqlength'}+1, 
			 $y1+$diff*$i, $gd->red);
	    }
	    else {
		$gd->im->filledRectangle(
			 $x1+($x2-$x1)*($end+$offset+$overhang)/$self->{'_seqlength'}, 
			 $y1+$diff*$i, 
			 $x1+($x2-$x1)*($end+$offset+$overhang)/$self->{'_seqlength'}+1,
			 $y1+$diff*$i+5, $gd->blue);
	    }
	}
    }

    print "</map>";

    ##############################################################
    my ($prevBeg, $prevEnd);
    my ($nextBeg, $nextEnd);
    if ($beg > 1) {
	$prevEnd = $beg;
	$prevBeg = $beg - 100 + 1;
	if ($prevBeg <= 0) { $prevBeg = 1; }
    }
    if ($beg + 100 - 1 < $hitNum) {
	$nextBeg = $beg + 100 - 1;
	$nextEnd = $nextBeg + 100 - 1;
	if ($nextEnd > $hitNum) { $nextEnd = $hitNum; }
    }
    my $prevBegEnzyme = $enzymeBYindex{$prevBeg};
    my $prevEndEnzyme = $enzymeBYindex{$prevEnd};
    my $nextBegEnzyme = $enzymeBYindex{$nextBeg};
    my $nextEndEnzyme = $enzymeBYindex{$nextEnd};

    ##############################################################

    my $enzymeType = param('cut') || 'all';
    
    print p, center(b("Restriction Enzyme: ".
			  font({-color=>'red'},
			       $enzymeType)));
    if ($prevBeg || $nextBeg) {
	my $end = $beg + 100 - 1;
	if ($end > $hitNum) { $end = $hitNum; }
	my $begEnzyme = $enzymeBYindex{$beg};
	my $endEnzyme = $enzymeBYindex{$end};
	print center(b("( $begEnzyme to $endEnzyme )"));
    }
    print p, center(b(font({-color=>'red'}, 
			   "Click on any enzyme name to view its fragment sizes"))),p;
    
    ##############################################################

    my $gifDir = $configPath->tmpDir;
    my $gifHttp = $configUrl->dictyBaseHtmlTmp;
    my $gifname = "RESTmap.".$id.".".$beg.".gif";

    open (OUT, ">".$gifDir.$gifname)||die "Can't create ".$gifDir.$gifname.":$!";
    binmode OUT;
    print OUT $gd->im->gif;
    close OUT;
    print center(img({-src=>$gifHttp.$gifname,
		      -usemap=>'#RestGif'}));

    ##############################################################

    if ($prevBeg) { 
	print p, center("View previous page of Restriction Map, ".a({-href=>url."?id=$id&beg=$prevBeg&seqlength=".$self->{'_seqlength'}."&showNm=".$self->{'_showNm'}."&cut=".param('cut')}, "enzyme $prevBegEnzyme to $prevEndEnzyme"));
    }
    if ($nextBeg) {
	print p, center("View next page of Restriction Map, ".a({-href=>url."?id=$id&beg=$nextBeg&seqlength=".$self->{'_seqlength'}."&showNm=".$self->{'_showNm'}."&cut=".param('cut')}, "enzyme $nextBegEnzyme to $nextEndEnzyme"));
    }

    ##############################################################
}

########################################################################
sub createLegendGraph {
########################################################################
    my ($self) = @_;

    my $gifDir = $configPath->tmpDir;
    my $gifHttp = $configUrl->dictyBaseHtmlTmp;
    my $gifname = "RESTlegend.gif";

    if (! -e "$gifDir.$gifname") {

	my $gd = CreateGD->new(width=>$width,
			       height=>180);

	### Put a blue frame around the graph 
	$gd->im->rectangle(0, 0, $width-1, 179, $gd->blue);

	### key
	$gd->im->string(gdMediumBoldFont, 15, 5, "Key:", $gd->red);
	$gd->im->filledRectangle(15, 30, 17, 36, $gd->red);
	$gd->im->string(gdMediumBoldFont, 15, 27, " = Recognition sequence in Watson (5'-->3') strand", $gd->black);
	$gd->im->filledRectangle(15, 53, 17, 60, $gd->blue);
	$gd->im->string(gdMediumBoldFont, 15, 50, " = Recognition sequence in Crick (3'-->5') strand", $gd->black);
	
	### note
	$gd->im->string(gdMediumBoldFont, 15, 75, "Notes:", $gd->red);
	$gd->im->string(gdMediumBoldFont, 15, 97, "1. Each enzyme name is listed on the right of the map. Click on the name to retrieve", $gd->black);
	$gd->im->string(gdMediumBoldFont, 15, 112, "   the fragment size details.", $gd->black);
	$gd->im->string(gdMediumBoldFont, 15, 127, "2. ", $gd->black);
	$gd->im->string(gdMediumBoldFont, 35, 127, "Green", $gd->darkGreen);
	$gd->im->string(gdMediumBoldFont, 70, 127, " enzyme name   = 3' overhang", $gd->black);
	$gd->im->string(gdMediumBoldFont, 35, 142, "Magenta", $gd->magenta);
	$gd->im->string(gdMediumBoldFont, 84, 142, " enzyme name = 5' overhang", $gd->black);
	$gd->im->string(gdMediumBoldFont, 35, 157, "Orange", $gd->orange);
	$gd->im->string(gdMediumBoldFont, 77, 157, " enzyme name  = blunt end", $gd->black);
   
	open (OUT, ">".$gifDir.$gifname)||die "Can't create ".$gifDir.$gifname.":$!";
	binmode OUT;
	print OUT $gd->im->gif;
	close OUT;
    }
    print p, center(img({-src=>$gifHttp.$gifname}));
}

########################################################################
sub printEnzymes {
########################################################################
    my ($self, $notCutEnzymeArrayRef) = @_;

    if (!@$notCutEnzymeArrayRef) { return; }

    if (param('cut') =~ /^Six-base/i) {
	print p, "Six-base cutter enzymes that do not cut:", p;
    }
    elsif (param('cut') =~ /^3/i) {
	print p, "3' overhang enzymes that do not cut:", p;
    }
    elsif (param('cut') =~ /^5/i) {
	print p, "5' overhang enzymes that do not cut:", p;
    }
    elsif (param('cut') =~ /^blunt/i) {
	print p, "Blunt end enzymes that do not cut:", p;
    }
    else {
	print p, "Enzymes that do not cut:",p;
    }
    print "<pre>";
    my ($rows, $row, $count);
    my %found;
    foreach my $enzyme (@$notCutEnzymeArrayRef) {
	if ($found{$enzyme}) { next; }
	$count++;
	if ($count < 10) {
	    print $enzyme, "\t";
	}
	else {
	    print $enzyme, "\n";
	    $count = 1;
	}
	$found{$enzyme}++;
    }
    print "</pre>";
}

########################################################################
sub displayMatchList {
########################################################################
    my ($self, $matchList) = @_;

    &printStartPage($self->database, $self->title, $self->help);

    print "Available sequence names that begin with ".b($self->{'_query'}.":"),p;

    my @seqname = split(/\:/, $matchList);
    my $list;
    foreach my $seqname (@seqname) {
	my $linkNm = $seqname;
	$linkNm =~ s/ \(GB\)$//i;
	$list .= li(a({-href=>url."?seqname=$linkNm&cut=".param('cut')}, $seqname));
    }
    print ul($list),p;
    
    &printEndPage;
}

########################################################################
sub displayFragments {
########################################################################
    my ($self) = @_;
    
    my $showNm = uc(param('showNm')) || 'Unknown Sequence';
    $self->{'_title'} = "Fragment Sizes for $showNm Cut with ".param('enzyme');
    
    &printStartPage($self->database, $self->title, $self->help);

    my $datafile = $configPath->tmpDir."restriction.out.".param('id').".tmp";

    open(IN, $datafile) || die "Can't open '$datafile' for reading:$!\n";
    
    my $restEnzyme = param('enzyme');

    my $i = 0, 
    my $cutNo = 0;
    my $enzyme;
    my ($offset, $overhang, $pattern);
    my (@cutNum, @Num1, @Num2);
    while(<IN>) {
	chomp;
	if (/^\>\>$restEnzyme\: ([0-9]+) (-?[0-9]+) (.+)$/) {
	    $offset = $1;
	    $overhang = $2;
	    $pattern = $3;
	}
	if ($pattern && /^\>[A-Za-z0-9\_\-]+\:\[([0-9]+)\,([0-9]+)\]$/) {
	    my $num1 = $1;
	    my $num2 = $2;
	    if ($num1 < $num2) {
		$Num1[$i] = $num1;
		$Num2[$i] = $num2;
		$cutNum[$cutNo++] = $num1 + $offset;
		$i++;
	    }
	    else {
		my $found;
		for (my $j = 0; $j < @Num1; $j++) {
		    if ($num1 == $Num2[$j] || $num2 == $Num1[$j]) {
			$found++;
			last;
		    }
		}
		if (!$found) {
		    $cutNum[$cutNo++] = $num2 + $offset + $overhang;
		}
	    }
	}
	if (/^\>\>([^\:]+)\:/) {
	    $enzyme = $1;
	}
	if ($pattern && $enzyme ne $restEnzyme) { last; }
    }
    close(IN);

    ############################################################

    print table({-align=>'center'},
		Tr(th({-align=>'right'},
		      font({-color=>'red'},
			   "offset (bp)")." : ").
		   td({-align=>'left'}, 
		      $offset)).
		Tr(th({-align=>'right'},
		      font({-color=>'red'},
			   "overhang (bp)")." : ").
		   td({-align=>'left'}, 
		      $overhang)).
		Tr(th({-align=>'right'},
		      font({-color=>'red'},
			   "recognition sequence")." : ").
		   td({-align=>'left'},
		      $pattern)));

    ############################################################

    my @sortedCutNum = sort { $a <=> $b }(@cutNum);
    unshift(@sortedCutNum, '0');
    push(@sortedCutNum, param('seqlength'));

    my $i;
    my $rows;
    my $preCut;
    my @SIZE;
    foreach my $cut (@sortedCutNum) {
	$i++;
	if ($i != 1) {
	    my $cutSize = $cut - $preCut;
            push(@SIZE, $cutSize);
	    $rows .= Tr({-align=>'center'},
			td(br).
			td($cutSize));
	}
	$rows .= Tr({-align=>'center'},
		     td($cut).
		     td(br));
        $preCut = $cut;
    }

    print p, table({-align=>'center',
		    -border=>'3'},
		   Tr({-align=>'center'},
		      th("Cut Site").
		      th("Fragment Size (bp)")).
		   $rows);


    ############################################################
 
    my $rows;
    foreach my $size (reverse( sort{ $a <=> $b }(@SIZE) )) {
	$rows .= Tr({-align=>'center'},
		    td($size));
    }
    print p, table({-align=>'center'},
		   Tr({-align=>'center'},
		      th("Sorted Fragment Size (bp)")).
		   $rows);

    print p, center(b(a({-href=>url."?id=".param('id')."&showNm=$showNm&seqlength=".param('seqlength')."&cut=".param('cut')}, 
			"Restriction Map for $showNm"))),p;

    &printEndPage;
}

########################################################################
sub setPatFile {
########################################################################
    my ($self) = @_;

    my $patfile;
    if (param('cut') =~ /^3/) {
	$patfile = $configPath->dataDir."rest_enzymes.3";
    }
    elsif (param('cut') =~ /^5/) {
	$patfile = $configPath->dataDir."rest_enzymes.5";
    }
    elsif (param('cut') =~ /^blunt/) {
	$patfile = $configPath->dataDir."rest_enzymes.blunt";
    }
    elsif (param('cut') =~ /^cut once/) {
	$patfile = $configPath->dataDir."rest_enzymes";
    }
    elsif (param('cut') =~ /^cut twice/) {
	$patfile = $configPath->dataDir."rest_enzymes";
    }
    elsif (param('cut') =~ /^Six-base cutters/) {
	$patfile = $configPath->dataDir."rest_enzymes.6base";
    }
    else {
	$patfile = $configPath->dataDir."rest_enzymes";
    }
    return $patfile;
}

########################################################################
sub setVariables {
########################################################################
    my ($self) = @_;

    my $scriptPath = $ENV{'DOCUMENT_ROOT'};
    # $scriptPath =~ s/\/html\/$//;
    $scriptPath =~ s/\/html$//;
    $scriptPath .= $ENV{'SCRIPT_NAME'};
    $scriptPath =~ s/[^\/]+$//;  
    $scan4matches = $scriptPath."scan_for_matches";

    $width = 620;
    $y1 = 30;
    $x1 = 10;
    $x2 = $width - 60;
    $diff = 12;

    $self->{'_showNm'} = param('showNm');
    $self->{'_seqlength'} = param('seqlength');
    
    if (!$self->{'_query'}) { return; }

    $self->{'_query'} =~ s/ \(GB\)$//i;

    my $seqObj = SeqParamTranslator->new(dbh=>$dbh,
					 query=>$self->{'_query'});
	
    if ($seqObj->error) {
	$self->err_report($seqObj->error);
	exit;
    }
	
    if ($seqObj->matchList) {
	$self->displayMatchList($seqObj->matchList);
	exit;
    }

    $self->{'_type'} = $seqObj->type;

    if ($seqObj->type =~ /^clone/i) {
	$self->{'_chr'} = $seqObj->chromosome;
	$self->{'_beg'} = $seqObj->start_coord;
	$self->{'_end'} = $seqObj->stop_coord;
    }
    elsif ($seqObj->featureObject) {
	my $featObj = $seqObj->featureObject;
	if ($self->{'_query'} =~ /[\*\%]/) {
	    $self->{'_query'} = $featObj->feature_name;
	}
	$self->{'_feature_name'} = $featObj->feature_name;
	$self->{'_chr'} = $featObj->chromosome;
	$self->{'_beg'} = $featObj->start_coord;
	$self->{'_end'} = $featObj->stop_coord;
	if ($featObj->strand =~ /^C/i) {
	    ($self->{'_beg'}, $self->{'_end'}) = 
		($self->{'_end'}, $self->{'_beg'});
	    $self->{'_rev'} = '-REV';
	}
    }
    else {
	# $self->{'_map'} = 'a2map';
    }

    $self->{'_map'} = 'a2map';
}


########################################################################
sub createTmpSeqFile {
########################################################################
    my ($self, $sequence) = @_;
    my $SEQTMP = $configPath->tmpDir."gcgseq.tmp.$$";  # Sequence temp file 
    open(OUT, ">$SEQTMP") || die "seqTools: Can't create tmp seq file:$!\n";
    if (param('seqtype') =~ /^protein/i) {
	print OUT "\!\!AA\_SEQUENCE\n";
    }
    else {
	print OUT "\!\!NA\_SEQUENCE\n";
    }
    print OUT "temp sequence ..\n\n";
    print OUT "$sequence\n";
    close (OUT);
}

########################################################################
sub getSequence {
########################################################################
    my ($self) = @_;

    $dbh = &ConnectToDatabase($self->database); 

    my $orf = param('name') || param('orf') || param('locus');

    if ($orf) {
	&DeleteUnwantedChar(\$orf);
	my $featObj = Feature->new(dbh=>$dbh,
				   feature_name=>$orf);
	if (!$featObj) {
	    print "The orf name is not found in our database.", br;
	    exit;
	}
	my $seqObj = Sequence->new(dbh=>$dbh,
				   feature_name=>$orf);
	
	if (!$seqObj || !$seqObj->getSequence) {
	    print "No sequence info available for ".param('name'), br;
	    exit;
	}
	$self->{'_sequence'} = $seqObj->getSequence;

    }
    elsif (param('chr')) {
	my $chr = param('chr');
	my $beg = param('beg');
	my $end = param('end');
	$beg -= param('flankl');
	$end += param('flankr');
	if ($beg <= 1) {
	    $beg = 1;
	}
	my $seqObj =  Sequence->new(dbh=>$dbh,
				    chrnum=>$chr,
				    beg=>$beg,
				    end=>$end,
				    rev=>param('rev'));
	if (!$seqObj || !$seqObj->getSequence) {
	    print "No sequence info available for chromosome ".param('chr'),br;
	}
	$self->{'_sequence'} = $seqObj->getSequence;
    }
    elsif (param('id')) {
	my $id = param('id');
	my $tmpDir = $configPath->tmpDir;
	my $seqtmp;
	if (-e "${tmpDir}gcgseq.tmp.2.$id") {
	    $seqtmp = $tmpDir."gcgseq.tmp.2.$id";
	}
	else {
	    $seqtmp = $tmpDir."gcgseq.tmp.$id";
	}
	open(TMP, "$seqtmp") || 
	    die "Can't open '${tmpDir}gcgseq.tmp.$id' for reading:$!";
	my $startSeq;
	while(<TMP>) {
	    if (/\.\.$/) {
		$startSeq = 1;
	    }
	    elsif ($startSeq == 1) {
		$self->{'_sequence'} .= $_;
	    }
	}
	close(TMP);
    }
}

########################################################################
sub err_report {
########################################################################
    my ($self, $err) = @_;

    &printStartPage($self->database, $self->title, $self->help);
   
    print b($err);
    
    &printEndPage;

    $dbh->disconnect;
    exit;
}


########################################################################
1; #####################################################################
########################################################################