[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.21, Tue Mar 27 17:57:49 2007 UTC revision 1.22, Thu Apr 12 18:35:28 2007 UTC
# Line 784  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 1151  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;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3