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

Annotation of /FigKernelPackages/ParalogResolution.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (view) (download) (as text)

1 : overbeek 1.1 # -*- perl -*-
2 :     ########################################################################
3 :     # Copyright (c) 2003-2009 University of Chicago and Fellowship
4 :     # for Interpretations of Genomes. All Rights Reserved.
5 :     #
6 :     # This file is part of the SEED Toolkit.
7 :     #
8 :     # The SEED Toolkit is free software. You can redistribute
9 :     # it and/or modify it under the terms of the SEED Toolkit
10 :     # Public License.
11 :     #
12 :     # You should have received a copy of the SEED Toolkit Public License
13 :     # along with this program; if not write to the University of Chicago
14 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
15 :     # Genomes at veronika@thefig.info or download a copy from
16 :     # http://www.theseed.org/LICENSE.TXT.
17 :     ########################################################################
18 :    
19 :     package ParalogResolution;
20 :    
21 :     use strict;
22 :     use FIG;
23 :    
24 :     use Data::Dumper;
25 :     use Carp;
26 : overbeek 1.2 use gjoalignment;
27 : overbeek 1.1
28 :     sub context_tags {
29 :     my($pegs,$dist) = @_;
30 :    
31 :     my $fig = new FIG;
32 :     my($i,%before,%after,%funcs);
33 :     foreach my $peg (@$pegs)
34 :     {
35 :     my $loc = $fig->feature_location($peg);
36 :     if ($loc)
37 :     {
38 :     my($contig,$beg,$end) = $fig->boundaries_of($loc);
39 :     if ($contig && $beg && $end)
40 :     {
41 :     my $min = &FIG::min($beg,$end) - $dist;
42 :     my $max = &FIG::max($beg,$end) + $dist;
43 :     my $feat;
44 :     ($feat,undef,undef) = $fig->genes_in_region(&FIG::genome_of($peg),$contig,$min,$max);
45 :     my(@before,@after,$seen_peg);
46 :     for ($i=0; ($i < @$feat); $i++)
47 :     {
48 :     my $fid = $feat->[$i];
49 :     if (&FIG::ftype($fid) eq 'peg')
50 :     {
51 :     if ($fid eq $peg)
52 :     {
53 :     $seen_peg = 1;
54 :     next;
55 :     }
56 :     my $func = $fig->function_of($fid,"",1);
57 :     $funcs{$func}++;
58 :    
59 :     if ((($beg < $end) && (! $seen_peg)) ||
60 :     (($beg > $end) && $seen_peg))
61 :     {
62 :     push(@before,[$fid,$func]);
63 :     }
64 :     else
65 :     {
66 :     push(@after,[$fid,$func]);
67 :     }
68 :     }
69 :     }
70 :     $before{$peg} = [($beg < $end) ? @before : reverse @before];
71 :     $after{$peg} = [($beg < $end) ? @after : reverse @after];
72 :     }
73 :     }
74 :     }
75 :     my @hits = sort { $funcs{$b} <=> $funcs{$a} }
76 :     grep { $funcs{$_} > 1 }
77 :     keys(%funcs);
78 :    
79 :     my %val;
80 :     my %val2func;
81 :     for ($i=0; ($i < @hits); $i++)
82 :     {
83 :     $val{$hits[$i]} = $i+1;
84 :     $val2func{$i+1} = $hits[$i];
85 :     }
86 :     my $tags = {};
87 :     foreach my $peg (@$pegs)
88 :     {
89 :     my $left = &process_neigh($peg,\%before,\%val);
90 :     my $right = &process_neigh($peg,\%after,\%val);
91 :     $tags->{$peg} = $left . "::" . $right;
92 :     }
93 :     return ($tags,[map { [$_,$val2func{$_}] } sort { $a <=> $b } keys(%val2func)]);
94 :     }
95 :    
96 :     sub process_neigh {
97 :     my($peg,$neigh,$val) = @_;
98 :    
99 :     my @pieces = ();
100 :     my $tuples = $neigh->{$peg};
101 :     foreach my $tuple (@$tuples)
102 :     {
103 :     my($peg1,$func1) = @$tuple;
104 :     my $x = $val->{$func1};
105 :     if ($x)
106 :     {
107 :     push(@pieces,$x);
108 :     }
109 :     }
110 :     return join(" ",@pieces);
111 :     }
112 : golsen 1.5
113 :    
114 : overbeek 1.1 sub reference_set_for_paralogs {
115 : golsen 1.5 my( $genomes, $roles, $parms ) = @_;
116 : overbeek 1.4
117 :     $parms ||= {};
118 :    
119 :     my $keep_multifunctional = $parms->{'keep_multifunctional'};
120 : golsen 1.5 if (! defined($keep_multifunctional)) { $keep_multifunctional = 1 }
121 : overbeek 1.4
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 : overbeek 1.1
128 :     if (! $max_sc) { $max_sc = 1.0e-5 }
129 :     my %ref_genomes = map { $_ => 1 } @$genomes;
130 :     my %ref_roles = map { $_ => 1 } @$roles;
131 :     my $fig = new FIG;
132 :     my %pegs;
133 : golsen 1.5 foreach my $role ( @$roles )
134 : overbeek 1.1 {
135 :     foreach my $peg ($fig->prots_for_role($role))
136 :     {
137 : golsen 1.5 if ( $ref_genomes{ &FIG::genome_of( $peg ) } && ( ! $fig->screwed_up( $peg ) ) )
138 : overbeek 1.1 {
139 :     my $func = $fig->function_of($peg,"",1);
140 : golsen 1.5 # Original version could be confused by comments
141 :     if ( 0 )
142 : overbeek 1.1 {
143 : golsen 1.5 if ($ref_roles{$func} || ($keep_multifunctional && &role_in($func,$role)))
144 :     {
145 :     # print STDERR "got $peg\n";
146 :     $pegs{$peg} = $func;
147 :     }
148 : overbeek 1.1 }
149 : golsen 1.5 else
150 :     {
151 :     my $funcF = $func;
152 :     $funcF =~ s/ ##? .*$//;
153 :     $pegs{ $peg } = $func if ! ( $funcF =~ / \/ / );
154 :     }
155 :     }
156 :     }
157 : overbeek 1.1 }
158 :     my @seqs = map { my $peg = $_;
159 :     my $gs = $fig->genus_species(&FIG::genome_of($peg));
160 :     my $func_of_peg = $pegs{$peg};
161 :     [$peg,"$func_of_peg \[$gs\]",$fig->get_translation($peg)]
162 :     }
163 :     keys(%pegs);
164 :    
165 :     # Now add paralogs from given genomes
166 : golsen 1.5 #
167 :     # Original code does not add paralogs in genomes that do not already have a
168 :     # gene annotated with a requested role.
169 :    
170 : overbeek 1.1 my @to_add = ();
171 : golsen 1.5 if ( 0 )
172 :     {
173 :     foreach my $tuple (@seqs)
174 :     {
175 :     my($peg,$comment,$seq) = @$tuple;
176 :     my $genome = &FIG::genome_of($peg);
177 :     if (! -s "$FIG_Config::organisms/$genome/Features/peg/fasta.psq")
178 :     {
179 :     system "$FIG_Config::ext_bin/formatdb -p T -i $FIG_Config::organisms/$genome/Features/peg/fasta";
180 :     }
181 :     my @blastout = &FIG::blastitP($peg,$seq,"$FIG_Config::organisms/$genome/Features/peg/fasta",$max_sc,"-FF");
182 :    
183 :     foreach my $sim (@blastout)
184 :     {
185 :     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);
189 :     my $gs = $fig->genus_species(&FIG::genome_of($sim->id2));
190 :     push(@to_add,[$sim->id2,"$pegs{$sim->id2} \[$gs\]",$fig->get_translation($sim->id2)]);
191 :     }
192 :     }
193 :     }
194 :     }
195 :     else
196 : overbeek 1.1 {
197 : golsen 1.5 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 );
230 : overbeek 1.1
231 : golsen 1.5 ( @seqs > 1 ) ? &gjoalignment::align_with_muscle( \@seqs ) # align, or align & tree
232 :     : wantarray ? () : undef;
233 : overbeek 1.1 }
234 :    
235 : golsen 1.5
236 : overbeek 1.1 sub role_in {
237 :     my($func,$role) = @_;
238 :    
239 :     my @roles_of_function = &FIG::roles_of_function($func);
240 :     my $i;
241 :     for ($i=0; ($i < @roles_of_function) && ($role ne $roles_of_function[$i]); $i++) {}
242 :     return ($i < @roles_of_function);
243 :     }
244 :    
245 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3