#!/usr/bin/perl
package RestGraphPage;

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

use MAP::RestGraph;

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

my $dbh;
my $dblink; 
my $configUrl;
my $configPath;

#######################################################################
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'},
			  table(Tr(td({-colspan=>2}, 
				      b("Choose Restriction Enzymes:"))).
				Tr(td("Overhang: ").
				   td(popup_menu(-name=>'overhang',
						 -values=>["all", 
							   "3' overhang", 
							   "5' overhang",
							   "blunt end"]))).
				Tr(td("Length of recognition site: ").
				   td(popup_menu(-name=>'siteLengthSign',
						 -values=>["=", ">", ">=", 
							   "<", "<="])." ".
				      popup_menu(-name=>'siteLength',
						 -values=>['all',
							   '4', '5', '6', 
							   '7', '8', '9', 
							   '10']))).
				Tr(td("Cut number: ").
				   td(popup_menu(-name=>'cutNumberSign',
						 -values=>["=", ">", ">=", 
							   "<", "<="])." ".
				      popup_menu(-name=>'cutNumber',
						 -values=>['all', 
							   '1', '2', '3', 
							   '4', '5', '6', 
							   '7', '8', '9',
							   '10'])))).p.
      			  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'} = "RestGraph";

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

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

    ####### retrieve sequence file
    $self->retrieveSequence;

    ####### print subtitle 
    $self->printSubTitle;

    ####### create and show map
    my $restGraph = $self->createAndShowMap;


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

    $dbh->disconnect;

    &printEndPage;
}

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

    if (param('id')) { return; }
	
    if ($self->{'_query'}) {

	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;
		$self->{'_sequence'} .= $_;
	    }
	    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'});
	    $self->{'_sequence'} = $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'});
	    $self->{'_sequence'} = $seqObj->getSequence;
	}
    }

    $self->{'_sequence'} =~ s/[^A-Za-z]//g;

}

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

    my $rows;

    if ($self->{'_overhang'}) {
	
	$rows = Tr(td({-align=>'right'}, b("Overhang: ")).td(br).
		   td(b(font({-color=>'red'}, $self->{'_overhang'}))));
		   
    }
    if ($self->{'_siteLength'}) {

	my $siteLength = $self->{'_siteLength'};

	$siteLength =~ s/^([<>]=?)([0-9]+)/$1 $2/;
	
	$rows .= Tr(td({-align=>'right'}, b("Length of recognition site: ")).
		    td(br).
		    td(b(font({-color=>'red'}, $siteLength))));

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

	my $cutNumber = $self->{'_cutNumber'};

	$cutNumber =~ s/^([<>]=?)([0-9]+)/$1 $2/;

	$rows .= Tr(td({-align=>'right'}, b("Cut number: ")).td(br).
		    td(b(font({-color=>'red'}, $cutNumber))));
	
    }
    if (!$rows) { 
	
	$rows = Tr(td({-align=>'right'}, b("Restriction Enzyme: ")).
		   td(b(font({-color=>'red'}, 'All'))));
    
    }
    print table({-align=>'center'}, $rows),p;


}

####################################################################
sub createAndShowMap {
####################################################################
    my ($self) = @_;
    
    my $restGraph = MAP::RestGraph->new(-seq=>$self->{'_sequence'},
                                        -display_id=>$self->{'_query'},
                                        -gifDir=>$configPath->tmpDir,
                                        -gifUrlRoot=>$configUrl->dictyBaseHtmlTmp,
                                        -gifFileName=>"restriction.$$.gif",
                                        -gifLabel=>'dictyBase',
					-overhang=>$self->{'_overhang'},
					-siteLength=>$self->{'_siteLength'},
                                        -cutNumber=>$self->{'_cutNumber'});


    if (!@{$restGraph->selectedEnzymeArrayRef}) {

	print center(b("No enzyme found matching your search criteria")),p;

    }
    else {

	$restGraph->showGraph;

    }
    return $restGraph;

}

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

    if (!@$notCutEnzymeArrayRef) { return; }

    my $desc;

    if ($self->{'_overhang'}) {
	
	$desc = $self->{'_overhang'}." enzymes ";
    
	$desc =~ s/^blunt/Blunt/;

    }
    else {

	$desc = "Enzymes ";

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

	my $siteLength;

	if ($self->{'_siteLength'} =~ /^[0-9]+/) {

	    $siteLength = "= ".$self->{'_siteLength'};

	}
	else {

	    $siteLength = $self->{'_siteLength'};

	    $siteLength =~ s/^([><]=?)([0-9]+)$/$1 $2/;

	}
	$desc .= "with recognition site ".$siteLength." ";

    }

    print p, $desc." 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 setVariables {
########################################################################
    my ($self) = @_;

    $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';
	}
    }

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

    ######## for graph
    
    if (param('overhang') !~ /^all/i) {

	$self->{'_overhang'} = param('overhang');

    }
    
    if (param('siteLength') !~ /^all/i) {

	my $sign;

	if (param('siteLengthSign') ne "=") {
	    
	    $sign = param('siteLengthSign');

	}

	$self->{'_siteLength'} = $sign.param('siteLength');

	if ($self->{'_siteLength'} eq "<4") {

	    $self->err_report("The minimum length of recognition site is 4. Please choose another number and try again.");

	    exit;

	}

    }
    
    if (param('cutNumber') !~ /^all/i) {

	my $sign;

	if (param('cutNumberSign') ne "=") {
	    
	    $sign = param('cutNumberSign');

	}
	$self->{'_cutNumber'} = $sign.param('cutNumber');

	if ($self->{'_cutNumber'} eq "<1") {

	    $self->err_report("The minimum cut number is 1. Please choose another number and try again.");

	    exit;

	}

    }
    

}

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