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

Annotation of /FigKernelPackages/ParalogResolution.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.4 - (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 :    
113 :     sub reference_set_for_paralogs {
114 : overbeek 1.4 my($genomes,$roles,$parms) = @_;
115 :    
116 :     $parms ||= {};
117 :    
118 :     my $keep_multifunctional = $parms->{'keep_multifunctional'};
119 :     if (! defined($keep_multifunctional)) { $keep_multifunctional = 1 };
120 :    
121 :     my $max_sc = $parms->{'max_sc'};
122 :     if (! defined($max_sc)) { $max_sc = 1.0e-5 }
123 :    
124 :     my $min_cov = $parms->{'min_cov'};
125 :     if (! defined($min_cov)) { $min_cov = 0.7 }
126 : overbeek 1.1
127 :     if (! $max_sc) { $max_sc = 1.0e-5 }
128 :     my %ref_genomes = map { $_ => 1 } @$genomes;
129 :     my %ref_roles = map { $_ => 1 } @$roles;
130 :     my $fig = new FIG;
131 :     my %pegs;
132 :     foreach my $role (@$roles)
133 :     {
134 :     foreach my $peg ($fig->prots_for_role($role))
135 :     {
136 :     if ($ref_genomes{&FIG::genome_of($peg)} && (! $fig->screwed_up($peg)))
137 :     {
138 :     my $func = $fig->function_of($peg,"",1);
139 :     if ($ref_roles{$func} || ($keep_multifunctional && &role_in($func,$role)))
140 :     {
141 :     # print STDERR "got $peg\n";
142 :     $pegs{$peg} = $func;
143 :     }
144 :     }
145 :     }
146 :     }
147 :     my @seqs = map { my $peg = $_;
148 :     my $gs = $fig->genus_species(&FIG::genome_of($peg));
149 :     my $func_of_peg = $pegs{$peg};
150 :     [$peg,"$func_of_peg \[$gs\]",$fig->get_translation($peg)]
151 :     }
152 :     keys(%pegs);
153 :    
154 :     # Now add paralogs from given genomes
155 :     my @to_add = ();
156 :     foreach my $tuple (@seqs)
157 :     {
158 :     my($peg,$comment,$seq) = @$tuple;
159 :     my $genome = &FIG::genome_of($peg);
160 :     if (! -s "$FIG_Config::organisms/$genome/Features/peg/fasta.psq")
161 :     {
162 :     system "$FIG_Config::ext_bin/formatdb -p T -i $FIG_Config::organisms/$genome/Features/peg/fasta";
163 :     }
164 :     my @blastout = &FIG::blastitP($peg,$seq,"$FIG_Config::organisms/$genome/Features/peg/fasta",$max_sc,"-FF");
165 :    
166 :     foreach my $sim (@blastout)
167 :     {
168 : overbeek 1.4 if ((! $pegs{$sim->id2}) && $fig->is_real_feature($sim->id2) &&
169 :     ((($sim->e2 + 1 - $sim->b2) /$sim->ln2) >= $min_cov))
170 : overbeek 1.1 {
171 :     $pegs{$sim->id2} = $fig->function_of($sim->id2);
172 :     my $gs = $fig->genus_species(&FIG::genome_of($sim->id2));
173 :     push(@to_add,[$sim->id2,"$pegs{$sim->id2} \[$gs\]",$fig->get_translation($sim->id2)]);
174 :     }
175 :     }
176 :     }
177 :     push(@seqs,@to_add);
178 :     if (@seqs < 2) { return undef }
179 : overbeek 1.2 return &gjoalignment::align_with_muscle(\@seqs);
180 : overbeek 1.1 }
181 :    
182 :     sub role_in {
183 :     my($func,$role) = @_;
184 :    
185 :     my @roles_of_function = &FIG::roles_of_function($func);
186 :     my $i;
187 :     for ($i=0; ($i < @roles_of_function) && ($role ne $roles_of_function[$i]); $i++) {}
188 :     return ($i < @roles_of_function);
189 :     }
190 :    
191 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3