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

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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3