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

Diff of /FigKernelScripts/find_neighbors_using_figfams.pl

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

revision 1.5, Wed Jan 30 22:50:37 2008 UTC revision 1.9, Fri May 15 23:29:42 2009 UTC
# Line 23  Line 23 
23  my $fig = new FIG;  my $fig = new FIG;
24    
25  use FigFams;  use FigFams;
 my $figfams = new FigFams($fig);  
   
 print STDERR "$0 using figfams in $figfams->{dir}\n";  
   
26  use ToCall;  use ToCall;
27  use NewGenome;  use NewGenome;
28    
# Line 41  Line 37 
37    
38  my ($i, $fam_data, $to_call_dir, $N, $foundF);  my ($i, $fam_data, $to_call_dir, $N, $foundF);
39    
40  # First, let's determine which FigfamsData directory to use  #...First, let's determine which FigfamsData directory to use:
 $fam_data = "$FIG_Config::data/FigfamsData";  
41  for ($i=0; ($i < @ARGV) && ($ARGV[$i] !~ /^-fams=/); $i++) {}  for ($i=0; ($i < @ARGV) && ($ARGV[$i] !~ /^-fams=/); $i++) {}
42  if ($i < @ARGV) {  if ($i < @ARGV) {
43      if ($ARGV[$i] =~ /^-fams=(\S+)/) {      if ($ARGV[$i] =~ /^-fams=(\S+)/) {
# Line 56  Line 51 
51      }      }
52  }  }
53    
54    my $figfams = new FigFams($fig, $fam_data);
55    
56  # Now pick up the normal parameters  # Now pick up the normal parameters
57  (  (
58   ($to_call_dir = shift @ARGV) &&   ($to_call_dir = shift @ARGV) &&
59   ($N           = shift @ARGV) &&   ($N           = shift @ARGV) &&
60   ($foundF      = shift @ARGV) && open(FOUNDFAMS, ">>$foundF")   ($foundF      = shift @ARGV)
61  )  )
62      || die "\n   $usage\n\n";      || die "\n   $usage\n\n";
63    
# Line 71  Line 68 
68  # we are doing) you would use ToCall.pm.  # we are doing) you would use ToCall.pm.
69    
70  my $to_call;  my $to_call;
71  my $old;  my $keep_original_calls = 0;
   
72  if ((@ARGV > 0) && ($ARGV[0] =~ /^old/i)) {  if ((@ARGV > 0) && ($ARGV[0] =~ /^old/i)) {
73      $old = 1;      print STDERR qq(Re-annotating original calls\n) if $ENV{VERBOSE};
74        $keep_original_calls = 1;
75      $to_call = new ToCall($to_call_dir);      $to_call = new ToCall($to_call_dir);
76  }  }
77  else {  else {
78      $to_call = new NewGenome($to_call_dir);      $to_call = new NewGenome($to_call_dir);
79  }  }
80    
81    
82    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
83    #...Open auxiliary files (unless "Keep Original Calls" selected)
84    #-----------------------------------------------------------------------
85    open(FOUNDFAMS, ">>$foundF") || die "Could not append-open $foundF";
86    
87    my $called_by_file = $keep_original_calls ? qq(/dev/null)
88                                         : qq($to_call_dir/called_by);
89    open(CALLED_BY, ">>$called_by_file")
90        || die "Could not append-open $called_by_file";
91    
92    
93    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
94    #...Get the current ORF-calls if they exist;
95    #   else, try to call a set of candidate ORFs.
96    #-----------------------------------------------------------------------
97  my @orf_ids = $to_call->get_fids_for_type('orf');  my @orf_ids = $to_call->get_fids_for_type('orf');
98  if (not @orf_ids) {  if (not @orf_ids) {
99      $to_call->possible_orfs;      $to_call->possible_orfs();
100      @orf_ids = $to_call->get_fids_for_type('orf');      @orf_ids = $to_call->get_fids_for_type('orf');
101  }  }
102    
# Line 95  Line 108 
108      }      }
109  }  }
110    
111    
112  my @peg_ids = ();  my @peg_ids = ();
113  my %genomes_hit;  my %genomes_hit;
114  my %weight_of_hits;  my %weight_of_hits;
# Line 115  Line 129 
129              if $ENV{VERBOSE};              if $ENV{VERBOSE};
130    
131          my $fid;          my $fid;
132          if (! $old)          if (! $keep_original_calls)
133          {          {
134                my $annot = [ qq(RAST),
135                              qq($func\nCalled by "$self" using FIGfam $fam_id.)
136                              ];
137              my $orf = &NewGenome::ORF::new('NewGenome::ORF', $to_call, $orf_id);              my $orf = &NewGenome::ORF::new('NewGenome::ORF', $to_call, $orf_id);
138              $fid = $orf->promote_to_peg($sims, $func);              $fid      = $orf->promote_to_peg($sims, $func, $annot);
139    
140                if ($fid) {
141                    print CALLED_BY "$fid\t$self\n";
142                    print FOUNDFAMS "$fid\t$fam_id\t$func\n";
143          }          }
144          else              else {
145          {                  die "Could not promote $orf_id";
146                }
147            }
148            else {
149              $fid = $orf_id;              $fid = $orf_id;
150          }          }
151    
152          if ($fid) {          if ($fid) {
153              push @peg_ids, $fid;              push @peg_ids, $fid;
154    
             print FOUNDFAMS "$fid\t$fam_id\t$func\n";  
   
155              for ($i=0; ($i < $N) && ($i < @$sims); ++$i) {              for ($i=0; ($i < $N) && ($i < @$sims); ++$i) {
156                  ++$genomes_hit{&FIG::genome_of($sims->[$i]->id2)};                  ++$genomes_hit{&FIG::genome_of($sims->[$i]->id2)};
157                    next unless $sims->[$i]->ln2;
158                  $weight_of_hits{&FIG::genome_of($sims->[$i]->id2)} +=                  $weight_of_hits{&FIG::genome_of($sims->[$i]->id2)} +=
159                      ($sims->[$i]->bsc / $sims->[$i]->ln2);                      ($sims->[$i]->bsc / $sims->[$i]->ln2);
160              }              }
# Line 155  Line 178 
178  for ($i=0; ($i < $N) && ($i < @best); ++$i) {  for ($i=0; ($i < $N) && ($i < @best); ++$i) {
179      print "$best[$i]\t$genomes_hit{$best[$i]}\t$weight_of_hits{$best[$i]}\n";      print "$best[$i]\t$genomes_hit{$best[$i]}\t$weight_of_hits{$best[$i]}\n";
180  }  }
181    
182    $to_call->export_features();
183    
184    close(CALLED_BY);
185  close(FOUNDFAMS);  close(FOUNDFAMS);
 $to_call->export_features;  
186    
187  exit(0);  exit(0);

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3