[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.1, Mon Jun 22 19:35:59 2009 UTC revision 1.5, Wed Jul 1 00:08:52 2009 UTC
# Line 23  Line 23 
23    
24  use Data::Dumper;  use Data::Dumper;
25  use Carp;  use Carp;
26    use gjoalignment;
27    
28  sub context_tags {  sub context_tags {
29      my($pegs,$dist) = @_;      my($pegs,$dist) = @_;
# Line 109  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,$keep_multifunctional,$max_sc) = @_;      my( $genomes, $roles, $parms ) = @_;
116    
117        $parms ||= {};
118    
119        my $keep_multifunctional = $parms->{'keep_multifunctional'};
120        if (! defined($keep_multifunctional)) { $keep_multifunctional = 1 }
121    
122        my $max_sc               = $parms->{'max_sc'};
123        if (! defined($max_sc))  { $max_sc = 1.0e-5 }
124    
125        my $min_cov              = $parms->{'min_cov'};
126        if (! defined($min_cov)) { $min_cov = 0.7 }
127    
128      if (! $max_sc) { $max_sc = 1.0e-5 }      if (! $max_sc) { $max_sc = 1.0e-5 }
129      my %ref_genomes = map { $_ => 1 } @$genomes;      my %ref_genomes = map { $_ => 1 } @$genomes;
# Line 124  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 140  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 153  Line 182 
182    
183          foreach my $sim (@blastout)          foreach my $sim (@blastout)
184          {          {
185              if (! $pegs{$sim->id2})                  if ((! $pegs{$sim->id2}) && $fig->is_real_feature($sim->id2) &&
186                        ((($sim->e2 + 1 - $sim->b2) /$sim->ln2) >= $min_cov))
187              {              {
188                  $pegs{$sim->id2} = $fig->function_of($sim->id2);                  $pegs{$sim->id2} = $fig->function_of($sim->id2);
189                  my $gs = $fig->genus_species(&FIG::genome_of($sim->id2));                  my $gs = $fig->genus_species(&FIG::genome_of($sim->id2));
# Line 161  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 &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.1  
changed lines
  Added in v.1.5

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3