[Bio] / Kmers2 / build_data_directory.pl Repository:
ViewVC logotype

Diff of /Kmers2/build_data_directory.pl

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

revision 1.2, Sun Jan 20 14:44:16 2013 UTC revision 1.5, Mon Jul 22 20:00:15 2013 UTC
# Line 1  Line 1 
1  use strict;  use strict;
2  use Data::Dumper;  use Data::Dumper;
3  use SeedEnv;  use SeedEnv;
4    use FIG;
5    my $fig = new FIG;
6    
7  my $usage = "usage: build_data_directory DataDir K";  my $usage = "usage: build_data_directory DataDir K";
8  my($genus,$dataD,$k);  my($genus,$dataD,$k);
# Line 12  Line 14 
14  print STDERR "building $dataD for Kmers of size $k\n";  print STDERR "building $dataD for Kmers of size $k\n";
15  mkdir($dataD,0777);  mkdir($dataD,0777);
16    
17    my %seen;
18    my $existing = $ENV{'SAS_SERVER'};
19    $ENV{'SAS_SERVER'} = 'PUBSEED';
20    open(G,">$dataD/all.genomes") || die "could not open $dataD/all.genomes";
21    foreach $_ (`svr_all_genomes`)
22    {
23        if (($_ !~ /[pP]hage\b/) && ($_ !~ /[pP]lasmid\b/))
24        {
25            if (($_ =~ /^(\S+)\s(\S+).*\s(\d+\.\d+)$/) && ($2 ne "sp."))
26            {
27                my $genus = $1;
28                my $species = $2;
29                my $genome = $3;
30                if ($species !~ /^sp\.?$/i)
31                {
32                    my $md5 = $fig->genome_md5sum($genome);
33                    if ((! $seen{$md5}) && $fig->is_prokaryotic($genome))
34                    {
35                        $seen{$md5} = 1;
36                        print G $_;
37                    }
38                }
39            }
40        }
41    }
42    close(G);
43    $ENV{'SAS_SERVER'} = $existing;
44    
45  &build_function_index($dataD);  &build_function_index($dataD);
46  &build_otu_index($dataD);  &build_otu_index($dataD);
47  &build_reduced_kmers($dataD,$k);  &build_reduced_kmers($dataD,$k);
# Line 27  Line 57 
57      my %funcs;      my %funcs;
58      while (defined($_ = <ASSIGNMENTS>))      while (defined($_ = <ASSIGNMENTS>))
59      {      {
60          if ($_ =~ /^\S+\t(\S[^\#]*\S+)/)          if ($_ =~ /^\S+\t(\S.*\S+)/)
61          {          {
62              $funcs{$1}++;              my $stripped = &SeedUtils::strip_func($1);   #### CHECK THIS
63                $funcs{$stripped}++;
64          }          }
65      }      }
66      close(ASSIGNMENTS);      close(ASSIGNMENTS);
# Line 107  Line 138 
138              $to_fI{$2} = $1;              $to_fI{$2} = $1;
139          }          }
140      }      }
141    #
142    #   we take assignments from both the pubSEED and the ASEED.  ASEED functions
143    #   override PSEED functions
144    #
145      my %peg_to_fI;      my %peg_to_fI;
146      foreach $_ (`cut -f2 $dataD/all.genomes | svr_all_features peg | svr_function_of`)      &load_peg_to_fI($dataD,\%to_fI,'PUBSEED',\%peg_to_fI);
147      {      &load_peg_to_fI($dataD,\%to_fI,'SEED',\%peg_to_fI);
         if ($_ =~ /^(\S+)\t(\S.*\S)/)  
         {  
             my $peg = $1;  
             my $func = $2;  
             if ($to_fI{$func})  
             {  
                 $peg_to_fI{$peg} = $to_fI{$func};  
             }  
         }  
     }  
148    
149      open(RAW,"| sort -T . > $dataD/sorted.kmers") || die "could not open $dataD/sorted.kmers";      open(RAW,"| sort -T . > $dataD/sorted.kmers") || die "could not open $dataD/sorted.kmers";
150      foreach my $g (`cut -f2 $dataD/all.genomes`)      foreach my $g (`cut -f2 $dataD/all.genomes`)
# Line 167  Line 191 
191      close(REDUCED);      close(REDUCED);
192  }  }
193    
194    sub load_peg_to_fI {
195        my($dataD,$to_fI,$which_seed,$peg_to_fI) = @_;
196        my $existing = $ENV{'SAS_SERVER'};
197        $ENV{'SAS_SERVER'} = $which_seed;
198    
199        foreach $_ (`cut -f2 $dataD/all.genomes | svr_all_features peg | svr_function_of`)
200        {
201            if ($_ =~ /^(\S+)\t(\S.*\S)/)
202            {
203                my $peg = $1;
204                my $func = &SeedUtils::strip_func($2);
205                if ($to_fI->{$func})
206                {
207                    $peg_to_fI->{$peg} = $to_fI->{$func};
208                }
209            }
210        }
211        $ENV{'SAS_SERVER'} = $existing;
212    }
213    
214    
215  sub process_set {  sub process_set {
216      my($kmer,$set,$fh) = @_;      my($kmer,$set,$fh) = @_;
217    
# Line 177  Line 222 
222          $funcs{$fI}++;          $funcs{$fI}++;
223      }      }
224      my @tmp = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);      my @tmp = sort { $funcs{$b} <=> $funcs{$a} } keys(%funcs);
225      if ($tmp[0] && ($funcs{$tmp[0]} >= (0.8 * @$set)))      if ($tmp[0] && ($funcs{$tmp[0]} > (0.5 * @$set)))
226      {      {
227          my $best_fI = $tmp[0];          my $best_fI = $tmp[0];
228          my %otus;          my %otus;

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3