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

Diff of /FigKernelPackages/FIGO.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2, Thu Nov 23 16:15:30 2006 UTC revision 1.3, Sat Feb 17 17:14:11 2007 UTC
# Line 25  Line 25 
25  use Tracer;  use Tracer;
26  use Data::Dumper;  use Data::Dumper;
27  use FigFams;  use FigFams;
28    use gjoparseblast;
29    
30  sub new {  sub new {
31      my($class,$low_level) = @_;      my($class,$low_level) = @_;
# Line 41  Line 42 
42    
43      my $self = {};      my $self = {};
44      $self->{_fig} = $fig;      $self->{_fig} = $fig;
45        $self->{_tmp_dir} = $FIG_Config::temp;
46      return bless $self, $class;      return bless $self, $class;
47  }  }
48    
# Line 324  Line 326 
326      return map { &AnnotationO::new('AnnotationO',@$_) } $fig->feature_annotations($self->id,1);      return map { &AnnotationO::new('AnnotationO',@$_) } $fig->feature_annotations($self->id,1);
327  }  }
328    
329    sub possibly_truncated {
330        my($self) = @_;
331        my $figO = $self->{_figO};
332        my $fig  = $figO->{_fig};
333    
334        return $fig->possibly_truncated($self->id);
335    }
336    
337    sub possible_frameshift {
338        my($self) = @_;
339        my $figO = $self->{_figO};
340        my($tmp_dir) = $figO->{_tmp_dir};
341    
342        if (! $self->possibly_truncated)
343        {
344            my @sims = $self->sims( -max => 1, -cutoff => 1.0e-50);
345            if (my $sim = shift @sims)
346            {
347                my $peg2 = $sim->id2;
348                my $ln1  = $sim->ln1;
349                my $ln2  = $sim->ln2;
350                my $b2   = $sim->b2;
351                my $e2   = $sim->e2;
352                my $adjL = 100 + (($b2-1) * 3);
353                my $adjR = 100 + (($ln2 - $e2) * 3);
354                if ($ln2 > (1.2 * $ln1))
355                {
356                    my $loc = $self->location;
357                    if ($loc =~ /^(\S+)_(\d+)_(\d+)/)
358                    {
359                        my $contig = $1;
360                        my $beg    = $2;
361                        my $end = $3;
362                        my $contigO = new ContigO($figO,$self->genome,$contig);
363                        my $begA = &max(1,$beg - $adjL);
364                        my $endA = &min($end+$adjR,$contigO->contig_length);
365                        my $dna  = $contigO->dna_seq($begA,$endA);
366                        open(TMP,">$tmp_dir/tmp_dna") || die "couild not open tmp_dna";
367                        print TMP ">dna\n$dna\n";
368                        close(TMP);
369    
370                        my $peg2O = new FeatureO($figO,$peg2);
371                        my $prot  = $peg2O->prot_seq;
372                        open(TMP,">$tmp_dir/tmp_prot") || die "could not open tmp_prot";
373                        print TMP ">tmp_prot\n$prot\n";
374                        close(TMP);
375                        &run("formatdb -i $tmp_dir/tmp_dna -pF");
376                        open(BLAST,"blastall -i $tmp_dir/tmp_prot -d $tmp_dir/tmp_dna -p tblastn -FF -e 1.0e-50 |")
377                            || die "could not blast";
378    
379                        my $db_seq_out = &gjoparseblast::next_blast_subject(\*BLAST,1);
380                        my @hsps       = sort { $a->[0] <=> $b->[0] }
381                                         map { [$_->[9],$_->[10],$_->[12],$_->[13]] }
382                                         grep { $_->[1] < 1.0e-50 }
383                                         @{$db_seq_out->[6]};
384                        my @prot = map { [$_->[0],$_->[1]] } @hsps;
385                        my @dna  = map { [$_->[2],$_->[3]] } @hsps;
386                        if (&covers(\@prot,length($prot),3) && &covers(\@dna,3*length($prot),9))
387                        {
388                            return 1;
389                        }
390                    }
391                }
392            }
393        }
394        return 0;
395    }
396    
397    sub run {
398        my($cmd) = @_;
399        (system($cmd) == 0) || Confess("FAILED: $cmd");
400    }
401    
402    sub max {
403        my($x,$y) = @_;
404        return ($x < $y) ? $y : $x;
405    }
406    
407    sub min {
408        my($x,$y) = @_;
409        return ($x < $y) ? $x : $y;
410    }
411    
412    sub covers {
413        my($hsps,$ln,$diff) = @_;
414    
415        my $hsp1 = shift @$hsps;
416        my $hsp2;
417        while ($hsp1 && ($hsp2 = shift @$hsps) && ($hsp1 = &merge($hsp1,$hsp2,$diff))) {}
418        return ($hsp1 && (($hsp1->[1] - $hsp1->[0]) > (0.9 * $ln)));
419    }
420    
421    sub merge {
422        my($hsp1,$hsp2,$diff) = @_;
423    
424        my($b1,$e1) = @$hsp1;
425        my($b2,$e2) = @$hsp2;
426        return (($e2 > $e1) && (abs($b2-$e1) <= $diff)) ? [$b1,$e2] : undef;
427    }
428    
429  use Sim;  use Sim;
430  sub sims {  sub sims {
431      my($self,%args) = @_;      my($self,%args) = @_;

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3