[Bio] / FigKernelPackages / Observation.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/Observation.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2, Tue Jun 12 16:51:53 2007 UTC revision 1.25, Wed Jul 11 15:17:49 2007 UTC
# Line 1  Line 1 
1  package Observation;  package Observation;
2    
3    use lib '/vol/ontologies';
4    use DBMaster;
5    
6  require Exporter;  require Exporter;
7  @EXPORT_OK = qw(get_objects);  @EXPORT_OK = qw(get_objects);
8    
9    use FIG_Config;
10  use strict;  use strict;
11  use warnings;  #use warnings;
12    use HTML;
13    
14  1;  1;
15    
# Line 21  Line 26 
26    
27  The data can be used to display information for a given sequence feature (protein or other, but primarily information is computed for proteins).  The data can be used to display information for a given sequence feature (protein or other, but primarily information is computed for proteins).
28    
 Example:  
   
 use FIG;  
 use Observation;  
   
 my $fig = new FIG;  
 my $fid = "fig|83333.1.peg.3";  
   
 my $observations = Observation::get_objects($fid);  
 foreach my $observation (@$observations) {  
     print "ID: " . $fid . "\n";  
     print "Start: " . $observation->start() . "\n";  
     ...  
 }  
   
 B<return an array of objects>  
   
   
 print "$Observation->acc\n" prints the Accession number if present for the Observation  
   
29  =cut  =cut
30    
31  =head1 BACKGROUND  =head1 BACKGROUND
# Line 64  Line 49 
49    
50  The public methods this package provides are listed below:  The public methods this package provides are listed below:
51    
 =head3 acc()  
52    
53  A valid accession or remote ID (in the style of a db_xref) or a valid local ID (FID) in case this is supported.  =head3 context()
54    
55    Returns close or diverse for purposes of displaying genomic context
56    
57  =cut  =cut
58    
59  sub acc {  sub context {
60    my ($self) = @_;    my ($self) = @_;
61    
62    return $self->{acc};    return $self->{context};
63  }  }
64    
65  =head3 description()  =head3 rows()
   
 The description of the hit. Taken from the data or from the our Ontology database for some cases e.g. IPR or PFAM.  
66    
67  B<Please note:>  each row in a displayed table
 Either remoteid or description is required.  
68    
69  =cut  =cut
70    
71  sub description {  sub rows {
72    my ($self) = @_;    my ($self) = @_;
73    
74      return $self->{rows};
75    }
76    
77    =head3 acc()
78    
79    A valid accession or remote ID (in the style of a db_xref) or a valid local ID (FID) in case this is supported.
80    
81    =cut
82    
83    sub acc {
84      my ($self) = @_;
85    return $self->{acc};    return $self->{acc};
86  }  }
87    
# Line 100  Line 94 
94    
95  =over 9  =over 9
96    
97  =item sim (seq)  =item IDENTICAL (seq)
98    
99  =item bbh (seq)  =item SIM (seq)
100    
101  =item pch (fc)  =item BBH (seq)
102    
103  =item figfam (seq)  =item PCH (fc)
104    
105  =item ipr (dom)  =item FIGFAM (seq)
106    
107  =item cdd (dom)  =item IPR (dom)
108    
109  =item pfam (dom)  =item CDD (dom)
110    
111  =item signalp (dom)  =item PFAM (dom)
112    
113  =item cello (loc)  =item SIGNALP_CELLO_TMPRED (loc)
114    
115  =item tmhmm (loc)  =item PDB (seq)
116    
117  =item hmmtop (loc)  =item TMHMM (loc)
118    
119    =item HMMTOP (loc)
120    
121  =back  =back
122    
# Line 182  Line 178 
178    return $self->{stop};    return $self->{stop};
179  }  }
180    
181  =head3 evalue()  =head3 start()
182    
183  E-value or P-Value if present.  Start of hit in query sequence.
184    
185  =cut  =cut
186    
187  sub evalue {  sub qstart {
188    my ($self) = @_;    my ($self) = @_;
189    
190    return $self->{evalue};      return $self->{qstart};
191  }  }
192    
193  =head3 score()  =head3 qstop()
   
 Score if present.  
194    
195  B<Please note: >  End of the hit in query sequence.
 Either score or eval are required.  
196    
197  =cut  =cut
198    
199  sub score {  sub qstop {
200    my ($self) = @_;    my ($self) = @_;
201    
202    return $self->{score};      return $self->{qstop};
203  }  }
204    
205    =head3 hstart()
206    
207  =head3 display_method()  Start of hit in hit sequence.
208    
209  If available use the function specified here to display the "raw" observation.  =cut
 In the case of a BLAST alignment of fid1 and fid2 a cgi script  
 will be called to display the results of running the command "bl2seq fid1 fid2".  
210    
211  B<Please note> that URL linked to in display_method() is an external component and needs to added to the code for every class of evidence.  sub hstart {
212        my ($self) = @_;
213    
214        return $self->{hstart};
215    }
216    
217    =head3 end()
218    
219    End of the hit in hit sequence.
220    
221  =cut  =cut
222    
223  sub display_method {  sub hstop {
224    my ($self) = @_;    my ($self) = @_;
225    
226    # add code here      return $self->{hstop};
   
   return $self->{display_method};  
227  }  }
228    
229  =head3 rank()  =head3 qlength()
   
 Returns an integer from 1 - 10 indicating the importance of this observations.  
230    
231  Currently always returns 1.  length of the query sequence in similarities
232    
233  =cut  =cut
234    
235  sub rank {  sub qlength {
236    my ($self) = @_;    my ($self) = @_;
237    
238  #  return $self->{rank};      return $self->{qlength};
   
   return 1;  
239  }  }
240    
241  =head3 supports_annotation()  =head3 hlength()
   
 Does a this observation support the annotation of its feature?  
