[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.11, Thu Jun 21 21:15:23 2007 UTC revision 1.21, Thu Jun 28 22:08:05 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;  use HTML;
13    
14  1;  1;
# Line 118  Line 122 
122    
123  =item PFAM (dom)  =item PFAM (dom)
124    
125  =item SIGNALP (dom)  =item SIGNALP_CELLO_TMPRED (loc)
126    
127  =item  CELLO(loc)  =item PDB (seq)
128    
129  =item TMHMM (loc)  =item TMHMM (loc)
130    
# Line 287  Line 291 
291  }  }
292    
293    
294  =head3 display_method()  =head3 display()
   
 If available use the function specified here to display the "raw" observation.  
 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".  
295    
296  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.  will be different for each type
297    
298  =cut  =cut
299    
# Line 397  Line 397 
397  sub get_objects {  sub get_objects {
398      my ($self,$fid,$classes) = @_;      my ($self,$fid,$classes) = @_;
399    
   
400      my $objects = [];      my $objects = [];
401      my @matched_datasets=();      my @matched_datasets=();
402    
# Line 411  Line 410 
410          get_functional_coupling($fid,\@matched_datasets);          get_functional_coupling($fid,\@matched_datasets);
411      }      }
412      else{      else{
         #IPR,CDD,CELLO,PFAM,SIGNALP - attribute based  
413          my %domain_classes;          my %domain_classes;
414          my $identical_flag=0;          my $identical_flag=0;
415          my $pch_flag=0;          my $pch_flag=0;
416            my $location_flag = 0;
417          my $sims_flag=0;          my $sims_flag=0;
418            my $cluster_flag = 0;
419            my $pdb_flag = 0;
420          foreach my $class (@$classes){          foreach my $class (@$classes){
421              if($class =~ /(IPR|CDD|PFAM)/){              if($class =~ /(IPR|CDD|PFAM)/){
422                  $domain_classes{$class} = 1;                  $domain_classes{$class} = 1;
# Line 428  Line 429 
429              {              {
430                  $pch_flag = 1;                  $pch_flag = 1;
431              }              }
432                elsif ($class =~/(SIGNALP_CELLO_TMPRED)/)
433                {
434                    $location_flag = 1;
435                }
436              elsif ($class eq "SIM")              elsif ($class eq "SIM")
437              {              {
438                  $sims_flag = 1;                  $sims_flag = 1;
439              }              }
440                elsif ($class eq "CLUSTER")
441                {
442                    $cluster_flag = 1;
443                }
444                elsif ($class eq "PDB")
445                {
446                    $pdb_flag = 1;
447                }
448    
449          }          }
450    
451          if ($identical_flag ==1)          if ($identical_flag ==1)
# Line 450  Line 464 
464              get_sims_observations($fid,\@matched_datasets);              get_sims_observations($fid,\@matched_datasets);
465          }          }
466    
467          #add CELLO and SignalP later          if ($location_flag == 1)
468            {
469                get_attribute_based_location_observations($fid,\@matched_datasets);
470            }
471            if ($cluster_flag == 1)
472            {
473                get_cluster_observations($fid,\@matched_datasets);
474            }
475            if ($pdb_flag == 1)
476            {
477                get_pdb_observations($fid,\@matched_datasets);
478            }
479    
480    
481      }      }
482    
483      foreach my $dataset (@matched_datasets) {      foreach my $dataset (@matched_datasets) {
# Line 464  Line 491 
491          if ($dataset->{'class'} eq "IDENTICAL"){          if ($dataset->{'class'} eq "IDENTICAL"){
492              $object = Observation::Identical->new($dataset);              $object = Observation::Identical->new($dataset);
493          }          }
494            if ($dataset->{'class'} eq "SIGNALP_CELLO_TMPRED"){
495                $object = Observation::Location->new($dataset);
496            }
497          if ($dataset->{'class'} eq "SIM"){          if ($dataset->{'class'} eq "SIM"){
498              $object = Observation::Sims->new($dataset);              $object = Observation::Sims->new($dataset);
499          }          }
500            if ($dataset->{'class'} eq "CLUSTER"){
501                $object = Observation::Cluster->new($dataset);
502            }
503            if ($dataset->{'class'} eq "PDB"){
504                $object = Observation::PDB->new($dataset);
505            }
506    
507          push (@$objects, $object);          push (@$objects, $object);
508      }      }
509    
# Line 585  Line 622 
622      }      }
623  }  }
624    
625    sub get_attribute_based_location_observations{
626    
627        my ($fid,$datasets_ref) = (@_);
628        my $fig = new FIG;
629    
630        my $location_attributes = ['SignalP','CELLO','TMPRED'];
631    
632        my $dataset = {'type' => "loc", 'class' => 'SIGNALP_CELLO_TMPRED'};
633        foreach my $attr_ref ($fig->get_attributes($fid,$location_attributes)) {
634            my $key = @$attr_ref[1];
635            my @parts = split("::",$key);
636            my $sub_class = $parts[0];
637            my $sub_key = $parts[1];
638            my $value = @$attr_ref[2];
639            if($sub_class eq "SignalP"){
640                if($sub_key eq "cleavage_site"){
641                    my @value_parts = split(";",$value);
642                    $dataset->{'cleavage_prob'} = $value_parts[0];
643                    $dataset->{'cleavage_loc'} = $value_parts[1];
644                }
645                elsif($sub_key eq "signal_peptide"){
646                    $dataset->{'signal_peptide_score'} = $value;
647                }
648            }
649            elsif($sub_class eq "CELLO"){
650                $dataset->{'cello_location'} = $sub_key;
651                $dataset->{'cello_score'} = $value;
652            }
653            elsif($sub_class eq "TMPRED"){
654                my @value_parts = split(";",$value);
655                $dataset->{'tmpred_score'} = $value_parts[0];
656                $dataset->{'tmpred_locations'} = $value_parts[1];
657            }
658        }
659    
660        push (@{$datasets_ref} ,$dataset);
661    
662    }
663    
664    
665  =head3 get_attribute_based_evidence (internal)  =head3 get_attribute_based_evidence (internal)
666    
667  This method retrieves evidence from the attribute server  This method retrieves evidence from the attribute server
# Line 657  Line 734 
734      }      }
735  }  }
736    
737    =head3 get_pdb_observations() (internal)
738    
739    This methods sets the type and class for pdb observations
740    
741    =cut
742    
743    sub get_pdb_observations{
744        my ($fid,$datasets_ref) = (@_);
745    
746        my $fig = new FIG;
747    
748        print STDERR "get pdb obs called\n";
749        foreach my $attr_ref ($fig->get_attributes($fid,'PDB')) {
750    
751            my $key = @$attr_ref[1];
752            my($key1,$key2) =split("::",$key);
753            my $value = @$attr_ref[2];
754            my ($evalue,$location) = split(";",$value);
755    
756            if($evalue =~/(\d+)\.(\d+)/){
757                my $part2 = 1000 - $1;
758                my $part1 = $2/100;
759                $evalue = $part1."e-".$part2;
760            }
761    
762            my($start,$stop) =split("-",$location);
763    
764            my $url = @$attr_ref[3];
765            my $dataset = {'class' => 'PDB',
766                           'type' => 'seq' ,
767                           'acc' => $key2,
768                           'evalue' => $evalue,
769                           'start' => $start,
770                           'stop' => $stop
771                           };
772    
773            push (@{$datasets_ref} ,$dataset);
774        }
775    
776    }
777    
778    
779    
780    
781    =head3 get_cluster_observations() (internal)
782    
783    This methods sets the type and class for cluster observations
784    
785    =cut
786    
787    sub get_cluster_observations{
788        my ($fid,$datasets_ref) = (@_);
789    
790        my $dataset = {'class' => 'CLUSTER',
791                       'type' => 'fc'
792                       };
793        push (@{$datasets_ref} ,$dataset);
794    }
795    
796    
797  =head3 get_sims_observations() (internal)  =head3 get_sims_observations() (internal)
798    
799  This methods retrieves sims fills the internal data structures.  This methods retrieves sims fills the internal data structures.
# Line 965  Line 1102 
1102      return $self->{database};      return $self->{database};
1103  }  }
1104    
1105    ############################################################
1106    ############################################################
1107    package Observation::PDB;
1108    
1109    use base qw(Observation);
1110    
1111    sub new {
1112    
1113        my ($class,$dataset) = @_;
1114        my $self = $class->SUPER::new($dataset);
1115        $self->{acc} = $dataset->{'acc'};
1116        $self->{evalue} = $dataset->{'evalue'};
1117        $self->{start} = $dataset->{'start'};
1118        $self->{stop} = $dataset->{'stop'};
1119        bless($self,$class);
1120        return $self;
1121    }
1122    
1123    =head3 display()
1124    
1125    displays data stored in best_PDB attribute and in Ontology server for given PDB id
1126    
1127    =cut
1128    
1129    sub display{
1130        my ($self,$gd,$fid) = @_;
1131    
1132        my $dbmaster = DBMaster->new(-database =>'Ontology');
1133    
1134        print STDERR "PDB::display called\n";
1135    
1136        my $acc = $self->acc;
1137    
1138        print STDERR "acc:$acc\n";
1139        my ($pdb_description,$pdb_source,$pdb_ligand);
1140        my $pdb_objs = $dbmaster->pdb->get_objects( { 'id' => $acc } );
1141        if(!scalar(@$pdb_objs)){
1142            $pdb_description = "not available";
1143            $pdb_source = "not available";
1144            $pdb_ligand = "not available";
1145        }
1146        else{
1147            my $pdb_obj = $pdb_objs->[0];
1148            $pdb_description = $pdb_obj->description;
1149            $pdb_source = $pdb_obj->source;
1150            $pdb_ligand = $pdb_obj->ligand;
1151        }
1152    
1153        my $lines = [];
1154        my $line_data = [];
1155        my $line_config = { 'title' => "PDB hit for $fid",
1156                            'short_title' => "best PDB",
1157                            'basepair_offset' => '1' };
1158    
1159        my $fig = new FIG;
1160        my $seq = $fig->get_translation($fid);
1161        my $fid_stop = length($seq);
1162    
1163        my $fid_element_hash = {
1164            "title" => $fid,
1165            "start" => '1',
1166            "end" =>  $fid_stop,
1167            "color"=> '1',
1168            "zlayer" => '1'
1169            };
1170    
1171        push(@$line_data,$fid_element_hash);
1172    
1173        my $links_list = [];
1174        my $descriptions = [];
1175    
1176        my $name;
1177        $name = {"title" => 'id',
1178                 "value" => $acc};
1179        push(@$descriptions,$name);
1180    
1181        my $description;
1182        $description = {"title" => 'pdb description',
1183                        "value" => $pdb_description};
1184        push(@$descriptions,$description);
1185    
1186        my $score;
1187        $score = {"title" => "score",
1188                  "value" => $self->evalue};
1189        push(@$descriptions,$score);
1190    
1191        my $start_stop;
1192        my $start_stop_value = $self->start."_".$self->stop;
1193        $start_stop = {"title" => "start-stop",
1194                       "value" => $start_stop_value};
1195        push(@$descriptions,$start_stop);
1196    
1197        my $source;
1198        $source = {"title" => "source",
1199                  "value" => $pdb_source};
1200        push(@$descriptions,$source);
1201    
1202        my $ligand;
1203        $ligand = {"title" => "pdb ligand",
1204                   "value" => $pdb_ligand};
1205        push(@$descriptions,$ligand);
1206    
1207        my $link;
1208        my $link_url ="http://www.rcsb.org/pdb/explore/explore.do?structureId=".$acc;
1209    
1210        $link = {"link_title" => $acc,
1211                 "link" => $link_url};
1212        push(@$links_list,$link);
1213    
1214        my $pdb_element_hash = {
1215            "title" => "PDB homology",
1216            "start" => $self->start,
1217            "end" =>  $self->stop,
1218            "color"=> '6',
1219            "zlayer" => '3',
1220            "links_list" => $links_list,
1221            "description" => $descriptions};
1222    
1223        push(@$line_data,$pdb_element_hash);
1224        $gd->add_line($line_data, $line_config);
1225    
1226        return $gd;
1227    }
1228    
1229    1;
1230    
1231  ############################################################  ############################################################
1232  ############################################################  ############################################################
# Line 1124  Line 1385 
1385      my $links_list = [];      my $links_list = [];
1386      my $descriptions = [];      my $descriptions = [];
1387    
1388      my $description_function;      my $db_and_id = $thing->acc;
1389      $description_function = {"title" => $thing->class,      my ($db,$id) = split("::",$db_and_id);
                              "value" => $thing->acc};  
