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

Annotation of /FigKernelPackages/ParalogResolution.pm

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (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 :    
27 :     sub context_tags {
28 :     my($pegs,$dist) = @_;
29 :    
30 :     my $fig = new FIG;
31 :     my($i,%before,%after,%funcs);
32 :     foreach my $peg (@$pegs)
33 :     {
34 :     my $loc = $fig->feature_location($peg);
35 :     if ($loc)
36 :     {
37 :     my($contig,$beg,$end) = $fig->boundaries_of($loc);
38 :     if ($contig && $beg && $end)
39 :     {
40 :     my $min = &FIG::min($beg,$end) - $dist;
41 :     my $max = &FIG::max($beg,$end) + $dist;
42 :     my $feat;
43 :     ($feat,undef,undef) = $fig->genes_in_region(&FIG::genome_of($peg),$contig,$min,$max);
44 :     my(@before,@after,$seen_peg);
45 :     for ($i=0; ($i < @$feat); $i++)
46 :     {
47 :     my $fid = $feat->[$i];
48 :     if (&FIG::ftype($fid) eq 'peg')
49 :     {
50 :     if ($fid eq $peg)
51 :     {
52 :     $seen_peg = 1;
53 :     next;
54 :     }
55 :     my $func = $fig->function_of($fid,"",1);
56 :     $funcs{$func}++;
57 :    
58 :     if ((($beg < $end) && (! $seen_peg)) ||
59 :     (($beg > $end) && $seen_peg))
60 :     {
61 :     push(@before,[$fid,$func]);
62 :     }
63 :     else
64 :     {
65 :     push(@after,[$fid,$func]);
66 :     }
67 :     }
68 :     }
69 :     $before{$peg} = [($beg < $end) ? @before : reverse @before];
70 :     $after{$peg} = [($beg < $end) ? @after : reverse @after];
71 :     }
72 :     }
73 :     }
74 :     my @hits = sort { $funcs{$b} <=> $funcs{$a} }
75 :     grep { $funcs{$_} > 1 }
76 :     keys(%funcs);
77 :    
78 :     my %val;
79 :     my %val2func;
80 :     for ($i=0; ($i < @hits); $i++)
81 :     {
82 :     $val{$hits[$i]} = $i+1;
83 :     $val2func{$i+1} = $hits[$i];
84 :     }
85 :     my $tags = {};
86 :     foreach my $peg (@$pegs)
87 :     {
88 :     my $left = &process_neigh($peg,\%before,\%val);
89 :     my $right = &process_neigh($peg,\%after,\%val);
90 :     $tags->{$peg} = $left . "::" . $right;
91 :     }
92 :     return ($tags,[map { [$_,$val2func{$_}] } sort { $a <=> $b } keys(%val2func)]);
93 :     }
94 :    
95 :     sub process_neigh {
96 :     my($peg,$neigh,$val) = @_;
97 :    
98 :     my @pieces = ();
99 :     my $tuples = $neigh->{$peg};
100 :     foreach my $tuple (@$tuples)
101 :     {
102 :     my($peg1,$func1) = @$tuple;
103 :     my $x = $val->{$func1};
104 :     if ($x)
105 :     {
106 :     push(@pieces,$x);
107 :     }
108 :     }
109 :     return join(" ",@pieces);
110 :     }
111 :    
112 :     sub reference_set_for_paralogs {
113 :     my($genomes,$roles,$keep_multifunctional,$max_sc) = @_;
114 :    
115 :     if (! $max_sc) { $max_sc = 1.0e-5 }
116 :     my %ref_genomes = map { $_ => 1 } @$genomes;
117 :     my %ref_roles = map { $_ => 1 } @$roles;
118 :     my $fig = new FIG;
119 :     my %pegs;
120 :     foreach my $role (@$roles)
121 :     {
122 :     foreach my $peg ($fig->prots_for_role($role))
123 :     {
124 :     if ($ref_genomes{&FIG::genome_of($peg)} && (! $fig->screwed_up($peg)))
125 :     {
126 :     my $func = $fig->function_of($peg,"",1);
127 :     if ($ref_roles{$func} || ($keep_multifunctional && &role_in($func,$role)))
128 :     {
129 :     # print STDERR "got $peg\n";
130 :     $pegs{$peg} = $func;
131 :     }
132 :     }
133 :     }
134 :     }
135 :     my @seqs = map { my $peg = $_;
136 :     my $gs = $fig->genus_species(&FIG::genome_of($peg));
137 :     my $func_of_peg = $pegs{$peg};
138 :     [$peg,"$func_of_peg \[$gs\]",$fig->get_translation($peg)]
139 :     }
140 :     keys(%pegs);
141 :    
142 :     # Now add paralogs from given genomes
143 :     my @to_add = ();
144 :     foreach my $tuple (@seqs)
145 :     {
146 :     my($peg,$comment,$seq) = @$tuple;
147 :     my $genome = &FIG::genome_of($peg);
148 :     if (! -s "$FIG_Config::organisms/$genome/Features/peg/fasta.psq")
149 :     {
150 :     system "$FIG_Config::ext_bin/formatdb -p T -i $FIG_Config::organisms/$genome/Features/peg/fasta";
151 :     }
152 :     my @blastout = &FIG::blastitP($peg,$seq,"$FIG_Config::organisms/$genome/Features/peg/fasta",$max_sc,"-FF");
153 :    
154 :     foreach my $sim (@blastout)
155 :     {
156 :     if (! $pegs{$sim->id2})
157 :     {
158 :     $pegs{$sim->id2} = $fig->function_of($sim->id2);
159 :     my $gs = $fig->genus_species(&FIG::genome_of($sim->id2));
160 :     push(@to_add,[$sim->id2,"$pegs{$sim->id2} \[$gs\]",$fig->get_translation($sim->id2)]);
161 :     }
162 :     }
163 :     }
164 :     push(@seqs,@to_add);
165 :     if (@seqs < 2) { return undef }
166 :     return &align_with_muscle(\@seqs);
167 :     }
168 :    
169 :     sub role_in {
170 :     my($func,$role) = @_;
171 :    
172 :     my @roles_of_function = &FIG::roles_of_function($func);
173 :     my $i;
174 :     for ($i=0; ($i < @roles_of_function) && ($role ne $roles_of_function[$i]); $i++) {}
175 :     return ($i < @roles_of_function);
176 :     }
177 :    
178 :     1;

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3