242    
243  Returns  length of the hit sequence in similarities
244    
245  =over 3  =cut
246    
247  =item 10, if feature annotation is identical to $self->description  sub hlength {
248        my ($self) = @_;
249    
250  =item 1, Feature annotation is similar to $self->annotation; this is computed using FIG::SameFunc()      return $self->{hlength};
251    }
252    
253  =item undef  =head3 evalue()
254    
255  =back  E-value or P-Value if present.
256    
257  =cut  =cut
258    
259  sub supports_annotation {  sub evalue {
260    my ($self) = @_;    my ($self) = @_;
261    
262    # no code here so far    return $self->{evalue};
   
   return $self->{supports_annotation};  
263  }  }
264    
265  =head3 url()  =head3 score()
266    
267  URL describing the subject. In case of a BLAST hit against a sequence, this URL will lead to a page displaying the sequence record for the sequence. In case of an HMM hit, the URL will be to the URL description.  Score if present.
268    
269  =cut  =cut
270    
271  sub url {  sub score {
272    my ($self) = @_;    my ($self) = @_;
273      return $self->{score};
274    }
275    
276    =head3 display()
277    
278    my $url = get_url($self->type, $self->acc);  will be different for each type
279    
280    =cut
281    
282    sub display {
283    
284      die "Abstract Method Called\n";
285    
   return $url;  
286  }  }
287    
288  =head3 get_objects()  =head3 display_table()
289    
290  This is the B<REAL WORKHORSE> method of this Package.  will be different for each type
291    
292    =cut
293    
294    sub display_table {
295    
296  It will probably have to:    die "Abstract Table Method Called\n";
297    
 - get all sims for the feature  
 - get all bbhs for the feature  
 - copy information from sim to bbh (bbh have no match location etc)  
 - get pchs (difficult)  
 - get attributes (there is code for this that in get_attribute_based_observations  
 - get_attributes_based_observations returns an array of arrays of hashes like this"  
   
   my $datasets =  
      [  
        [ { name => 'acc', value => '1234' },  
         { name => 'from', value => '4' },  
         { name => 'to', value => '400' },  
         ....  
        ],  
        [ { name => 'acc', value => '456' },  
         { name => 'from', value => '1' },  
         { name => 'to', value => '100' },  
         ....  
        ],  
        ...  
      ];  
    return $datasets;  
298   }   }
299    
300  It will invoke the required calls to the SEED API to retrieve the information required.  =head3 get_objects()
301    
302    This is the B<REAL WORKHORSE> method of this Package.
303    
304  =cut  =cut
305    
306  sub get_objects {  sub get_objects {
307      my ($self,$fid) = @_;      my ($self,$fid,$scope) = @_;
308    
309      my $objects = [];      my $objects = [];
310    my @matched_datasets=();    my @matched_datasets=();
311    
312    # call function that fetches attribut based observations      # call function that fetches attribute based observations
313    # returns an array of arrays of hashes    # returns an array of arrays of hashes
   #  
   get_attribute_based_observations($fid,\@matched_datasets);  
314    
315    # read sims + bbh (enrich BBHs with sims coordindates etc)      if($scope){
316    # read pchs          get_cluster_observations($fid,\@matched_datasets,$scope);
317    # read figfam match data from 48hr directory (BobO knows how do do this!)      }
318    # what sources of evidence did I miss?      else{
319            my %domain_classes;
320            $domain_classes{'CDD'} = 1;
321            get_identical_proteins($fid,\@matched_datasets);
322            get_attribute_based_domain_observations($fid,\%domain_classes,\@matched_datasets);
323            get_sims_observations($fid,\@matched_datasets);
324            get_functional_coupling($fid,\@matched_datasets);
325            get_attribute_based_location_observations($fid,\@matched_datasets);
326            get_pdb_observations($fid,\@matched_datasets);
327        }
328    
329    foreach my $dataset (@matched_datasets) {    foreach my $dataset (@matched_datasets) {
330      my $object = $self->new();          my $object;
331      foreach my $attribute (@$dataset) {          if($dataset->{'type'} eq "dom"){
332        $object->{$attribute->{'name'}} = $attribute->{'value'};              $object = Observation::Domain->new($dataset);
333      }      }
334  #    $object->{$attribute->{'feature_id'}} = $attribute->{$fid};          if($dataset->{'class'} eq "PCH"){
335      push (@$objects, $object);              $object = Observation::FC->new($dataset);
336            }
337            if ($dataset->{'class'} eq "IDENTICAL"){
338                $object = Observation::Identical->new($dataset);
339            }
340            if ($dataset->{'class'} eq "SIGNALP_CELLO_TMPRED"){
341                $object = Observation::Location->new($dataset);
342            }
343            if ($dataset->{'class'} eq "SIM"){
344                $object = Observation::Sims->new($dataset);
345            }
346            if ($dataset->{'class'} eq "CLUSTER"){
347                $object = Observation::Cluster->new($dataset);
348            }
349            if ($dataset->{'class'} eq "PDB"){
350                $object = Observation::PDB->new($dataset);
351      }      }
352    
353            push (@$objects, $object);
354        }
355    
356    return $objects;    return $objects;
357    
358  }  }
359    
360  =head1 Internal Methods  =head1 Internal Methods
# Line 355  Line 365 
365    
366  =cut  =cut
367    
368    sub get_attribute_based_domain_observations{
369    
370  =head3 get_url (internal)      # we read a FIG ID and a reference to an array (of arrays of hashes, see above)
371        my ($fid,$domain_classes,$datasets_ref) = (@_);
 get_url() return a valid URL or undef for any observation.  
372    
373  URLs are constructed by looking at the Accession acc()  and  name()      my $fig = new FIG;
374    
375  Info from both attributes is combined with a table of base URLs stored in this function.      foreach my $attr_ref ($fig->get_attributes($fid)) {
376            my $key = @$attr_ref[1];
377            my @parts = split("::",$key);
378            my $class = $parts[0];
379    
380            if($domain_classes->{$parts[0]}){
381                my $val = @$attr_ref[2];
382                if($val =~/^(\d+\.\d+|0\.0);(\d+)-(\d+)/){
383                    my $raw_evalue = $1;
384                    my $from = $2;
385                    my $to = $3;
386                    my $evalue;
387                    if($raw_evalue =~/(\d+)\.(\d+)/){
388                        my $part2 = 1000 - $1;
389                        my $part1 = $2/100;
390                        $evalue = $part1."e-".$part2;
391                    }
392                    else{
393                        $evalue = "0.0";
394                    }
395    
396                    my $dataset = {'class' => $class,
397                                   'acc' => $key,
398                                   'type' => "dom" ,
399                                   'evalue' => $evalue,
400                                   'start' => $from,
401                                   'stop' => $to,
402                                   'fig_id' => $fid,
403                                   'score' => $raw_evalue
404                                   };
405    
406  =cut                  push (@{$datasets_ref} ,$dataset);
407                }
408            }
409        }
410    }
411    
412  sub get_url {  sub get_attribute_based_location_observations{
413    
414   my ($self) = @_;      my ($fid,$datasets_ref) = (@_);
415   my $url='';      my $fig = new FIG;
416    
417  # a hash with a URL for each observation; identified by name()      my $location_attributes = ['SignalP','CELLO','TMPRED'];
 #my $URL             => { 'PFAM' => "http://www.sanger.ac.uk/cgi-bin/Pfam/getacc?" ,\  
 #                       'IPR'    => "http://www.ebi.ac.uk/interpro/DisplayIproEntry?ac=" ,\  
 #                          'CDD' => "http://www.ncbi.nlm.nih.gov/Structure/cdd/cddsrv.cgi?uid=",\  
 #                       'PIR'    => "http://www.ncbi.nlm.nih.gov/Structure/cdd/cddsrv.cgi?uid=",\  
 #                       'FIGFAM' => '',\  
 #                          'sim'=> "http://www.theseed.org/linkin.cgi?id=",\  
 #                          'bbh'=> "http://www.theseed.org/linkin.cgi?id="  
 #};  
418    
419  # if (defined $URL{$self->name}) {      my $dataset = {'type' => "loc", 'class' => 'SIGNALP_CELLO_TMPRED','fig_id' => $fid};
420  #     $url = $URL{$self->name}.$self->acc;      foreach my $attr_ref ($fig->get_attributes($fid,$location_attributes)) {
421  #     return $url;          my $key = @$attr_ref[1];
422  # }          my @parts = split("::",$key);
423  # else          my $sub_class = $parts[0];
424       return undef;          my $sub_key = $parts[1];
425            my $value = @$attr_ref[2];
426            if($sub_class eq "SignalP"){
427                if($sub_key eq "cleavage_site"){
428                    my @value_parts = split(";",$value);
429                    $dataset->{'cleavage_prob'} = $value_parts[0];
430                    $dataset->{'cleavage_loc'} = $value_parts[1];
431                }
432                elsif($sub_key eq "signal_peptide"){
433                    $dataset->{'signal_peptide_score'} = $value;
434                }
435            }
436            elsif($sub_class eq "CELLO"){
437                $dataset->{'cello_location'} = $sub_key;
438                $dataset->{'cello_score'} = $value;
439            }
440            elsif($sub_class eq "TMPRED"){
441                my @value_parts = split(";",$value);
442                $dataset->{'tmpred_score'} = $value_parts[0];
443                $dataset->{'tmpred_locations'} = $value_parts[1];
444  }  }
445        }
446    
447        push (@{$datasets_ref} ,$dataset);
448    
449  =head3 get_display_method (internal)  }
450    
451  get_display_method() return a valid URL or undef for any observation.  =head3 get_pdb_observations() (internal)
452    
453  URLs are constructed by looking at the Accession acc()  and  name()  This methods sets the type and class for pdb observations
 and Info from both attributes is combined with a table of base URLs stored in this function.  