1390    
1391      push(@$descriptions,$description_function);      my $dbmaster = DBMaster->new(-database =>'Ontology');
1392    
1393        my ($name_title,$name_value,$description_title,$description_value);
1394        if($db eq "CDD"){
1395            my $cdd_objs = $dbmaster->cdd->get_objects( { 'id' => $id } );
1396            if(!scalar(@$cdd_objs)){
1397                $name_title = "name";
1398                $name_value = "not available";
1399                $description_title = "description";
1400                $description_value = "not available";
1401            }
1402            else{
1403                my $cdd_obj = $cdd_objs->[0];
1404                $name_title = "name";
1405                $name_value = $cdd_obj->term;
1406                $description_title = "description";
1407                $description_value = $cdd_obj->description;
1408            }
1409        }
1410    
1411        my $name;
1412        $name = {"title" => $name_title,
1413                 "value" => $name_value};
1414        push(@$descriptions,$name);
1415    
1416        my $description;
1417        $description = {"title" => $description_title,
1418                                 "value" => $description_value};
1419        push(@$descriptions,$description);
1420    
1421      my $score;      my $score;
1422      $score = {"title" => "score",      $score = {"title" => "score",
# Line 1136  Line 1424 
1424      push(@$descriptions,$score);      push(@$descriptions,$score);
1425    
1426      my $link_id;      my $link_id;
1427      if ($thing->acc =~/CDD::(\d+)/){      if ($thing->acc =~/\w+::(\d+)/){
1428          $link_id = $1;          $link_id = $1;
1429      }      }
1430    
1431      my $link;      my $link;
1432        my $link_url;
1433        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"}
1434        elsif($thing->class eq "PFAM"){$link_url = "http://www.sanger.ac.uk/cgi-bin/Pfam/getacc?$link_id"}
1435        else{$link_url = "NO_URL"}
1436    
1437      $link = {"link_title" => $thing->acc,      $link = {"link_title" => $thing->acc,
1438               "link" => "http://0-www.ncbi.nlm.nih.gov.library.vu.edu.au:80/Structure/cdd/cddsrv.cgi?uid=$link_id"};               "link" => $link_url};
1439      push(@$links_list,$link);      push(@$links_list,$link);
1440    
1441      my $element_hash = {      my $element_hash = {
# Line 1163  Line 1456 
1456    
1457  #########################################  #########################################
1458  #########################################  #########################################
1459    package Observation::Location;
1460    
1461    use base qw(Observation);
1462    
1463    sub new {
1464    
1465        my ($class,$dataset) = @_;
1466        my $self = $class->SUPER::new($dataset);
1467        $self->{cleavage_prob} = $dataset->{'cleavage_prob'};
1468        $self->{cleavage_loc} = $dataset->{'cleavage_loc'};
1469        $self->{signal_peptide_score} = $dataset->{'signal_peptide_score'};
1470        $self->{cello_location} = $dataset->{'cello_location'};
1471        $self->{cello_score} = $dataset->{'cello_score'};
1472        $self->{tmpred_score} = $dataset->{'tmpred_score'};
1473        $self->{tmpred_locations} = $dataset->{'tmpred_locations'};
1474    
1475        bless($self,$class);
1476        return $self;
1477    }
1478    
1479    sub display {
1480        my ($thing,$gd,$fid) = @_;
1481    
1482        my $fig= new FIG;
1483        my $length = length($fig->get_translation($fid));
1484    
1485        my $cleavage_prob;
1486        if($thing->cleavage_prob){$cleavage_prob = $thing->cleavage_prob;}
1487        my ($cleavage_loc_begin,$cleavage_loc_end) = split("-",$thing->cleavage_loc);
1488        my $signal_peptide_score = $thing->signal_peptide_score;
1489        my $cello_location = $thing->cello_location;
1490        my $cello_score = $thing->cello_score;
1491        my $tmpred_score = $thing->tmpred_score;
1492        my @tmpred_locations = split(",",$thing->tmpred_locations);
1493    
1494        my $lines = [];
1495        my $line_config = { 'title' => 'Localization Evidence',
1496                            'short_title' => 'Local',
1497                            'basepair_offset' => '1' };
1498    
1499        #color is
1500        my $color = "5";
1501    
1502        my $line_data = [];
1503    
1504        if($cello_location){
1505            my $cello_descriptions = [];
1506            my $description_cello_location = {"title" => 'Best Cello Location',
1507                                              "value" => $cello_location};
1508    
1509            push(@$cello_descriptions,$description_cello_location);
1510    
1511            my $description_cello_score = {"title" => 'Cello Score',
1512                                           "value" => $cello_score};
1513    
1514            push(@$cello_descriptions,$description_cello_score);
1515    
1516            my $element_hash = {
1517                "title" => "CELLO",
1518                "start" => "1",
1519                "end" =>  $length + 1,
1520                "color"=> $color,
1521                "type" => 'box',
1522                "zlayer" => '2',
1523                "description" => $cello_descriptions};
1524    
1525            push(@$line_data,$element_hash);
1526        }
1527    
1528        my $color = "6";
1529        #if(0){
1530        if($tmpred_score){
1531            foreach my $tmpred (@tmpred_locations){
1532                my $descriptions = [];
1533                my ($begin,$end) =split("-",$tmpred);
1534                my $description_tmpred_score = {"title" => 'TMPRED score',
1535                                 "value" => $tmpred_score};
1536    
1537                push(@$descriptions,$description_tmpred_score);
1538    
1539                my $element_hash = {
1540                "title" => "transmembrane location",
1541                "start" => $begin + 1,
1542                "end" =>  $end + 1,
1543                "color"=> $color,
1544                "zlayer" => '5',
1545                "type" => 'smallbox',
1546                "description" => $descriptions};
1547    
1548                push(@$line_data,$element_hash);
1549            }
1550        }
1551    
1552        my $color = "1";
1553        if($signal_peptide_score){
1554            my $descriptions = [];
1555            my $description_signal_peptide_score = {"title" => 'signal peptide score',
1556                                                    "value" => $signal_peptide_score};
1557    
1558            push(@$descriptions,$description_signal_peptide_score);
1559    
1560            my $description_cleavage_prob = {"title" => 'cleavage site probability',
1561                                             "value" => $cleavage_prob};
1562    
1563            push(@$descriptions,$description_cleavage_prob);
1564    
1565            my $element_hash = {
1566                "title" => "SignalP",
1567                "start" => $cleavage_loc_begin - 2,
1568                "end" =>  $cleavage_loc_end + 3,
1569                "type" => 'bigbox',
1570                "color"=> $color,
1571                "zlayer" => '10',
1572                "description" => $descriptions};
1573    
1574            push(@$line_data,$element_hash);
1575        }
1576    
1577        $gd->add_line($line_data, $line_config);
1578    
1579        return ($gd);
1580    
1581    }
1582    
1583    sub cleavage_loc {
1584      my ($self) = @_;
1585    
1586      return $self->{cleavage_loc};
1587    }
1588    
1589    sub cleavage_prob {
1590      my ($self) = @_;
1591    
1592      return $self->{cleavage_prob};
1593    }
1594    
1595    sub signal_peptide_score {
1596      my ($self) = @_;
1597    
1598      return $self->{signal_peptide_score};
1599    }
1600    
1601    sub tmpred_score {
1602      my ($self) = @_;
1603    
1604      return $self->{tmpred_score};
1605    }
1606    
1607    sub tmpred_locations {
1608      my ($self) = @_;
1609    
1610      return $self->{tmpred_locations};
1611    }
1612    
1613    sub cello_location {
1614      my ($self) = @_;
1615    
1616      return $self->{cello_location};
1617    }
1618    
1619    sub cello_score {
1620      my ($self) = @_;
1621    
1622      return $self->{cello_score};
1623    }
1624    
1625    
1626    #########################################
1627    #########################################
1628  package Observation::Sims;  package Observation::Sims;
1629    
1630  use base qw(Observation);  use base qw(Observation);
# Line 1294  Line 1756 
1756  }  }
1757    
1758  sub html_enc { $_ = $_[0]; s/\&/&amp;/g; s/\>/&gt;/g; s/\</&lt;/g; $_ }  sub html_enc { $_ = $_[0]; s/\&/&amp;/g; s/\>/&gt;/g; s/\</&lt;/g; $_ }
1759    
1760    
1761    
1762    ############################
1763    package Observation::Cluster;
1764    
1765    use base qw(Observation);
1766    
1767    sub new {
1768    
1769        my ($class,$dataset) = @_;
1770        my $self = $class->SUPER::new($dataset);
1771    
1772        bless($self,$class);
1773        return $self;
1774    }
1775    
1776    sub display {
1777        my ($self,$gd, $fid) = @_;
1778    
1779        my $fig = new FIG;
1780        my $all_regions = [];
1781    
1782        #get the organism genome
1783        my $target_genome = $fig->genome_of($fid);
1784    
1785        # get location of the gene
1786        my $data = $fig->feature_location($fid);
1787        my ($contig, $beg, $end);
1788    
1789        if ($data =~ /(.*)_(\d+)_(\d+)$/){
1790            $contig = $1;
1791            $beg = $2;
1792            $end = $3;
1793        }
1794    
1795        my ($region_start, $region_end);
1796        if ($beg < $end)
1797        {
1798            $region_start = $beg - 4000;
1799            $region_end = $end+4000;
1800        }
1801        else
1802        {
1803            $region_start = $end-4000;
1804            $region_end = $beg+4000;
1805        }
1806    
1807        # call genes in region
1808        my ($target_gene_features, $reg_beg, $reg_end) = $fig->genes_in_region($target_genome, $contig, $region_start, $region_end);
1809        push(@$all_regions,$target_gene_features);
1810        my (@start_array_region);
1811        push (@start_array_region, $region_start);
1812    
1813        my %all_genes;
1814        my %all_genomes;
1815        foreach my $feature (@$target_gene_features){ $all_genes{$feature} = 1;}
1816        my $compare_regions_flag = 1; # set it for compare regions view (0 -> no view, 1-> yes view)
1817        my $functional_coupling_flag = 0; # set functional coupling for view (0 -> no view, 1-> yes view)
1818    
1819        if ($functional_coupling_flag == 1)
1820        {
1821            my @coup = grep { $_->[1]} $fig->coupling_and_evidence($fid,5000,1e-10,4,1);
1822    
1823            my $coup_count = 0;
1824    
1825            foreach my $pair (@{$coup[0]->[2]}) {
1826                #   last if ($coup_count > 10);
1827                my ($peg1,$peg2) = @$pair;
1828    
1829                my $location = $fig->feature_location($peg1);
1830                my ($pair_contig,$pair_beg,$pair_end,$pair_region_start,$pair_region_stop,$pair_genome);
1831                if($location =~/(.*)_(\d+)_(\d+)$/){
1832                    $pair_contig = $1;
1833                    $pair_beg = $2;
1834                    $pair_end = $3;
1835                    if ($pair_beg < $pair_end)
1836                    {
1837                        $pair_region_start = $pair_beg - 4000;
1838                        $pair_region_stop = $pair_end+4000;
1839                    }
1840                    else
1841                    {
1842                        $pair_region_start = $pair_end-4000;
1843                        $pair_region_stop = $pair_beg+4000;
1844                    }
1845    
1846                    push (@start_array_region, $pair_region_start);
1847    
1848                    $pair_genome = $fig->genome_of($peg1);
1849                    $all_genomes{$pair_genome} = 1;
1850                    my ($pair_features) = $fig->genes_in_region($pair_genome, $pair_contig, $pair_region_start, $pair_region_stop);
1851                    push(@$all_regions,$pair_features);
1852                    foreach my $pair_feature (@$pair_features){ $all_genes{$pair_feature} = 1;}
1853                }
1854                $coup_count++;
1855            }
1856        }
1857    
1858        if ($compare_regions_flag)
1859        {
1860            # make a hash of genomes that are phylogenetically close
1861            #my $close_threshold = ".26";
1862            #my @genomes = $fig->genomes('complete');
1863            #my %close_genomes = ();
1864            #foreach my $compared_genome (@genomes)
1865            #{
1866            #    my $dist = $fig->crude_estimate_of_distance($target_genome,$compared_genome);
1867            #    #$close_genomes{$compared_genome} = $dist;
1868            #    if ($dist <= $close_threshold)
1869            #    {
1870            #       $all_genomes{$compared_genome} = 1;
1871            #    }
1872            #}
1873            $all_genomes{"216592.1"} = 1;
1874            $all_genomes{"79967.1"} = 1;
1875            $all_genomes{"199310.1"} = 1;
1876            $all_genomes{"216593.1"} = 1;
1877            $all_genomes{"155864.1"} = 1;
1878            $all_genomes{"83334.1"} = 1;
1879            $all_genomes{"316407.3"} = 1;
1880    
1881            foreach my $comp_genome (keys %all_genomes){
1882                my $return = $fig->bbh_list($comp_genome,[$fid]);
1883                my $feature_list = $return->{$fid};
1884                foreach my $peg1 (@$feature_list){
1885                    my $location = $fig->feature_location($peg1);
1886                    my ($pair_contig,$pair_beg,$pair_end,$pair_region_start,$pair_region_stop,$pair_genome);
1887                    if($location =~/(.*)_(\d+)_(\d+)$/){
1888                        $pair_contig = $1;
1889                        $pair_beg = $2;
1890                        $pair_end = $3;
1891                        if ($pair_beg < $pair_end)
1892                        {
1893                            $pair_region_start = $pair_beg - 4000;
1894                            $pair_region_stop = $pair_end + 4000;
1895                            print STDERR "begFIG: $peg1, START:$pair_region_start, END: $pair_region_stop";
1896                        }
1897                        else
1898                        {
1899                            $pair_region_start = $pair_end-4000;
1900                            $pair_region_stop = $pair_beg+4000;
1901                            print STDERR "endFIG: $peg1, START:$pair_region_start, END: $pair_region_stop";
1902                        }
1903    
1904                        push (@start_array_region, $pair_region_start);
1905    
1906                        $pair_genome = $fig->genome_of($peg1);
1907                        $all_genomes{$pair_genome} = 1;
1908                        my ($pair_features) = $fig->genes_in_region($pair_genome, $pair_contig, $pair_region_start, $pair_region_stop);
1909                        push(@$all_regions,$pair_features);
1910                        foreach my $pair_feature (@$pair_features){ $all_genes{$pair_feature} = 1;}
1911                    }
1912                }
1913            }
1914        }
1915    
1916        # get the PCH to each of the genes
1917        my $pch_sets = [];
1918        my %pch_already;
1919        foreach my $gene_peg (keys %all_genes)
1920        {
1921            if ($pch_already{$gene_peg}){next;};
1922            my $gene_set = [$gene_peg];
1923            foreach my $pch_peg ($fig->in_pch_pin_with($gene_peg)) {
1924                $pch_peg =~ s/,.*$//;
1925                my $pch_genome = $fig->genome_of($pch_peg);
1926                if ( ($gene_peg ne $pch_peg) && ($all_genomes{$pch_genome})) {
1927                    push(@$gene_set,$pch_peg);
1928                    $pch_already{$pch_peg}=1;
1929                }
1930                $pch_already{$gene_peg}=1;
1931            }
1932            push(@$pch_sets,$gene_set);
1933        }
1934    
1935        #create a rank of the pch's
1936        my %pch_set_rank;
1937        my $order = 0;
1938        foreach my $set (@$pch_sets){
1939            my $count = scalar(@$set);
1940            $pch_set_rank{$order} = $count;
1941            $order++;
1942        }
1943    
1944        my %peg_rank;
1945        my $counter =  1;
1946        foreach my $pch_order (sort {$pch_set_rank{$b} <=> $pch_set_rank{$a}} keys %pch_set_rank){
1947            my $good_set = @$pch_sets[$pch_order];
1948            my $flag_set = 0;
1949            if (scalar (@$good_set) > 1)
1950            {
1951                foreach my $peg (@$good_set){
1952                    if ((!$peg_rank{$peg})){
1953                        $peg_rank{$peg} = $counter;
1954                        $flag_set = 1;
1955                    }
1956                }
1957                $counter++ if ($flag_set == 1);
1958            }
1959            else
1960            {
1961                foreach my $peg (@$good_set){
1962                    $peg_rank{$peg} = 100;
1963                }
1964            }
1965        }
1966    
1967    
1968    #    my $bbh_sets = [];
1969    #    my %already;
1970    #    foreach my $gene_key (keys(%all_genes)){
1971    #       if($already{$gene_key}){next;}
1972    #       my $gene_set = [$gene_key];
1973    #
1974    #       my $gene_key_genome = $fig->genome_of($gene_key);
1975    #
1976    #       foreach my $genome_key (keys(%all_genomes)){
1977    #           #next if ($gene_key_genome eq $genome_key);
1978    #           my $return = $fig->bbh_list($genome_key,[$gene_key]);
1979    #
1980    #           my $feature_list = $return->{$gene_key};
1981    #           foreach my $fl (@$feature_list){
1982    #               push(@$gene_set,$fl);
1983    #           }
1984    #       }
1985    #       $already{$gene_key} = 1;
1986    #       push(@$bbh_sets,$gene_set);
1987    #    }
1988    #
1989    #    my %bbh_set_rank;
1990    #    my $order = 0;
1991    #    foreach my $set (@$bbh_sets){
1992    #       my $count = scalar(@$set);
1993    #       $bbh_set_rank{$order} = $count;
1994    #       $order++;
1995    #    }
1996    #
1997    #    my %peg_rank;
1998    #    my $counter =  1;
1999    #    foreach my $bbh_order (sort {$bbh_set_rank{$b} <=> $bbh_set_rank{$a}} keys %bbh_set_rank){
2000    #       my $good_set = @$bbh_sets[$bbh_order];
2001    #       my $flag_set = 0;
2002    #       if (scalar (@$good_set) > 1)
2003    #       {
2004    #           foreach my $peg (@$good_set){
2005    #               if ((!$peg_rank{$peg})){
2006    #                   $peg_rank{$peg} = $counter;
2007    #                   $flag_set = 1;
2008    #               }
2009    #           }
2010    #           $counter++ if ($flag_set == 1);
2011    #       }
2012    #       else
2013    #       {
2014    #           foreach my $peg (@$good_set){
2015    #               $peg_rank{$peg} = 100;
2016    #           }
2017    #       }
2018    #    }
2019    
2020        foreach my $region (@$all_regions){
2021            my $sample_peg = @$region[0];
2022            my $region_genome = $fig->genome_of($sample_peg);
2023            my $region_gs = $fig->genus_species($region_genome);
2024            my $abbrev_name = $fig->abbrev($region_gs);
2025            my $line_config = { 'title' => $region_gs,
2026                                'short_title' => $abbrev_name,
2027                                'basepair_offset' => '0'
2028                                };
2029    
2030            my $offset = shift @start_array_region;
2031    
2032            my $line_data = [];
2033            foreach my $fid1 (@$region){
2034                my $element_hash;
2035                my $links_list = [];
2036                my $descriptions = [];
2037    
2038                my $color = $peg_rank{$fid1};
2039    
2040                # get subsystem information
2041                my $function = $fig->function_of($fid1);
2042                my $url_link = "http://seed-viewer.theseed.org/index.cgi?action=ShowAnnotation&prot=".$fid1;
2043    
2044                my $link;
2045                $link = {"link_title" => $fid1,
2046                         "link" => $url_link};
2047                push(@$links_list,$link);
2048    
2049                my @subsystems = $fig->peg_to_subsystems($fid1);
2050                foreach my $subsystem (@subsystems){
2051                    my $link;
2052                    $link = {"link" => "http://seed-viewer.theseed.org/index.cgi?action=ShowSubsystem&subsystem_name=$subsystem",
2053                             "link_title" => $subsystem};
2054                    push(@$links_list,$link);
2055                }
2056    
2057                my $description_function;
2058                $description_function = {"title" => "function",
2059                                         "value" => $function};
2060                push(@$descriptions,$description_function);
2061    
2062                my $description_ss;
2063                my $ss_string = join (",", @subsystems);
2064                $description_ss = {"title" => "subsystems",
2065                                   "value" => $ss_string};
2066                push(@$descriptions,$description_ss);
2067    
2068    
2069                my $fid_location = $fig->feature_location($fid1);
2070                if($fid_location =~/(.*)_(\d+)_(\d+)$/){
2071                    my($start,$stop);
2072                    $start = $2 - $offset;
2073                    $stop = $3 - $offset;
2074                    $element_hash = {
2075                        "title" => $fid1,
2076                        "start" => $start,
2077                        "end" =>  $stop,
2078                        "type"=> 'arrow',
2079                        "color"=> $color,
2080                        "zlayer" => "2",
2081                        "links_list" => $links_list,
2082                        "description" => $descriptions
2083                    };
2084                    push(@$line_data,$element_hash);
2085                }
2086            }
2087            $gd->add_line($line_data, $line_config);
2088        }
2089        return $gd;
2090    }
2091    
2092    

Legend:
Removed from v.1.11  
changed lines
  Added in v.1.21

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3