[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.26, Wed Sep 26 04:03:48 2007 UTC revision 1.27, Thu Oct 4 00:58:56 2007 UTC
# Line 1204  Line 1204 
1204  sub possible_frameshift {  sub possible_frameshift {
1205      my($self) = @_;      my($self) = @_;
1206      my $figO = $self->{_figO};      my $figO = $self->{_figO};
1207        my $fig = $figO->{_fig};
1208      my($tmp_dir) = $figO->{_tmp_dir};      my($tmp_dir) = $figO->{_tmp_dir};
1209    
1210      #...Skip tests and return '0' if truncated...      #...Skip tests and return '0' if truncated...
1211      if (! $self->possibly_truncated)      if (! $self->possibly_truncated)
1212      {      {
1213          #...Get best precomputed BLAST hit if E-value < 1.0e-50:          #...Get best precomputed BLAST hit if E-value < 1.0e-20:
1214          my @sims = $self->sims( -max => 1, -cutoff => 1.0e-50);          my @sims = $self->sims( -max => 5, -cutoff => 1.0e-20);
1215            while ((@sims > 0) && $fig->possibly_truncated($sims[0]->id2)) { shift @sims }
1216    
1217          #...If a sim was returned:          #...If a sim was returned:
1218          if (my $sim = shift @sims)          if (my $sim = shift @sims)
# Line 1226  Line 1228 
1228              my $adjL = 100 + (($b2-1) * 3);              my $adjL = 100 + (($b2-1) * 3);
1229              my $adjR = 100 + (($ln2 - $e2) * 3);              my $adjR = 100 + (($ln2 - $e2) * 3);
1230    
1231                if ($ENV{DEBUG}) { print STDERR "adjL = $adjL adjR = $adjR ln1 = $ln1 peg2 = $peg2 ln2 = $ln2\n" }
1232              #...If hit is more than 20% longer than query:              #...If hit is more than 20% longer than query:
1233              if ($ln2 > (1.2 * $ln1))              if ($ln2 > (1.2 * $ln1))
1234              {              {
# Line 1241  Line 1244 
1244                      my $contigO = new ContigO($figO, $self->genome->id, $contig);                      my $contigO = new ContigO($figO, $self->genome->id, $contig);
1245    
1246                      #...Extract DNA subsequence, including guard regions:                      #...Extract DNA subsequence, including guard regions:
1247                      my $begA = &max(1, $beg - $adjL);                      my($begA,$endA,$dna);
1248                      my $endA = &min($end+$adjR, $contigO->contig_length);                      if ($beg < $end)
1249                      my $dna  = $contigO->dna_seq($begA,$endA);                      {
1250                            $begA = &max(1, $beg - $adjL);
1251                            $endA = &min($end+$adjR, $contigO->contig_length);
1252                            $dna  = $contigO->dna_seq($begA,$endA);
1253                        }
1254                        else
1255                        {
1256                            $endA = &max(1, $beg - $adjL);
1257                            $begA = &min($end+$adjR, $contigO->contig_length);
1258                            $dna  = $contigO->dna_seq($begA,$endA);
1259                        }
1260    
1261                      if (defined($dna) && (length($dna) > 90))                      if (defined($dna) && (length($dna) > 90))
1262                      {                      {
1263                          #...Open tmp-file and write FASTA containing DNA subregion to be BLASTed:                          #...Open tmp-file and write FASTA containing DNA subregion to be BLASTed:
# Line 1265  Line 1279 
1279                              #...Build BLAST nucleotide database for extracted DNA region,                              #...Build BLAST nucleotide database for extracted DNA region,
1280                              #   and TBLASTN $peg2 against the DNA:                              #   and TBLASTN $peg2 against the DNA:
1281                              &run("formatdb -i $tmp_dir/tmp_dna -pF");                              &run("formatdb -i $tmp_dir/tmp_dna -pF");
1282                              open(BLAST,"blastall -i $tmp_dir/tmp_prot -d $tmp_dir/tmp_dna -p tblastn -FF -e 1.0e-50 |")                              open(BLAST,"blastall -i $tmp_dir/tmp_prot -d $tmp_dir/tmp_dna -p tblastn -FF -e 1.0e-20 |")
1283                                  || die "could not blast";                                  || die "could not blast";
1284    
1285                              #...Parse the TBLASTN output; find and sort HSPs by left boundary:                              #...Parse the TBLASTN output; find and sort HSPs by left boundary:
1286                              my $db_seq_out = &gjoparseblast::next_blast_subject(\*BLAST,1);                              my $db_seq_out = &gjoparseblast::next_blast_subject(\*BLAST,1);
1287                                if ($ENV{DEBUG}) { print STDERR &Dumper(['blast output',$db_seq_out]) }
1288                              my @hsps       = sort { $a->[0] <=> $b->[0] }                              my @hsps       = sort { $a->[0] <=> $b->[0] }
1289                                                map { [$_->[9], $_->[10], $_->[12], $_->[13]] }                                                map { [$_->[9], $_->[10], $_->[12], $_->[13]] }
1290                                                grep { $_->[1] < 1.0e-50 }                                                grep { $_->[1] < 1.0e-20 }
1291                                                @ { $db_seq_out->[6] };                                                @ { $db_seq_out->[6] };
1292    
1293                              #...Extract HSP boundary pairs:                              #...Extract HSP boundary pairs:
1294                              my @prot = map { [$_->[0], $_->[1]] } @hsps;                              my @prot = map { [$_->[0], $_->[1]] } @hsps;
1295                              my @dna  = map { [$_->[2], $_->[3]] } @hsps;                              my @dna  = map { [$_->[2], $_->[3]] } @hsps;
1296                                if ($ENV{DEBUG}) { print STDERR &Dumper(\@prot,\@dna) }
1297    
1298                              #...If the "cover" of the HSPs covers more than 90% of $peg2 w gaps < 3 AA,                              #...If the "cover" of the HSPs covers more than 90% of $peg2 w gaps < 3 AA,
1299                              #   and the "cover" of the HPSs cover more than 90% of the extracted DNA                              #   and the "cover" of the HPSs cover more than 90% of the extracted DNA
# Line 1400  Line 1416 
1416  sub covers {  sub covers {
1417      my($hsps,$ln,$diff,$must_shift) = @_;      my($hsps,$ln,$diff,$must_shift) = @_;
1418    
1419        if ($ENV{DEBUG}) { print STDERR &Dumper(['hsps',$hsps,'ln',$ln,'diff',$diff,'must_shift',$must_shift]) }
1420      my $hsp1 = shift @$hsps;      my $hsp1 = shift @$hsps;
1421      my $hsp2;      my $hsp2;
1422      my $merged = 0;      my $merged = 0;
1423      while ($hsp1 && ($hsp2 = shift @$hsps) &&      while ($hsp1 && ($hsp2 = shift @$hsps) &&
1424             ($must_shift ? &diff_frames($hsp1,$hsp2) : 1) &&             ($must_shift ? &diff_frames($hsp1,$hsp2) : 1) &&
1425             ($hsp1 = &merge($hsp1,$hsp2,$diff))) { $merged = 1 }             ($hsp1 = &merge($hsp1,$hsp2,$diff)))
1426        {
1427            $merged = 1;
1428            if ($ENV{DEBUG}) { print STDERR &Dumper(['merged',$hsp1]) }
1429        }
1430      return ($merged && $hsp1 && (($hsp1->[1] - $hsp1->[0]) > (0.9 * $ln)));      return ($merged && $hsp1 && (($hsp1->[1] - $hsp1->[0]) > (0.9 * $ln)));
1431  }  }
1432    
1433  sub diff_frames {  sub diff_frames {
1434      my($hsp1,$hsp2) = @_;      my($hsp1,$hsp2) = @_;
1435      return (($hsp1->[0] % 3) != ($hsp2->[0] % 3));      return ((($hsp1->[1]+1) % 3) != ($hsp2->[0] % 3));
1436  }  }
1437    
1438    
# Line 1429  Line 1450 
1450    
1451      my($b1,$e1) = @$hsp1;      my($b1,$e1) = @$hsp1;
1452      my($b2,$e2) = @$hsp2;      my($b2,$e2) = @$hsp2;
1453      return (($e2 > $e1) && (abs($b2-$e1) <= $diff)) ? [$b1,$e2] : undef;      return (($e2 > $e1) && (($b2-$e1) <= $diff)) ? [$b1,$e2] : undef;
1454  }  }
1455    
1456    

Legend:
Removed from v.1.26  
changed lines
  Added in v.1.27

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3