454    
455  =cut  =cut
456    
457  sub get_display_method {  sub get_pdb_observations{
458        my ($fid,$datasets_ref) = (@_);
459    
460   my ($self) = @_;      my $fig = new FIG;
461    
462  # a hash with a URL for each observation; identified by name()      foreach my $attr_ref ($fig->get_attributes($fid,'PDB')) {
 #my $URL             => { 'sim'=> "http://www.theseed.org/featalign.cgi?id1=",\  
 #                        'bbh'=> "http://www.theseed.org/featalign.cgi?id1="  
 # };  
463    
464  #if (defined $URL{$self->name}) {          my $key = @$attr_ref[1];
465  #     $url = $URL{$self->name}.$self->feature_id."&id2=".$self->acc;          my($key1,$key2) =split("::",$key);
466  #     return $url;          my $value = @$attr_ref[2];
467  # }          my ($evalue,$location) = split(";",$value);
468  # else  
469       return undef;          if($evalue =~/(\d+)\.(\d+)/){
470                my $part2 = 1000 - $1;
471                my $part1 = $2/100;
472                $evalue = $part1."e-".$part2;
473            }
474    
475            my($start,$stop) =split("-",$location);
476    
477            my $url = @$attr_ref[3];
478            my $dataset = {'class' => 'PDB',
479                           'type' => 'seq' ,
480                           'acc' => $key2,
481                           'evalue' => $evalue,
482                           'start' => $start,
483                           'stop' => $stop,
484                           'fig_id' => $fid
485                           };
486    
487            push (@{$datasets_ref} ,$dataset);
488        }
489  }  }
490    
491  =head3 get_attribute_based_evidence (internal)  =head3 get_cluster_observations() (internal)
492    
493  This method retrieves evidence from the attribute server  This methods sets the type and class for cluster observations
494    
495  =cut  =cut
496    
497  sub get_attribute_based_observations{  sub get_cluster_observations{
498        my ($fid,$datasets_ref,$scope) = (@_);
499    
500        my $dataset = {'class' => 'CLUSTER',
501                       'type' => 'fc',
502                       'context' => $scope,
503                       'fig_id' => $fid
504                       };
505        push (@{$datasets_ref} ,$dataset);
506    }
507    
     # we read a FIG ID and a reference to an array (of arrays of hashes, see above)  
     my ($fid,$datasets_ref) = (@_);  
508    
509      my $_myfig = new FIG;  =head3 get_sims_observations() (internal)
510    
511      foreach my $attr_ref ($_myfig->get_attributes($fid)) {  This methods retrieves sims fills the internal data structures.
512    
513          # convert the ref into a string for easier handling  =cut
         my ($string) = "@$attr_ref";  
514    
515  #       print "S:$string\n";  sub get_sims_observations{
         my ($key,$val) = ( $string =~ /\S+\s(\S+)\s(\S+)/);  
516    
517          # THIS SHOULD BE DONE ANOTHER WAY FM->TD      my ($fid,$datasets_ref) = (@_);
518          # we need to do the right thing for each type, ie no evalue for CELLO and no coordinates, but a score, etc      my $fig = new FIG;
519          # as fas as possible this should be configured so that the type of observation and the regexp are      my @sims= $fig->nsims($fid,100,1e-20,"all");
520          # stored somewhere for easy expansion      my ($dataset);
521          #      foreach my $sim (@sims){
522            my $hit = $sim->[1];
523            my $percent = $sim->[2];
524            my $evalue = $sim->[10];
525            my $qfrom = $sim->[6];
526            my $qto = $sim->[7];
527            my $hfrom = $sim->[8];
528            my $hto = $sim->[9];
529            my $qlength = $sim->[12];
530            my $hlength = $sim->[13];
531            my $db = get_database($hit);
532            my $func = $fig->function_of($hit);
533            my $organism = $fig->org_of($hit);
534    
535            $dataset = {'class' => 'SIM',
536                        'acc' => $hit,
537                        'identity' => $percent,
538                        'type' => 'seq',
539                        'evalue' => $evalue,
540                        'qstart' => $qfrom,
541                        'qstop' => $qto,
542                        'hstart' => $hfrom,
543                        'hstop' => $hto,
544                        'database' => $db,
545                        'organism' => $organism,
546                        'function' => $func,
547                        'qlength' => $qlength,
548                        'hlength' => $hlength,
549                        'fig_id' => $fid
550                        };
551    
552            push (@{$datasets_ref} ,$dataset);
553        }
554    }
555    
556    =head3 get_database (internal)
557    This method gets the database association from the sequence id
558    
559    =cut
560    
561          if (($key =~ /PFAM::/) || ( $key =~ /IPR::/) || ( $key =~ /CDD::/) ) {  sub get_database{
562        my ($id) = (@_);
563    
564              # some keys are composite CDD::1233244 or PFAM:PF1233      my ($db);
565        if ($id =~ /^fig\|/)              { $db = "FIG" }
566        elsif ($id =~ /^gi\|/)            { $db = "NCBI" }
567        elsif ($id =~ /^^[NXYZA]P_/)      { $db = "RefSeq" }
568        elsif ($id =~ /^sp\|/)            { $db = "SwissProt" }
569        elsif ($id =~ /^uni\|/)           { $db = "UniProt" }
570        elsif ($id =~ /^tigr\|/)          { $db = "TIGR" }
571        elsif ($id =~ /^pir\|/)           { $db = "PIR" }
572        elsif ($id =~ /^kegg\|/)          { $db = "KEGG" }
573        elsif ($id =~ /^tr\|/)            { $db = "TrEMBL" }
574        elsif ($id =~ /^eric\|/)          { $db = "ASAP" }
575        elsif ($id =~ /^img\|/)           { $db = "JGI" }
576    
577        return ($db);
578    
             if ( $key =~ /::/ ) {  
                 my ($firstkey,$restkey) = ( $key =~ /([a-zA-Z0-9]+)::(.*)/);  
                 $val=$restkey.";".$val;  
                 $key=$firstkey;  
579              }              }
580    
             my ($acc,$raw_evalue, $from,$to) = ($val =~ /(\S+);(\S+);(\d+)-(\d+)/ );  
581    
582              my $evalue= 255;  =head3 get_identical_proteins() (internal)
             if (defined $raw_evalue) { # some of the tool do not give us an evalue  
583    
584                  my ($k,$expo) = ( $raw_evalue =~ /(\d+).(\d+)/);  This methods retrieves sims fills the internal data structures.
                 my ($new_k, $new_exp);  
585    
586                  #  =cut
587                  #  THIS DOES NOT WORK PROPERLY  
588                  #  sub get_identical_proteins{
589                  if($raw_evalue =~/(\d+).(\d+)/){  
590        my ($fid,$datasets_ref) = (@_);
591        my $fig = new FIG;
592        my $funcs_ref;
593    
594  #                   $new_exp = (1000+$expo);      my @maps_to = grep { $_ ne $fid and $_ !~ /^xxx/ } map { $_->[0] } $fig->mapped_prot_ids($fid);
         #           $new_k = $k / 100;  
595    
596        foreach my $id (@maps_to) {
597            my ($tmp, $who);
598            if (($id ne $fid) && ($tmp = $fig->function_of($id))) {
599                $who = &get_database($id);
600                push(@$funcs_ref, [$id,$who,$tmp]);
601                  }                  }
                 $evalue = "0.01"#new_k."e-".$new_exp;  
602              }              }
603    
604              # unroll it all into an array of hashes      my ($dataset);
605              # this needs to be done differently for different types of observations      my $dataset = {'class' => 'IDENTICAL',
606              my $dataset = [ { name => 'class', value => $key },                     'type' => 'seq',
607                              { name => 'acc' , value => $acc},                     'fig_id' => $fid,
608                              { name => 'type', value => "dom"} , # this clearly needs to be done properly FM->TD                     'rows' => $funcs_ref
609                              { name => 'evalue', value => $evalue },                     };
                             { name => 'start', value => $from},  
                             { name => 'stop' , value => $to}  
                             ];  
610    
611              push (@{$datasets_ref} ,$dataset);              push (@{$datasets_ref} ,$dataset);
612          }  
613      }  
614  }  }
615    
616  =head3 get_sims_and_bbhs() (internal)  =head3 get_functional_coupling() (internal)
617    
618  This methods retrieves sims and also BBHs and fills the internal data structures.  This methods retrieves the functional coupling of a protein given a peg ID
619    
620  =cut  =cut
621    
622  #     sub get_sims_and_bbhs{  sub get_functional_coupling{
623    
624  #       # blast m8 output format      my ($fid,$datasets_ref) = (@_);
625  #       # id1, id2, %ident, align len, mismatches, gaps, q.start, q.stop, s. start, s.stop, eval, bit      my $fig = new FIG;
626        my @funcs = ();
 #       my $Sims=();  
 #       @sims_src = $fig->sims($fid,80,500,"fig",0);  
 #       print "found $#sims_src SIMs\n";  
 #       foreach $sims (@sims_src) {  
 #           my ($sims_string) = "@$sims";  
 # #       print "$sims_string\n";  
 #           my ($rfid,$start,$stop,$eval) = ( $sims_string =~ /\S+\s+(\S+)\s+\S+\s\S+\s+(\S+)\s+(\S+)\s+  
 #                                             \S+\s+\S+\s+\S+\s+\S+\s+(\S+)+.*/);  
 # #       print "ID: $rfid, E:$eval, Start:$start stop:$stop\n";  
 #           $Sims{$rfid}{'eval'}=$eval;  
 #           $Sims{$rfid}{'start'}=$start;  
 #           $Sims{$rfid}{'stop'}=$stop;  
 #           print "$rfid $Sims{$rfid}{'eval'}\n";  
 #       }  
   
 #       # BBHs  
 #       my $BBHs=();  
   
 #       @bbhs_src = $fig->bbhs($fid,1.0e-10);  
 #       print "found $#bbhs_src BBHs\n";  
 #       foreach $bbh (@bbhs_src) {  
 #           #print "@$bbh\n";  
 #           my ($bbh_string) = "@$bbh";  
 #           my ($rfid,$eval,$score) = ( $bbh_string =~ /(\S+)\s(\S+)\s(\S+)/);  
 #           #print "ID: $rfid, E:$eval, S:$score\n";  
 #           $BBHs{$rfid}{'eval'}=$eval;  
 #           $BBHs{$rfid}{'score'}=$score;  
 # #print "$rfid $BBHs{$rfid}{'eval'}\n";  
 #       }  
