[Bio] / FigWebServices / diffsF.cgi Repository:
ViewVC logotype

Annotation of /FigWebServices/diffsF.cgi

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (view) (download)

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 use GenoGraphics;
20 :    
21 :     use FIG;
22 :     my $fig = new FIG;
23 :    
24 : overbeek 1.4 use FigFam;
25 :     my $default_dir = "$FIG_Config::data/FigfamsData";
26 :    
27 : overbeek 1.1 use HTML;
28 :    
29 :     use CGI;
30 :     my $cgi = new CGI;
31 :    
32 :     if (0)
33 :     {
34 :     my $VAR1;
35 :     eval(join("",`cat /tmp/diffsF_parms`));
36 :     $cgi = $VAR1;
37 :     # print STDERR &Dumper($cgi);
38 :     }
39 :    
40 :     if (0)
41 :     {
42 :     print $cgi->header;
43 :     my @params = $cgi->param;
44 :     print "<pre>\n";
45 :     foreach $_ (@params)
46 :     {
47 :     print "$_\t:",join(",",$cgi->param($_)),":\n";
48 :     }
49 :    
50 :     if (0)
51 :     {
52 :     if (open(TMP,">/tmp/diffsF_parms"))
53 :     {
54 :     print TMP &Dumper($cgi);
55 :     close(TMP);
56 :     }
57 :     }
58 :     exit;
59 :     }
60 :    
61 :     use constant FID => 0;
62 :     use constant LOCUS => 1;
63 :     use constant CONTIG => 2;
64 :     use constant START => 3;
65 :     use constant STOP => 4;
66 :     use constant LEN => 5;
67 :     use constant STRAND => 6;
68 :     use constant TYPE => 7;
69 :     use constant TAXID => 8;
70 :    
71 :     my $html = [];
72 :     unshift @$html, "<TITLE>Feature Differences</TITLE>\n";
73 : overbeek 1.4 $org1 = $cgi->param('org1');
74 :     $org2 = $cgi->param('org2');
75 :     $gene = $cgi->param('gene');
76 :     $align1 = $cgi->param('align1');
77 :     $align2 = $cgi->param('align2');
78 : overbeek 1.1
79 :     if (!((-d $org1) || (-d $org2)))
80 :     {
81 : overbeek 1.4 push(@$html,$cgi->h1("bad parms: gene=$gene align1=$align1 align2=$align2 org1=$org1 org2=$org2"));
82 : overbeek 1.1 }
83 :     elsif ($gene = $cgi->param('gene'))
84 :     {
85 :     &show_region($fig,$cgi,$html,$org1,$gene);
86 :     }
87 : overbeek 1.4 elsif (($gene1 = $cgi->param('align1')) && ($gene2 = $cgi->param('align2')))
88 :     {
89 :     &show_alignment($fig, $cgi, $html, $gene1, $gene2, $org1, $org2);
90 :     }
91 : overbeek 1.1 else
92 :     {
93 : overbeek 1.4 &list_diffs($fig, $cgi, $html, $org1, $org2);
94 : overbeek 1.1 }
95 :     &HTML::show_page($cgi,$html);
96 :    
97 :     sub show_region {
98 :     my($fig,$cgi,$html,$org,$gene) = @_;
99 :     my @feat = &load_tbl2($org,$gene);
100 :     for ($i=0; ($i < @feat) && ($feat[$i]->[FID] ne $gene); $i++) {}
101 :     if ($i < @feat)
102 :     {
103 :     $first = ($i > 4) ? $i-5 : 0;
104 :     $last = ($i < (@feat - 5)) ? $i+5 : @feat - 1;
105 :     my $gg = [];
106 :     my $genes = [];
107 :     $min = 10000000;
108 :     $max = 0;
109 :     while ($first <= $last)
110 :     {
111 :     $f = $feat[$first];
112 :     $beg = $f->[START];
113 :     $end = $f->[STOP];
114 :     $min = &FIG::min($min,$beg);
115 :     $min = &FIG::min($min,$end);
116 :     $max = &FIG::max($max,$beg);
117 :     $max = &FIG::max($max,$end);
118 :     $fid = $f->[FID];
119 : overbeek 1.4 # $fid =~ /\.([a-z]+\.\d+)$/;
120 :     # $info = $1;
121 : overbeek 1.1 push(@$genes,[ &FIG::min($beg,$end),
122 :     &FIG::max($beg,$end),
123 :     ($beg < $end) ? "rightArrow" : "leftArrow",
124 :     ($fid eq $gene) ? "red" : "blue",
125 :     "",
126 :     "",
127 :     "$fid: $beg to $end",
128 :     ""
129 :     ] );
130 :     $first++;
131 :     }
132 :     my $map = ['Region',0,$max+1-$min,&decr_coords($genes,$min)];
133 :     push(@$gg,$map);
134 :     push( @$html, @{ &GenoGraphics::render( $gg, 1000, 4, 1 ) } );
135 :     }
136 :     }
137 :    
138 :     sub decr_coords {
139 :     my($genes,$min) = @_;
140 :     my($gene);
141 :    
142 :     foreach $gene (@$genes)
143 :     {
144 :     $gene->[0] -= $min;
145 :     $gene->[1] -= $min;
146 :     }
147 :     return $genes;
148 :     }
149 :    
150 :     sub list_diffs {
151 : overbeek 1.3 my ($fig, $cgi, $html, $org1, $org2) = @_;
152 :    
153 : overbeek 1.1 $org1T = &load_tbl1($org1);
154 :     $org2T = &load_tbl1($org2);
155 : overbeek 1.3
156 :     my @both;
157 :     my @just1;
158 :     foreach $key (sort { by_locus($org1T->{$a}, $org1T->{$b}) } keys %$org1T)
159 : overbeek 1.1 {
160 :     $type1 = $org1T->{$key}->[TYPE];
161 :     $locus1 = $org1T->{$key}->[LOCUS];
162 :    
163 :     if ($org2T->{$key})
164 :     {
165 :     $type2 = $org2T->{$key}->[TYPE];
166 :     $locus2 = $org2T->{$key}->[LOCUS];
167 :     if (($type1 ne $type2) || ($locus1 ne $locus2))
168 :     {
169 :     $fid1 = $org1T->{$key}->[FID];
170 :     $fid2 = $org2T->{$key}->[FID];
171 :     $link1 = &link($org1,$fid1);
172 :     $link2 = &link($org2,$fid2);
173 : overbeek 1.2
174 :     if (($len1 = $org1T->{$key}->[LEN]) < ($len2 = $org2T->{$key}->[LEN])) {
175 :     $changed = qq(lengthened);
176 :     }
177 :     else {
178 :     $changed = qq(shortened);
179 :     }
180 : overbeek 1.3
181 : overbeek 1.6 push @both, [ $changed, $fid1, $link1, $len1, $fid2, $link2, $len2, ($len2 - $len1) ];
182 : overbeek 1.1 }
183 :     }
184 :     else
185 :     {
186 :     $fid1 = $org1T->{$key}->[FID];
187 :     $link1 = &link($org1,$fid1);
188 : overbeek 1.6 push @just1, [ $fid1, $link1, $org1T->{$key}->[LEN] ];
189 : overbeek 1.1 }
190 :     }
191 : overbeek 1.3
192 :     my @just2;
193 : overbeek 1.1 foreach $key (sort { by_locus($org2T->{$a},$org2T->{$b}) } keys %$org2T)
194 :     {
195 :     if (! $org1T->{$key})
196 :     {
197 :     $type2 = $org2T->{$key}->[TYPE];
198 :     $locus2 = $org2T->{$key}->[LOCUS];
199 :     $fid2 = $org2T->{$key}->[FID];
200 :     $link2 = &link($org2,$fid2);
201 : overbeek 1.6 push @just2, [ $fid2, $link2, $org2T->{$key}->[LEN] ];
202 : overbeek 1.1 }
203 :     }
204 : overbeek 1.3
205 : overbeek 1.4 push @$html, map {
206 : overbeek 1.6 $align = qq(<A HREF="diffsF.cgi?align1=$_->[1]&align2=$_->[4]&org1=$org1&org2=$org2">align</A>);
207 :     $cgi->h2("$align $_->[0]: $_->[2] ($_->[3] bp) => $_->[5] ($_->[6] bp), diff=$_->[7]")
208 : overbeek 1.4 }
209 :     sort {
210 : overbeek 1.6 abs($b->[7]) <=> abs($a->[7])
211 : overbeek 1.4 } @both;
212 :    
213 :    
214 :     push @$html, map {
215 : overbeek 1.6 $cgi->h2("just2: $_->[1] ($_->[2] bp)")
216 : overbeek 1.4 }
217 :     sort {
218 : overbeek 1.6 $b->[2] <=> $a->[2]
219 : overbeek 1.4 } @just2;
220 :    
221 :    
222 :     push @$html, map {
223 : overbeek 1.6 $cgi->h2("just1: $_->[1] ($_->[2] bp)")
224 : overbeek 1.4 }
225 :     sort {
226 : overbeek 1.6 $b->[2] <=> $a->[2]
227 : overbeek 1.4 } @just1;
228 :     }
229 :    
230 :     sub show_alignment {
231 :     my ($fig, $cgi, $html, $gene1, $gene2, $org1, $org2) = @_;
232 :    
233 :     my %neighbors = map { m/^(\d+\.\d+)/o ? ($1 => 1) : () } `cat $org2/neighbors`;
234 :    
235 :     my (@tmp, $tmp);
236 :     if (@tmp = grep { m/^(\S+)/o && ($1 eq $gene2) } `cat $org2/found`) {
237 :     if ($tmp[0] =~ m/^\S+\t(FIG\d{6})/o) {
238 :     $figfam = new FigFam($fig, $1, $default_dir);
239 : overbeek 1.6 my $fam_id = $figfam->family_id();
240 :     push @$html, $cgi->h2("$gene2 is in Family $fam_id");
241 :     # my @members = grep { $neighbors{$fig->genome_of($_)} } $figfam->list_members();
242 :     my @members = $figfam->list_members();
243 :     @members || push @$html, $cgi->h2("No members in $fam_id");
244 :    
245 :     my $tmp_seqs = "$FIG_Config::temp/tmp_seqs_$$";
246 :     my $tmp_fid = "$FIG_Config::temp/tmp_fid_$$";
247 :     open(TMP_SEQS, ">$tmp_seqs")
248 : overbeek 1.4 || ((push @$html, $cgi->h2("Could not write-open $tmp_seqs")) && return);
249 :    
250 :     foreach my $id (@members) {
251 :     my $seq = $fig->get_translation($id);
252 : overbeek 1.6 if ($seq) {
253 :     # push @$html, $cgi->h2("Got seq for $id");
254 :     &FIG::display_id_and_seq($id, \$seq, \*TMP_SEQS);
255 :     }
256 :     else {
257 :     push @$html, $cgi->h2("No seq for $id");
258 :     }
259 : overbeek 1.4 }
260 : overbeek 1.6
261 :     open(TMP_FID, ">$tmp_fid")
262 :     || ((push @$html, $cgi->h2("Could not write-open $tmp_fid")) && return);
263 :     print TMP_FID "$gene1\n";
264 :     $tmp = `$FIG_Config::bin/pull_fasta_entries $org1/Features/peg/fasta < $tmp_fid`;
265 :     $tmp =~ s/^>fig/>old/;
266 :     print TMP_SEQS $tmp;
267 :    
268 :     open(TMP_FID, ">$tmp_fid")
269 :     || ((push @$html, $cgi->h2("Could not write-open $tmp_fid")) && return);
270 :     print TMP_FID "$gene2\n";
271 :     $tmp = `$FIG_Config::bin/pull_fasta_entries $org2/Features/peg/fasta < $tmp_fid`;
272 :     $tmp =~ s/^>fig/>new/;
273 :     print TMP_SEQS $tmp;
274 :    
275 :     close(TMP_SEQS)
276 : overbeek 1.4 || ((push @$html, $cgi->h2("Could not close $tmp_seqs")) && return);
277 :    
278 : overbeek 1.6 system("$FIG_Config::ext_bin/clustalw -infile=$tmp_seqs -align -outorder=aligned > /dev/null")
279 : overbeek 1.4 && ((push @$html, $cgi->h2("clustalw failed on file $tmp_seqs")) && return);
280 :    
281 :     push @$html, qq(<pre>\n);
282 : overbeek 1.6 @tmp = `cat $tmp_seqs.aln`;
283 :     push @$html, (@tmp ? (@tmp) : ("No lines read from $tmp_seqs.aln"));
284 : overbeek 1.4 push @$html, qq(</pre>\n);
285 : overbeek 1.6
286 :     system("rm -f $tmp_seqs* $tmp_fid")
287 :     && ((push @$html, $cgi->h2("Could not remove temporary files")) && return);
288 : overbeek 1.4 }
289 :     else {
290 :     $tmp[0] =~ s/\t/ /go;
291 :     push @$html, $cgi->h2("Could not parse family from entry: $tmp[0]");
292 :     }
293 :     }
294 :     else {
295 :     push @$html, $cgi->h2("No family found for gene2=$gene2");
296 :     }
297 : overbeek 1.1 }
298 :    
299 :     sub load_tbl1
300 :     {
301 :     my ($dir) = @_;
302 :     my ($entry, $id, $locus, $contig, $beg, $end, $len, $strand, $taxid, $type);
303 :    
304 :     open(TBL, "cat $dir/Features/*/tbl |") || die "Could not open $dir";
305 :    
306 :     my $tbl = {};
307 :     while (defined($entry = <TBL>))
308 :     {
309 :     chomp $entry;
310 :    
311 :     ($id, $locus) = split /\t/, $entry;
312 :     $id =~ m/^[^\|]+\|(\d+\.\d+)\.([^\.]+)/;
313 :     ($taxid, $type) = ($1, $2);
314 :    
315 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
316 :     && defined($contig) && $contig
317 :     && defined($beg) && $beg
318 :     && defined($end) && $end
319 :     && defined($len) && $len
320 :     && defined($strand) && $strand
321 :     )
322 :     {
323 :     $tbl->{"$contig\t$strand$end"} = [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $taxid ];
324 :     }
325 :     else
326 :     {
327 :     warn "INVALID ENTRY:\t$entry\n";
328 :     }
329 :     }
330 :    
331 :     return $tbl;
332 :     }
333 :    
334 :     sub load_tbl2
335 :     {
336 :     my ($dir,$gene) = @_;
337 :    
338 : overbeek 1.3 my ($entry, $id, $locus, $contig, $beg, $end, $len, $strand, $taxid, $type, $keep);
339 : overbeek 1.1
340 :     open(TBL, "cat $dir/Features/*/tbl |") || die "Could not open $dir";
341 :    
342 :     my $tbl = {};
343 :     while (defined($entry = <TBL>))
344 :     {
345 :     chomp $entry;
346 :    
347 :     ($id, $locus) = split /\t/, $entry;
348 :     $id =~ m/^[^\|]+\|(\d+\.\d+)\.([^\.]+)/;
349 :     ($taxid, $type) = ($1, $2);
350 :    
351 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
352 :     && defined($contig) && $contig
353 :     && defined($beg) && $beg
354 :     && defined($end) && $end
355 :     && defined($len) && $len
356 :     && defined($strand) && $strand
357 :     )
358 :     {
359 : overbeek 1.3 push(@{$tbl->{$contig}}, [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $taxid ]);
360 : overbeek 1.1 if ($id eq $gene) { $keep = $contig }
361 :     }
362 :     else
363 :     {
364 :     warn "INVALID ENTRY:\t$entry\n";
365 :     }
366 :     }
367 : overbeek 1.5 $keep || die "no contig for gene=$gene";
368 : overbeek 1.1 return sort { by_locus($a,$b) } @{$tbl->{$keep}};
369 :     }
370 :    
371 :     sub from_locus
372 :     {
373 :     my ($locus) = @_;
374 :     my ($contig, $beg, $end);
375 :    
376 :     if ( ($locus =~ m/^([^,]+)_(\d+)_\d+/) && ($contig = $1) && ($beg = $2)
377 :     && ($locus =~ m/[^,]+_\d+_(\d+)$/) && ($end = $1)
378 :     )
379 :     {
380 :     return ($contig, $beg, $end, (1+abs($end-$beg)), (($beg < $end) ? '+' : '-'));
381 :     }
382 :     else
383 :     {
384 :     die "Invalid locus $locus";
385 :     }
386 :    
387 :     return ();
388 :     }
389 :    
390 :     sub by_locus
391 :     {
392 :     my ($a, $b) = @_;
393 :    
394 :     my (undef, undef, $A_contig, $A_beg, $A_end, $A_len, $A_strand) = @$a;
395 :     my (undef, undef, $B_contig, $B_beg, $B_end, $B_len, $B_strand) = @$b;
396 :    
397 :     return ( ($A_contig cmp $B_contig)
398 :     || (&FIG::min($A_beg, $A_end) <=> &FIG::min($B_beg, $B_end))
399 :     || ($B_len <=> $A_len)
400 :     || ($A_strand cmp $B_strand)
401 :     );
402 :     }
403 :    
404 :     sub link {
405 :     my($org,$fid) = @_;
406 :    
407 :     return "<a href=diffsF.cgi?org1=$org&gene=$fid>$fid</a>";
408 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3