[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.1, Thu May 31 02:08:28 2007 UTC revision 1.24, Tue Jul 10 20:11:38 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";  
   
 $fig->get_evidence($fid);  
 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 58  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 94  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 176  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    my $url = get_url($self->type, $self->acc);  =head3 display()
277    
278    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  It will probably have to:  sub display_table {
295    
296      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 349  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  =head3 get_display_method (internal)      push (@{$datasets_ref} ,$dataset);
448    
449  get_display_method() return a valid URL or undef for any observation.  }
450    
451    =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 $fig = new FIG;
461    
462   my ($self) = @_;      foreach my $attr_ref ($fig->get_attributes($fid,'PDB')) {
463    
464  # a hash with a URL for each observation; identified by name()          my $key = @$attr_ref[1];
465  #my $URL             => { 'sim'=> "http://www.theseed.org/featalign.cgi?id1=",\          my($key1,$key2) =split("::",$key);
466  #                        'bbh'=> "http://www.theseed.org/featalign.cgi?id1="          my $value = @$attr_ref[2];
467  # };          my ($evalue,$location) = split(";",$value);
468    
469            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  #if (defined $URL{$self->name}) {          push (@{$datasets_ref} ,$dataset);
488  #     $url = $URL{$self->name}.$self->feature_id."&id2=".$self->acc;      }
 #     return $url;  
 # }  
 # else  
      return undef;  
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      # we read a FIG ID and a reference to an array (of arrays of hashes, see above)      my $dataset = {'class' => 'CLUSTER',
501      my ($fid,$datasets_ref) = (@_);                     'type' => 'fc',
502                       'context' => $scope,
503                       'fig_id' => $fid
504                       };
505        push (@{$datasets_ref} ,$dataset);
506    }
507    
     my $_myfig = new FIG;  
