[Bio] / FigKernelScripts / initialize_alignments_for_subsystems.pl Repository:
ViewVC logotype

Annotation of /FigKernelScripts/initialize_alignments_for_subsystems.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : overbeek 1.2 # -*- perl -*-
2 :     #
3 :     # Copyright (c) 2003-2006 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 : overbeek 1.1
20 :     use FIG;
21 :     my $fig = new FIG;
22 :    
23 :     # usage: initialize_alignments_for_subsystems FileOfSub [no argument means "all"]
24 :    
25 :     my @subsys;
26 :     if (@ARGV == 1)
27 :     {
28 :     @subsys = map { ($_ =~ /(\S.*\S)/) ? $1 : () } `cat $ARGV[0]`;
29 :     }
30 :     else
31 :     {
32 :     @subsys = $fig->all_subsystems;
33 :     }
34 :    
35 :     foreach my $subsys (@subsys)
36 :     {
37 :     print STDERR "processing $subsys\n";
38 :     my @roles = $fig->subsystem_to_roles($subsys);
39 :     for (my $i = 0; ($i < @roles); $i++)
40 :     {
41 :     print STDERR " $roles[$i]\n";
42 :     &align_column($fig,\@roles,$i+1,$subsys);
43 :     }
44 :     }
45 :    
46 :     sub align_column {
47 :     my($fig,$roles,$colN,$subsys) = @_;
48 :    
49 :     &check_index("$FIG_Config::data/Subsystems/$name/Alignments",$roles);
50 :     if (($role = &which_role_for_column($colN,$roles)) &&
51 :     ((@pegs = &seqs_to_align($fig,$role,$subsys)) > 1))
52 :     {
53 :     my $tmpF = "/tmp/seqs.fasta.$$";
54 :     open(TMP,">$tmpF") || die "could not open $tmpF";
55 :    
56 :     foreach $peg (@pegs)
57 :     {
58 :     if ($pseq = $fig->get_translation($peg))
59 :     {
60 :     $pseq =~ s/[uU]/x/g;
61 :     print TMP ">$peg\n$pseq\n";
62 :     }
63 :     }
64 :     close(TMP);
65 :    
66 :     my $dir = "$FIG_Config::data/Subsystems/$subsys/Alignments/$colN";
67 :    
68 :     if (-d $dir)
69 :     {
70 :     system "rm -rf \"$dir\"";
71 :     }
72 :    
73 :     &FIG::run("$FIG_Config::bin/split_and_trim_sequences \"$dir/split_info\" < $tmpF");
74 :    
75 :     if (-s "$dir/split_info/set.sizes")
76 :     {
77 :     open(SZ,"<$dir/split_info/set.sizes") || die " could not open $dir/split_info/set.sizes";
78 :     while (defined($_ = <SZ>))
79 :     {
80 :     if (($_ =~ /^(\d+)\t(\d+)/) && ($2 > 3))
81 :     {
82 :     my $n = $1;
83 :     &FIG::run("$FIG_Config::bin/make_phob_from_seqs \"$dir/$n\" < \"$dir/split_info\"/$n");
84 :     }
85 :     }
86 :     close(SZ);
87 :     &update_index("$FIG_Config::data/Subsystems/$subsys/Alignments/index",$colN,$role);
88 :     }
89 :     else
90 :     {
91 :     system("rm -rf \"$dir\"");
92 :     }
93 :     }
94 :     }
95 :    
96 :    
97 :    
98 :     sub check_index {
99 :     my($alignments,$roles) = @_;
100 :    
101 :     if (-s "$alignments/index")
102 :     {
103 :     my $ok = 1;
104 :     foreach $_ (`cat \"$alignments/index\"`)
105 :     {
106 :     $ok = $ok && (($_ =~ /^(\d+)\t(\S.*\S)/) && ($roles->[$1 - 1] eq $2));
107 :     }
108 :     if (! $ok)
109 :     {
110 :     system "rm -rf \"$alignments\"";
111 :     return 0;
112 :     }
113 :     return 1;
114 :     }
115 :     else
116 :     {
117 :     system "rm -rf \"$alignments\"";
118 :     }
119 :     return 0;
120 :     }
121 :    
122 :     sub update_index {
123 :     my($file,$colN,$role) = @_;
124 :    
125 :     my @lines = ();
126 :     if (-s $file)
127 :     {
128 :     @lines = grep { $_ !~ /^$colN\t/ } `cat $file`;
129 :     }
130 :     push(@lines,"$colN\t$role\n");
131 :     open(TMP,">$file") || die "could not open $file";
132 :     foreach $_ (@lines)
133 :     {
134 :     print TMP $_;
135 :     }
136 :     close(TMP);
137 :     }
138 :    
139 :     sub which_role_for_column {
140 :     my($col,$roles) = @_;
141 :     my($i);
142 :    
143 :     if (($col =~ /^(\d+)/) && ($1 <= @$roles))
144 :     {
145 :     return $roles->[$1-1];
146 :     }
147 :     return undef;
148 :     }
149 :    
150 :     sub seqs_to_align {
151 :     my($fig,$role,$subsys) = @_;
152 :     my($genome,$peg,%existing_function);
153 :    
154 :     my @seqs = ();
155 :     my @genomes = $fig->subsystem_genomes($subsys);
156 :     foreach $genome (map { $_->[0] } @{$fig->subsystem_genomes($subsys)})
157 :     {
158 :     foreach $peg ($fig->pegs_in_subsystem_cell($subsys,$genome,$role))
159 :     {
160 :     push(@seqs,$peg);
161 :     my $func = $fig->function_of($peg);
162 :     if (! &FIG::hypo($func))
163 :     {
164 :     my @roles_of_func = $fig->roles_of_function($func);
165 :     for (my $i=0; ($i < @roles_of_func) && ($roles_of_func[$i] ne $role); $i++) {}
166 :     if ($i < @roles_of_func)
167 :     {
168 :     $existing_function{$func} = 1;
169 :     }
170 :     }
171 :     }
172 :     }
173 :    
174 :     my $new = 0;
175 :     my @similar = &get_homologs($fig,\@seqs,1.0e-20);
176 :     while (($new < 100) && ($peg = shift @similar))
177 :     {
178 :     push(@seqs,$peg);
179 :     my $func = $fig->function_of($peg);
180 :     if (! $existing_function{$func})
181 :     {
182 :     $new++;
183 :     }
184 :     }
185 :     return @seqs;
186 :     }
187 :    
188 :     sub get_homologs {
189 :     my($fig,$pegs,$cutoff) = @_;
190 :     my($peg,$sim,$id2);
191 :    
192 :     my @homologs = ();
193 :     my %got = map { $_ => 1 } @$pegs;
194 :     my %new;
195 :    
196 :     foreach $peg (@$pegs)
197 :     {
198 :     foreach $sim ($fig->sims($peg,300,$cutoff,"fig"))
199 :     {
200 :     $id2 = $sim->id2;
201 :     if ((! $got{$id2}) && ((! $new{$id2}) || ($new{$id2} > $sim->psc)))
202 :     {
203 :     $new{$id2} = $sim->psc;
204 :     }
205 :     }
206 :     }
207 :     @homologs = sort { $new{$a} <=> $new{$b} } keys(%new);
208 :     return @homologs;
209 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3