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

Diff of /FigKernelPackages/FIGV.pm

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

revision 1.54, Thu May 29 23:00:34 2008 UTC revision 1.55, Tue Jun 24 16:25:54 2008 UTC
# Line 28  Line 28 
28  use Data::Dumper;  use Data::Dumper;
29  use vars qw($AUTOLOAD);  use vars qw($AUTOLOAD);
30  use DB_File;  use DB_File;
31    use File::Basename;
32  use FileHandle;  use FileHandle;
33    
34  sub new {  sub new {
# Line 844  Line 845 
845      }      }
846      else      else
847      {      {
848          &load_tbl($self);          $self->load_feature_indexes();
849          my $tblH = $self->{_tbl};          #
850            # Use the recno index if exists.
851            #
852    
853          my $maxV = 0;          my $maxV = 0;
854          my $minV = 1000000000;          my $minV = 1000000000;
855          my $genes = [];          my $genes = [];
856    
857            my $recnos = $self->{_feat_recno};
858    
859            if (ref($recnos))
860            {
861                while (my($ftype, $list) = each (%$recnos))
862                {
863                    #
864                    # Look up start/end of this contig in the btree index.
865                    #
866    
867                    my $inf = $self->{_feat_btree}->{$ftype}->{$contig};
868                    my($istart, $iend);
869                    if ($inf)
870                    {
871                        ($istart, $iend) = split(/$;/, $inf);
872                    }
873                    else
874                    {
875                        $istart = 0;
876                        $iend = $#$list;
877                    }
878    
879                    for (my $idx = $istart; $idx <= $iend; $idx++)
880                    {
881                        my($fid, $fcontig, $fbeg, $fend, $fstrand) = split(/$;/, $list->[$idx]);
882                        if ($contig eq $fcontig and &overlaps($beg, $end, $fbeg, $fend))
883                        {
884                            $minV = &FIG::min($minV,$fbeg,$fend);
885                            $maxV = &FIG::max($maxV,$fbeg,$fend);
886                            push(@$genes,$fid);
887                        }
888                    }
889                }
890            }
891            else
892            {
893                &load_tbl($self);
894                my $tblH = $self->{_tbl};
895          while ( my($fid,$tuple) = each %$tblH)          while ( my($fid,$tuple) = each %$tblH)
896          {          {
897              if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))              if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))
# Line 867  Line 910 
910                  }                  }
911              }              }
912          }          }
913            }
914          return ($genes,$minV,$maxV);          return ($genes,$minV,$maxV);
915      }      }
916  }  }
# Line 892  Line 936 
936      }      }
937      else      else
938      {      {
939          &load_tbl($self);          $self->load_feature_indexes();
940    
941          my %contigs;          my %contigs;
942            my $btrees = $self->{_feat_btree};
943    
944            if (ref($btrees))
945            {
946                while (my($ftype, $btree) = each (%$btrees))
947                {
948                    map { $contigs{$_} = 1 } grep { !/^fig/ } keys %$btree ;
949                }
950            }
951            else
952            {
953                &load_tbl($self);
954          my $tblH = $self->{_tbl};          my $tblH = $self->{_tbl};
955          while ( my($fid,$tuple) = each %$tblH)          while ( my($fid,$tuple) = each %$tblH)
956          {          {
# Line 902  Line 959 
959                  $contigs{$1} = 1;                  $contigs{$1} = 1;
960              }              }
961          }          }
962            }
963          return keys(%contigs);          return keys(%contigs);
964      }      }
965  }  }
# Line 919  Line 977 
977      }      }
978      else      else
979      {      {
980            $self->load_feature_indexes();
981    
982            my %contigs;
983            my $btrees = $self->{_feat_btree};
984    
985            if (ref($btrees))
986            {
987                my $btree = $btrees->{$type};
988                return sort { &FIG::by_fig_id($a, $b) }
989                    grep { /^fig/ } keys %$btree;
990            }
991            else
992            {
993          &load_tbl($self);          &load_tbl($self);
994          my $tblH = $self->{_tbl};          my $tblH = $self->{_tbl};
995          return sort { &FIG::by_fig_id($a,$b) }          return sort { &FIG::by_fig_id($a,$b) }
# Line 926  Line 997 
997                 keys(%$tblH);                 keys(%$tblH);
998      }      }
999  }  }
1000    }
1001    
1002  sub all_features_detailed_fast {  sub all_features_detailed_fast {
1003      my($self,$genome, $regmin, $regmax, $contig) = @_;      my($self,$genome, $regmin, $regmax, $contig) = @_;
# Line 940  Line 1012 
1012      }      }
1013      else      else
1014      {      {
1015            $self->load_feature_indexes();
1016            my $feat_details = [];
1017            my $recnos = $self->{_feat_recno};
1018    
1019            if (ref($recnos))
1020            {
1021                while (my($ftype, $list) = each (%$recnos))
1022                {
1023                    #
1024                    # Look up start/end of this contig in the btree index.
1025                    #
1026    
1027                    my $inf = $self->{_feat_btree}->{$ftype}->{$contig};
1028                    my($istart, $iend);
1029                    if ($inf)
1030                    {
1031                        ($istart, $iend) = split(/$;/, $inf);
1032                    }
1033                    else
1034                    {
1035                        $istart = 0;
1036                        $iend = $#$list;
1037                    }
1038    
1039                    for (my $idx = $istart; $idx <= $iend; $idx++)
1040                    {
1041                        my($fid, $fcontig, $fbeg, $fend, $fstrand) = split(/$;/, $list->[$idx]);
1042    
1043                        if (not defined($regmin) or not defined($regmax) or not defined($contig) or
1044                            (($contig eq $fcontig) and
1045                             ($fbeg < $regmin and $regmin < $fend) or ($fbeg < $regmax and $regmax < $fend) or ($fbeg > $regmin and $fend < $regmax)))
1046                        {
1047                            my($loc, $index, @aliases) = split(/$;/, $self->{_feat_btree}->{$ftype}->{$fid});
1048    
1049                            push(@$feat_details,[$fid, $loc, join(",", @aliases), $ftype, $fbeg, $fend, $self->function_of($fid),'master','']);
1050                        }
1051                    }
1052                }
1053            }
1054            else
1055            {
1056          &load_tbl($self);          &load_tbl($self);
1057          my $tblH = $self->{_tbl};          my $tblH = $self->{_tbl};
         my $feat_details = [];  
1058          while ( my($fid,$tuple) = each %$tblH)          while ( my($fid,$tuple) = each %$tblH)
1059          {          {
1060              if ($fid =~ /^fig\|\d+\.\d+\.(\S+)\.\d+/)              if ($fid =~ /^fig\|\d+\.\d+\.(\S+)\.\d+/)
# Line 973  Line 1085 
1085                  }                  }
1086              }              }
1087          }          }
1088            }
1089          return $feat_details;          return $feat_details;
1090      }      }
1091  }  }
# Line 1035  Line 1148 
1148      my $newG    = $self->{_genome};      my $newG    = $self->{_genome};
1149      my $newGdir = $self->{_orgdir};      my $newGdir = $self->{_orgdir};
1150    
1151      if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))      if (($fid =~ /^fig\|(\d+\.\d+)\.([^.]+)/) && ($1 eq $newG))
1152        {
1153            my $ftype = $2;
1154    
1155            $self->load_feature_indexes();
1156    
1157            my $btree = $self->{_feat_btree}->{$ftype};
1158            if ($btree)
1159            {
1160                my($loc, $idx, @aliases) = split(/$;/, $btree->{$fid});
1161                return wantarray ? split(/,/, $loc) : $loc;
1162            }
1163            else
1164      {      {
1165          &load_tbl($self);          &load_tbl($self);
1166          if (my $x = $self->{_tbl}->{$fid})          if (my $x = $self->{_tbl}->{$fid})
# Line 1054  Line 1179 
1179              return undef;              return undef;
1180          }          }
1181      }      }
1182        }
1183      else      else
1184      {      {
1185          return scalar $fig->feature_location($fid);          return scalar $fig->feature_location($fid);
# Line 1215  Line 1341 
1341    
1342      if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))      if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1343      {      {
1344          &load_pseq($self);          &load_feature_indexes($self);
1345          return $self->{_pseq}->{$peg};  
1346            my $out = $self->{_feat}->{peg}->{$peg};
1347            if (defined($out))
1348            {
1349                return $out;
1350            }
1351    
1352            #
1353            # If we have a blast-formatted fasta, use fastacmd to
1354            # do the lookup, and cache the output for later use.
1355            #
1356    
1357            if ($self->{_feat_fasta}->{peg})
1358            {
1359                my $id = "gnl|$peg";
1360                my $cmd = "$FIG_Config::ext_bin/fastacmd -d $self->{_feat_fasta}->{peg} -s '$id'";
1361                open(P, "$cmd|") or die "get_translation: cmd failed with $!: $cmd";
1362                $_ = <P>;
1363                my $out;
1364                while (<P>)
1365                {
1366                    s/\s+//g;
1367                    $out .= $_;
1368                }
1369                close(P);
1370                $self->{_feat}->{$peg} = $out;
1371                return $out;
1372            }
1373            else
1374            {
1375                return $self->{_feat}->{$peg};
1376            }
1377      }      }
1378      else      else
1379      {      {
# Line 1234  Line 1391 
1391    
1392      if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))      if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1393      {      {
1394          &load_pseq($self);          my $t = $self->get_translation($peg);
1395          return length($self->{_pseq}->{$peg});          return length($t);
1396      }      }
1397      else      else
1398      {      {
# Line 1253  Line 1410 
1410    
1411      if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))      if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1412      {      {
1413          &load_pseq($self);          return $self->translation_length($peg) > 0 ? 1 : 0;
         return length($self->{_pseq}->{$peg}) > 0 ? 1 : 0;;  
1414      }      }
1415      else      else
1416      {      {
# Line 1298  Line 1454 
1454          return $fig->pegs_of($genome);          return $fig->pegs_of($genome);
1455      }      }
1456    
1457        $self->load_feature_indexes();
1458    
1459        my $t = $self->{_feat_btree}->{peg};
1460        if ($t)
1461        {
1462            return keys %$t;
1463        }
1464        else
1465        {
1466            warn "Reverting to load_tbl for $newG\n";
1467      $self->load_tbl();      $self->load_tbl();
1468      return grep { /\.peg\./ } keys %{$self->{_tbl}};      return grep { /\.peg\./ } keys %{$self->{_tbl}};
1469  }  }
1470    }
1471    
1472    
1473  sub rnas_of  sub rnas_of
# Line 1313  Line 1480 
1480    
1481      if ($genome ne $newG)      if ($genome ne $newG)
1482      {      {
1483          return $fig->pegs_of($genome);          return $fig->rna_of($genome);
1484      }      }
1485    
1486        $self->load_feature_indexes();
1487    
1488        my $t = $self->{_feat_btree}->{rna};
1489        if ($t)
1490        {
1491            return keys %$t;
1492        }
1493        else
1494        {
1495            warn "Reverting to load_tbl for $newG\n";
1496      $self->load_tbl();      $self->load_tbl();
1497      return grep { /\.rna\./ } keys %{$self->{_tbl}};      return grep { /\.rna\./ } keys %{$self->{_tbl}};
1498  }  }
1499    
1500    }
1501    
1502  sub is_virtual_feature  sub is_virtual_feature
1503  {  {
1504      my($self, $peg) = @_;      my($self, $peg) = @_;
# Line 1330  Line 1509 
1509    
1510      if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))      if (($peg =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1511      {      {
         &load_pseq($self);  
1512          return 1;          return 1;
1513      }      }
1514      else      else
# Line 1382  Line 1560 
1560      return @vals;      return @vals;
1561  }  }
1562    
1563    sub load_feature_indexes
1564    {
1565        my($self) = @_;
1566    
1567        if ($self->{_feat}) { return };
1568    
1569        my $newGdir = $self->{_orgdir};
1570    
1571        for my $fdir (<$newGdir/Features/*>)
1572        {
1573            my $ftype = basename($fdir);
1574    
1575            #
1576            # If we have a tbl.btree, tie that for our use.
1577            #
1578    
1579            my $tbl_idx = {};
1580            my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1581            if ($tie)
1582            {
1583                $self->{_feat_tie}->{$ftype} = $tie;
1584                $self->{_feat_btree}->{$ftype} = $tbl_idx;
1585            }
1586    
1587            my $tbl_list = [];
1588            my $ltie = tie @$tbl_list, 'DB_File', "$fdir/tbl.recno", O_RDONLY, 0666, $DB_RECNO;
1589            if ($tie)
1590            {
1591                $self->{_feat_ltie}->{$ftype} = $ltie;
1592                $self->{_feat_recno}->{$ftype} = $tbl_list;
1593            }
1594    
1595            #
1596            # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1597            #
1598    
1599            my $pseq     = {};
1600    
1601            if (-f "$fdir/fasta.norm.phr")
1602            {
1603                $self->{_feat_fasta}->{$ftype} = "$fdir/fasta.norm";
1604            }
1605            else
1606            {
1607                #
1608                # Otherwise, we need to load the data.
1609                #
1610    
1611                if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1612                {
1613                    local $/ = "\n>";
1614                    my $x;
1615                    while (defined($x = <FASTA>))
1616                    {
1617                        chomp $x;
1618                        if ($x =~ />?(\S+)[^\n]*\n(.*)/s)
1619                        {
1620                            my $peg = $1;
1621                            my $seq = $2;
1622                            $seq =~ s/\s//gs;
1623                            $pseq->{$peg} = $seq;
1624                        }
1625                    }
1626                    close(FASTA);
1627                }
1628            }
1629            $self->{_feat}->{$ftype} = $pseq;
1630        }
1631    }
1632    
1633  sub load_pseq {  sub load_pseq {
1634      my($self) = @_;      my($self) = @_;
1635    
1636      if ($self->{_pseq}) { return };      if ($self->{_pseq}) { return };
1637    
1638      my $newGdir = $self->{_orgdir};      my $newGdir = $self->{_orgdir};
1639        my $fdir = "$newGdir/Features/peg";
1640    
1641        #
1642        # If we have a tbl.btree, tie that for our use.
1643        #
1644    
1645        my $tbl_idx = {};
1646        my $tie = tie %$tbl_idx, 'DB_File', "$fdir/tbl.btree", O_RDONLY, 0666, $DB_BTREE;
1647        if ($tie)
1648        {
1649            $self->{_tbl_tie} = $tie;
1650            $self->{_tbl_btree} = $tbl_idx;
1651        }
1652    
1653        #
1654        # If we have fasta.norm.phr, set _pseq_fasta to the fasta file to use with fastacmd.
1655        #
1656    
1657      my $pseq     = {};      my $pseq     = {};
1658    
1659        if (-f "$fdir/fasta.norm.phr")
1660        {
1661            $self->{_pseq_fasta} = "$fdir/fasta.norm";
1662        }
1663        else
1664        {
1665            #
1666            # Otherwise, we need to load the data.
1667            #
1668    
1669      if (open(FASTA,"<$newGdir/Features/peg/fasta"))      if (open(FASTA,"<$newGdir/Features/peg/fasta"))
1670      {      {
1671          local $/ = "\n>";          local $/ = "\n>";
# Line 1406  Line 1683 
1683          }          }
1684          close(FASTA);          close(FASTA);
1685      }      }
1686        }
1687      $self->{_pseq} = $pseq;      $self->{_pseq} = $pseq;
1688  }  }
1689    
# Line 1490  Line 1768 
1768      my $newGdir = $self->{_orgdir};      my $newGdir = $self->{_orgdir};
1769    
1770      my $functions  = {};      my $functions  = {};
1771    
1772        my $bfile = "$newGdir/assigned_functions.btree";
1773        my $tie;
1774        if (-f $bfile)
1775        {
1776            $tie = tie %$functions, 'DB_File', $bfile, O_RDONLY, 0666, $DB_BTREE;
1777        }
1778    
1779        if (!$tie)
1780        {
1781      foreach my $x (`cat $newGdir/*functions`)      foreach my $x (`cat $newGdir/*functions`)
1782      {      {
1783          if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))          if (($x =~ /^(fig\|(\d+\.\d+)\.\S+)\t(\S[^\t]*\S)/) && ($2 eq $newG))
# Line 1497  Line 1785 
1785              $functions->{$1} = $3;              $functions->{$1} = $3;
1786          }          }
1787      }      }
1788        }
1789    
1790      $self->{_functions} = $functions;      $self->{_functions} = $functions;
1791  }  }
1792    

Legend:
Removed from v.1.54  
changed lines
  Added in v.1.55

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3