508    
509      foreach my $attr_ref ($_myfig->get_attributes($fid)) {  =head3 get_sims_observations() (internal)
510    
511          # convert the ref into a string for easier handling  This methods retrieves sims fills the internal data structures.
         my ($string) = "@$attr_ref";  
512    
513  #       print "S:$string\n";  =cut
         my ($key,$val) = ( $string =~ /\S+\s(\S+)\s(\S+)/);  
514    
515          # THIS SHOULD BE DONE ANOTHER WAY FM->TD  sub get_sims_observations{
         # we need to do the right thing for each type, ie no evalue for CELLO and no coordinates, but a score, etc  
         # as fas as possible this should be configured so that the type of observation and the regexp are  
         # stored somewhere for easy expansion  
         #  
516    
517          if (($key =~ /PFAM::/) || ( $key =~ /IPR::/) || ( $key =~ /CDD::/) ) {      my ($fid,$datasets_ref) = (@_);
518        my $fig = new FIG;
519        my @sims= $fig->nsims($fid,100,1e-20,"all");
520        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              # some keys are composite CDD::1233244 or PFAM:PF1233  =head3 get_database (internal)
557    This method gets the database association from the sequence id
558    
559    =cut
560    
561    sub get_database{
562        my ($id) = (@_);
563    
564        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{
                 if($raw_evalue =~/(\d+).(\d+)/){  
589    
590  #                   $new_exp = (1000+$expo);      my ($fid,$datasets_ref) = (@_);
591          #           $new_k = $k / 100;      my $fig = new FIG;
592        my $funcs_ref;
593    
594        my @maps_to = grep { $_ ne $fid and $_ !~ /^xxx/ } map { $_->[0] } $fig->mapped_prot_ids($fid);
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 533  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 565  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_table()
1309    
1310    If available use the function specified here to display the "raw" observation.
1311    This code will display a table for the similarities protein
1312    
1313    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.
1314    
1315    =cut
1316    
1317    sub display_table {
1318        my ($self,$dataset) = @_;
1319    
1320        my $data = [];
1321        my $count = 0;
1322        my $content;
1323        my $fig = new FIG;
1324        my $cgi = new CGI;
1325        foreach my $thing (@$dataset) {
1326            my $single_domain = [];
1327            next if ($thing->class ne "SIM");
1328            $count++;
1329    
1330            my $id = $thing->acc;
1331    
1332            # add the subsystem information
1333            my @in_sub  = $fig->peg_to_subsystems($id);
1334            my $in_sub;
1335    
1336            if (@in_sub > 0) {
1337                $in_sub = @in_sub;
1338    
1339                # RAE: add a javascript popup with all the subsystems
1340                my $ss_list=join "<br>", map { my $g = $_; $g =~ s/\_/ /g; $_ = $g } sort {$a cmp $b} @in_sub;
1341                $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);
1342            } else {
1343                $in_sub = "&nbsp;";
1344            }
1345    
1346            # add evidence code with tool tip
1347            my $ev_codes=" &nbsp; ";
1348            my @ev_codes = "";
1349            if ($id =~ /^fig\|\d+\.\d+\.peg\.\d+$/) {
1350                my @codes = grep { $_->[1] =~ /^evidence_code/i } $fig->get_attributes($id);
1351                @ev_codes = ();
1352                foreach my $code (@codes) {
1353                    my $pretty_code = $code->[2];
1354                    if ($pretty_code =~ /;/) {
1355                        my ($cd, $ss) = split(";", $code->[2]);
1356                        $ss =~ s/_/ /g;
1357                        $pretty_code = $cd;# . " in " . $ss;
1358                    }
1359                    push(@ev_codes, $pretty_code);
1360                }
1361            }
1362    
1363            if (scalar(@ev_codes) && $ev_codes[0]) {
1364                my $ev_code_help=join("<br />", map {&HTML::evidence_codes_explain($_)} @ev_codes);
1365                $ev_codes = $cgi->a(
1366                                    {
1367                                        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));
1368            }
1369    
1370            # add the aliases
1371            my $aliases = undef;
1372            $aliases = &html_enc( join( ", ", $fig->feature_aliases($id) ) );
1373            $aliases = &HTML::set_prot_links( $cgi, $aliases );
1374            $aliases ||= "&nbsp;";
1375    
1376            my $iden    = $thing->identity;
1377            my $ln1     = $thing->qlength;
1378            my $ln2     = $thing->hlength;
1379            my $b1      = $thing->qstart;
1380            my $e1      = $thing->qstop;
1381            my $b2      = $thing->hstart;
1382            my $e2      = $thing->hstop;
1383            my $d1      = abs($e1 - $b1) + 1;
1384            my $d2      = abs($e2 - $b2) + 1;
1385            my $reg1    = "$b1-$e1 (<b>$d1/$ln1</b>)";
1386            my $reg2    = "$b2-$e2 (<b>$d2/$ln2</b>)";
1387    
1388    
1389            push(@$single_domain,$thing->database);
1390            push(@$single_domain,&HTML::set_prot_links($cgi,$thing->acc));
1391            push(@$single_domain,$thing->evalue);
1392            push(@$single_domain,"$iden\%");
1393            push(@$single_domain,$reg1);
1394            push(@$single_domain,$reg2);
1395            push(@$single_domain,$in_sub);
1396            push(@$single_domain,$ev_codes);
1397            push(@$single_domain,$thing->organism);
1398            push(@$single_domain,$thing->function);
1399            push(@$single_domain,$aliases);
1400            push(@$data,$single_domain);
1401        }
1402    
1403        if ($count >0){
1404            $content = $data;
1405        }
1406        else
1407        {
1408            $content = "<p>This PEG does not have any similarities</p>";
1409        }
1410        return ($content);
1411    }
1412    
1413    sub html_enc { $_ = $_[0]; s/\&/&amp;/g; s/\>/&gt;/g; s/\</&lt;/g; $_ }
1414    
1415    
1416    
1417    ############################
1418    package Observation::Cluster;
1419    
1420    use base qw(Observation);
1421    
1422    sub new {
1423    
1424        my ($class,$dataset) = @_;
1425        my $self = $class->SUPER::new($dataset);
1426        $self->{context} = $dataset->{'context'};
1427        bless($self,$class);
1428        return $self;
1429    }
1430    
1431    sub display {
1432        my ($self,$gd) = @_;
1433    
1434        my $fid = $self->fig_id;
1435        my $compare_or_coupling = $self->context;
1436        my $gd_window_size = $gd->window_size;
1437        my $fig = new FIG;
1438        my $all_regions = [];
1439    
1440        #get the organism genome
1441        my $target_genome = $fig->genome_of($fid);
1442    
1443        # get location of the gene
1444        my $data = $fig->feature_location($fid);
1445        my ($contig, $beg, $end);
1446        my %reverse_flag;
1447    
1448        if ($data =~ /(.*)_(\d+)_(\d+)$/){
1449            $contig = $1;
1450            $beg = $2;
1451            $end = $3;
1452        }
1453    
1454        my $offset;
1455        my ($region_start, $region_end);
1456        if ($beg < $end)
1457        {
1458            $region_start = $beg - 4000;
1459            $region_end = $end+4000;
1460            $offset = ($2+(($3-$2)/2))-($gd_window_size/2);
1461        }
1462        else
1463        {
1464            $region_start = $end-4000;
1465            $region_end = $beg+4000;
1466            $offset = ($3+(($2-$3)/2))-($gd_window_size/2);
1467            $reverse_flag{$target_genome} = 1;
1468        }
1469    
1470        # call genes in region
1471        my ($target_gene_features, $reg_beg, $reg_end) = $fig->genes_in_region($target_genome, $contig, $region_start, $region_end);
1472        push(@$all_regions,$target_gene_features);
1473        my (@start_array_region);
1474        push (@start_array_region, $offset);
1475    
1476        my %all_genes;
1477        my %all_genomes;
1478        foreach my $feature (@$target_gene_features){ $all_genes{$feature} = 1;}
1479    
1480        if ($compare_or_coupling eq "diverse")
1481        {
1482            my @coup = grep { $_->[1]} $fig->coupling_and_evidence($fid,5000,1e-10,4,1);
1483    
1484            my $coup_count = 0;
1485    
1486            foreach my $pair (@{$coup[0]->[2]}) {
1487                #   last if ($coup_count > 10);
1488                my ($peg1,$peg2) = @$pair;
1489    
1490                my ($pair_contig,$pair_beg,$pair_end,$pair_region_start,$pair_region_stop,$pair_genome);
1491                $pair_genome = $fig->genome_of($peg1);
1492    
1493                my $location = $fig->feature_location($peg1);
1494                if($location =~/(.*)_(\d+)_(\d+)$/){
1495                    $pair_contig = $1;
1496                    $pair_beg = $2;
1497                    $pair_end = $3;
1498                    if ($pair_beg < $pair_end)
1499                    {
1500                        $pair_region_start = $pair_beg - 4000;
1501                        $pair_region_stop = $pair_end+4000;
1502                        $offset = ($2+(($3-$2)/2))-($gd_window_size/2);
1503                    }
1504                    else
1505                    {
1506                        $pair_region_start = $pair_end-4000;
1507                        $pair_region_stop = $pair_beg+4000;
1508                        $offset = ($3+(($2-$3)/2))-($gd_window_size/2);
1509                        $reverse_flag{$pair_genome} = 1;
1510                    }
1511    
1512                    push (@start_array_region, $offset);
1513    
1514                    $all_genomes{$pair_genome} = 1;
1515                    my ($pair_features) = $fig->genes_in_region($pair_genome, $pair_contig, $pair_region_start, $pair_region_stop);
1516                    push(@$all_regions,$pair_features);
1517                    foreach my $pair_feature (@$pair_features){ $all_genes{$pair_feature} = 1;}
1518                }
1519                $coup_count++;
1520            }
1521        }
1522    
1523        elsif ($compare_or_coupling eq "close")
1524        {
1525            # make a hash of genomes that are phylogenetically close
1526            #my $close_threshold = ".26";
1527            #my @genomes = $fig->genomes('complete');
1528            #my %close_genomes = ();
1529            #foreach my $compared_genome (@genomes)
1530            #{
1531            #    my $dist = $fig->crude_estimate_of_distance($target_genome,$compared_genome);
1532            #    #$close_genomes{$compared_genome} = $dist;
1533            #    if ($dist <= $close_threshold)
1534            #    {
1535            #       $all_genomes{$compared_genome} = 1;
1536            #    }
1537            #}
1538            $all_genomes{"216592.1"} = 1;
1539            $all_genomes{"79967.1"} = 1;
1540            $all_genomes{"199310.1"} = 1;
1541            $all_genomes{"216593.1"} = 1;
1542            $all_genomes{"155864.1"} = 1;
1543            $all_genomes{"83334.1"} = 1;
1544            $all_genomes{"316407.3"} = 1;
1545    
1546            foreach my $comp_genome (keys %all_genomes){
1547                my $return = $fig->bbh_list($comp_genome,[$fid]);
1548                my $feature_list = $return->{$fid};
1549                foreach my $peg1 (@$feature_list){
1550                    my $location = $fig->feature_location($peg1);
1551                    my ($pair_contig,$pair_beg,$pair_end,$pair_region_start,$pair_region_stop,$pair_genome);
1552                    $pair_genome = $fig->genome_of($peg1);
1553    
1554                    if($location =~/(.*)_(\d+)_(\d+)$/){
1555                        $pair_contig = $1;
1556                        $pair_beg = $2;
1557                        $pair_end = $3;
1558                        if ($pair_beg < $pair_end)
1559                        {
1560                            $pair_region_start = $pair_beg - 4000;
1561                            $pair_region_stop = $pair_end + 4000;
1562                            $offset = ($2+(($3-$2)/2))-($gd_window_size/2);
1563                        }
1564                        else
1565                        {
1566                            $pair_region_start = $pair_end-4000;
1567                            $pair_region_stop = $pair_beg+4000;
1568                            $offset = ($3+(($2-$3)/2))-($gd_window_size/2);
1569                            $reverse_flag{$pair_genome} = 1;
1570                        }
1571    
1572                        push (@start_array_region, $offset);
1573                        $all_genomes{$pair_genome} = 1;
1574                        my ($pair_features) = $fig->genes_in_region($pair_genome, $pair_contig, $pair_region_start, $pair_region_stop);
1575                        push(@$all_regions,$pair_features);
1576                        foreach my $pair_feature (@$pair_features){ $all_genes{$pair_feature} = 1;}
1577                    }
1578                }
1579            }
1580        }
1581    
1582        # get the PCH to each of the genes
1583        my $pch_sets = [];
1584        my %pch_already;
1585        foreach my $gene_peg (keys %all_genes)
1586        {
1587            if ($pch_already{$gene_peg}){next;};
1588            my $gene_set = [$gene_peg];
1589            foreach my $pch_peg ($fig->in_pch_pin_with($gene_peg)) {
1590                $pch_peg =~ s/,.*$//;
1591                my $pch_genome = $fig->genome_of($pch_peg);
1592                if ( ($gene_peg ne $pch_peg) && ($all_genomes{$pch_genome})) {
1593                    push(@$gene_set,$pch_peg);
1594                    $pch_already{$pch_peg}=1;
1595                }
1596                $pch_already{$gene_peg}=1;
1597            }
1598            push(@$pch_sets,$gene_set);
1599        }
1600    
1601        #create a rank of the pch's
1602        my %pch_set_rank;
1603        my $order = 0;
1604        foreach my $set (@$pch_sets){
1605            my $count = scalar(@$set);
1606            $pch_set_rank{$order} = $count;
1607            $order++;
1608        }
1609    
1610        my %peg_rank;
1611        my $counter =  1;
1612        foreach my $pch_order (sort {$pch_set_rank{$b} <=> $pch_set_rank{$a}} keys %pch_set_rank){
1613            my $good_set = @$pch_sets[$pch_order];
1614            my $flag_set = 0;
1615            if (scalar (@$good_set) > 1)
1616            {
1617                foreach my $peg (@$good_set){
1618                    if ((!$peg_rank{$peg})){
1619                        $peg_rank{$peg} = $counter;
1620                        $flag_set = 1;
1621                    }
1622                }
1623                $counter++ if ($flag_set == 1);
1624            }
1625            else
1626            {
1627                foreach my $peg (@$good_set){
1628                    $peg_rank{$peg} = 100;
1629                }
1630            }
1631        }
1632    
1633    
1634    #    my $bbh_sets = [];
1635    #    my %already;
1636    #    foreach my $gene_key (keys(%all_genes)){
1637    #       if($already{$gene_key}){next;}
1638    #       my $gene_set = [$gene_key];
1639    #
1640    #       my $gene_key_genome = $fig->genome_of($gene_key);
1641    #
1642    #       foreach my $genome_key (keys(%all_genomes)){
1643    #           #next if ($gene_key_genome eq $genome_key);
1644    #           my $return = $fig->bbh_list($genome_key,[$gene_key]);
1645    #
1646    #           my $feature_list = $return->{$gene_key};
1647    #           foreach my $fl (@$feature_list){
1648    #               push(@$gene_set,$fl);
1649    #           }
1650    #       }
1651    #       $already{$gene_key} = 1;
1652    #       push(@$bbh_sets,$gene_set);
1653    #    }
1654    #
1655    #    my %bbh_set_rank;
1656    #    my $order = 0;
1657    #    foreach my $set (@$bbh_sets){
1658    #       my $count = scalar(@$set);
1659    #       $bbh_set_rank{$order} = $count;
1660    #       $order++;
1661    #    }
1662    #
1663    #    my %peg_rank;
1664    #    my $counter =  1;
1665    #    foreach my $bbh_order (sort {$bbh_set_rank{$b} <=> $bbh_set_rank{$a}} keys %bbh_set_rank){
1666    #       my $good_set = @$bbh_sets[$bbh_order];
1667    #       my $flag_set = 0;
1668    #       if (scalar (@$good_set) > 1)
1669    #       {
1670    #           foreach my $peg (@$good_set){
1671    #               if ((!$peg_rank{$peg})){
1672    #                   $peg_rank{$peg} = $counter;
1673    #                   $flag_set = 1;
1674    #               }
1675    #           }
1676    #           $counter++ if ($flag_set == 1);
1677    #       }
1678    #       else
1679    #       {
1680    #           foreach my $peg (@$good_set){
1681    #               $peg_rank{$peg} = 100;
1682    #           }
1683    #       }
1684    #    }
1685    
1686        foreach my $region (@$all_regions){
1687            my $sample_peg = @$region[0];
1688            my $region_genome = $fig->genome_of($sample_peg);
1689            my $region_gs = $fig->genus_species($region_genome);
1690            my $abbrev_name = $fig->abbrev($region_gs);
1691            my $line_config = { 'title' => $region_gs,
1692                                'short_title' => $abbrev_name,
1693                                'basepair_offset' => '0'
1694                                };
1695    
1696            my $offsetting = shift @start_array_region;
1697    
1698            my $line_data = [];
1699            foreach my $fid1 (@$region){
1700                my $element_hash;
1701                my $links_list = [];
1702                my $descriptions = [];
1703    
1704                my $color = $peg_rank{$fid1};
1705    
1706                # get subsystem information
1707                my $function = $fig->function_of($fid1);
1708                my $url_link = "http://seed-viewer.theseed.org/index.cgi?action=ShowAnnotation&prot=".$fid1;
1709    
1710                my $link;
1711                $link = {"link_title" => $fid1,
1712                         "link" => $url_link};
1713                push(@$links_list,$link);
1714    
1715                my @subsystems = $fig->peg_to_subsystems($fid1);
1716                foreach my $subsystem (@subsystems){
1717                    my $link;
1718                    $link = {"link" => "http://seed-viewer.theseed.org/index.cgi?action=ShowSubsystem&subsystem_name=$subsystem",
1719                             "link_title" => $subsystem};
1720                    push(@$links_list,$link);
1721                }
1722    
1723                my $description_function;
1724                $description_function = {"title" => "function",
1725                                         "value" => $function};
1726                push(@$descriptions,$description_function);
1727    
1728                my $description_ss;
1729                my $ss_string = join (",", @subsystems);
1730                $description_ss = {"title" => "subsystems",
1731                                   "value" => $ss_string};
1732                push(@$descriptions,$description_ss);
1733    
1734    
1735                my $fid_location = $fig->feature_location($fid1);
1736                if($fid_location =~/(.*)_(\d+)_(\d+)$/){
1737                    my($start,$stop);
1738                    $start = $2 - $offsetting;
1739                    $stop = $3 - $offsetting;
1740    
1741                    if (defined($reverse_flag{$region_genome})){
1742                        $start = $gd_window_size - $start;
1743                        $stop = $gd_window_size - $stop;
1744                    }
1745    
1746                    $element_hash = {
1747                        "title" => $fid1,
1748                        "start" => $start,
1749                        "end" =>  $stop,
1750                        "type"=> 'arrow',
1751                        "color"=> $color,
1752                        "zlayer" => "2",
1753                        "links_list" => $links_list,
1754                        "description" => $descriptions
1755                    };
1756                    push(@$line_data,$element_hash);
1757                }
1758            }
1759            $gd->add_line($line_data, $line_config);
1760        }
1761        return $gd;
1762    }
1763    
1764    

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.24

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3