[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.20, Sun Mar 25 14:21:41 2007 UTC revision 1.22, Thu Apr 12 18:35:28 2007 UTC
# Line 117  Line 117 
117      return bless $self, $class;      return bless $self, $class;
118  }  }
119    
120    sub function_of {
121        my($self,$id) = @_;
122    
123        my $fig  = $self->{_fig};
124        my $func = $fig->function_of($id);
125    
126        return ($func ? $func : "");
127    }
128    
129  =head3 genomes  =head3 genomes
130    
# Line 777  Line 784 
784    
785  =item RETURNS:  =item RETURNS:
786    
787  The TAxon-ID for the "GenomeO" object containg the feature.  The Taxon-ID for the "GenomeO" object containing the feature.
788    
789  =back  =back
790    
# Line 1144  Line 1151 
1151      my $figO = $self->{_figO};      my $figO = $self->{_figO};
1152      my($tmp_dir) = $figO->{_tmp_dir};      my($tmp_dir) = $figO->{_tmp_dir};
1153    
1154        #...Skip tests and return '0' if truncated...
1155      if (! $self->possibly_truncated)      if (! $self->possibly_truncated)
1156      {      {
1157            #...Get best precomputed BLAST hit if E-value < 1.0e-50:
1158          my @sims = $self->sims( -max => 1, -cutoff => 1.0e-50);          my @sims = $self->sims( -max => 1, -cutoff => 1.0e-50);
1159    
1160            #...If a sim was returned:
1161          if (my $sim = shift @sims)          if (my $sim = shift @sims)
1162          {          {
1163                #...Get best hit FID and boundaries:
1164              my $peg2 = $sim->id2;              my $peg2 = $sim->id2;
1165              my $ln1  = $sim->ln1;              my $ln1  = $sim->ln1;
1166              my $ln2  = $sim->ln2;              my $ln2  = $sim->ln2;
1167              my $b2   = $sim->b2;              my $b2   = $sim->b2;
1168              my $e2   = $sim->e2;              my $e2   = $sim->e2;
1169    
1170                #...Convert from AA to BP, and pad out w/ 100 bp guard region:
1171              my $adjL = 100 + (($b2-1) * 3);              my $adjL = 100 + (($b2-1) * 3);
1172              my $adjR = 100 + (($ln2 - $e2) * 3);              my $adjR = 100 + (($ln2 - $e2) * 3);
1173    
1174                #...If hit is more than 20% longer than query:
1175              if ($ln2 > (1.2 * $ln1))              if ($ln2 > (1.2 * $ln1))
1176              {              {
1177                    #...Get and parse query location:
1178                  my $loc = $self->location;                  my $loc = $self->location;
1179                  if ($loc =~ /^(\S+)_(\d+)_(\d+)/)                  if ($loc =~ /^(\S+)_(\d+)_(\d+)/)
1180                  {                  {
1181                      my $contig = $1;                      my $contig = $1;
1182                      my $beg    = $2;                      my $beg    = $2;
1183                      my $end = $3;                      my $end = $3;
1184    
1185                        #...Create new ContigO object:
1186                      my $contigO = new ContigO($figO,$self->genome->id,$contig);                      my $contigO = new ContigO($figO,$self->genome->id,$contig);
1187    
1188                        #...Extract DNA subsequence, including guard regions:
1189                      my $begA = &max(1,$beg - $adjL);                      my $begA = &max(1,$beg - $adjL);
1190                      my $endA = &min($end+$adjR,$contigO->contig_length);                      my $endA = &min($end+$adjR,$contigO->contig_length);
1191                      my $dna  = $contigO->dna_seq($begA,$endA);                      my $dna  = $contigO->dna_seq($begA,$endA);
1192                      open(TMP,">$tmp_dir/tmp_dna") || die "couild not open tmp_dna";  
1193                        #...Open tmp-file and write FASTA containing DNA subregion to be BLASTed:
1194                        open( TMP, ">$tmp_dir/tmp_dna") || die "could not open tmp_dna";
1195                      print TMP ">dna\n$dna\n";                      print TMP ">dna\n$dna\n";
1196                      close(TMP);                      close(TMP);
1197    
1198                      my $peg2O = new FeatureO($figO,$peg2);                      #...Create new FeatureO object corresponding tp $peg2:
1199                      my $prot  = $peg2O->prot_seq;                      my $pegO2 = new FeatureO($figO,$peg2);
1200    
1201                        #...Fetch its translation, and print to tmp FASTA file for BLASTing:
1202                        my $prot  = $pegO2->prot_seq;
1203                      open(TMP,">$tmp_dir/tmp_prot") || die "could not open tmp_prot";                      open(TMP,">$tmp_dir/tmp_prot") || die "could not open tmp_prot";
1204                      print TMP ">tmp_prot\n$prot\n";                      print TMP ">tmp_prot\n$prot\n";
1205                      close(TMP);                      close(TMP);
1206    
1207                        #...Build BLAST nucleotide database for extracted DNA region,
1208                        #   and TBLASTN $peg2 against the DNA:
1209                      &run("formatdb -i $tmp_dir/tmp_dna -pF");                      &run("formatdb -i $tmp_dir/tmp_dna -pF");
1210                      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-50 |")
1211                          || die "could not blast";                          || die "could not blast";
1212    
1213                        #...Parse the TBLASTN output; find and sort HSPs by left boundary:
1214                      my $db_seq_out = &gjoparseblast::next_blast_subject(\*BLAST,1);                      my $db_seq_out = &gjoparseblast::next_blast_subject(\*BLAST,1);
1215                      my @hsps       = sort { $a->[0] <=> $b->[0] }                      my @hsps       = sort { $a->[0] <=> $b->[0] }
1216                                       map { [$_->[9],$_->[10],$_->[12],$_->[13]] }                                       map { [$_->[9],$_->[10],$_->[12],$_->[13]] }
1217                                       grep { $_->[1] < 1.0e-50 }                                       grep { $_->[1] < 1.0e-50 }
1218                                       @{$db_seq_out->[6]};                                       @{$db_seq_out->[6]};
1219    
1220                        #...Extract HSP boundary pairs:
1221                      my @prot = map { [$_->[0],$_->[1]] } @hsps;                      my @prot = map { [$_->[0],$_->[1]] } @hsps;
1222                      my @dna  = map { [$_->[2],$_->[3]] } @hsps;                      my @dna  = map { [$_->[2],$_->[3]] } @hsps;
1223    
1224                        #...If the "cover" of the HSPs covers more than 90% of $peg2 w gaps < 3 AA,
1225                        #   and the "cover" of the HPSs cover more than 90% of the extracted DNA
1226                        #   w/ gaps < 9 bp (but not a multiple of 3), suspect a possible frameshift:
1227                      if (&covers(\@prot,length($prot),3,0) && &covers(\@dna,3*length($prot),9,1))                      if (&covers(\@prot,length($prot),3,0) && &covers(\@dna,3*length($prot),9,1))
1228                      {                      {
1229                          return 1;                          return 1;
# Line 1374  Line 1410 
1410      my $fig  = $figO->{_fig};      my $fig  = $figO->{_fig};
1411    
1412      my $cutoff = $args{-cutoff} ? $args{-cutoff} : 1.0e-5;      my $cutoff = $args{-cutoff} ? $args{-cutoff} : 1.0e-5;
1413      my $all    = $args{-all}    ? $args{-all}    : "fig";      my $all    = $args{-all}    ? 'all'          : "fig";
1414      my $max    = $args{-max}    ? $args{-max}    : 10000;      my $max    = $args{-max}    ? $args{-max}    : 10000;
1415    
1416      my @sims = $fig->sims($self->id,$max,$cutoff,$all);      my @sims = $fig->sims($self->id,$max,$cutoff,$all);
# Line 1383  Line 1419 
1419          my $peg1 = FeatureO->new($figO, $sims[0]->[0]);          my $peg1 = FeatureO->new($figO, $sims[0]->[0]);
1420    
1421          foreach my $sim (@sims) {          foreach my $sim (@sims) {
1422              $sim->[0] = $peg1;  #           $sim->[0] = $peg1;
1423              $sim->[1] = FeatureO->new($figO, $sim->[1]);  #           $sim->[1] = FeatureO->new($figO, $sim->[1]);
1424          }          }
1425      }      }
1426    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3