[Bio] / FigKernelPackages / ParalogResolution.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/ParalogResolution.pm

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

revision 1.4, Sat Jun 27 18:55:29 2009 UTC revision 1.5, Wed Jul 1 00:08:52 2009 UTC
# Line 110  Line 110 
110      return join(" ",@pieces);      return join(" ",@pieces);
111  }  }
112    
113    
114  sub reference_set_for_paralogs {  sub reference_set_for_paralogs {
115      my($genomes,$roles,$parms) = @_;      my($genomes,$roles,$parms) = @_;
116    
117      $parms ||= {};      $parms ||= {};
118    
119      my $keep_multifunctional = $parms->{'keep_multifunctional'};      my $keep_multifunctional = $parms->{'keep_multifunctional'};
120      if (! defined($keep_multifunctional)) { $keep_multifunctional = 1 };      if (! defined($keep_multifunctional)) { $keep_multifunctional = 1 }
121    
122      my $max_sc               = $parms->{'max_sc'};      my $max_sc               = $parms->{'max_sc'};
123      if (! defined($max_sc))  { $max_sc = 1.0e-5 }      if (! defined($max_sc))  { $max_sc = 1.0e-5 }
# Line 136  Line 137 
137              if ($ref_genomes{&FIG::genome_of($peg)} && (! $fig->screwed_up($peg)))              if ($ref_genomes{&FIG::genome_of($peg)} && (! $fig->screwed_up($peg)))
138              {              {
139                  my $func = $fig->function_of($peg,"",1);                  my $func = $fig->function_of($peg,"",1);
140                    #  Original version could be confused by comments
141                    if ( 0 )
142                    {
143                  if ($ref_roles{$func} || ($keep_multifunctional && &role_in($func,$role)))                  if ($ref_roles{$func} || ($keep_multifunctional && &role_in($func,$role)))
144                  {                  {
145  #                   print STDERR "got $peg\n";  #                   print STDERR "got $peg\n";
146                      $pegs{$peg} = $func;                      $pegs{$peg} = $func;
147                  }                  }
148              }              }
149                    else
150                    {
151                        my $funcF = $func;
152                        $funcF =~ s/ ##? .*$//;
153                        $pegs{ $peg } = $func if ! ( $funcF =~ / \/ / );
154                    }
155                }
156          }          }
157      }      }
158      my @seqs = map { my $peg = $_;      my @seqs = map { my $peg = $_;
# Line 152  Line 163 
163                     keys(%pegs);                     keys(%pegs);
164    
165  #   Now add paralogs from given genomes  #   Now add paralogs from given genomes
166    #
167    #  Original code does not add paralogs in genomes that do not already have a
168    #  gene annotated with a requested role.
169    
170      my @to_add = ();      my @to_add = ();
171        if ( 0 )
172        {
173      foreach my $tuple (@seqs)      foreach my $tuple (@seqs)
174      {      {
175          my($peg,$comment,$seq) = @$tuple;          my($peg,$comment,$seq) = @$tuple;
# Line 174  Line 191 
191              }              }
192          }          }
193      }      }
194        }
195        else
196        {
197            foreach my $genome ( @$genomes )
198            {
199                #  Verify the blast db
200                my $fasta = "$FIG_Config::organisms/$genome/Features/peg/fasta";
201                if ( ! -s "$fasta.psq" )
202                {
203                    next if ! -s $fasta;
204                    system "$FIG_Config::ext_bin/formatdb -p T -i $fasta";
205                    next if ! -s "$fasta.psq";
206                }
207    
208                my $gs = $fig->genus_species( $genome );
209                my %seq_done;
210                foreach my $tuple ( @seqs )
211                {
212                    my( $peg, $comment, $seq ) = @$tuple;
213                    next if $seq_done{ $seq }++;
214                    foreach my $sim ( &FIG::blastitP( $peg, $seq, $fasta, $max_sc, "-FF" ) )
215                    {
216                        my $id2 = $sim->id2;
217                        if ( ( ! $pegs{ $id2 } )
218                          && $fig->is_real_feature( $id2 )
219                          && ( ( ( $sim->e2 + 1 - $sim->b2 ) / $sim->ln2 ) >= $min_cov )
220                           )
221                        {
222                            $pegs{ $id2 } = $fig->function_of( $id2 ) || 'undefined';
223                            push( @to_add, [ $id2, "$pegs{$id2} \[$gs\]", $fig->get_translation($id2) ] );
224                        }
225                    }
226                }
227            }
228        }
229      push(@seqs,@to_add);      push(@seqs,@to_add);
230      if (@seqs < 2) { return undef }  
231      return &gjoalignment::align_with_muscle(\@seqs);      ( @seqs > 1 ) ? &gjoalignment::align_with_muscle( \@seqs )  # align, or align & tree
232                      : wantarray ? () : undef;
233  }  }
234    
235    
236  sub role_in {  sub role_in {
237      my($func,$role) = @_;      my($func,$role) = @_;
238    

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3