627    
628  #     }      # initialize some variables
629        my($sc,$neigh);
630    
631        # set default parameters for coupling and evidence
632        my ($bound,$sim_cutoff,$coupling_cutoff) = (5000, 1.0e-10, 4);
633    
634        # get the fc data
635        my @fc_data = $fig->coupling_and_evidence($fid,$bound,$sim_cutoff,$coupling_cutoff,1);
636    
637        # retrieve data
638        my @rows = map { ($sc,$neigh) = @$_;
639                         [$sc,$neigh,scalar $fig->function_of($neigh)]
640                      } @fc_data;
641    
642        my ($dataset);
643        my $dataset = {'class' => 'PCH',
644                       'type' => 'fc',
645                       'fig_id' => $fid,
646                       'rows' => \@rows
647                       };
648    
649        push (@{$datasets_ref} ,$dataset);
650    
651    }
652    
653  =head3 new (internal)  =head3 new (internal)
654    
# Line 539  Line 657 
657  =cut  =cut
658    
659  sub new {  sub new {
660    my ($self) = @_;    my ($class,$dataset) = @_;
661    
662    $self = { acc => '',    my $self = { class => $dataset->{'class'},
663              description => '',                 type => $dataset->{'type'},
664              class => '',                 fig_id => $dataset->{'fig_id'},
665              type => '',                 score => $dataset->{'score'},
             start => '',  
             stop => '',  
             evalue => '',  
             score => '',  
             display_method => '',  
             feature_id => '',  
             rank => '',  
             supports_annotation => ''  
666            };            };
667    
668    bless($self, 'Observation');    bless($self,$class);
669    
670    return $self;    return $self;
671  }  }
672    
673    =head3 identity (internal)
674    
675    Returns the % identity of the similar sequence
676    
677    =cut
678    
679    sub identity {
680        my ($self) = @_;
681    
682        return $self->{identity};
683    }
684    
685    =head3 fig_id (internal)
686    
687    =cut
688    
689    sub fig_id {
690      my ($self) = @_;
691      return $self->{fig_id};
692    }
693    
694  =head3 feature_id (internal)  =head3 feature_id (internal)
695    
 Returns the ID  of the feature these Observations belong to.  
