[Bio] / FigKernelScripts / load_features.pl Repository:
ViewVC logotype

Diff of /FigKernelScripts/load_features.pl

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

revision 1.5, Fri Dec 26 13:01:58 2003 UTC revision 1.20, Mon Dec 5 18:56:37 2005 UTC
# Line 1  Line 1 
1  # -*- perl -*-  # -*- perl -*-
2    #
3    # Copyright (c) 2003-2006 University of Chicago and Fellowship
4    # for Interpretations of Genomes. All Rights Reserved.
5    #
6    # This file is part of the SEED Toolkit.
7    #
8    # The SEED Toolkit is free software. You can redistribute
9    # it and/or modify it under the terms of the SEED Toolkit
10    # Public License.
11    #
12    # You should have received a copy of the SEED Toolkit Public License
13    # along with this program; if not write to the University of Chicago
14    # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15    # Genomes at veronika@thefig.info or download a copy from
16    # http://www.theseed.org/LICENSE.TXT.
17    #
18    
19    
20  use strict;  use strict;
21  use FIG;  use FIG;
22  my $fig = new FIG;  my $fig = new FIG;
23    
24  use DBrtns;  use Tracer;
   
25    
26    Trace("Preparing to load features.") if T(2);
27    my ($mode, @genomes) = FIG::parse_genome_args(@ARGV);
28  my $temp_dir = "$FIG_Config::temp";  my $temp_dir = "$FIG_Config::temp";
29  my($organisms_dir) = "$FIG_Config::organisms";  my $organisms_dir = "$FIG_Config::organisms";
30    
31  my($genome,@types,$type,$id,$loc,@aliases,$aliases,$contig);  my($genome,@types,$type,$id,$loc,@aliases,$aliases,$contig);
32    
33  # usage: load_features [G1 G2 G3 ... ]  # usage: load_features [G1 G2 G3 ... ]
34    
35  open(REL,">$temp_dir/tmpfeat$$") || die "could not open $temp_dir/tmpfeat$$";  Open(\*REL, ">$temp_dir/tmpfeat$$");
36    Open(\*ALIAS, "| sort -T $temp_dir -u > $temp_dir/tmpalias$$");
37    
 my $dbf = $fig->{_dbf};  
