[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.20, Wed Jun 27 22:14:01 2007 UTC revision 1.21, Thu Jun 28 22:08:05 2007 UTC
# Line 1800  Line 1800 
1800      }      }
1801      else      else
1802      {      {
1803          $region_end = $end+4000;          $region_start = $end-4000;
1804          $region_start = $beg-4000;          $region_end = $beg+4000;
1805      }      }
1806    
1807      # call genes in region      # call genes in region
# Line 1813  Line 1813 
1813      my %all_genes;      my %all_genes;
1814      my %all_genomes;      my %all_genomes;
1815      foreach my $feature (@$target_gene_features){ $all_genes{$feature} = 1;}      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);      my @coup = grep { $_->[1]} $fig->coupling_and_evidence($fid,5000,1e-10,4,1);
1822    
1823      my $coup_count = 0;      my $coup_count = 0;
1824    
1825      foreach my $pair (@{$coup[0]->[2]}) {      foreach my $pair (@{$coup[0]->[2]}) {
1826          last if ($coup_count > 10);              #   last if ($coup_count > 10);
1827          my ($peg1,$peg2) = @$pair;          my ($peg1,$peg2) = @$pair;
1828    
1829          my $location = $fig->feature_location($peg1);          my $location = $fig->feature_location($peg1);
# Line 1835  Line 1839 
1839              }              }
1840              else              else
1841              {              {
1842                  $pair_region_stop = $pair_end+4000;                      $pair_region_start = $pair_end-4000;
1843                  $pair_region_start = $pair_beg-4000;                      $pair_region_stop = $pair_beg+4000;
1844              }              }
1845    
1846              push (@start_array_region, $pair_region_start);              push (@start_array_region, $pair_region_start);
# Line 1849  Line 1853 
1853          }          }
1854          $coup_count++;          $coup_count++;
1855      }      }
1856        }
1857    
1858      my $bbh_sets = [];      if ($compare_regions_flag)
1859      my %already;      {
1860      foreach my $gene_key (keys(%all_genes)){          # make a hash of genomes that are phylogenetically close
1861          if($already{$gene_key}){next;}          #my $close_threshold = ".26";
1862          my $gene_set = [$gene_key];          #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          my $gene_key_genome = $fig->genome_of($gene_key);                      push (@start_array_region, $pair_region_start);
1905    
1906          foreach my $genome_key (keys(%all_genomes)){                      $pair_genome = $fig->genome_of($peg1);
1907              #next if ($gene_key_genome eq $genome_key);                      $all_genomes{$pair_genome} = 1;
1908              my $return = $fig->bbh_list($genome_key,[$gene_key]);                      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              my $feature_list = $return->{$gene_key};      # get the PCH to each of the genes
1917              foreach my $fl (@$feature_list){      my $pch_sets = [];
1918                  push(@$gene_set,$fl);      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          $already{$gene_key} = 1;          push(@$pch_sets,$gene_set);
         push(@$bbh_sets,$gene_set);  
1933      }      }
1934    
1935      my %bbh_set_rank;      #create a rank of the pch's
1936        my %pch_set_rank;
1937      my $order = 0;      my $order = 0;
1938      foreach my $set (@$bbh_sets){      foreach my $set (@$pch_sets){
1939          my $count = scalar(@$set);          my $count = scalar(@$set);
1940          $bbh_set_rank{$order} = $count;          $pch_set_rank{$order} = $count;
1941          $order++;          $order++;
1942      }      }
1943    
1944      my %peg_rank;      my %peg_rank;
1945      my $counter =  1;      my $counter =  1;
1946      foreach my $bbh_order (sort {$bbh_set_rank{$b} <=> $bbh_set_rank{$a}} keys %bbh_set_rank){      foreach my $pch_order (sort {$pch_set_rank{$b} <=> $pch_set_rank{$a}} keys %pch_set_rank){
1947          my $good_set = @$bbh_sets[$bbh_order];          my $good_set = @$pch_sets[$pch_order];
1948          my $flag_set = 0;          my $flag_set = 0;
1949          if (scalar (@$good_set) > 1)          if (scalar (@$good_set) > 1)
1950          {          {
# Line 1902  Line 1964 
1964          }          }
1965      }      }
1966    
1967      open (FH, ">$FIG_Config::temp/good_sets.txt");  
1968      foreach my $pr (sort {$peg_rank{$a} <=> $peg_rank{$b}} keys(%peg_rank)){ print FH "rank:$peg_rank{$pr}\tpr:$pr\n";}  #    my $bbh_sets = [];
1969      close (FH);  #    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){      foreach my $region (@$all_regions){
2021          my $sample_peg = @$region[0];          my $sample_peg = @$region[0];
# Line 1925  Line 2036 
2036              my $descriptions = [];              my $descriptions = [];
2037    
2038              my $color = $peg_rank{$fid1};              my $color = $peg_rank{$fid1};
             if ($color == 1) {  
                 print STDERR "PEG: $fid1, RANK: $color";  
             }  
2039    
2040              # get subsystem information              # get subsystem information
2041              my $function = $fig->function_of($fid1);              my $function = $fig->function_of($fid1);
# Line 1961  Line 2069 
2069              my $fid_location = $fig->feature_location($fid1);              my $fid_location = $fig->feature_location($fid1);
2070              if($fid_location =~/(.*)_(\d+)_(\d+)$/){              if($fid_location =~/(.*)_(\d+)_(\d+)$/){
2071                  my($start,$stop);                  my($start,$stop);
2072                  if ($2 < $3){$start = $2; $stop = $3;}                  $start = $2 - $offset;
2073                  else{$stop = $2; $start = $3;}                  $stop = $3 - $offset;
                 $start = $start - $offset;  
                 $stop = $stop - $offset;  
2074                  $element_hash = {                  $element_hash = {
2075                      "title" => $fid1,                      "title" => $fid1,
2076                      "start" => $start,                      "start" => $start,

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3