[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.11, Thu Nov 5 19:45:28 2009 UTC
# Line 20  Line 20 
20  use warnings;  use warnings;
21    
22  use FIG;  use FIG;
23  my $fig = new FIG;  my $fig = FIG->new();
   
 use FigFams;  
 my $figfams = new FigFams($fig);  
   
 print STDERR "$0 using figfams in $figfams->{dir}\n";  
24    
25    use FFs;
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    $fam_data = $fig->get_figfams_data($fam_data);
55    my $figfams = new FFs($fam_data, $fig);
56    
57  # Now pick up the normal parameters  # Now pick up the normal parameters
58  (  (
59   ($to_call_dir = shift @ARGV) &&   ($to_call_dir = shift @ARGV) &&
60   ($N           = shift @ARGV) &&   ($N           = shift @ARGV) &&
61   ($foundF      = shift @ARGV) && open(FOUNDFAMS, ">>$foundF")   ($foundF      = shift @ARGV)
62  )  )
63      || die "\n   $usage\n\n";      || die "\n   $usage\n\n";
64    
# Line 71  Line 69 
69  # we are doing) you would use ToCall.pm.  # we are doing) you would use ToCall.pm.
70    
71  my $to_call;  my $to_call;
72  my $old;  my $keep_original_calls = 0;
   
73  if ((@ARGV > 0) && ($ARGV[0] =~ /^old/i)) {  if ((@ARGV > 0) && ($ARGV[0] =~ /^old/i)) {
74      $old = 1;      print STDERR qq(Re-annotating original calls\n) if $ENV{VERBOSE};
75      $to_call = new ToCall($to_call_dir);      $keep_original_calls = 1;
76        $to_call = ToCall->new($to_call_dir);
77  }  }
78  else {  else {
79      $to_call = new NewGenome($to_call_dir);      $to_call = NewGenome->new($to_call_dir);
80  }  }
81    
82    
83    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
84    #...Open auxiliary files (unless "Keep Original Calls" selected)
85    #-----------------------------------------------------------------------
86    open(FOUNDFAMS, ">>$foundF") || die "Could not append-open $foundF";
87    
88    my $called_by_file = 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 114  Line 128 
128              , " sims)\t$func\n"              , " sims)\t$func\n"
129              if $ENV{VERBOSE};              if $ENV{VERBOSE};
130    
131          my $fid;          my $annot = [ qq(RAST),
132          if (! $old)                        qq($func\nCalled by "$self" using FIGfam $fam_id.)
133          {                        ];
134              my $orf = &NewGenome::ORF::new('NewGenome::ORF', $to_call, $orf_id);  
135              $fid = $orf->promote_to_peg($sims, $func);          my $orf;
136          }          if ($keep_original_calls) {
137          else              $orf = &ToCall::PEG::new(   'ToCall::PEG',    $to_call, $orf_id, $seq_of{$orf_id});
138          {          }
139              $fid = $orf_id;          else {
140                $orf = &NewGenome::ORF::new('NewGenome::ORF', $to_call, $orf_id);
141          }          }
142    
143          if ($fid) {          my $fid  = $orf->promote_to_peg($sims, $func, $annot);
             push @peg_ids, $fid;  
144    
145            if ($fid) {
146                print CALLED_BY "$fid\t$self\n";
147              print FOUNDFAMS "$fid\t$fam_id\t$func\n";              print FOUNDFAMS "$fid\t$fam_id\t$func\n";
148            }
149            else {
150                die "Could not promote $orf_id";
151            }
152    
153            if ($fid) {
154                push @peg_ids, $fid;
155    
156              for ($i=0; ($i < $N) && ($i < @$sims); ++$i) {              for ($i=0; ($i < $N) && ($i < @$sims); ++$i) {
157                  ++$genomes_hit{&FIG::genome_of($sims->[$i]->id2)};                  my $org_2 = &FIG::genome_of($sims->[$i]->id2);
158                  $weight_of_hits{&FIG::genome_of($sims->[$i]->id2)} +=                  if (not defined($genomes_hit{$org_2})) {
159                      ($sims->[$i]->bsc / $sims->[$i]->ln2);                      $genomes_hit{$org_2}    = 0;
160                        $weight_of_hits{$org_2} = 0;
161                    }
162    
163                    ++$genomes_hit{$org_2};
164    
165                    if ($sims->[$i]->ln2) {
166                        $weight_of_hits{$org_2} +=
167                            ($sims->[$i]->bsc() / $sims->[$i]->ln2());
168                    }
169              }              }
170          }          }
171          else {          else {
# Line 153  Line 185 
185    
186  my @best = sort { $weight_of_hits{$b} <=> $weight_of_hits{$a} } keys(%genomes_hit);  my @best = sort { $weight_of_hits{$b} <=> $weight_of_hits{$a} } keys(%genomes_hit);
187  for ($i=0; ($i < $N) && ($i < @best); ++$i) {  for ($i=0; ($i < $N) && ($i < @best); ++$i) {
188      print "$best[$i]\t$genomes_hit{$best[$i]}\t$weight_of_hits{$best[$i]}\n";      print STDOUT (join(qq(\t), ($best[$i],
189                                    $genomes_hit{$best[$i]},
190                                    $weight_of_hits{$best[$i]},
191                                    $fig->genus_species($best[$i])
192                                    )
193                           ),
194                      qq(\n)
195                      );
196  }  }
197    close(CALLED_BY);
198  close(FOUNDFAMS);  close(FOUNDFAMS);
199  $to_call->export_features;  
200    $to_call->export_features();
201    
202  exit(0);  exit(0);

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3