[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.19, Mon Mar 19 19:49:25 2007 UTC revision 1.23, Mon May 7 00:43:26 2007 UTC
# Line 72  Line 72 
72  use SproutFIG;  use SproutFIG;
73  use Tracer;  use Tracer;
74  use Data::Dumper;  use Data::Dumper;
75    use Carp;
76  use FigFams;  use FigFams;
77  use gjoparseblast;  use gjoparseblast;
78    
# Line 117  Line 118 
118      return bless $self, $class;      return bless $self, $class;
119  }  }
120    
121    sub function_of {
122        my($self,$id) = @_;
123    
124        my $fig  = $self->{_fig};
125        my $func = $fig->function_of($id);
126    
127        return ($func ? $func : "");
128    }
129    
130  =head3 genomes  =head3 genomes
131    
# Line 697  Line 705 
705  package FeatureO;  package FeatureO;
706  ########################################################################  ########################################################################
707  use Data::Dumper;  use Data::Dumper;
708    use Carp;
709    
710  =head1 FeatureO  =head1 FeatureO
711    
# Line 777  Line 786 
786    
787  =item RETURNS:  =item RETURNS:
788    
789  The TAxon-ID for the "GenomeO" object containg the feature.  The Taxon-ID for the "GenomeO" object containing the feature.
790    
791  =back  =back
792    
# Line 1144  Line 1153 
1153      my $figO = $self->{_figO};      my $figO = $self->{_figO};
1154      my($tmp_dir) = $figO->{_tmp_dir};      my($tmp_dir) = $figO->{_tmp_dir};
1155    
1156        #...Skip tests and return '0' if truncated...
1157      if (! $self->possibly_truncated)      if (! $self->possibly_truncated)
1158      {      {
1159            #...Get best precomputed BLAST hit if E-value < 1.0e-50:
1160          my @sims = $self->sims( -max => 1, -cutoff => 1.0e-50);          my @sims = $self->sims( -max => 1, -cutoff => 1.0e-50);
1161    
1162            #...If a sim was returned:
1163          if (my $sim = shift @sims)          if (my $sim = shift @sims)
1164          {          {
1165                #...Get best hit FID and boundaries:
1166              my $peg2 = $sim->id2;              my $peg2 = $sim->id2;
1167              my $ln1  = $sim->ln1;              my $ln1  = $sim->ln1;
1168              my $ln2  = $sim->ln2;              my $ln2  = $sim->ln2;
1169              my $b2   = $sim->b2;              my $b2   = $sim->b2;
1170              my $e2   = $sim->e2;              my $e2   = $sim->e2;
1171    
1172                #...Convert from AA to BP, and pad out w/ 100 bp guard region:
1173              my $adjL = 100 + (($b2-1) * 3);              my $adjL = 100 + (($b2-1) * 3);
1174              my $adjR = 100 + (($ln2 - $e2) * 3);              my $adjR = 100 + (($ln2 - $e2) * 3);
1175    
1176                #...If hit is more than 20% longer than query:
1177              if ($ln2 > (1.2 * $ln1))              if ($ln2 > (1.2 * $ln1))
1178              {              {
1179                    #...Get and parse query location:
1180                  my $loc = $self->location;                  my $loc = $self->location;
1181                  if ($loc =~ /^(\S+)_(\d+)_(\d+)/)                  if ($loc =~ /^(\S+)_(\d+)_(\d+)/)
1182                  {                  {
1183                      my $contig = $1;                      my $contig = $1;
1184                      my $beg    = $2;                      my $beg    = $2;
1185                      my $end = $3;                      my $end = $3;
1186    
1187                        #...Create new ContigO object:
1188                      my $contigO = new ContigO($figO,$self->genome->id,$contig);                      my $contigO = new ContigO($figO,$self->genome->id,$contig);
1189    
1190                        #...Extract DNA subsequence, including guard regions:
1191                      my $begA = &max(1,$beg - $adjL);                      my $begA = &max(1,$beg - $adjL);
1192                      my $endA = &min($end+$adjR,$contigO->contig_length);                      my $endA = &min($end+$adjR,$contigO->contig_length);
1193                      my $dna  = $contigO->dna_seq($begA,$endA);                      my $dna  = $contigO->dna_seq($begA,$endA);
1194                      open(TMP,">$tmp_dir/tmp_dna") || die "couild not open tmp_dna";                      if (defined($dna) && (length($dna) > 90))
1195                        {
1196                            #...Open tmp-file and write FASTA containing DNA subregion to be BLASTed:
1197                            open( TMP, ">$tmp_dir/tmp_dna") || die "could not open tmp_dna";
1198                      print TMP ">dna\n$dna\n";                      print TMP ">dna\n$dna\n";
1199                      close(TMP);                      close(TMP);
1200    
1201                      my $peg2O = new FeatureO($figO,$peg2);                          #...Create new FeatureO object corresponding tp $peg2:
1202                      my $prot  = $peg2O->prot_seq;                          my $pegO2 = new FeatureO($figO,$peg2);
1203    
1204                            #...Fetch its translation, and print to tmp FASTA file for BLASTing:
1205                            my $prot  = $pegO2->prot_seq;
1206                            if (defined($prot) && (length($prot) > 30))
1207                            {
1208                      open(TMP,">$tmp_dir/tmp_prot") || die "could not open tmp_prot";                      open(TMP,">$tmp_dir/tmp_prot") || die "could not open tmp_prot";
1209                      print TMP ">tmp_prot\n$prot\n";                      print TMP ">tmp_prot\n$prot\n";
1210                      close(TMP);                      close(TMP);
1211    
1212                                #...Build BLAST nucleotide database for extracted DNA region,
1213                                #   and TBLASTN $peg2 against the DNA:
1214                      &run("formatdb -i $tmp_dir/tmp_dna -pF");                      &run("formatdb -i $tmp_dir/tmp_dna -pF");
1215                      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 |")
1216                          || die "could not blast";                          || die "could not blast";
1217    
1218                                #...Parse the TBLASTN output; find and sort HSPs by left boundary:
1219                      my $db_seq_out = &gjoparseblast::next_blast_subject(\*BLAST,1);                      my $db_seq_out = &gjoparseblast::next_blast_subject(\*BLAST,1);
1220                      my @hsps       = sort { $a->[0] <=> $b->[0] }                      my @hsps       = sort { $a->[0] <=> $b->[0] }
1221                                       map { [$_->[9],$_->[10],$_->[12],$_->[13]] }                                       map { [$_->[9],$_->[10],$_->[12],$_->[13]] }
1222                                       grep { $_->[1] < 1.0e-50 }                                       grep { $_->[1] < 1.0e-50 }
1223                                       @{$db_seq_out->[6]};                                       @{$db_seq_out->[6]};
1224    
1225                                #...Extract HSP boundary pairs:
1226                      my @prot = map { [$_->[0],$_->[1]] } @hsps;                      my @prot = map { [$_->[0],$_->[1]] } @hsps;
1227                      my @dna  = map { [$_->[2],$_->[3]] } @hsps;                      my @dna  = map { [$_->[2],$_->[3]] } @hsps;
1228    
1229                                #...If the "cover" of the HSPs covers more than 90% of $peg2 w gaps < 3 AA,
1230                                #   and the "cover" of the HPSs cover more than 90% of the extracted DNA
1231                                #   w/ gaps < 9 bp (but not a multiple of 3), suspect a possible frameshift:
1232                      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))
1233                      {                      {
1234                          return 1;                          return 1;
# Line 1196  Line 1237 
1237              }              }
1238          }          }
1239      }      }
1240            }
1241        }
1242      return 0;      return 0;
1243  }  }
1244    
# Line 1226  Line 1269 
1269    
1270  sub run {  sub run {
1271      my($cmd) = @_;      my($cmd) = @_;
1272      (system($cmd) == 0) || Confess("FAILED: $cmd");      (system($cmd) == 0) || confess("FAILED: $cmd");
1273  }  }
1274    
1275    
# Line 1374  Line 1417 
1417      my $fig  = $figO->{_fig};      my $fig  = $figO->{_fig};
1418    
1419      my $cutoff = $args{-cutoff} ? $args{-cutoff} : 1.0e-5;      my $cutoff = $args{-cutoff} ? $args{-cutoff} : 1.0e-5;
1420      my $all    = $args{-all}    ? $args{-all}    : "fig";      my $all    = $args{-all}    ? 'all'          : "fig";
1421      my $max    = $args{-max}    ? $args{-max}    : 10000;      my $max    = $args{-max}    ? $args{-max}    : 10000;
1422    
1423      return $fig->sims($self->id,$max,$cutoff,$all);      my @sims = $fig->sims($self->id,$max,$cutoff,$all);
1424    
1425        if (@sims) {
1426            my $peg1 = FeatureO->new($figO, $sims[0]->[0]);
1427    
1428            foreach my $sim (@sims) {
1429    #           $sim->[0] = $peg1;
1430    #           $sim->[1] = FeatureO->new($figO, $sim->[1]);
1431            }
1432        }
1433    
1434        return @sims;
1435  }  }
1436    
1437    

Legend:
Removed from v.1.19  
changed lines
  Added in v.1.23

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3