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

Diff of /FigKernelPackages/SeedV.pm

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

revision 1.5, Thu Aug 19 19:39:41 2010 UTC revision 1.6, Fri Aug 27 19:23:51 2010 UTC
# Line 474  Line 474 
474      return $self->{_contig_len_cache}->{$contig};      return $self->{_contig_len_cache}->{$contig};
475  }  }
476    
477    sub contig_entrypoints
478    {
479        my($self) = @_;
480        $self->load_tbl();
481        return $self->{_contig_entrypoints};
482    }
483    
484  sub contigs_of  sub contigs_of
485  {  {
486      my ($self) = @_;      my ($self) = @_;
# Line 687  Line 694 
694              my $tblH = $self->{_tbl};              my $tblH = $self->{_tbl};
695              while ( my($fid,$tuple) = each %$tblH)              while ( my($fid,$tuple) = each %$tblH)
696              {              {
697                    next if $self->is_deleted_fid($fid);
698                  if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))                  if (($tuple->[0]->[0] =~ /^(\S+)_(\d+)_\d+$/) && ($1 eq $contig))
699                  {                  {
700                      my $beg1 = $2;                      my $beg1 = $2;
# Line 996  Line 1004 
1004    
1005      if (($feature_id =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))      if (($feature_id =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))
1006      {      {
1007          warn "add_annotation on non-seedv fid\n";          warn "add_annotation on non-seedv fid '$feature_id'\n";
1008          return 0;          return 0;
1009      }      }
1010    
# Line 1021  Line 1029 
1029          my $ann = $self->{_ann};          my $ann = $self->{_ann};
1030          push(@{$ann->{$feature_id}}, [$feature_id, $time_made, $user, $annotation . "\n"]);          push(@{$ann->{$feature_id}}, [$feature_id, $time_made, $user, $annotation . "\n"]);
1031      }      }
1032        else
1033        {
1034            warn "Cannot write $file: $!";
1035        }
1036      return 0;      return 0;
1037  }  }
1038    
# Line 1445  Line 1457 
1457      my $newGdir = $self->{_orgdir};      my $newGdir = $self->{_orgdir};
1458      my $tbl     = {};      my $tbl     = {};
1459    
1460        my $contig_entries = [];
1461        my $cur_contig;
1462    
1463      for my $tbl_file (<$newGdir/Features/*/tbl>)      for my $tbl_file (<$newGdir/Features/*/tbl>)
1464      {      {
1465            my($type) = $tbl_file =~ m,/([^/]+)/tbl$,;
1466            print "Load $type\n";
1467          if (open(my $fh, "<", $tbl_file))          if (open(my $fh, "<", $tbl_file))
1468          {          {
1469              while (defined(my $x = <$fh>))              while (defined(my $x = <$fh>))
# Line 1456  Line 1473 
1473                  {                  {
1474                      my $fid = $1;                      my $fid = $1;
1475                      my $loc = [split(/,/,$2)];                      my $loc = [split(/,/,$2)];
1476                        my($contig) = $loc->[0] =~ /^(.*)_\d+_\d+$/;
1477                      my $aliases = $4 ? [split(/\t/,$4)] : [];                      my $aliases = $4 ? [split(/\t/,$4)] : [];
1478                      if (! $self->is_deleted_fid($fid))                      if (! $self->is_deleted_fid($fid))
1479                      {                      {
1480                          $tbl->{$fid} = [$loc,$aliases];                          $tbl->{$fid} = [$loc,$aliases];
1481                            if ($type eq 'peg' && $contig ne $cur_contig)
1482                            {
1483                                push(@$contig_entries, [$contig, $fid]);
1484                                $cur_contig = $contig;
1485                            }
1486                      }                      }
1487                  }                  }
1488                  else {                  else {
# Line 1474  Line 1497 
1497          }          }
1498      }      }
1499      print STDERR ("Loaded ", (scalar keys %$tbl), " features from $newGdir\n") if $ENV{FIG_VERBOSE};      print STDERR ("Loaded ", (scalar keys %$tbl), " features from $newGdir\n") if $ENV{FIG_VERBOSE};
1500    
1501        $self->{_contig_entrypoints} = $contig_entries;
1502    
1503      $self->{_tbl} = $tbl;      $self->{_tbl} = $tbl;
1504  }  }
1505    
# Line 1586  Line 1612 
1612      }      }
1613  }  }
1614    
1615  sub is_deleted_fid {  sub is_deleted_fid
1616    {
1617        my($self, $fid) = @_;
1618        my $newG    = $self->{_genome};
1619    
1620        if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 eq $newG))
1621        {
1622            my $delH;
1623            $delH = $self->{_deleted_features};
1624            if (! $delH)
1625            {
1626                $delH = $self->load_deleted_features_hash();
1627            }
1628            return $delH->{$fid} ? 1 : 0;
1629        }
1630      return 0;      return 0;
1631  }  }
1632    
# Line 1599  Line 1638 
1638    
1639      my $dir = $self->{_orgdir};      my $dir = $self->{_orgdir};
1640      my $index = {};      my $index = {};
1641        my $rev_index = {};
1642    
1643      for my $cfile (<$dir/CorrToReferenceGenomes/*>)      for my $cfile (<$dir/CorrToReferenceGenomes/*>)
1644      {      {
# Line 1610  Line 1650 
1650                  {                  {
1651                      my $ent = CorrTableEntry->new($_);                      my $ent = CorrTableEntry->new($_);
1652                      push(@{$index->{$ent->id1}}, $ent);                      push(@{$index->{$ent->id1}}, $ent);
1653                        push(@{$rev_index->{$ent->id2}}, $ent);
1654                  }                  }
1655                  close($fh);                  close($fh);
1656              }              }
# Line 1620  Line 1661 
1661          }          }
1662      }      }
1663      $self->{correspondence_index} = $index;      $self->{correspondence_index} = $index;
1664        $self->{correspondence_index_rev} = $rev_index;
1665  }  }
1666    
1667  sub get_correspondences  sub get_correspondences
# Line 1631  Line 1673 
1673      return $self->{correspondence_index}->{$id};      return $self->{correspondence_index}->{$id};
1674  }  }
1675    
1676    sub get_correspondences_rev
1677    {
1678        my($self, $id) = @_;
1679    
1680        $self->load_correspondences();
1681    
1682        return $self->{correspondence_index_rev}->{$id};
1683    }
1684    
1685  sub get_best_correspondences  sub get_best_correspondences
1686  {  {
1687      my($self, $id, $cutoff) = @_;      my($self, $id, $cutoff) = @_;
# Line 1642  Line 1693 
1693      for my $hit (sort { SeedUtils::genome_of($a->id2) <=> SeedUtils::genome_of($b->id2) or      for my $hit (sort { SeedUtils::genome_of($a->id2) <=> SeedUtils::genome_of($b->id2) or
1694                          $a->psc <=> $b->psc }                          $a->psc <=> $b->psc }
1695                   grep { (!defined($cutoff) or $_->psc < $cutoff) and                   grep { (!defined($cutoff) or $_->psc < $cutoff) and
1696                          ($_->hitinfo eq '<=>' or $_->func1 eq $_->func2) }                          ($_->hitinfo eq '<=>' or SeedUtils::strip_func_comment($_->func1) eq SeedUtils::strip_func_comment($_->func2)) }
1697                   @$corrs)                   @$corrs)
1698      {      {
1699          if (defined($lg) && $lg ne SeedUtils::genome_of($hit->id2))          if (defined($lg) && $lg ne SeedUtils::genome_of($hit->id2))
# Line 1685  Line 1736 
1736                                                            int($center - $width / 2),                                                            int($center - $width / 2),
1737                                                            int($center + $width / 2));                                                            int($center + $width / 2));
1738    
1739      my $left_extent = $center - $minV;  #    my $left_extent = $center - $minV;
1740      my $right_extent = $maxV - $center;  #    my $right_extent = $maxV - $center;
1741        my $left_extent = int($width / 2);
1742        my $right_extent = int($width / 2);
1743    
1744      #      #
1745      # Determine other context. Get locations for the pin, then find the      # Determine other context. Get locations for the pin, then find the
# Line 1761  Line 1814 
1814          my @row;          my @row;
1815          for my $peg (keys %$region)          for my $peg (keys %$region)
1816          {          {
1817                my $func = defined($funcs->{$peg}) ? $funcs->{$peg} : "";
1818              my($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($region->{$peg});              my($contig, $min, $max, $dir) = &SeedUtils::boundaries_of($region->{$peg});
1819              my $ent = [$peg, $contig, $min, $max, $dir, $funcs->{$peg}, $row];              my $ent = [$peg, $contig, $min, $max, $dir, $func, $row];
1820              $cref->{$peg} = $ent;              $cref->{$peg} = $ent;
1821              push(@row, $ent);              push(@row, $ent);
1822          }          }
# Line 1774  Line 1828 
1828      return \@context, $cref, $genome_names;      return \@context, $cref, $genome_names;
1829  }  }
1830    
1831    sub add_feature {
1832        my( $self, $user, $genome, $type, $location, $aliases, $sequence) = @_;
1833    
1834        my $newG    = $self->{_genome};
1835        my $newGdir = $self->{_orgdir};
1836    
1837        if ($genome ne $newG) {
1838            return undef;
1839        }
1840    
1841        # perform sanity checks
1842        unless ($genome && $type && $location && $sequence) {
1843            print STDERR "SEED error: add_feature failed due to missing parameter\n";
1844            return undef;
1845        }
1846    
1847        if ( $type !~ /^[0-9A-Za-z_]+$/ )
1848        {
1849            print STDERR "SEED error: add_feature failed due to bad type: $type\n";
1850            return undef;
1851        }
1852    
1853        if ( length ( $location ) > 5000 )
1854        {
1855            print STDERR "SEED error: add_feature failed because location is over 5000 char:\n";
1856            print STDERR "$location\n";
1857            return undef;
1858        }
1859    
1860        my @loc  = split( /,/, $location );
1861        my @loc2 = grep { $_->[0] && $_->[1] && $_->[2] }
1862                   map  { [ $_ =~ m/^(.+)_(\d+)_(\d+)$/ ] }
1863                   @loc;
1864    
1865        if ( (@loc2 == 0) || ( @loc != @loc2 ) )
1866        {
1867            print STDERR "SEED error: add_feature failed because location is missing or malformed:\n";
1868            print STDERR "$location\n";
1869            return undef;
1870        }
1871    
1872        if ( my @bad_names = grep { length( $_->[0] ) > 96 } @loc2 )
1873        {
1874            print STDERR "SEED error: add_feature failed because location contains a contig name of over 96 char:\n";
1875            print STDERR join( ", ", @bad_names ) . "\n";
1876            return undef;
1877        }
1878    
1879        # create type directory if it does not exist
1880        unless (-d "$newGdir/Features/$type") {
1881            &SeedUtils::verify_dir("$newGdir/Features/$type");
1882            (-d "$newGdir/Features/$type")
1883                || die qq(Feature directory \'$newGdir/Features/$type\' does not exist, and could not be created);
1884    
1885            open(TMP, ">", "$newGdir/Features/$type/tbl")
1886                || die "Could not create empty $newGdir/Features/$type/tbl: $!";
1887            close(TMP);
1888    
1889            open(TMP, ">", "$newGdir/Features/$type/fasta")
1890                || die "Could not create empty $newGdir/Features/$type/fasta: $!";
1891            close(TMP);
1892        }
1893    
1894        # create an id
1895        my $id = "fig|$genome.$type.";
1896        my $feature_dir = "$newGdir/Features/$type";
1897        my $file = "$feature_dir/tbl";
1898        if (-f $file) {
1899            unless (open(FILE, "<$file")) {
1900                print STDERR "SEED error: could not open tbl file: $@\n";
1901                return undef;
1902            }
1903            my $entry;
1904            my $max = 0;
1905            while (defined($entry = <FILE>)) {
1906                chomp $entry;
1907                if ($entry =~ /^fig\|$genome\.$type\.(\d+)/) {
1908                    my $curr_id = $1;
1909                    if ($curr_id > $max) {
1910                        $max = $curr_id;
1911                    }
1912                }
1913                else {
1914                    confess qq(Could not parse $type tbl entry: $entry);
1915                }
1916            }
1917            close FILE;
1918            ++$max;
1919            $id .= $max;
1920        } else {
1921            $id .= "1";
1922        }
1923    
1924        # append to tbl file
1925        unless (open(FILE, ">>$file")) {
1926            print STDERR "SEED error: could not open tbl file: $@\n";
1927            return undef;
1928        }
1929    
1930        $aliases =~ s/,/\t/g;
1931        print FILE "$id\t$location\t$aliases\n";
1932        close FILE;
1933        chmod(0777,$file);
1934    
1935        my $typeH = $self->{_features}->{$type};
1936        if ($typeH)
1937        {
1938            $typeH->{$id} = 1;
1939        }
1940    
1941        my $tbl = $self->{_tbl};
1942        if ($tbl)
1943        {
1944            $tbl->{$id} = [[@loc],[split(/\t/,$aliases)]];
1945        }
1946    
1947        # append to fasta file
1948        $sequence =~ s/\s//g;
1949    
1950        $file = "$feature_dir/fasta";
1951        unless (open(FILE, ">>$file")) {
1952            print STDERR "SEED error: could not open fasta file: $@\n";
1953            return undef;
1954        }
1955        print FILE ">$id\n$sequence\n";
1956        close FILE;
1957        chmod(0777,$file);
1958    
1959        # append to called_by file
1960        $file = "$newGdir/called_by";
1961        unless (open(FILE, ">>$file")) {
1962            print STDERR "SEED error: could not open called_by file: $@\n";
1963            return undef;
1964        }
1965        print FILE "$id\tmanual: $user\n";
1966        close FILE;
1967        chmod(0777,$file);
1968    
1969        #
1970        # If a btree was created for this, clear the ref to it and delete the files since
1971        # they are now invalid.
1972        #
1973    
1974        my $tie = $self->{_feat_tie}->{$type};
1975        my $btree = $self->{_feat_btree}->{$type};
1976    
1977        if ($tie)
1978        {
1979            untie $tie;
1980            delete $self->{$_}->{$type} for qw(_feat_tie _feat_btree _feat_ltie _feat_recno _feat_fasta);
1981    
1982            unlink("$feature_dir/tbl.btree");
1983            unlink("$feature_dir/tbl.recno");
1984        }
1985    
1986        if (-f "$feature_dir/fasta.norm.phr")
1987        {
1988            unlink(<$feature_dir/fasta.norm.*>);
1989        }
1990    
1991    
1992        # declare success
1993        return $id;
1994    }
1995    
1996    sub delete_feature {
1997        my($self,$user,$fid) = @_;
1998    
1999        my $newG    = $self->{_genome};
2000        my $newGdir = $self->{_orgdir};
2001    
2002        if (($fid =~ /^fig\|(\d+\.\d+)/) && ($1 ne $newG))
2003        {
2004            warn "delete_feature on non-seedv fid\n";
2005            return 0;
2006        }
2007    
2008        if (open(DEL,">>$newGdir/deleted.fids"))
2009        {
2010            print DEL "$fid\t$user\n";
2011            close(DEL);
2012            $self->load_deleted_features_hash();
2013        }
2014        else
2015        {
2016            carp "could not open $newGdir/deleted.fids: failed to delete $fid (user=$user)\n";
2017        }
2018    }
2019    
2020    sub load_deleted_features_hash {
2021        my($self) = @_;
2022    
2023        my $newGdir = $self->{_orgdir};
2024        my $deletedH = {};
2025        if (open(DEL,"<$newGdir/deleted.fids"))
2026        {
2027            local $/ = "\n";
2028            while (<DEL>)
2029            {
2030                if ($_ =~ /^(\S+)/)
2031                {
2032                    $deletedH->{$1} = 1;
2033                }
2034            }
2035            close(DEL);
2036        }
2037        $self->{_deleted_features} = $deletedH;
2038    
2039        return $deletedH;
2040    }
2041    
2042    #
2043    # List of feature types that exist in this genome.
2044    #
2045    sub feature_types
2046    {
2047        my($self) = @_;
2048    
2049        my $newGdir = $self->{_orgdir};
2050    
2051        my %sort_prefs = (peg => 2,
2052                          rna => 1);
2053    
2054        if (opendir(my $dh, "$newGdir/Features"))
2055        {
2056            my @types = grep { ! /^\./ && -d "$newGdir/Features/$_" } readdir($dh);
2057            closedir($dh);
2058            return sort { $sort_prefs{$b} <=> $sort_prefs{$a} or $a cmp $b } @types;
2059        }
2060        return ();
2061    }
2062    
2063  1;  1;

Legend:
Removed from v.1.5  
changed lines
  Added in v.1.6

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3