[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.1 - (view) (download) (as text)

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3