#!/usr/bin/perl
package sdevCreateSeqPage;

#######################################################################
##### Author : Mira Kaloper	
##### Date   : May 13 2002
##### Description : This package contains all necessary methods for 
#####               creating and displaing sequences . 
#####              
#######################################################################
use strict;

# use Bio::SeqIO;
use Bio::Seq;
use DBI;
use CGI qw/:all/;

use lib "/usr/local/dicty/www_dictybase/db/lib/common";
use Login qw (ConnectToDatabase);
use lib "/usr/local/dicty/www_dictybase/db/lib/dictyBase";
use dictyBaseCentralMod qw(:formatPage :getInfo);

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

use lib "/usr/local/dicty/www_dictybase/db/lib/dictyBase/Objects";
use Feature;
use ConfigURLdictyBase;
use Chromosome;
use Sequence;
use Subfeature;

$| = 1;

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

my $dbh;
my $dblink; 
my $configUrl;
my $locusObj;
my $featObj;


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

	$self = {};

	bless $self;

        $self->{'_database'} = $args{'database'}; 
      	#$self->{'_dbh'} = $args{'dbh'};
	$self->{'_help'} = $args{'help'};
	$self->{'_title'} = $args{'title'};
        $self->{'_feature'} = $args{'feature'};
        $self->{'_user'} = $args{'user'};
        $self->{'_seqtype'} = $args{'seqtype'};
        $self->{'_new_start'} = $args{'start'};
        $self->{'_new_stop'} = $args{'stop'};
        $self->{'_map'} = $args{'map'};

        #mira
        $self->{'_subFeatureObj'} = $args{'subFeatureObj'}; 

        #$self->{'_subFeatureObj'}->test;
        #my @array = $self->{'_subFeatureObj'}->getExonsArray;
        #print join(' ', @array), br;

        $dblink = ($self->{'_database'} =~ m/dev/i) ? "dictyBaseDEV" : "dictyBase"; 

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


        $self->setVariables;
        $self->setSubFeature;

        #create sequences        
        $self->createGenomicDNASequence;
        $self->createGenomic1KSequence; 
        $self->createCodingDNASequence;
        $self->createProteinSequence;

    	return $self;
}

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

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

}

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

    	my ($self) = @_;

        $configUrl = ConfigURLdictyBase->new;

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


        if ($self->{'_seqtype'} eq 'genomicDNA') {
           $self->printGenomicDNA;
           exit;
        }

        if ($self->{'_seqtype'} eq 'genomic1KDNA') {
           $self->printGenomic1K;
           exit;
        }


       if ($self->{'_seqtype'} eq 'codingDNA') {
           $self->printCodingDNA;
           exit;
        }


       if ($self->{'_seqtype'} eq 'protein') {
           $self->printProtein;
           exit;
        }


       if ($self->{'_map'}) {
           $self->printGCGGraphics;
           exit;
        }

        &printEndPage;

}


#######################################################################
sub setSubFeature {
#######################################################################
  my $self = shift;

  my $feature_no = $self->{'_feature_no'};

  my %subfeatureHash = Subfeature->GetStartStopTypeHash(dbh=>$dbh,
                                                     feature_no=>$feature_no);

  my @array = keys (%subfeatureHash);

  my $hashSize = $#array + 1;

  my ($startDiff, $stopDiff);

  $startDiff = $self->{'_start'} - $self->{'_new_start'}
                        if ($self->{'_new_start'});

  $stopDiff = $self->{'_new_stop'} - $self->{'_stop'}
                        if ($self->{'_new_stop'});

  my $crick = 1;
  if ($self->{'_strand'} eq 'C') {
     $startDiff = $startDiff * (-1);
     $stopDiff = $stopDiff * (-1);
     $crick = -1;
  }


  if (param('Exon1')) {
    $self->getSubFeatureFromParam;
  } elsif ($self->{'_subFeatureObj'}) {
    $self->getSubFeatureFromObject;
  } else {
    #this will not execute - it will get values from param instead
    $self->getSubFeatureFromDatabase;
  }

 
}


#######################################################################
sub getSubFeatureFromObject {
#######################################################################
  my $self = shift;

  my @array = $self->{'_subFeatureObj'}->getExonsArray;

  foreach my $element (@array) {

  my ($heading, $start, $stop, $subfeature_no) = split(/:/,$element);

      $self->{"$heading"} = "$start:$stop:$subfeature_no";

  }

  @array = $self->{'_subFeatureObj'}->getIntronsArray;

  foreach my $element (@array) {

  my ($heading, $start, $stop, $subfeature_no) = split(/:/,$element);

      $self->{"$heading"} = "$start:$stop:$subfeature_no";

  }


}

#######################################################################
sub getSubFeatureFromDatabase {
#######################################################################
  my $self = shift;

  my $feature_no = $self->{'_feature_no'};

  my @subfeatureArray = Subfeature->GetStartStopArray(dbh=>$dbh,
                                                     feature_no=>$feature_no);

  my $arraySize = $#subfeatureArray + 1;
  my $heading;
  my $i = 1;

  foreach my $element (@subfeatureArray) {
    my ($start, $stop, $type, $subfeature_no ) = split(/\|/,$element);

    if ($type eq "Exon") {
       $heading = $type . $i;
    } else {
       $heading = $type . $i++;
    }


    #set exon/intron info
    $self->{"_$heading"} = "$start:$stop:$subfeature_no";
  }




}