38    
39    if ($mode eq 'all') {
40    
41        # Here we extract external aliases from the peg.synonyms table, when they can be inferred
42  my @genomes;      # accurately.
43  if (@ARGV == 0)          Trace("Extracting external aliases from the peg.synonyms table.") if T(2);
44        open(\*SYN, "<$FIG_Config::global/peg.synonyms");
45        while (defined($_ = <SYN>))
46  {  {
47      $dbf->drop_table( tbl => "features" );                  chop;
48      if ($FIG_Config::dbms eq "Pg")                  my($x,$y) = split(/\t/,$_);
49                    my @ids = map { $_ =~ /^([^,]+),(\d+)/; [$1,$2] } ($x,split(/;/,$y));
50                    my @fig = ();
51                    my(@nonfig) = ();
52                    foreach $_ (@ids)
53      {      {
54          $dbf->create_table( tbl  => "features",                          if ($_->[0] =~ /^fig\|/)
55                              flds => "id varchar(32) UNIQUE NOT NULL, type varchar(16),genome varchar(16),"  .                          {
56                                      "location varchar(5000),"  .                                  push(@fig,$_);
                                     "contig varchar(96), minloc INTEGER, maxloc INTEGER,"  .  
                                     "aliases TEXT, PRIMARY KEY ( id )"  
                             );  
57      }      }
58      elsif ($FIG_Config::dbms eq "mysql")                          else
59      {      {
60          $dbf->create_table( tbl  => "features",                                  push(@nonfig,$_);
                             flds => "id varchar(32) UNIQUE NOT NULL, type varchar(16),genome varchar(16),"  .  
                                     "location TEXT,"  .  
                                     "contig varchar(96), minloc INTEGER, maxloc INTEGER,"  .  
                                     "aliases TEXT, PRIMARY KEY ( id )"  
                             );  
61      }      }
   
     @genomes = $fig->genomes;  
62  }  }
63  else  
64                    foreach $x (@fig)
65  {  {
66      @genomes = @ARGV;                          my($peg,$peg_ln) = @$x;
67      foreach $genome (@genomes)                          my $genome = &FIG::genome_of($peg);
68                            foreach $_ (@nonfig)
69                            {
70                                    if ((@fig == 1) || ($peg_ln == $_->[1]))
71      {      {
72          $dbf->SQL("DELETE FROM features WHERE ( genome = \'$genome\' )");                                          print ALIAS "$peg\t$_->[0]\t$genome\n";
73                                            Trace("Alias record $peg, $_->[0] for $genome.") if T(4);
74      }      }
75  }  }
76                    }
77        }
78        close(SYN);
79    }
80    
81  foreach $genome (@genomes)  foreach $genome (@genomes)
82  {  {
83            Trace("Processing $genome.") if T(3);
84      opendir(FEAT,"$organisms_dir/$genome/Features")      opendir(FEAT,"$organisms_dir/$genome/Features")
85          || die "could not open $genome/Features";          || die "could not open $genome/Features";
86      @types = grep { $_ =~ /^[a-zA-Z]+$/ } readdir(FEAT);      @types = grep { $_ =~ /^[a-zA-Z]+$/ } readdir(FEAT);
# Line 66  Line 91 
91          if ((-s "$organisms_dir/$genome/Features/$type/tbl") &&          if ((-s "$organisms_dir/$genome/Features/$type/tbl") &&
92              open(TBL,"<$organisms_dir/$genome/Features/$type/tbl"))              open(TBL,"<$organisms_dir/$genome/Features/$type/tbl"))
93          {          {
94  #           print STDERR "loading $genome/Features/$type/tbl\n";              Trace("Loading $genome/Features/$type/tbl") if T(4);
95              while (defined($_ = <TBL>))              my @tbl = <TBL>;
96                close(TBL);
97                my %seen;
98    
99                while ($_ = pop @tbl)
100              {              {
101                  chop;                  chop;
102                  ($id,$loc,@aliases) = split(/\t/,$_);                  ($id,$loc,@aliases) = split(/\t/,$_);
103                  if ($id)  
104                    if ($id && (! $seen{$id}))
105                  {                  {
106                        $seen{$id} = 1;
107                      my($minloc,$maxloc);                      my($minloc,$maxloc);
108                      if ($loc)                      if ($loc)
109                      {                      {
# Line 92  Line 123 
123                      if (@aliases > 0)                      if (@aliases > 0)
124                      {                      {
125                          $aliases = join(",",grep(/\S/,@aliases));                          $aliases = join(",",grep(/\S/,@aliases));
126                            my $alias;
127                            foreach $alias (@aliases)
128                            {
129                                if ($alias =~ /^([NXYZA]P_|gi\||sp\|\tr\||kegg\||uni\|)/)
130                                {
131    
132                                    print ALIAS "$id\t$alias\t$genome\tOVERRIDE\n";
133                                    Trace("$id override alias $alias for $genome") if T(4);
134                                }
135                            }
136                      }                      }
137                      else                      else
138                      {                      {
# Line 99  Line 140 
140                      }                      }
141                      $minloc = (! $minloc) ? 0 : $minloc;                      $minloc = (! $minloc) ? 0 : $minloc;
142                      $maxloc = (! $maxloc) ? 0 : $maxloc;                      $maxloc = (! $maxloc) ? 0 : $maxloc;
143                      if ((length($loc) < 5000) && (length($contig) < 96) && (length($id) < 32))                      if ((length($loc) < 5000) && (length($contig) < 96) && (length($id) < 32) && ($id =~ /(\d+)$/))
144                      {                      {
145                          print REL "$id\t$type\t$genome\t$loc\t$contig\t$minloc\t$maxloc\t$aliases\n";                          print REL "$id\t$1\t$type\t$genome\t$loc\t$contig\t$minloc\t$maxloc\t$aliases\n";
146                      }                      }
147                  }                  }
148              }              }
             close(TBL);  
149          }          }
150      }      }
151  }  }
152  close(REL);  close(REL);
153    close(ALIAS);
154    Open(\*ALIASIN, "<$temp_dir/tmpalias$$");
155    Open(\*ALIASOUT, ">$temp_dir/tmpalias$$.1");
156    Trace("Parsing alias file.") if T(2);
157    $_ = <ALIASIN>;
158    while ($_ && ($_ =~ /^(\S+)/))
159    {
160        my @aliases = ();
161        my $curr = $1;
162        while ($_ && ($_ =~ /^(\S+)\t(\S+)(\t(\S+))?/) && ($1 eq $curr))
163        {
164                    push(@aliases,[$2,$3 ? 1 : 0]);
165                    $_ = <ALIASIN>;
166        }
167        my $x;
168        my $genome = &FIG::genome_of($curr);
169        foreach $x (@aliases)
170        {
171            if ($x->[1])
172            {
173                print ALIASOUT "$curr\t$x->[0]\t$genome\n";
174            }
175            else
176            {
177                my $i;
178                for ($i=0; ($i < @aliases) && ((! $aliases[$i]->[1]) || (! &same_class($x->[0],$aliases[$i]->[0]))); $i++) {}
179                if ($i == @aliases)
180                {
181                    print ALIASOUT "$curr\t$x->[0]\t$genome\n";
182                }
183            }
184        }
185    }
186    close(ALIASIN);
187    close(ALIASOUT);
188    unlink("$temp_dir/tmpalias$$");
189    
190  $dbf->load_table( tbl => "features",  $fig->reload_table($mode, 'features',
191                    file => "$temp_dir/tmpfeat$$" );                                     "id varchar(32), idN INTEGER, type varchar(16),genome varchar(16),"  .
192                                        "location TEXT,"  .
193                                        "contig varchar(96), minloc INTEGER, maxloc INTEGER,"  .
194                                        "aliases TEXT",
195                                            { features_id_ix => "id", features_org_ix => "genome",
196                                              features_type_ix => "type", features_beg_ix => "genome, contig, minloc" },
197                                            "$temp_dir/tmpfeat$$", \@genomes);
198    unlink("$temp_dir/tmpfeat$$");
199    
200  if (@ARGV == 0)  $fig->reload_table($mode, 'ext_alias',
201  {                                          "id varchar(32), alias varchar(32), genome varchar(16)",
202      $dbf->create_index( idx  => "features_org_ix",                                          { ext_alias_alias_ix => "alias", ext_alias_genome_ix => "genome",
203                          tbl  => "features",                                            ext_alias_ix_id => "id" },
204                          type => "btree",                                          "$temp_dir/tmpalias$$.1", \@genomes );
                         flds => "genome" );  
     $dbf->create_index( idx  => "features_type_ix",  
                         type => "btree",  
                         tbl  => "features",  
                         flds => "type" );  
     $dbf->create_index( idx  => "features_beg_ix",  
                         type => "btree",  
                         tbl  => "features",  
                         flds => "genome,contig,minloc" );  
205    
206      $dbf->vacuum_it("features")  unlink("$temp_dir/tmpalias$$.1");
207    Trace("Features loaded.") if T(2);
208    
209    sub same_class {
210        my($x,$y) = @_;
211    
212        my $class1 = &classA($x);
213        my $class2 = &classA($y);
214        return ($class1 && ($class1 eq $class2));
215    }
216    
217    sub classA {
218        my($alias) = @_;
219    
220        if ($alias =~ /^([^\|]+)\|/)
221        {
222                    return $1;
223        }
224        elsif ($alias =~ /^[NXYZA]P_[0-9\.]+$/)
225        {
226                    return "refseq";
227        }
228        else
229        {
230                    return "";
231        }
232  }  }
 unlink("$temp_dir/tmpfeat$$");  

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3