[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.37, Tue Sep 4 18:34:13 2007 UTC revision 1.38, Mon Sep 10 15:10:04 2007 UTC
# Line 438  Line 438 
438      my ($observation, $fid) = @_;      my ($observation, $fid) = @_;
439      my $fig = new FIG;      my $fig = new FIG;
440      my %families;      my %families;
441      my @sims= $fig->nsims($fid,20000,10,"all");      my @sims= $fig->nsims($fid,20000,10,"fig");
442    
443      foreach my $sim (@sims){      foreach my $sim (@sims){
444          next if ($sim->[1] !~ /fig\|/);          next if ($sim->[1] !~ /fig\|/);
445          my $genome = $fig->genome_of($sim->[1]);          my $genome = $fig->genome_of($sim->[1]);
446          my $taxonomy = $fig->taxonomy_of($fig->genome_of($sim->[1]));          my $taxonomy = $fig->taxonomy_of($fig->genome_of($sim->[1]));
447          my $parent_tax = "Root";          my $parent_tax = "Root";
448            my @currLineage = ($parent_tax);
449          foreach my $tax (split(/\; /, $taxonomy)){          foreach my $tax (split(/\; /, $taxonomy)){
450              push (@{$families{children}{$parent_tax}}, $tax);              push (@{$families{children}{$parent_tax}}, $tax);
451                push (@currLineage, $tax);
452              $families{parent}{$tax} = $parent_tax;              $families{parent}{$tax} = $parent_tax;
453                $families{lineage}{$tax} = join(";", @currLineage);
454              $parent_tax = $tax;              $parent_tax = $tax;
455          }          }
456      }      }
# Line 642  Line 645 
645    
646      my ($fid,$datasets_ref) = (@_);      my ($fid,$datasets_ref) = (@_);
647      my $fig = new FIG;      my $fig = new FIG;
648      my @sims= $fig->nsims($fid,500,1e-20,"all");      my @sims= $fig->nsims($fid,500,10,"fig");
649      my ($dataset);      my ($dataset);
650    
651      my %id_list;      my %id_list;
# Line 2149  Line 2152 
2152  }  }
2153    
2154  sub display {  sub display {
2155      my ($self,$gd) = @_;      my ($self,$gd,$selected_taxonomies) = @_;
2156    
2157      my $fid = $self->fig_id;      my $fid = $self->fig_id;
2158      my $compare_or_coupling = $self->context;      my $compare_or_coupling = $self->context;
2159      my $gd_window_size = $gd->window_size;      my $gd_window_size = $gd->window_size;
2160      my $fig = new FIG;      my $fig = new FIG;
2161      my $all_regions = [];      my $all_regions = [];
2162        my $gene_associations={};
2163    
2164      #get the organism genome      #get the organism genome
2165      my $target_genome = $fig->genome_of($fid);      my $target_genome = $fig->genome_of($fid);
2166        $gene_associations->{$fid}->{"organism"} = $target_genome;
2167        $gene_associations->{$fid}->{"main_gene"} = $fid;
2168        $gene_associations->{$fid}->{"reverse_flag"} = 0;
2169    
2170      # get location of the gene      # get location of the gene
2171      my $data = $fig->feature_location($fid);      my $data = $fig->feature_location($fid);
# Line 2185  Line 2192 
2192          $region_end = $beg+4000;          $region_end = $beg+4000;
2193          $offset = ($3+(($2-$3)/2))-($gd_window_size/2);          $offset = ($3+(($2-$3)/2))-($gd_window_size/2);
2194          $reverse_flag{$target_genome} = $fid;          $reverse_flag{$target_genome} = $fid;
2195            $gene_associations->{$fid}->{"reverse_flag"} = 1;
2196      }      }
2197    
2198      # call genes in region      # call genes in region
# Line 2195  Line 2203 
2203    
2204      my %all_genes;      my %all_genes;
2205      my %all_genomes;      my %all_genomes;
2206      foreach my $feature (@$target_gene_features){ $all_genes{$feature} = $fid;}      foreach my $feature (@$target_gene_features){ $all_genes{$feature} = $fid; $gene_associations->{$feature}->{"main_gene"}=$fid;}
2207    
2208      if ($compare_or_coupling eq "diverse")      if ($compare_or_coupling eq "diverse")
2209      {      {
# Line 2219  Line 2227 
2227                  {                  {
2228                      $pair_region_start = $pair_beg - 4000;                      $pair_region_start = $pair_beg - 4000;
2229                      $pair_region_stop = $pair_end+4000;                      $pair_region_stop = $pair_end+4000;
2230                      $offset = ($2+(($3-$2)/2))-($gd_window_size/2);                      $offset = ($pair_beg+(($pair_end-$pair_beg)/2))-($gd_window_size/2);
2231                  }                  }
2232                  else                  else
2233                  {                  {
2234                      $pair_region_start = $pair_end-4000;                      $pair_region_start = $pair_end-4000;
2235                      $pair_region_stop = $pair_beg+4000;                      $pair_region_stop = $pair_beg+4000;
2236                      $offset = ($3+(($2-$3)/2))-($gd_window_size/2);                      $offset = ($pair_end+(($pair_beg-$pair_end)/2))-($gd_window_size/2);
2237                      $reverse_flag{$pair_genome} = $peg1;                      $reverse_flag{$pair_genome} = $peg1;
2238                  }                  }
2239    
# Line 2241  Line 2249 
2249      }      }
2250      elsif ($compare_or_coupling eq "sims"){      elsif ($compare_or_coupling eq "sims"){
2251          # get the selected boxes          # get the selected boxes
2252          my @selected_taxonomoy = ("Deltaproteobacteria", "Vibrionales", "Viridiplantae");          #my @selected_taxonomy = ("Streptococcaceae", "Enterobacteriales");
2253            my @selected_taxonomy = @$selected_taxonomies;
2254    
2255          # get the similarities and store only the ones that match the lineages selected          # get the similarities and store only the ones that match the lineages selected
2256          my @selected_sims;          my @selected_sims;
2257          my @sims= $fig->nsims($fid,20000,10,"all");          my @sims= $fig->nsims($fid,20000,10,"fig");
2258    
2259            if (@selected_taxonomy > 0){
2260          foreach my $sim (@sims){          foreach my $sim (@sims){
2261              next if ($sim->[1] !~ /fig\|/);              next if ($sim->[1] !~ /fig\|/);
2262              my $genome = $fig->genome_of($sim->[1]);              my $genome = $fig->genome_of($sim->[1]);
# Line 2259  Line 2269 
2269              my %saw;              my %saw;
2270              @selected_sims = grep(!$saw{$_}++, @selected_sims);              @selected_sims = grep(!$saw{$_}++, @selected_sims);
2271          }          }
2272            }
2273    
2274          # get the gene context for the sorted matches          # get the gene context for the sorted matches
2275          foreach my $sim_fid(@selected_sims){          foreach my $sim_fid(@selected_sims){
2276              #get the organism genome              #get the organism genome
2277              my $sim_genome = $fig->genome_of($sim_fid);              my $sim_genome = $fig->genome_of($sim_fid);
2278                $gene_associations->{$sim_fid}->{"organism"} = $sim_genome;
2279                $gene_associations->{$sim_fid}->{"main_gene"} = $sim_fid;
2280                $gene_associations->{$sim_fid}->{"reverse_flag"} = 0;
2281    
2282              # get location of the gene              # get location of the gene
2283              my $data = $fig->feature_location($sim_fid);              my $data = $fig->feature_location($sim_fid);
2284              my ($contig, $beg, $end);              my ($contig, $beg, $end);
             my %reverse_flag;  
2285    
2286              if ($data =~ /(.*)_(\d+)_(\d+)$/){              if ($data =~ /(.*)_(\d+)_(\d+)$/){
2287                  $contig = $1;                  $contig = $1;
# Line 2282  Line 2295 
2295              {              {
2296                  $region_start = $beg - 4000;                  $region_start = $beg - 4000;
2297                  $region_end = $end+4000;                  $region_end = $end+4000;
2298                  $offset = ($2+(($3-$2)/2))-($gd_window_size/2);                  $offset = ($beg+(($end-$beg)/2))-($gd_window_size/2);
2299              }              }
2300              else              else
2301              {              {
2302                  $region_start = $end-4000;                  $region_start = $end-4000;
2303                  $region_end = $beg+4000;                  $region_end = $beg+4000;
2304                  $offset = ($3+(($2-$3)/2))-($gd_window_size/2);                  $offset = ($end+(($beg-$end)/2))-($gd_window_size/2);
2305                  $reverse_flag{$target_genome} = $sim_fid;                  $reverse_flag{$sim_genome} = $sim_fid;
2306                    $gene_associations->{$sim_fid}->{"reverse_flag"} = 1;
2307              }              }
2308    
2309              # call genes in region              # call genes in region
2310              my ($sim_gene_features, $reg_beg, $reg_end) = $fig->genes_in_region($sim_genome, $contig, $region_start, $region_end);              my ($sim_gene_features, $reg_beg, $reg_end) = $fig->genes_in_region($sim_genome, $contig, $region_start, $region_end);
2311              push(@$all_regions,$sim_gene_features);              push(@$all_regions,$sim_gene_features);
             my (@start_array_region);  
2312              push (@start_array_region, $offset);              push (@start_array_region, $offset);
2313                foreach my $feature (@$sim_gene_features){ $all_genes{$feature} = $sim_fid;$gene_associations->{$feature}->{"main_gene"}=$sim_fid;}
2314              my %all_genes;              $all_genomes{$sim_genome} = 1;
             my %all_genomes;  
             foreach my $feature (@$sim_gene_features){ $all_genes{$feature} = $sim_fid;}  
         }  
2315      }      }
2316    
     elsif ($compare_or_coupling eq "close")  
     {  
         # make a hash of genomes that are phylogenetically close  
         #my $close_threshold = ".26";  
         #my @genomes = $fig->genomes('complete');  
         #my %close_genomes = ();  
         #foreach my $compared_genome (@genomes)  
         #{  
         #    my $dist = $fig->crude_estimate_of_distance($target_genome,$compared_genome);  
         #    #$close_genomes{$compared_genome} = $dist;  
         #    if ($dist <= $close_threshold)  
         #    {  
         #       $all_genomes{$compared_genome} = 1;  
         #    }  
         #}  
         $all_genomes{"216592.1"} = 1;  
         $all_genomes{"79967.1"} = 1;  
         $all_genomes{"199310.1"} = 1;  
         $all_genomes{"216593.1"} = 1;  
         $all_genomes{"155864.1"} = 1;  
         $all_genomes{"83334.1"} = 1;  
         $all_genomes{"316407.3"} = 1;  
   
         foreach my $comp_genome (keys %all_genomes){  
             my $return = $fig->bbh_list($comp_genome,[$fid]);  
             my $feature_list = $return->{$fid};  
             foreach my $peg1 (@$feature_list){  
                 my $location = $fig->feature_location($peg1);  
                 my ($pair_contig,$pair_beg,$pair_end,$pair_region_start,$pair_region_stop,$pair_genome);  
                 $pair_genome = $fig->genome_of($peg1);  
   
                 if($location =~/(.*)_(\d+)_(\d+)$/){  
                     $pair_contig = $1;  
                     $pair_beg = $2;  
                     $pair_end = $3;  
                     if ($pair_beg < $pair_end)  
                     {  
                         $pair_region_start = $pair_beg - 4000;  
                         $pair_region_stop = $pair_end + 4000;  
                         $offset = ($2+(($3-$2)/2))-($gd_window_size/2);  
                     }  
                     else  
                     {  
                         $pair_region_start = $pair_end-4000;  
                         $pair_region_stop = $pair_beg+4000;  
                         $offset = ($3+(($2-$3)/2))-($gd_window_size/2);  
                         $reverse_flag{$pair_genome} = $peg1;  
2317                      }                      }
2318    
2319                      push (@start_array_region, $offset);      # cluster the genes
2320                      $all_genomes{$pair_genome} = 1;      my @all_pegs = keys %all_genes;
2321                      my ($pair_features) = $fig->genes_in_region($pair_genome, $pair_contig, $pair_region_start, $pair_region_stop);      my $color_sets = &cluster_genes($fig,\@all_pegs,$fid);
                     push(@$all_regions,$pair_features);  
                     foreach my $pair_feature (@$pair_features){ $all_genes{$pair_feature} = $peg1;}  
                 }  
             }  
         }  
     }  
   
     # get the PCH to each of the genes  
     my $pch_sets = [];  
     my %pch_already;  
     foreach my $gene_peg (keys %all_genes)  
     {  
         if ($pch_already{$gene_peg}){(next);};  
         my $gene_set = [$gene_peg];  
         foreach my $pch_peg ($fig->in_pch_pin_with($gene_peg)) {  
             $pch_peg =~ s/,.*$//;  
             my $pch_genome = $fig->genome_of($pch_peg);  
             if ( ($gene_peg ne $pch_peg) && ($all_genomes{$pch_genome})) {  
                 push(@$gene_set,$pch_peg);  
                 $pch_already{$pch_peg}=1;  
             }  
             $pch_already{$gene_peg}=1;  
         }  
         push(@$pch_sets,$gene_set);  
     }  
   
     #create a rank of the pch's  
     my %pch_set_rank;  
     my $order = 0;  
     foreach my $set (@$pch_sets){  
         my $count = scalar(@$set);  
         $pch_set_rank{$order} = $count;  
         $order++;  
     }  
   
     my %peg_rank;  
     my $counter =  1;  
     foreach my $pch_order (sort {$pch_set_rank{$b} <=> $pch_set_rank{$a}} keys %pch_set_rank){  
         my $good_set = @$pch_sets[$pch_order];  
         my $flag_set = 0;  
         if (scalar (@$good_set) > 1)  
         {  
             foreach my $peg (@$good_set){  
                 if ((!$peg_rank{$peg})){  
                     $peg_rank{$peg} = $counter;  
                     $flag_set = 1;  
                 }  
             }  
             $counter++ if ($flag_set == 1);  
         }  
         else  
         {  
             foreach my $peg (@$good_set){  
                 $peg_rank{$peg} = "20";  
             }  
         }  
     }  
   
   
 #    my $bbh_sets = [];  
 #    my %already;  
 #    foreach my $gene_key (keys(%all_genes)){  
 #       if($already{$gene_key}){(next);}  
 #       my $gene_set = [$gene_key];  
 #  
 #       my $gene_key_genome = $fig->genome_of($gene_key);  
 #  
 #       foreach my $genome_key (keys(%all_genomes)){  
 #           #(next) if ($gene_key_genome eq $genome_key);  
 #           my $return = $fig->bbh_list($genome_key,[$gene_key]);  
 #  
 #           my $feature_list = $return->{$gene_key};  
 #           foreach my $fl (@$feature_list){  
 #               push(@$gene_set,$fl);  
 #           }  
 #       }  
 #       $already{$gene_key} = 1;  
 #       push(@$bbh_sets,$gene_set);  
 #    }  
 #  
 #    my %bbh_set_rank;  
 #    my $order = 0;  
 #    foreach my $set (@$bbh_sets){  
 #       my $count = scalar(@$set);  
 #       $bbh_set_rank{$order} = $count;  
 #       $order++;  
 #    }  
 #  
 #    my %peg_rank;  
 #    my $counter =  1;  
 #    foreach my $bbh_order (sort {$bbh_set_rank{$b} <=> $bbh_set_rank{$a}} keys %bbh_set_rank){  
 #       my $good_set = @$bbh_sets[$bbh_order];  
 #       my $flag_set = 0;  
 #       if (scalar (@$good_set) > 1)  
 #       {  
 #           foreach my $peg (@$good_set){  
 #               if ((!$peg_rank{$peg})){  
 #                   $peg_rank{$peg} = $counter;  
 #                   $flag_set = 1;  
 #               }  
 #           }  
 #           $counter++ if ($flag_set == 1);  
 #       }  
 #       else  
 #       {  
 #           foreach my $peg (@$good_set){  
 #               $peg_rank{$peg} = "20";  
 #           }  
 #       }  
 #    }  
2322    
2323      foreach my $region (@$all_regions){      foreach my $region (@$all_regions){
2324          my $sample_peg = @$region[0];          my $sample_peg = @$region[0];
# Line 2481  Line 2334 
2334    
2335          my $second_line_config = { 'title' => "$region_gs",          my $second_line_config = { 'title' => "$region_gs",
2336                                     'short_title' => "",                                     'short_title' => "",
2337                                     'basepair_offset' => '0'                                     'basepair_offset' => '0',
2338                                       'no_middle_line' => '1'
2339                                     };                                     };
2340    
2341          my $line_data = [];          my $line_data = [];
# Line 2498  Line 2352 
2352              my $links_list = [];              my $links_list = [];
2353              my $descriptions = [];              my $descriptions = [];
2354    
2355              my $color = $peg_rank{$fid1};              my $color = $color_sets->{$fid1};
2356    
2357              # get subsystem information              # get subsystem information
2358              my $function = $fig->function_of($fid1);              my $function = $fig->function_of($fid1);
# Line 2566  Line 2420 
2420                  # if there is an overlap, put into second line                  # if there is an overlap, put into second line
2421                  if ($second_line_flag == 1){ push(@$second_line_data,$element_hash); $prev_second_flag = 1;}                  if ($second_line_flag == 1){ push(@$second_line_data,$element_hash); $prev_second_flag = 1;}
2422                  else{ push(@$line_data,$element_hash); $prev_second_flag = 0;}                  else{ push(@$line_data,$element_hash); $prev_second_flag = 0;}
   
2423              }              }
2424          }          }
2425          $gd->add_line($line_data, $line_config);          $gd->add_line($line_data, $line_config);
# Line 2575  Line 2428 
2428      return $gd;      return $gd;
2429  }  }
2430    
2431    sub cluster_genes {
2432        my($fig,$all_pegs,$peg) = @_;
2433        my(%seen,$i,$j,$k,$x,$cluster,$conn,$pegI,$red_set);
2434    
2435        my @color_sets = ();
2436    
2437        $conn = &get_connections_by_similarity($fig,$all_pegs);
2438    
2439        for ($i=0; ($i < @$all_pegs); $i++) {
2440            if ($all_pegs->[$i] eq $peg) { $pegI = $i }
2441            if (! $seen{$i}) {
2442                $cluster = [$i];
2443                $seen{$i} = 1;
2444                for ($j=0; ($j < @$cluster); $j++) {
2445                    $x = $conn->{$cluster->[$j]};
2446                    foreach $k (@$x) {
2447                        if (! $seen{$k}) {
2448                            push(@$cluster,$k);
2449                            $seen{$k} = 1;
2450                        }
2451                    }
2452                }
2453    
2454                if ((@$cluster > 1) || ($cluster->[0] eq $pegI)) {
2455                    push(@color_sets,$cluster);
2456                }
2457            }
2458        }
2459        for ($i=0; ($i < @color_sets) && (! &in($pegI,$color_sets[$i])); $i++) {}
2460        $red_set = $color_sets[$i];
2461        splice(@color_sets,$i,1);
2462        @color_sets = sort { @$b <=> @$a } @color_sets;
2463        unshift(@color_sets,$red_set);
2464    
2465        my $color_sets = {};
2466        for ($i=0; ($i < @color_sets); $i++) {
2467            foreach $x (@{$color_sets[$i]}) {
2468                $color_sets->{$all_pegs->[$x]} = $i;
2469            }
2470        }
2471        return $color_sets;
2472    }
2473    
2474    sub get_connections_by_similarity {
2475        my($fig,$all_pegs) = @_;
2476        my($i,$j,$tmp,$peg,%pos_of);
2477        my($sim,%conn,$x,$y);
2478    
2479        for ($i=0; ($i < @$all_pegs); $i++) {
2480            $tmp = $fig->maps_to_id($all_pegs->[$i]);
2481            push(@{$pos_of{$tmp}},$i);
2482            if ($tmp ne $all_pegs->[$i]) {
2483                push(@{$pos_of{$all_pegs->[$i]}},$i);
2484            }
2485        }
2486    
2487        foreach $y (keys(%pos_of)) {
2488            $x = $pos_of{$y};
2489            for ($i=0; ($i < @$x); $i++) {
2490                for ($j=$i+1; ($j < @$x); $j++) {
2491                    push(@{$conn{$x->[$i]}},$x->[$j]);
2492                    push(@{$conn{$x->[$j]}},$x->[$i]);
2493                }
2494            }
2495        }
2496    
2497        for ($i=0; ($i < @$all_pegs); $i++) {
2498            foreach $sim ($fig->nsims($all_pegs->[$i],500,10,"raw")) {
2499                if (defined($x = $pos_of{$sim->id2})) {
2500                    foreach $y (@$x) {
2501                        push(@{$conn{$i}},$y);
2502                    }
2503                }
2504            }
2505        }
2506        return \%conn;
2507    }
2508    
2509    sub in {
2510        my($x,$xL) = @_;
2511        my($i);
2512    
2513        for ($i=0; ($i < @$xL) && ($x != $xL->[$i]); $i++) {}
2514        return ($i < @$xL);
2515    }

Legend:
Removed from v.1.37  
changed lines
  Added in v.1.38

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3