######################################################################
sub getSubFeatureFromParam {
###################################################################### 
  my $self = shift;

  my $i = 1;
  while (param("Exon$i")) {
    #set exon/intron info
    $self->{"_Exon$i"} = param("Exon$i");
    $i++;
  }

  my $i = 1;
  while (param("Intron$i")) {
    #set exon/intron info
    $self->{"_Intron$i"} = param("Intron$i");
    $i++;
  }

}


#######################################################################
sub getChrSequence {
#######################################################################
  my $self = shift;

  my $chromObject;

  eval {
     $chromObject =  Chromosome->new(dbh=>$dbh,
                                     chromosome=>$self->{'_chr'});
  };

  if ($@){
     print "Could not create chromosome object:  $@\n" ;
     $dbh->disconnect;
     exit;
  }


  if (!$chromObject) {
     print "Could not create object for ", $self->{'_chr'}, br;
     return;
  }

  my $sequence = $chromObject->chr_seq;

  if (!$sequence) {
     print "Sequence was not found in the database";
     exit;
  }

  $self->{'_sequence'} = $sequence; 

}



#######################################################################
sub createGenomicDNASequence {
#######################################################################
  my $self = shift;


  my $start = $self->{'_new_start'};
  my $stop =  $self->{'_new_stop'};

  my $rev;
  if ($self->{'_strand'} eq 'C') {
     ($start,$stop) = ($stop,$start);
     $rev = 1;
  }

  my $sequenceObj = Sequence->new(dbh=>$dbh,
                             chrnum=>$self->{'_chr'},
                             beg=>$start,
                             end=>$stop,  
                             rev=>$rev);


  $self->{'_genomicDNA'} = $sequenceObj->getSequence;

}

#######################################################################
sub getGenomicDNASequence {
#######################################################################
  my $self = shift;

  return $self->{'_genomicDNA'};


}


#######################################################################
sub createGenomic1KSequence {
#######################################################################
  my $self = shift;


  my $start = $self->{'_new_start'};
  my $stop =  $self->{'_new_stop'};

  my $rev;
  if ($self->{'_strand'} eq 'C') {
     ($start,$stop) = ($stop,$start);
     $rev = 1;
  }

  $start = $start - 1000;
  $stop = $stop + 1000;

  if ($start <=0) {
     $start = 1;
  }


  my $chrObj = Chromosome->new(dbh=>$dbh,
                               chromosome=>$self->{'_chr'});

  if (!$chrObj) {
     print "Could not create Chromosome object for : ", $self->{'_chr'}, br;
     exit;
  }


  my $len = $chrObj->physical_length;

  if ($stop > $len ) {
     $stop = $len;
  }


  my $sequenceObj = Sequence->new(dbh=>$dbh,
                             chrnum=>$self->{'_chr'},
                             beg=>$start,
                             end=>$stop,
                             rev=>$rev);


  $self->{'_genomic1KDNA'} = $sequenceObj->getSequence;

}


#######################################################################
sub getGenomic1KDNASequence {
#######################################################################
  my $self = shift;

  return $self->{'_genomic1KDNA'};


}

#######################################################################
sub createCodingDNASequence {
#######################################################################
  my $self = shift;

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

  my $sequence = $self->{'_genomicDNA'};


  my $new_start = $self->{'_new_start'};
  my $new_stop = $self->{'_new_stop'};

  if ($self->{'_strand'} eq 'C')  {
     ($new_start, $new_stop) = ($new_stop, $new_start);
  }


  my $seqObj = Bio::Seq->new('-seq'=>$sequence,
                            '-desc'=>'',
                            '-display_id'=>$self->{'_feature'},
                            '-accession_number'=>'none',
                            '-alphabet'=>'');

  #print "I am after Bio sequence", br;

  #remove introns from sequence....
  my $new_seq;
  my $beg = 1;
  my $i = 1;
  my $heading = '_Intron' . $i;

  while ($self->{$heading}) {
    my ($start,$stop,$no) = split(/:/,$self->{$heading});

    #print $start, br;
    #print $stop, br;

    my $newObj = $seqObj->trunc($beg,$start-1);

    $beg = $stop+1;
    #$beg = $new_start - $self->{'_start'};

    $new_seq .= $newObj->seq();

    $i++;
    $heading = '_Intron' . $i;
  }

  my $newObj = $seqObj->trunc($beg,$new_stop-$new_start+1);
 
  $new_seq .= $newObj->seq();

  $self->{'_codingDNA'} = $new_seq;

 
}


#######################################################################
sub getCodingDNASequence {
#######################################################################
  my $self = shift;

  return $self->{'_codingDNA'};

}


