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

View of /FigKernelPackages/ComparedRegions.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (download) (as text) (annotate)
Fri Jun 15 21:24:49 2007 UTC (12 years, 9 months ago) by paczian
Branch: MAIN
initial version

package ComparedRegions;

1;

use strict;
use warnings;

use FIG;
use FIGV;

sub get_compared_regions {
  my ($params) = @_;

  # check parameters
  my $fig = $params->{fig} || return undef;
  my $peg = $params->{id} || return undef;
  my $is_sprout = $params->{is_sprout} || 0;
  my $region_size = $params->{region_size} || 16000;
  my $number_of_genomes = $params->{number_of_genomes} || 5;

  # initialize data variable
  my $data = [];

  # get the n pegs closest to the one passed
  my @closest_pegs = &get_closest_pegs($fig, $peg, $is_sprout, $number_of_genomes);
  unshift(@closest_pegs, $peg);
  
  # iterate over the returned pegs
  foreach my $peg (@closest_pegs) { 
    my $loc = $fig->feature_location($peg);
    my ($contig,$beg,$end) = $fig->boundaries_of($loc);
    if ($contig && $beg && $end) {
      my $mid = int(($beg + $end) / 2);
      my $min = int($mid - ($region_size / 2));
      my $max = int($mid + ($region_size / 2));
      my $features = [];
      my $feature_ids;
      ($feature_ids, undef, undef) = $fig->genes_in_region($fig->genome_of($peg),$contig,$min,$max);
      foreach my $fid (@$feature_ids) {
	my $floc = $fig->feature_location($fid);
	my ($contig1,$beg1,$end1) = $fig->boundaries_of($fig->feature_location($fid));
	$beg1 = &in_bounds($min,$max,$beg1);
	$end1 = &in_bounds($min,$max,$end1);
	push(@$features, { 'start' => $beg1,
			   'stop' => $end1,
			   'id' => $fid });
      }
      push(@$data, { 'features' => $features,
		     'contig' => $contig,
		     'offset' => $mid - 8000 });
    }
  }

  # return the data
  return $data;
}

# returns the n closest pegs, sorted by taxonomy
sub get_closest_pegs {
  my ($fig, $peg, $is_sprout, $n) = @_;

  my($id2,$d,$peg2,$i);

  my @closest;
  if ($is_sprout) {
    @closest = map { $_->[0] } sort { $a->[1] <=> $b->[1] } $fig->bbhs($peg, 1.0e-10);
  } else {
    @closest = map { $id2 = $_->id2; ($id2 =~ /^fig\|/) ? $id2 : () } $fig->sims($peg,&FIG::max(20,$n*4),1.0e-20,"fig",&FIG::max(20,$n*4));
  }

  if (@closest >= ($n-1)) { 
    $#closest = $n-2 ;
  }
  my %closest = map { $_ => 1 } @closest;
  
  # there are dragons flying around...
  my @pinned_to = grep { ($_ ne $peg) && (! $closest{$_}) } $fig->in_pch_pin_with($peg);
  my $g1 = $fig->genome_of($peg);
  @pinned_to = map {$_->[1] } sort { $a->[0] <=> $b->[0] } map { $peg2 = $_; $d = $fig->crude_estimate_of_distance($g1,$fig->genome_of($peg2)); [$d,$peg2] } @pinned_to;
  
  if (@closest == ($n-1)) {
    $#closest = ($n - 2) - &FIG::min(scalar @pinned_to,int($n/2));
    for ($i=0; ($i < @pinned_to) && (@closest < ($n-1)); $i++) {
      if (! $closest{$pinned_to[$i]}) {
	$closest{$pinned_to[$i]} = 1;
	push(@closest,$pinned_to[$i]);
      }
    }
  }

  if ($fig->possibly_truncated($peg)) {
    push(@closest, &possible_extensions($fig, $peg, \@closest));
  }
  @closest = $fig->sort_fids_by_taxonomy(@closest);

  return @closest;

}

sub in_bounds {
    my($min,$max,$x) = @_;

    if     ($x < $min)     { return $min }
    elsif  ($x > $max)     { return $max }
    else                   { return $x   }
}

sub possible_extensions {
  my($fig, $peg,$closest_pegs) = @_;
  my($g,$sim,$id2,$peg1,%poss);
  
  $g = &FIG::genome_of($peg);
  
  foreach $peg1 (@$closest_pegs) {
    if ($g ne &genome_of($peg1)) {
      foreach $sim ($fig->sims($peg1,500,1.0e-5,"all")) {
	$id2 = $sim->id2;
	if (($id2 ne $peg) && ($id2 =~ /^fig\|$g\./) && $fig->possibly_truncated($id2)) {
	  $poss{$id2} = 1;
	}
      }
    }
  }
  return keys(%poss);
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3