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

Annotation of /FigKernelScripts/svr_compare_feature_tables.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : gdpusch 1.1 # -*- perl -*-
2 :    
3 :     #
4 :     # This is a SAS Component.
5 :     #
6 :    
7 :     #
8 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
9 :     # for Interpretations of Genomes. All Rights Reserved.
10 :     #
11 :     # This file is part of the SEED Toolkit.
12 :     #
13 :     # The SEED Toolkit is free software. You can redistribute
14 :     # it and/or modify it under the terms of the SEED Toolkit
15 :     # Public License.
16 :     #
17 :     # You should have received a copy of the SEED Toolkit Public License
18 :     # along with this program; if not write to the University of Chicago
19 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
20 :     # Genomes at veronika@thefig.info or download a copy from
21 :     # http://www.theseed.org/LICENSE.TXT.
22 :     #
23 :    
24 : olson 1.2 #
25 :     # Input files:
26 :     #
27 :     # ID seed-format-location function
28 :     #
29 :    
30 : gdpusch 1.4 # usage: svr_compare_feature_tables old_features.tab new_fatures.tab [summary.yaml] > comparison.tab 2> summary.txt
31 : gdpusch 1.1
32 :     use strict;
33 :     use warnings;
34 :    
35 :     use SeedUtils;
36 :     use Data::Dumper;
37 : olson 1.2 use YAML::Any;
38 : gdpusch 1.1 # use Carp;
39 :    
40 :     $0 =~ m/([^\/]+)$/;
41 :     my $self = $1;
42 : olson 1.2 my $usage = "$self old_features.tab new_fatures.tab [summary.yaml] \> comparison.tab 2\> summary.txt";
43 : gdpusch 1.1
44 : gdpusch 1.4 my $old_tab_file;
45 :     (($old_tab_file = shift) && (-f $old_tab_file))
46 :     || die "Could not find old_tab_file $old_tab_file\n\n\tusage: $usage\n\n";
47 :    
48 :     my $new_tab_file;
49 :     (($new_tab_file = shift) && (-f $new_tab_file))
50 :     || die "Could not find new_tab_file $new_tab_file\n\n\tusage: $usage\n\n";
51 : gdpusch 1.1
52 : olson 1.2 my $summary_yaml;
53 : gdpusch 1.4 if (@ARGV) {
54 : olson 1.2 $summary_yaml = shift;
55 :     }
56 :    
57 : gdpusch 1.4 my ($old_tbl, $old_num_pegs) = &load_tbl($old_tab_file);
58 :     my ($new_tbl, $new_num_pegs) = &load_tbl($new_tab_file);
59 : gdpusch 1.1
60 :    
61 : gdpusch 1.4 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 ENTRY => 8;
70 :     use constant FUNC => 9;
71 :     use constant ALT_IDS => 10;
72 : gdpusch 1.1
73 :    
74 :     my $identical = 0;
75 :     my $same_stop = 0;
76 :     my $differ = 0;
77 :     my $short = 0;
78 :     my $long = 0;
79 :     my $added = 0;
80 :     my $lost = 0;
81 :    
82 :    
83 :     my (%keys, @keys);
84 :     foreach my $key ((keys %$old_tbl), (keys %$new_tbl)) {
85 :     $keys{$key} = 1;
86 :     }
87 :     @keys = sort { &by_key($a,$b) } (keys %keys);
88 : olson 1.2
89 : gdpusch 1.1 print STDERR (q(Num keys = ), (scalar @keys), qq(\n\n)) if $ENV{VERBOSE};
90 :    
91 : gdpusch 1.4 print STDOUT (q(#), join(qq(\t), qw(Comparison Old_ID New_ID Old_Length New_Length Length_Diff Old_Loc New_Loc Old_Function New_Function Old_Alt_IDs New_Alt_IDs)), qq(\n));
92 :    
93 : gdpusch 1.1 foreach my $key (sort { &by_key($a,$b) } @keys) {
94 :     my $case = q();
95 :    
96 :     my $old_fid = q();
97 :     my $old_func = q();
98 :     my $old_loc = q();
99 :     my $old_len = 0;
100 : gdpusch 1.4 my $old_alt = q();
101 : gdpusch 1.1 if (defined($old_tbl->{$key})) {
102 :     $old_fid = $old_tbl->{$key}->[FID];
103 :     $old_func = $old_tbl->{$key}->[FUNC];
104 :     $old_loc = $old_tbl->{$key}->[LOCUS];
105 :     $old_len = $old_tbl->{$key}->[LEN];
106 : gdpusch 1.4 $old_alt = $old_tbl->{$key}->[ALT_IDS];
107 : gdpusch 1.1 }
108 :    
109 :     my $new_fid = q();
110 :     my $new_func = q();
111 :     my $new_loc = q();
112 :     my $new_len = 0;
113 : gdpusch 1.4 my $new_alt = q();
114 : gdpusch 1.1 if (defined($new_tbl->{$key})) {
115 :     $new_fid = $new_tbl->{$key}->[FID];
116 :     $new_func = $new_tbl->{$key}->[FUNC];
117 :     $new_loc = $new_tbl->{$key}->[LOCUS];
118 :     $new_len = $new_tbl->{$key}->[LEN];
119 : gdpusch 1.4 $new_alt = $new_tbl->{$key}->[ALT_IDS];
120 : gdpusch 1.1 }
121 :    
122 :     if (defined($old_tbl->{$key})) {
123 :     if (not defined($new_tbl->{$key})) {
124 :     $case = q(lost);
125 :    
126 :     ++$lost;
127 :     ++$differ;
128 :     die Dumper($old_tbl->{$key}) unless $old_len;
129 :     }
130 :     else {
131 :     ++$same_stop;
132 :     if ($old_len == $new_len) {
133 :     $case = q(ident);
134 :    
135 :     ++$identical;
136 :     }
137 :     elsif ($old_len > $new_len) {
138 :     $case = q(short);
139 :    
140 :     ++$short;
141 :     ++$differ;
142 :     }
143 :     elsif ($old_len < $new_len) {
144 :     $case = q(long);
145 :    
146 :     ++$long;
147 :     ++$differ;
148 :     }
149 :     else {
150 :     die "Could not handle $key";
151 :     }
152 :     }
153 :     }
154 :     else {
155 :     $case = q(added);
156 :    
157 :     ++$added;
158 :     ++$differ;
159 :     }
160 :     my $diff = $new_len - $old_len;
161 :    
162 : gdpusch 1.4 print STDOUT (join(qq(\t), ($case, $old_fid, $new_fid, $old_len, $new_len, $diff, $old_loc, $new_loc, $old_func, $new_func, $old_alt, $new_alt)), qq(\n));
163 :     }
164 : olson 1.2
165 :     if (defined($summary_yaml))
166 :     {
167 :     if (open(my $fh, ">", $summary_yaml))
168 :     {
169 :     &write_summary_yaml($fh, $old_num_pegs, $new_num_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost);
170 :     }
171 :     else
172 :     {
173 :     die "Error opening $summary_yaml for writing: $!";
174 :     }
175 :     }
176 :     else
177 :     {
178 :     &write_summary($old_num_pegs, $new_num_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost);
179 : gdpusch 1.1 }
180 :    
181 :     exit(0);
182 :    
183 :    
184 :    
185 :     sub load_tbl
186 :     {
187 :     my ($file) = @_;
188 :     my ($tbl, $num_pegs);
189 : gdpusch 1.4 my ($key, $entry, $id, $locus, $func, $rest, $contig, $beg, $end, $len, $strand, $type);
190 : gdpusch 1.1
191 :     open(TBL, "<$file") || die "Could not read-open \'$file\'";
192 :     while (defined($entry = <TBL>))
193 :     {
194 : gdpusch 1.4 next if ($entry =~ m/^\#/);
195 :    
196 : gdpusch 1.1 chomp $entry;
197 : olson 1.5 my @fields = split /\t/, $entry, -1;
198 : gdpusch 1.4 $id = shift @fields;
199 :     $func = pop @fields;
200 :     $locus = pop @fields;
201 :    
202 :     $rest = join(q(,), (grep { $_ } @fields) );
203 : gdpusch 1.1
204 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
205 :     && defined($contig) && $contig
206 :     && defined($beg) && $beg
207 :     && defined($end) && $end
208 :     && defined($len) && $len
209 :     && defined($strand) && $strand
210 :     )
211 :     {
212 :     $key = "$contig\t$strand$end";
213 :    
214 :     $type = q(peg); #...Until such time as RNAs are handled properly...
215 :     # if (($type eq 'peg') || ($type eq 'orf')) {
216 :     ++$num_pegs;
217 :     # }
218 :     # else {
219 :     # warn "Unknown feature type: $entry\n";
220 :     # }
221 :    
222 :     if (not defined($tbl->{$key})) {
223 : gdpusch 1.4 $tbl->{$key} = [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $entry, ($func || q()), $rest ];
224 : gdpusch 1.1 }
225 :     else {
226 :     warn "Skipping same-STOP TBL entry for $file, $key:\n"
227 :     , "$tbl->{$key}->[ENTRY]\n$entry\n\n";
228 :     }
229 :     }
230 :     else {
231 :     warn "INVALID ENTRY:\t$entry\n";
232 :     }
233 :     }
234 :    
235 :     return ($tbl, $num_pegs);
236 :     }
237 :    
238 :     sub from_locus {
239 :     my ($locus) = @_;
240 :    
241 :     if ($locus) {
242 :     my ($contig, $left, $right, $dir) = SeedUtils::boundaries_of($locus);
243 :     my ($beg, $end) = ($dir eq q(+)) ? ($left, $right) : ($right, $left);
244 :    
245 : gdpusch 1.4 if ($contig && $left && $right && $dir) {
246 :     return ($contig, $beg, $end, (1 + abs($right - $left)), $dir);
247 :     }
248 :     else {
249 :     die "Invalid locus $locus";
250 :     }
251 : gdpusch 1.1 }
252 :     else {
253 : gdpusch 1.4 die "Missing locus";
254 : gdpusch 1.1 }
255 :    
256 :     return ();
257 :     }
258 :    
259 :     sub by_locus {
260 :     my ($x, $a, $b) = @_;
261 :    
262 :     my (undef, undef, $A_contig, $A_beg, $A_end, $A_len, $A_strand) = $x->[$a];
263 :     my (undef, undef, $B_contig, $B_beg, $B_end, $B_len, $B_strand) = $x->[$b];
264 :    
265 :     return ( ($A_contig cmp $B_contig)
266 :     || (&FIG::min($A_beg, $A_end) <=> &FIG::min($B_beg, $B_end))
267 :     || ($B_len <=> $A_len)
268 :     || ($A_strand cmp $B_strand)
269 :     );
270 :     }
271 :     sub by_key {
272 :     my ($x, $y) = @_;
273 :    
274 :     my ($X_contig, $X_strand, $X_end) = ($x =~ m/^([^\t]+)\t([+-])(\d+)/o);
275 :     my ($Y_contig, $Y_strand, $Y_end) = ($y =~ m/^([^\t]+)\t([+-])(\d+)/o);
276 :    
277 :     return ( ($X_contig cmp $Y_contig)
278 :     || ($X_end <=> $Y_end)
279 :     || ($X_strand cmp $Y_strand)
280 :     );
281 :     }
282 :    
283 :    
284 :     sub write_summary {
285 :     my ($old_pegs, $new_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost) = @_;
286 :    
287 :     print STDERR "old_num = $old_pegs PEGs\n";
288 :     print STDERR "new_num = $new_pegs PEGs\n\n";
289 :    
290 :     print STDERR ' Num %_Old %_New', qq(\n);
291 :     printf STDERR "same_stop = %4u %5.2f%% %5.2f%%\n"
292 :     , $same_stop, 100*$same_stop/$old_pegs, 100*$same_stop/$new_pegs;
293 :    
294 :    
295 :     printf STDERR "added = %4u %5.2f%% %5.2f%%\n"
296 :     , $added, 100*$added/$old_pegs, 100*$added/$new_pegs;
297 :    
298 :     printf STDERR "lost = %4u %5.2f%% %5.2f%%\n\n"
299 :     , $lost, 100*$lost/$old_pegs, 100*$lost/$new_pegs;
300 :    
301 :    
302 :     print STDERR ' Num %_Old %_New %_Common', qq(\n);
303 :     printf STDERR "identical = %4u %5.2f%% %5.2f%% %5.2f%%\n"
304 :     , $identical, 100*$identical/$old_pegs, 100*$identical/$new_pegs, 100*$identical/$same_stop;;
305 :    
306 :     printf STDERR "differ = %4u %5.2f%% %5.2f%% %5.2f%% %s\n"
307 :     , $differ, 100*$differ/$old_pegs, 100*$differ/$new_pegs, 100*$differ/$same_stop, q((Includes Lost and Added));
308 :    
309 :     printf STDERR "short = %4u %5.2f%% %5.2f%% %5.2f%%\n"
310 :     , $short, 100*$short/$old_pegs, 100*$short/$new_pegs, 100*$short/$same_stop;
311 :    
312 :     printf STDERR "long = %4u %5.2f%% %5.2f%% %5.2f%%\n\n"
313 :     , $long, 100*$long/$old_pegs, 100*$long/$new_pegs, 100*$long/$same_stop;
314 :    
315 :    
316 :     return 1;
317 :     }
318 : gdpusch 1.4
319 :    
320 : olson 1.2 sub write_summary_yaml {
321 :     my ($fh, $old_pegs, $new_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost) = @_;
322 :    
323 :     my $dat = {
324 :     old_num => $old_pegs,
325 :     new_num => $new_pegs,
326 :     };
327 :    
328 :     for my $what (qw(same_stop added lost identical differ short long))
329 :     {
330 :     my $val = eval "\$$what";
331 :     $dat->{$what} = $val;
332 :     $dat->{"${what}_pct_old"} = 100 * $val / $old_pegs;
333 :     $dat->{"${what}_pct_new"} = 100 * $val / $new_pegs;
334 :     }
335 : olson 1.3 for my $what (qw(identical differ short long))
336 :     {
337 :     my $val = eval "\$$what";
338 :     $dat->{"${what}_pct_common"} = 100 * $val / $same_stop;
339 :     }
340 : olson 1.2
341 :     print $fh Dump($dat);
342 :     return 1;
343 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3