696    
697  =cut  =cut
698    
# Line 571  Line 701 
701    
702    return $self->{feature_id};    return $self->{feature_id};
703  }  }
704    
705    =head3 id (internal)
706    
707    Returns the ID  of the identical sequence
708    
709    =cut
710    
711    sub id {
712        my ($self) = @_;
713    
714        return $self->{id};
715    }
716    
717    =head3 organism (internal)
718    
719    Returns the organism  of the identical sequence
720    
721    =cut
722    
723    sub organism {
724        my ($self) = @_;
725    
726        return $self->{organism};
727    }
728    
729    =head3 function (internal)
730    
731    Returns the function of the identical sequence
732    
733    =cut
734    
735    sub function {
736        my ($self) = @_;
737    
738        return $self->{function};
739    }
740    
741    =head3 database (internal)
742    
743    Returns the database of the identical sequence
744    
745    =cut
746    
747    sub database {
748        my ($self) = @_;
749    
750        return $self->{database};
751    }
752    
753    sub score {
754      my ($self) = @_;
755    
756      return $self->{score};
757    }
758    
759    ############################################################
760    ############################################################
761    package Observation::PDB;
762    
763    use base qw(Observation);
764    
765    sub new {
766    
767        my ($class,$dataset) = @_;
768        my $self = $class->SUPER::new($dataset);
769        $self->{acc} = $dataset->{'acc'};
770        $self->{evalue} = $dataset->{'evalue'};
771        $self->{start} = $dataset->{'start'};
772        $self->{stop} = $dataset->{'stop'};
773        bless($self,$class);
774        return $self;
775    }
776    
777    =head3 display()
778    
779    displays data stored in best_PDB attribute and in Ontology server for given PDB id
780    
781    =cut
782    
783    sub display{
784        my ($self,$gd) = @_;
785    
786        my $fid = $self->fig_id;
787        my $dbmaster = DBMaster->new(-database =>'Ontology');
788    
789        my $acc = $self->acc;
790    
791        print STDERR "acc:$acc\n";
792        my ($pdb_description,$pdb_source,$pdb_ligand);
793        my $pdb_objs = $dbmaster->pdb->get_objects( { 'id' => $acc } );
794        if(!scalar(@$pdb_objs)){
795            $pdb_description = "not available";
796            $pdb_source = "not available";
797            $pdb_ligand = "not available";
798        }
799        else{
800            my $pdb_obj = $pdb_objs->[0];
801            $pdb_description = $pdb_obj->description;
802            $pdb_source = $pdb_obj->source;
803            $pdb_ligand = $pdb_obj->ligand;
804        }
805    
806        my $lines = [];
807        my $line_data = [];
808        my $line_config = { 'title' => "PDB hit for $fid",
809                            'short_title' => "best PDB",
810                            'basepair_offset' => '1' };
811    
812        my $fig = new FIG;
813        my $seq = $fig->get_translation($fid);
814        my $fid_stop = length($seq);
815    
816        my $fid_element_hash = {
817            "title" => $fid,
818            "start" => '1',
819            "end" =>  $fid_stop,
820            "color"=> '1',
821            "zlayer" => '1'
822            };
823    
824        push(@$line_data,$fid_element_hash);
825    
826        my $links_list = [];
827        my $descriptions = [];
828    
829        my $name;
830        $name = {"title" => 'id',
831                 "value" => $acc};
832        push(@$descriptions,$name);
833    
834        my $description;
835        $description = {"title" => 'pdb description',
836                        "value" => $pdb_description};
837        push(@$descriptions,$description);
838    
839        my $score;
840        $score = {"title" => "score",
841                  "value" => $self->evalue};
842        push(@$descriptions,$score);
843    
844        my $start_stop;
845        my $start_stop_value = $self->start."_".$self->stop;
846        $start_stop = {"title" => "start-stop",
847                       "value" => $start_stop_value};
848        push(@$descriptions,$start_stop);
849    
850        my $source;
851        $source = {"title" => "source",
852                  "value" => $pdb_source};
853        push(@$descriptions,$source);
854    
855        my $ligand;
856        $ligand = {"title" => "pdb ligand",
857                   "value" => $pdb_ligand};
858        push(@$descriptions,$ligand);
859    
860        my $link;
861        my $link_url ="http://www.rcsb.org/pdb/explore/explore.do?structureId=".$acc;
862    
863        $link = {"link_title" => $acc,
864                 "link" => $link_url};
865        push(@$links_list,$link);
866    
867        my $pdb_element_hash = {
868            "title" => "PDB homology",
869            "start" => $self->start,
870            "end" =>  $self->stop,
871            "color"=> '6',
872            "zlayer" => '3',
873            "links_list" => $links_list,
874            "description" => $descriptions};
875    
876        push(@$line_data,$pdb_element_hash);
877        $gd->add_line($line_data, $line_config);
878    
879        return $gd;
880    }
881    
882    1;
883    
884    ############################################################
885    ############################################################
886    package Observation::Identical;
887    
888    use base qw(Observation);
889    
890    sub new {
891    
892        my ($class,$dataset) = @_;
893        my $self = $class->SUPER::new($dataset);
894        $self->{rows} = $dataset->{'rows'};
895    
896        bless($self,$class);
897        return $self;
898    }
899    
900    =head3 display_table()
901    
902    If available use the function specified here to display the "raw" observation.
903    This code will display a table for the identical protein
904    
905    
906    B<Please note> that URL linked to in display_method() is an external component and needs to added to the code for every class of evi
907    dence.
908    
909    =cut
910    
911    
912    sub display_table{
913        my ($self) = @_;
914    
915        my $fig = new FIG;
916        my $fid = $self->fig_id;
917        my $rows = $self->rows;
918        my $cgi = new CGI;
919        my $all_domains = [];
920        my $count_identical = 0;
921        my $content;
922        foreach my $row (@$rows) {
923            my $id = $row->[0];
924            my $who = $row->[1];
925            my $assignment = $row->[2];
926            my $organism = $fig->org_of($fid);
927            my $single_domain = [];
928            push(@$single_domain,$who);
929            push(@$single_domain,&HTML::set_prot_links($cgi,$id));
930            push(@$single_domain,$organism);
931            push(@$single_domain,$assignment);
932            push(@$all_domains,$single_domain);
933            $count_identical++;
934        }
935    
936        if ($count_identical >0){
937            $content = $all_domains;
938        }
939        else{
940            $content = "<p>This PEG does not have any essentially identical proteins</p>";
941        }
942        return ($content);
943    }
944    
945    1;
946    
947    #########################################
948    #########################################
949    package Observation::FC;
950    1;
951    
952    use base qw(Observation);
953    
954    sub new {
955    
956        my ($class,$dataset) = @_;
957        my $self = $class->SUPER::new($dataset);
958        $self->{rows} = $dataset->{'rows'};
959    
960        bless($self,$class);
961        return $self;
962    }
963    
964    =head3 display_table()
965    
966    If available use the function specified here to display the "raw" observation.
967    This code will display a table for the identical protein
968    
969    
970    B<Please note> that URL linked to in display_method() is an external component and needs to added to the code for every class of evi
971    dence.
972    
973    =cut
974    
975    sub display_table {
976    
977        my ($self,$dataset) = @_;
978        my $fid = $self->fig_id;
979        my $rows = $self->rows;
980        my $cgi = new CGI;
981        my $functional_data = [];
982        my $count = 0;
983        my $content;
984    
985        foreach my $row (@$rows) {
986            my $single_domain = [];
987            $count++;
988    
989            # construct the score link
990            my $score = $row->[0];
991            my $toid = $row->[1];
992            my $link = $cgi->url(-relative => 1) . "?user=master&request=show_coupling_evidence&prot=$fid&to=$toid&SPROUT=";
993            my $sc_link = "<a href=$link>$score</a>";
994    
995            push(@$single_domain,$sc_link);
996            push(@$single_domain,$row->[1]);
997            push(@$single_domain,$row->[2]);
998            push(@$functional_data,$single_domain);
999        }
1000    
1001        if ($count >0){
1002            $content = $functional_data;
1003        }
1004        else
1005        {
1006            $content = "<p>This PEG does not have any functional coupling</p>";
1007        }
1008        return ($content);
1009    }
1010    
1011    
1012    #########################################
1013    #########################################
1014    package Observation::Domain;
1015    
1016    use base qw(Observation);
1017    
1018    sub new {
1019    
1020        my ($class,$dataset) = @_;
1021        my $self = $class->SUPER::new($dataset);
1022        $self->{evalue} = $dataset->{'evalue'};
1023        $self->{acc} = $dataset->{'acc'};
1024        $self->{start} = $dataset->{'start'};
1025        $self->{stop} = $dataset->{'stop'};
1026    
1027        bless($self,$class);
1028        return $self;
1029    }
1030    
1031    sub display {
1032        my ($thing,$gd) = @_;
1033        my $lines = [];
1034        my $line_config = { 'title' => $thing->acc,
1035                            'short_title' => $thing->type,
1036                            'basepair_offset' => '1' };
1037        my $color = "4";
1038    
1039        my $line_data = [];
1040        my $links_list = [];
1041        my $descriptions = [];
1042    
1043        my $db_and_id = $thing->acc;
1044        my ($db,$id) = split("::",$db_and_id);
1045    
1046        my $dbmaster = DBMaster->new(-database =>'Ontology');
1047    
1048        my ($name_title,$name_value,$description_title,$description_value);
1049        if($db eq "CDD"){
1050            my $cdd_objs = $dbmaster->cdd->get_objects( { 'id' => $id } );
1051            if(!scalar(@$cdd_objs)){
1052                $name_title = "name";
1053                $name_value = "not available";
1054                $description_title = "description";
1055                $description_value = "not available";
1056            }
1057            else{
1058                my $cdd_obj = $cdd_objs->[0];
1059                $name_title = "name";
1060                $name_value = $cdd_obj->term;
1061                $description_title = "description";
1062                $description_value = $cdd_obj->description;
1063            }
1064        }
1065    
1066        my $name;
1067        $name = {"title" => $name_title,
1068                 "value" => $name_value};
1069        push(@$descriptions,$name);
1070    
1071        my $description;
1072        $description = {"title" => $description_title,
1073                                 "value" => $description_value};
1074        push(@$descriptions,$description);
1075    
1076        my $score;
1077        $score = {"title" => "score",
1078                  "value" => $thing->evalue};
1079        push(@$descriptions,$score);
1080    
1081        my $link_id;
1082        if ($thing->acc =~/\w+::(\d+)/){
1083            $link_id = $1;
1084        }
1085    
1086        my $link;
1087        my $link_url;
1088        if ($thing->class eq "CDD"){$link_url = "http://0-www.ncbi.nlm.nih.gov.library.vu.edu.au:80/Structure/cdd/cddsrv.cgi?uid=$link_id"}
1089        elsif($thing->class eq "PFAM"){$link_url = "http://www.sanger.ac.uk/cgi-bin/Pfam/getacc?$link_id"}
1090        else{$link_url = "NO_URL"}
1091    
1092        $link = {"link_title" => $thing->acc,
1093                 "link" => $link_url};
1094        push(@$links_list,$link);
1095    
1096        my $element_hash = {
1097            "title" => $thing->type,
1098            "start" => $thing->start,
1099            "end" =>  $thing->stop,
1100            "color"=> $color,
1101            "zlayer" => '2',
1102            "links_list" => $links_list,
1103            "description" => $descriptions};
1104    
1105        push(@$line_data,$element_hash);
1106        $gd->add_line($line_data, $line_config);
1107    
1108        return $gd;
1109    
1110    }
1111    
1112    #########################################
1113    #########################################
1114    package Observation::Location;
1115    
1116    use base qw(Observation);
1117    
1118    sub new {
1119    
1120        my ($class,$dataset) = @_;
1121        my $self = $class->SUPER::new($dataset);
1122        $self->{cleavage_prob} = $dataset->{'cleavage_prob'};
1123        $self->{cleavage_loc} = $dataset->{'cleavage_loc'};
1124        $self->{signal_peptide_score} = $dataset->{'signal_peptide_score'};
1125        $self->{cello_location} = $dataset->{'cello_location'};
1126        $self->{cello_score} = $dataset->{'cello_score'};
1127        $self->{tmpred_score} = $dataset->{'tmpred_score'};
1128        $self->{tmpred_locations} = $dataset->{'tmpred_locations'};
1129    
1130        bless($self,$class);
1131        return $self;
1132    }
1133    
1134    sub display {
1135        my ($thing,$gd) = @_;
1136    
1137        my $fid = $thing->fig_id;
1138        my $fig= new FIG;
1139        my $length = length($fig->get_translation($fid));
1140    
1141        my $cleavage_prob;
1142        if($thing->cleavage_prob){$cleavage_prob = $thing->cleavage_prob;}
1143        my ($cleavage_loc_begin,$cleavage_loc_end) = split("-",$thing->cleavage_loc);
1144        my $signal_peptide_score = $thing->signal_peptide_score;
1145        my $cello_location = $thing->cello_location;
1146        my $cello_score = $thing->cello_score;
1147        my $tmpred_score = $thing->tmpred_score;
1148        my @tmpred_locations = split(",",$thing->tmpred_locations);
1149    
1150        my $lines = [];
1151        my $line_config = { 'title' => 'Localization Evidence',
1152                            'short_title' => 'Local',
1153                            'basepair_offset' => '1' };
1154    
1155        #color is
1156        my $color = "5";
1157    
1158        my $line_data = [];
1159    
1160        if($cello_location){
1161            my $cello_descriptions = [];
1162            my $description_cello_location = {"title" => 'Best Cello Location',
1163                                              "value" => $cello_location};
1164    
1165            push(@$cello_descriptions,$description_cello_location);
1166    
1167            my $description_cello_score = {"title" => 'Cello Score',
1168                                           "value" => $cello_score};
1169    
1170            push(@$cello_descriptions,$description_cello_score);
1171    
1172            my $element_hash = {
1173                "title" => "CELLO",
1174                "start" => "1",
1175                "end" =>  $length + 1,
1176                "color"=> $color,
1177                "type" => 'box',
1178                "zlayer" => '2',
1179                "description" => $cello_descriptions};
1180    
1181            push(@$line_data,$element_hash);
1182        }
1183    
1184        my $color = "6";
1185        if($tmpred_score){
1186            foreach my $tmpred (@tmpred_locations){
1187                my $descriptions = [];
1188                my ($begin,$end) =split("-",$tmpred);
1189                my $description_tmpred_score = {"title" => 'TMPRED score',
1190                                 "value" => $tmpred_score};
1191    
1192                push(@$descriptions,$description_tmpred_score);
1193    
1194                my $element_hash = {
1195                "title" => "transmembrane location",
1196                "start" => $begin + 1,
1197                "end" =>  $end + 1,
1198                "color"=> $color,
1199                "zlayer" => '5',
1200                "type" => 'smallbox',
1201                "description" => $descriptions};
1202    
1203                push(@$line_data,$element_hash);
1204            }
1205        }
1206    
1207        my $color = "1";
1208        if($signal_peptide_score){
1209            my $descriptions = [];
1210            my $description_signal_peptide_score = {"title" => 'signal peptide score',
1211                                                    "value" => $signal_peptide_score};
1212    
1213            push(@$descriptions,$description_signal_peptide_score);
1214    
1215            my $description_cleavage_prob = {"title" => 'cleavage site probability',
1216                                             "value" => $cleavage_prob};
1217    
1218            push(@$descriptions,$description_cleavage_prob);
1219    
1220            my $element_hash = {
1221                "title" => "SignalP",
1222                "start" => $cleavage_loc_begin - 2,
1223                "end" =>  $cleavage_loc_end + 3,
1224                "type" => 'bigbox',
1225                "color"=> $color,
1226                "zlayer" => '10',
1227                "description" => $descriptions};
1228    
1229            push(@$line_data,$element_hash);
1230        }
1231    
1232        $gd->add_line($line_data, $line_config);
1233    
1234        return ($gd);
1235    
1236    }
1237    
1238    sub cleavage_loc {
1239      my ($self) = @_;
1240    
1241      return $self->{cleavage_loc};
1242    }
1243    
1244    sub cleavage_prob {
1245      my ($self) = @_;
1246    
1247      return $self->{cleavage_prob};
1248    }
1249    
1250    sub signal_peptide_score {
1251      my ($self) = @_;
1252    
1253      return $self->{signal_peptide_score};
1254    }
1255    
1256    sub tmpred_score {
1257      my ($self) = @_;
1258    
1259      return $self->{tmpred_score};
1260    }
1261    
1262    sub tmpred_locations {
1263      my ($self) = @_;
1264    
1265      return $self->{tmpred_locations};
1266    }
1267    
1268    sub cello_location {
1269      my ($self) = @_;
1270    
1271      return $self->{cello_location};
1272    }
1273    
1274    sub cello_score {
1275      my ($self) = @_;
1276    
1277      return $self->{cello_score};
1278    }
1279    
1280    
1281    #########################################
1282    #########################################
1283    package Observation::Sims;
1284    
1285    use base qw(Observation);
1286    
1287    sub new {
1288    
1289        my ($class,$dataset) = @_;
1290        my $self = $class->SUPER::new($dataset);
1291        $self->{identity} = $dataset->{'identity'};
1292        $self->{acc} = $dataset->{'acc'};
1293        $self->{evalue} = $dataset->{'evalue'};
1294        $self->{qstart} = $dataset->{'qstart'};
1295        $self->{qstop} = $dataset->{'qstop'};
1296        $self->{hstart} = $dataset->{'hstart'};
1297        $self->{hstop} = $dataset->{'hstop'};
1298        $self->{database} = $dataset->{'database'};
1299        $self->{organism} = $dataset->{'organism'};
1300        $self->{function} = $dataset->{'function'};
1301        $self->{qlength} = $dataset->{'qlength'};
1302        $self->{hlength} = $dataset->{'hlength'};
1303    
1304        bless($self,$class);
1305        return $self;
1306    }
1307    
1308    =head3 display()
1309    
1310    If available use the function specified here to display a graphical observation.
1311    This code will display a graphical view of the similarities using the genome drawer object
1312    
1313    =cut
1314    
1315    sub display {
1316        my ($self,$gd) = @_;
1317    
1318        my $fig = new FIG;
1319        my $peg = $self->acc;
1320    
1321        my $organism = $self->organism;
1322        my $function = $self->function;
1323        my $abbrev_name = $fig->abbrev($organism);
1324        my $align_start = $self->qstart;
1325        my $align_stop = $self->qstop;
1326        my $hit_start = $self->hstart;
1327        my $hit_stop = $self->hstop;
1328    
1329        my $line_config = { 'title' => "$organism",
1330                            'short_title' => "$abbrev_name",
1331                            'basepair_offset' => '0'
1332                            };
1333    
1334        my $line_data = [];
1335    
1336        my $element_hash;
1337        my $links_list = [];
1338        my $descriptions = [];
1339    
1340        # get subsystem information
1341        my $url_link = "http://seed-viewer.theseed.org/index.cgi?action=ShowAnnotation&prot=".$peg;
1342    
1343        my $link;
1344        $link = {"link_title" => $peg,
1345                 "link" => $url_link};
1346        push(@$links_list,$link);
1347    
1348        my @subsystems = $fig->peg_to_subsystems($peg);
1349        foreach my $subsystem (@subsystems){
1350            my $link;
1351            $link = {"link" => "http://seed-viewer.theseed.org/index.cgi?action=ShowSubsystem&subsystem_name=$subsystem",
1352                     "link_title" => $subsystem};
1353            push(@$links_list,$link);
1354        }
1355    
1356        my $description_function;
1357        $description_function = {"title" => "function",
1358                                 "value" => $function};
1359        push(@$descriptions,$description_function);
1360    
1361        my $description_ss;
1362        my $ss_string = join (",", @subsystems);
1363        $description_ss = {"title" => "subsystems",
1364                           "value" => $ss_string};
1365        push(@$descriptions,$description_ss);
1366    
1367        my $description_loc;
1368        $description_loc = {"title" => "location start",
1369                            "value" => $hit_start};
1370        push(@$descriptions, $description_loc);
1371    
1372        $description_loc = {"title" => "location stop",
1373                            "value" => $hit_stop};
1374        push(@$descriptions, $description_loc);
1375    
1376        my $color;
1377        my $evalue = $self->evalue;
1378        while ($evalue =~ /-0/)
1379        {
1380            my ($chunk1, $chunk2) = split(/-/, $evalue);
1381            $chunk2 = substr($chunk2,1);
1382            $evalue = $chunk1 . "-" . $chunk2;
1383        }
1384    
1385        $color = 1;
1386    
1387        my $description_eval = {"title" => "E-Value",
1388                                "value" => $evalue};
1389        push(@$descriptions, $description_eval);
1390    
1391        my $identity = $self->identity;
1392        my $description_identity = {"title" => "Identity",
1393                                    "value" => $identity};
1394        push(@$descriptions, $description_identity);
1395    
1396        $element_hash = {
1397            "title" => $peg,
1398            "start" => $align_start,
1399            "end" =>  $align_stop,
1400            "type"=> 'box',
1401            "color"=> $color,
1402            "zlayer" => "2",
1403            "links_list" => $links_list,
1404            "description" => $descriptions
1405            };
1406        push(@$line_data,$element_hash);
1407        $gd->add_line($line_data, $line_config);
1408    
1409        return ($gd);
1410    
1411    }
1412    
1413    =head3 display_table()
1414    
1415    If available use the function specified here to display the "raw" observation.
1416    This code will display a table for the similarities protein
1417    
1418    B<Please note> that URL linked to in display_method() is an external component and needs to added to the code for every class of evidence.
1419    
1420    =cut
1421    
1422    sub display_table {
1423        my ($self,$dataset) = @_;
1424    
1425        my $data = [];
1426        my $count = 0;
1427        my $content;
1428        my $fig = new FIG;
1429        my $cgi = new CGI;
1430        foreach my $thing (@$dataset) {
1431            my $single_domain = [];
1432            next if ($thing->class ne "SIM");
1433            $count++;
1434    
1435            my $id = $thing->acc;
1436    
1437            # add the subsystem information
1438            my @in_sub  = $fig->peg_to_subsystems($id);
1439            my $in_sub;
1440    
1441            if (@in_sub > 0) {
1442                $in_sub = @in_sub;
1443    
1444                # RAE: add a javascript popup with all the subsystems
1445                my $ss_list=join "<br>", map { my $g = $_; $g =~ s/\_/ /g; $_ = $g } sort {$a cmp $b} @in_sub;
1446                $in_sub = $cgi->a( {id=>"subsystems", onMouseover=>"javascript:if(!this.tooltip) this.tooltip=new Popup_Tooltip(this, 'Subsystems', '$ss_list', ''); this.tooltip.addHandler(); return false;"}, $in_sub);
1447            } else {
1448                $in_sub = "&nbsp;";
1449            }
1450    
1451            # add evidence code with tool tip
1452            my $ev_codes=" &nbsp; ";
1453            my @ev_codes = "";
1454            if ($id =~ /^fig\|\d+\.\d+\.peg\.\d+$/) {
1455                my @codes = grep { $_->[1] =~ /^evidence_code/i } $fig->get_attributes($id);
1456                @ev_codes = ();
1457                foreach my $code (@codes) {
1458                    my $pretty_code = $code->[2];
1459                    if ($pretty_code =~ /;/) {
1460                        my ($cd, $ss) = split(";", $code->[2]);
1461                        $ss =~ s/_/ /g;
1462                        $pretty_code = $cd;# . " in " . $ss;
1463                    }
1464                    push(@ev_codes, $pretty_code);
1465                }
1466            }
1467    
1468            if (scalar(@ev_codes) && $ev_codes[0]) {
1469                my $ev_code_help=join("<br />", map {&HTML::evidence_codes_explain($_)} @ev_codes);
1470                $ev_codes = $cgi->a(
1471                                    {
1472                                        id=>"evidence_codes", onMouseover=>"javascript:if(!this.tooltip) this.tooltip=new Popup_Tooltip(this, 'Evidence Codes', '$ev_code_help', ''); this.tooltip.addHandler(); return false;"}, join("<br />", @ev_codes));
1473            }
1474    
1475            # add the aliases
1476            my $aliases = undef;
1477            $aliases = &html_enc( join( ", ", $fig->feature_aliases($id) ) );
1478            $aliases = &HTML::set_prot_links( $cgi, $aliases );
1479            $aliases ||= "&nbsp;";
1480    
1481            my $iden    = $thing->identity;
1482            my $ln1     = $thing->qlength;
1483            my $ln2     = $thing->hlength;
1484            my $b1      = $thing->qstart;
1485            my $e1      = $thing->qstop;
1486            my $b2      = $thing->hstart;
1487            my $e2      = $thing->hstop;
1488            my $d1      = abs($e1 - $b1) + 1;
1489            my $d2      = abs($e2 - $b2) + 1;
1490            my $reg1    = "$b1-$e1 (<b>$d1/$ln1</b>)";
1491            my $reg2    = "$b2-$e2 (<b>$d2/$ln2</b>)";
1492    
1493    
1494            push(@$single_domain,$thing->database);
1495            push(@$single_domain,&HTML::set_prot_links($cgi,$thing->acc));
1496            push(@$single_domain,$thing->evalue);
1497            push(@$single_domain,"$iden\%");
1498            push(@$single_domain,$reg1);
1499            push(@$single_domain,$reg2);
1500            push(@$single_domain,$in_sub);
1501            push(@$single_domain,$ev_codes);
1502            push(@$single_domain,$thing->organism);
1503            push(@$single_domain,$thing->function);
1504            push(@$single_domain,$aliases);
1505            push(@$data,$single_domain);
1506        }
1507    
1508        if ($count >0){
1509            $content = $data;
1510        }
1511        else
1512        {
1513            $content = "<p>This PEG does not have any similarities</p>";
1514        }
1515        return ($content);
1516    }
1517    
1518    sub html_enc { $_ = $_[0]; s/\&/&amp;/g; s/\>/&gt;/g; s/\</&lt;/g; $_ }
1519    
1520    
1521    
1522    ############################
1523    package Observation::Cluster;
1524    
1525    use base qw(Observation);
1526    
1527    sub new {
1528    
1529        my ($class,$dataset) = @_;
1530        my $self = $class->SUPER::new($dataset);
1531        $self->{context} = $dataset->{'context'};
1532        bless($self,$class);
1533        return $self;
1534    }
1535    
1536    sub display {
1537        my ($self,$gd) = @_;
1538    
1539        my $fid = $self->fig_id;
1540        my $compare_or_coupling = $self->context;
1541        my $gd_window_size = $gd->window_size;
1542        my $fig = new FIG;
1543        my $all_regions = [];
1544    
1545        #get the organism genome
1546        my $target_genome = $fig->genome_of($fid);
1547    
1548        # get location of the gene
1549        my $data = $fig->feature_location($fid);
1550        my ($contig, $beg, $end);
1551        my %reverse_flag;
1552    
1553        if ($data =~ /(.*)_(\d+)_(\d+)$/){
1554            $contig = $1;
1555            $beg = $2;
1556            $end = $3;
1557        }
1558    
1559        my $offset;
1560        my ($region_start, $region_end);
1561        if ($beg < $end)
1562        {
1563            $region_start = $beg - 4000;
1564            $region_end = $end+4000;
1565            $offset = ($2+(($3-$2)/2))-($gd_window_size/2);
1566        }
1567        else
1568        {
1569            $region_start = $end-4000;
1570            $region_end = $beg+4000;
1571            $offset = ($3+(($2-$3)/2))-($gd_window_size/2);
1572            $reverse_flag{$target_genome} = $fid;
1573        }
1574    
1575        # call genes in region
1576        my ($target_gene_features, $reg_beg, $reg_end) = $fig->genes_in_region($target_genome, $contig, $region_start, $region_end);
1577        push(@$all_regions,$target_gene_features);
1578        my (@start_array_region);
1579        push (@start_array_region, $offset);
1580    
1581        my %all_genes;
1582        my %all_genomes;
1583        foreach my $feature (@$target_gene_features){ $all_genes{$feature} = $fid;}
1584    
1585        if ($compare_or_coupling eq "diverse")
1586        {
1587            my @coup = grep { $_->[1]} $fig->coupling_and_evidence($fid,5000,1e-10,4,1);
1588    
1589            my $coup_count = 0;
1590    
1591            foreach my $pair (@{$coup[0]->[2]}) {
1592                #   last if ($coup_count > 10);
1593                my ($peg1,$peg2) = @$pair;
1594    
1595                my ($pair_contig,$pair_beg,$pair_end,$pair_region_start,$pair_region_stop,$pair_genome);
1596                $pair_genome = $fig->genome_of($peg1);
1597    
1598                my $location = $fig->feature_location($peg1);
1599                if($location =~/(.*)_(\d+)_(\d+)$/){
1600                    $pair_contig = $1;
1601                    $pair_beg = $2;
1602                    $pair_end = $3;
1603                    if ($pair_beg < $pair_end)
1604                    {
1605                        $pair_region_start = $pair_beg - 4000;
1606                        $pair_region_stop = $pair_end+4000;
1607                        $offset = ($2+(($3-$2)/2))-($gd_window_size/2);
1608                    }
1609                    else
1610                    {
1611                        $pair_region_start = $pair_end-4000;
1612                        $pair_region_stop = $pair_beg+4000;
1613                        $offset = ($3+(($2-$3)/2))-($gd_window_size/2);
1614                        $reverse_flag{$pair_genome} = $peg1;
1615                    }
1616    
1617                    push (@start_array_region, $offset);
1618    
1619                    $all_genomes{$pair_genome} = 1;
1620                    my ($pair_features) = $fig->genes_in_region($pair_genome, $pair_contig, $pair_region_start, $pair_region_stop);
1621                    push(@$all_regions,$pair_features);
1622                    foreach my $pair_feature (@$pair_features){ $all_genes{$pair_feature} = $peg1;}
1623                }
1624                $coup_count++;
1625            }
1626        }
1627    
1628        elsif ($compare_or_coupling eq "close")
1629        {
1630            # make a hash of genomes that are phylogenetically close
1631            #my $close_threshold = ".26";
1632            #my @genomes = $fig->genomes('complete');
1633            #my %close_genomes = ();
1634            #foreach my $compared_genome (@genomes)
1635            #{
1636            #    my $dist = $fig->crude_estimate_of_distance($target_genome,$compared_genome);
1637            #    #$close_genomes{$compared_genome} = $dist;
1638            #    if ($dist <= $close_threshold)
1639            #    {
1640            #       $all_genomes{$compared_genome} = 1;
1641            #    }
1642            #}
1643            $all_genomes{"216592.1"} = 1;
1644            $all_genomes{"79967.1"} = 1;
1645            $all_genomes{"199310.1"} = 1;
1646            $all_genomes{"216593.1"} = 1;
1647            $all_genomes{"155864.1"} = 1;
1648            $all_genomes{"83334.1"} = 1;
1649            $all_genomes{"316407.3"} = 1;
1650    
1651            foreach my $comp_genome (keys %all_genomes){
1652                my $return = $fig->bbh_list($comp_genome,[$fid]);
1653                my $feature_list = $return->{$fid};
1654                foreach my $peg1 (@$feature_list){
1655                    my $location = $fig->feature_location($peg1);
1656                    my ($pair_contig,$pair_beg,$pair_end,$pair_region_start,$pair_region_stop,$pair_genome);
1657                    $pair_genome = $fig->genome_of($peg1);
1658    
1659                    if($location =~/(.*)_(\d+)_(\d+)$/){
1660                        $pair_contig = $1;
1661                        $pair_beg = $2;
1662                        $pair_end = $3;
1663                        if ($pair_beg < $pair_end)
1664                        {
1665                            $pair_region_start = $pair_beg - 4000;
1666                            $pair_region_stop = $pair_end + 4000;
1667                            $offset = ($2+(($3-$2)/2))-($gd_window_size/2);
1668                        }
1669                        else
1670                        {
1671                            $pair_region_start = $pair_end-4000;
1672                            $pair_region_stop = $pair_beg+4000;
1673                            $offset = ($3+(($2-$3)/2))-($gd_window_size/2);
1674                            $reverse_flag{$pair_genome} = $peg1;
1675                        }
1676    
1677                        push (@start_array_region, $offset);
1678                        $all_genomes{$pair_genome} = 1;
1679                        my ($pair_features) = $fig->genes_in_region($pair_genome, $pair_contig, $pair_region_start, $pair_region_stop);
1680                        push(@$all_regions,$pair_features);
1681                        foreach my $pair_feature (@$pair_features){ $all_genes{$pair_feature} = $peg1;}
1682                    }
1683                }
1684            }
1685        }
1686    
1687        # get the PCH to each of the genes
1688        my $pch_sets = [];
1689        my %pch_already;
1690        foreach my $gene_peg (keys %all_genes)
1691        {
1692            if ($pch_already{$gene_peg}){next;};
1693            my $gene_set = [$gene_peg];
1694            foreach my $pch_peg ($fig->in_pch_pin_with($gene_peg)) {
1695                $pch_peg =~ s/,.*$//;
1696                my $pch_genome = $fig->genome_of($pch_peg);
1697                if ( ($gene_peg ne $pch_peg) && ($all_genomes{$pch_genome})) {
1698                    push(@$gene_set,$pch_peg);
1699                    $pch_already{$pch_peg}=1;
1700                }
1701                $pch_already{$gene_peg}=1;
1702            }
1703            push(@$pch_sets,$gene_set);
1704        }
1705    
1706        #create a rank of the pch's
1707        my %pch_set_rank;
1708        my $order = 0;
1709        foreach my $set (@$pch_sets){
1710            my $count = scalar(@$set);
1711            $pch_set_rank{$order} = $count;
1712            $order++;
1713        }
1714    
1715        my %peg_rank;
1716        my $counter =  1;
1717        foreach my $pch_order (sort {$pch_set_rank{$b} <=> $pch_set_rank{$a}} keys %pch_set_rank){
1718            my $good_set = @$pch_sets[$pch_order];
1719            my $flag_set = 0;
1720            if (scalar (@$good_set) > 1)
1721            {
1722                foreach my $peg (@$good_set){
1723                    if ((!$peg_rank{$peg})){
1724                        $peg_rank{$peg} = $counter;
1725                        $flag_set = 1;
1726                    }
1727                }
1728                $counter++ if ($flag_set == 1);
1729            }
1730            else
1731            {
1732                foreach my $peg (@$good_set){
1733                    $peg_rank{$peg} = 100;
1734                }
1735            }
1736        }
1737    
1738    
1739    #    my $bbh_sets = [];
1740    #    my %already;
1741    #    foreach my $gene_key (keys(%all_genes)){
1742    #       if($already{$gene_key}){next;}
1743    #       my $gene_set = [$gene_key];
1744    #
1745    #       my $gene_key_genome = $fig->genome_of($gene_key);
1746    #
1747    #       foreach my $genome_key (keys(%all_genomes)){
1748    #           #next if ($gene_key_genome eq $genome_key);
1749    #           my $return = $fig->bbh_list($genome_key,[$gene_key]);
1750    #
1751    #           my $feature_list = $return->{$gene_key};
1752    #           foreach my $fl (@$feature_list){
1753    #               push(@$gene_set,$fl);
1754    #           }
1755    #       }
1756    #       $already{$gene_key} = 1;
1757    #       push(@$bbh_sets,$gene_set);
1758    #    }
1759    #
1760    #    my %bbh_set_rank;
1761    #    my $order = 0;
1762    #    foreach my $set (@$bbh_sets){
1763    #       my $count = scalar(@$set);
1764    #       $bbh_set_rank{$order} = $count;
1765    #       $order++;
1766    #    }
1767    #
1768    #    my %peg_rank;
1769    #    my $counter =  1;
1770    #    foreach my $bbh_order (sort {$bbh_set_rank{$b} <=> $bbh_set_rank{$a}} keys %bbh_set_rank){
1771    #       my $good_set = @$bbh_sets[$bbh_order];
1772    #       my $flag_set = 0;
1773    #       if (scalar (@$good_set) > 1)
1774    #       {
1775    #           foreach my $peg (@$good_set){
1776    #               if ((!$peg_rank{$peg})){
1777    #                   $peg_rank{$peg} = $counter;
1778    #                   $flag_set = 1;
1779    #               }
1780    #           }
1781    #           $counter++ if ($flag_set == 1);
1782    #       }
1783    #       else
1784    #       {
1785    #           foreach my $peg (@$good_set){
1786    #               $peg_rank{$peg} = 100;
1787    #           }
1788    #       }
1789    #    }
1790    
1791        foreach my $region (@$all_regions){
1792            my $sample_peg = @$region[0];
1793            my $region_genome = $fig->genome_of($sample_peg);
1794            my $region_gs = $fig->genus_species($region_genome);
1795            my $abbrev_name = $fig->abbrev($region_gs);
1796            my $line_config = { 'title' => $region_gs,
1797                                'short_title' => $abbrev_name,
1798                                'basepair_offset' => '0'
1799                                };
1800    
1801            my $offsetting = shift @start_array_region;
1802    
1803            my $second_line_config = { 'title' => "$region_gs",
1804                                       'short_title' => "",
1805                                       'basepair_offset' => '0'
1806                                       };
1807    
1808            my $line_data = [];
1809            my $second_line_data = [];
1810    
1811            # initialize variables to check for overlap in genes
1812            my ($prev_start, $prev_stop, $prev_fig, $second_line_flag);
1813            my $major_line_flag = 0;
1814            my $prev_second_flag = 0;
1815    
1816            foreach my $fid1 (@$region){
1817                $second_line_flag = 0;
1818                my $element_hash;
1819                my $links_list = [];
1820                my $descriptions = [];
1821    
1822                my $color = $peg_rank{$fid1};
1823    
1824                # get subsystem information
1825                my $function = $fig->function_of($fid1);
1826                my $url_link = "http://seed-viewer.theseed.org/index.cgi?action=ShowAnnotation&prot=".$fid1;
1827    
1828                my $link;
1829                $link = {"link_title" => $fid1,
1830                         "link" => $url_link};
1831                push(@$links_list,$link);
1832    
1833                my @subsystems = $fig->peg_to_subsystems($fid1);
1834                foreach my $subsystem (@subsystems){
1835                    my $link;
1836                    $link = {"link" => "http://seed-viewer.theseed.org/index.cgi?action=ShowSubsystem&subsystem_name=$subsystem",
1837                             "link_title" => $subsystem};
1838                    push(@$links_list,$link);
1839                }
1840    
1841                my $description_function;
1842                $description_function = {"title" => "function",
1843                                         "value" => $function};
1844                push(@$descriptions,$description_function);
1845    
1846                my $description_ss;
1847                my $ss_string = join (",", @subsystems);
1848                $description_ss = {"title" => "subsystems",
1849                                   "value" => $ss_string};
1850                push(@$descriptions,$description_ss);
1851    
1852    
1853                my $fid_location = $fig->feature_location($fid1);
1854                if($fid_location =~/(.*)_(\d+)_(\d+)$/){
1855                    my($start,$stop);
1856                    $start = $2 - $offsetting;
1857                    $stop = $3 - $offsetting;
1858    
1859                    if ( (($prev_start) && ($prev_stop) ) &&
1860                         ( ($start < $prev_start) || ($start < $prev_stop) ||
1861                           ($stop < $prev_start) || ($stop < $prev_stop) )){
1862                        if (($second_line_flag == 0) && ($prev_second_flag == 0)) {
1863                            $second_line_flag = 1;
1864                            $major_line_flag = 1;
1865                        }
1866                    }
1867                    $prev_start = $start;
1868                    $prev_stop = $stop;
1869                    $prev_fig = $fid1;
1870    
1871                    if ((defined($reverse_flag{$region_genome})) && ($reverse_flag{$region_genome} eq $all_genes{$fid1})){
1872                        print STDERR "$region_genome - $fid1";
1873                        $start = $gd_window_size - $start;
1874                        $stop = $gd_window_size - $stop;
1875                    }
1876    
1877                    $element_hash = {
1878                        "title" => $fid1,
1879                        "start" => $start,
1880                        "end" =>  $stop,
1881                        "type"=> 'arrow',
1882                        "color"=> $color,
1883                        "zlayer" => "2",
1884                        "links_list" => $links_list,
1885                        "description" => $descriptions
1886                    };
1887    
1888                    # if there is an overlap, put into second line
1889                    if ($second_line_flag == 1){ push(@$second_line_data,$element_hash); $prev_second_flag = 1;}
1890                    else{ push(@$line_data,$element_hash); $prev_second_flag = 0;}
1891    
1892                }
1893            }
1894            $gd->add_line($line_data, $line_config);
1895            $gd->add_line($second_line_data, $second_line_config) if ($major_line_flag == 1);
1896        }
1897        return $gd;
1898    }
1899    
1900    

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.25

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3