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

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3