#######################################################################
sub createProteinSequence {
#######################################################################

  my $self = shift;

  my $sequence = $self->{'_codingDNA'};


  my $sequenceObj = Sequence->new(sequence=>$sequence,
                                  trans=>'1');



  # protein translation...
  $sequence = $sequenceObj->getSequence;
  $self->{'_protein'} = $sequence;

}


#######################################################################
sub getProteinSequence {
#######################################################################

  my $self = shift;

  return $self->{'_protein'};

}


#######################################################################
sub printGenomicDNA {
#######################################################################
  my $self = shift;

  my $sequence = $self->{'_genomicDNA'};
  
  my $seqObj = Sequence->new(sequence=>$sequence);
 
  $sequence = $seqObj->gcgFormatedSequence;

  print br, h4("Genomic DNA of Region:");

  print "<pre>";

  print $sequence;

  print "</pre>";

  
}


########################################################################
sub printCodingDNA {
########################################################################
  my $self = shift;

  my $sequence = $self->{'_codingDNA'};
  
  my $seqObj = Sequence->new(sequence=>$sequence);

  $sequence = $seqObj->gcgFormatedSequence;

  print br, h4("Coding DNA:");

  print "<pre>";

  print $sequence;

  print "</pre>";

}


########################################################################
sub printProtein {
########################################################################

  my $self = shift;


  my $sequence = $self->{'_protein'};
  
  my $seqObj = Sequence->new(sequence=>$sequence);

  $sequence = $seqObj->gcgFormatedSequence;

  print br, h4("Protein Translation:");

  print "<pre>";

  print $sequence;

  print "</pre>";

}



########################################################################
sub printGenomic1K {
########################################################################
  my $self = shift;


  my $sequence = $self->{'_genomic1KDNA'};

  my $seqObj = Sequence->new(sequence=>$sequence);

  $sequence = $seqObj->gcgFormatedSequence;

  print hr;
  print h4("Genomic +/- 1K"), br;
  print br;


  print "<pre>";

  print $sequence;

  print "</pre>";

  print hr;

}


#####################################################################
sub checkSequenceChange {
#####################################################################
  my $self = shift;

  my $error;

  #check if first 3 characters are "ATG"
  my $startDNA = substr($self->{'_codingDNA'},0,3);


  if ($startDNA ne "ATG") {
     $error++;
     print font({-size=>"+1",-color=>"red"},"The coding DNA sequence doesn't start with ATG!\n"), br;
  }


  #check if last 3 characters are: "TAA" or "TGA" or "TAG"
  my $stopDNA = substr($self->{'_codingDNA'},-3,3);


  if (($stopDNA ne "TAA") && ($stopDNA ne "TGA") && ($stopDNA ne "TAG")) {
     $error++;
     print font({-size=>"+1",-color=>"red"},"The coding sequence must end with TAA or TGA or TAG!\n"), br;
  }


  return $error;


}




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

    my $feature = $self->{'_feature'};

    &DeleteUnwantedChar(\$feature);

    if ($feature) {
        
        my $featureObj = Feature->new(dbh=>$dbh, feature_name=>$feature); 


	if (!$featureObj) {
	    $self->err_report("There is no feature associated with $feature in database.");
            exit;
	}
   
        $self->{'_feature_no'} = $featureObj->feature_no; 
        $self->{'_chrnum'} = $featureObj->chromosome;
        $self->{'_beg'} = $featureObj->start_coord;
        $self->{'_end'} = $featureObj->stop_coord;
        $self->{'_locus_name'} = $featureObj->locus_name; 
        $self->{'_chr'} = $featureObj->chromosome;
        $self->{'_start'} = $featureObj->start_coord;
        $self->{'_stop'} = $featureObj->stop_coord;
        $self->{'_strand'} = $featureObj->strand;
    }

}


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

    my $rev;
    my $start = $self->{'_new_start'};
    my $stop = $self->{'_new_stop'};


    ### o = open reading frames; n = no translation
    $self->{'_men'} = "-MEN=o";

    if ($self->{'_strand'} eq 'C') {
       $rev = '-REV'; 
       ($start,$stop) = ($stop,$start); 
    }

    my $gcg = GCG->new(map=>$self->{'_map'},
                       seqname=>$self->{'_feature'},
                       chrnum=>$self->{'_chrnum'},
                       beg=>$start,
                       end=>$stop,
                       rev=>$rev,
                       men=>$self->{'_men'});
  
    my $seqfile = $gcg->getSequence;

    print "<pre>";

    system("/usr/bin/cat $seqfile");

    print "</pre>";


    print "<hr> new call <BR> ";
   


    my $sequence = $self->{'_codingDNA'};

    print $sequence, "<BR>";

    my $gcgGenomic = GCG->new(map=>$self->{'_map'},
                              #seqname=>$self->{'_feature'},
                              sequence=>$sequence,
                              men=>$self->{'_men'});
 
    my $seqfileNew = $gcgGenomic->getSequence;
  
    print "<BR> $seqfileNew <br>";

    print "<pre>";

    system("/usr/bin/cat $seqfileNew");

    print "</pre>";


}



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

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

    $dbh->disconnect;
    exit;
}


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