[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.1, Fri Jun 1 02:14:30 2007 UTC revision 1.8, Wed Feb 13 23:04:18 2008 UTC
# Line 23  Line 23 
23  my $fig = new FIG;  my $fig = new FIG;
24    
25  use FigFams;  use FigFams;
 my $figfams = new FigFams;  
   
26  use ToCall;  use ToCall;
27  use NewGenome;  use NewGenome;
28    
# Line 39  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 54  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 69  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;
72    
73  if ((@ARGV > 0) && ($ARGV[0] =~ /^old/i)) {  if ((@ARGV > 0) && ($ARGV[0] =~ /^old/i)) {
74        $old = 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    open(FOUNDFAMS, ">>$foundF") || die "Could not append-open $foundF";
83    
84    my $called_by_file = "$to_call_dir/called_by";
85    open(CALLED_BY, ">>$called_by_file")
86        || die "Could not append-open $called_by_file";
87    
88    
89    #+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
90    #...Get the current ORF-calls if they exist;
91    #   else, try to call a set of candidate ORFs.
92    #-----------------------------------------------------------------------
93  my @orf_ids = $to_call->get_fids_for_type('orf');  my @orf_ids = $to_call->get_fids_for_type('orf');
94  if (not @orf_ids) {  if (not @orf_ids) {
95      $to_call->possible_orfs;      $to_call->possible_orfs();
96      @orf_ids = $to_call->get_fids_for_type('orf');      @orf_ids = $to_call->get_fids_for_type('orf');
97  }  }
98    
99  my (%seq_of, %len_of, $orf_len);  my (%seq_of, %len_of, $orf_len);
100  for my $fid (@orf_ids) {  for my $fid (@orf_ids) {
101      if ($orf_len = $to_call->get_feature_length($fid) > 900) {      if (($orf_len = $to_call->get_feature_length($fid)) > 900) {
102          $len_of{$fid} = $orf_len;          $len_of{$fid} = $orf_len;
103          $seq_of{$fid} = $to_call->get_feature_sequence($fid);          $seq_of{$fid} = $to_call->get_feature_sequence($fid);
104      }      }
105  }  }
106    
107    
108  my @peg_ids = ();  my @peg_ids = ();
109  my %genomes_hit;  my %genomes_hit;
110  my %weight_of_hits;  my %weight_of_hits;
111  foreach my $orf_id (sort { $len_of{$b} <=> $len_of{$a} } keys %len_of) {  foreach my $orf_id (sort { $len_of{$b} <=> $len_of{$a} } keys %len_of) {
     my ($fam, $sims) = $figfams->place_in_family($seq_of{$orf_id});  
     if (@$sims >= $N) {  
   
         my $fid = $to_call->promote_to_peg($orf_id);  
         push @peg_ids, $fid;  
112    
113        print STDERR "\nChecking $orf_id ($len_of{$orf_id})\n" if $ENV{VERBOSE};
114        my ($fam, $sims) = $figfams->place_in_family($seq_of{$orf_id});
115        if (not defined($fam)) {
116            print STDERR "\tCould not place in family --- skipping\n" if $ENV{VERBOSE};
117            next;
118        }
119        else {
120          my $fam_id = $fam->family_id();          my $fam_id = $fam->family_id();
121          my $func   = $fam->family_function();          my $func   = $fam->family_function();
122            print STDERR "\tplaced in $fam_id ("
123                , (scalar @$sims)
124                , " sims)\t$func\n"
125                if $ENV{VERBOSE};
126    
127            my $fid;
128            if (! $old)
129            {
130                my $annot = [ qq(RAST),
131                              qq($func\nCalled by "$self" using FIGfam $fam_id.)
132                              ];
133                my $orf   = &NewGenome::ORF::new('NewGenome::ORF', $to_call, $orf_id);
134                $fid      = $orf->promote_to_peg($sims, $func, $annot);
135    
136                if ($fid) {
137                    print CALLED_BY "$fid\t$self\n";
138          print FOUNDFAMS "$fid\t$fam_id\t$func\n";          print FOUNDFAMS "$fid\t$fam_id\t$func\n";
139                }
140                else {
141                    die "Could not promote $orf_id";
142                }
143            }
144    
145            if ($fid) {
146                push @peg_ids, $fid;
147    
148          for ($i=0; ($i < $N) && ($i < @$sims); ++$i) {          for ($i=0; ($i < $N) && ($i < @$sims); ++$i) {
149              ++$genomes_hit{&FIG::genome_of($sims->[$i]->id2)};              ++$genomes_hit{&FIG::genome_of($sims->[$i]->id2)};
150              $weight_of_hits{&FIG::genome_of($sims->[$i]->id2)} += $sims->[$i]->bsc;                  next unless $sims->[$i]->ln2;
151                    $weight_of_hits{&FIG::genome_of($sims->[$i]->id2)} +=
152                        ($sims->[$i]->bsc / $sims->[$i]->ln2);
153          }          }
154      }      }
155            else {
156                die "Could not promote $orf_id";
157            }
158        }
159    
160      last if (@peg_ids >= $N);      last if (@peg_ids >= $N);
161  }  }
162    
# Line 123  Line 171 
171  for ($i=0; ($i < $N) && ($i < @best); ++$i) {  for ($i=0; ($i < $N) && ($i < @best); ++$i) {
172      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";
173  }  }
174    
175    $to_call->export_features();
176    
177    close(CALLED_BY);
178  close(FOUNDFAMS);  close(FOUNDFAMS);
 $to_call->export_features;  
179    
180  exit(0);  exit(0);

Legend:
Removed from v.1.1  
changed lines
  Added in v.1.8

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3