[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.1 - (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 :     # usage: svr_compare_feature_tables old_features.tab new_fatures.tab > comparison.tab 2> summary.txt
25 :    
26 :     use strict;
27 :     use warnings;
28 :    
29 :     use SeedUtils;
30 :     use Data::Dumper;
31 :     # use Carp;
32 :    
33 :     $0 =~ m/([^\/]+)$/;
34 :     my $self = $1;
35 :     my $usage = "$self old_features.tab new_fatures.tab \> comparison.tab 2\> summary.txt";
36 :    
37 :     my $old_3col_file;
38 :     (($old_3col_file = shift) && (-f $old_3col_file))
39 :     || die "Could not find old_3col_file $old_3col_file\n\n\tusage: $usage\n\n";
40 :    
41 :     my $new_3col_file;
42 :     (($new_3col_file = shift) && (-f $new_3col_file))
43 :     || die "Could not find new_3col_file $new_3col_file\n\n\tusage: $usage\n\n";
44 :    
45 :     my ($old_tbl, $old_num_pegs) = &load_tbl($old_3col_file);
46 :     my ($new_tbl, $new_num_pegs) = &load_tbl($new_3col_file);
47 :    
48 :    
49 :     use constant FID => 0;
50 :     use constant LOCUS => 1;
51 :     use constant CONTIG => 2;
52 :     use constant START => 3;
53 :     use constant STOP => 4;
54 :     use constant LEN => 5;
55 :     use constant STRAND => 6;
56 :     use constant TYPE => 7;
57 :     use constant ENTRY => 8;
58 :     use constant FUNC => 9;
59 :    
60 :    
61 :     my $identical = 0;
62 :     my $same_stop = 0;
63 :     my $differ = 0;
64 :     my $short = 0;
65 :     my $long = 0;
66 :     my $added = 0;
67 :     my $lost = 0;
68 :    
69 :    
70 :     my (%keys, @keys);
71 :     foreach my $key ((keys %$old_tbl), (keys %$new_tbl)) {
72 :     $keys{$key} = 1;
73 :     }
74 :     @keys = sort { &by_key($a,$b) } (keys %keys);
75 :     print STDERR (q(Num keys = ), (scalar @keys), qq(\n\n)) if $ENV{VERBOSE};
76 :    
77 :    
78 :     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)), qq(\n));
79 :     foreach my $key (sort { &by_key($a,$b) } @keys) {
80 :     my $case = q();
81 :    
82 :     my $old_fid = q();
83 :     my $old_func = q();
84 :     my $old_loc = q();
85 :     my $old_len = 0;
86 :     if (defined($old_tbl->{$key})) {
87 :     $old_fid = $old_tbl->{$key}->[FID];
88 :     $old_func = $old_tbl->{$key}->[FUNC];
89 :     $old_loc = $old_tbl->{$key}->[LOCUS];
90 :     $old_len = $old_tbl->{$key}->[LEN];
91 :     }
92 :    
93 :     my $new_fid = q();
94 :     my $new_func = q();
95 :     my $new_loc = q();
96 :     my $new_len = 0;
97 :     if (defined($new_tbl->{$key})) {
98 :     $new_fid = $new_tbl->{$key}->[FID];
99 :     $new_func = $new_tbl->{$key}->[FUNC];
100 :     $new_loc = $new_tbl->{$key}->[LOCUS];
101 :     $new_len = $new_tbl->{$key}->[LEN];
102 :     }
103 :    
104 :     if (defined($old_tbl->{$key})) {
105 :     if (not defined($new_tbl->{$key})) {
106 :     $case = q(lost);
107 :    
108 :     ++$lost;
109 :     ++$differ;
110 :     die Dumper($old_tbl->{$key}) unless $old_len;
111 :     }
112 :     else {
113 :     ++$same_stop;
114 :     if ($old_len == $new_len) {
115 :     $case = q(ident);
116 :    
117 :     ++$identical;
118 :     }
119 :     elsif ($old_len > $new_len) {
120 :     $case = q(short);
121 :    
122 :     ++$short;
123 :     ++$differ;
124 :     }
125 :     elsif ($old_len < $new_len) {
126 :     $case = q(long);
127 :    
128 :     ++$long;
129 :     ++$differ;
130 :     }
131 :     else {
132 :     die "Could not handle $key";
133 :     }
134 :     }
135 :     }
136 :     else {
137 :     $case = q(added);
138 :    
139 :     ++$added;
140 :     ++$differ;
141 :     }
142 :     my $diff = $new_len - $old_len;
143 :    
144 :     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));
145 :     }
146 :    
147 :     &write_summary($old_num_pegs, $new_num_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost);
148 :    
149 :     exit(0);
150 :    
151 :    
152 :    
153 :     sub load_tbl
154 :     {
155 :     my ($file) = @_;
156 :     my ($tbl, $num_pegs);
157 :     my ($key, $entry, $id, $locus, $func, $contig, $beg, $end, $len, $strand, $type);
158 :    
159 :     open(TBL, "<$file") || die "Could not read-open \'$file\'";
160 :     while (defined($entry = <TBL>))
161 :     {
162 :     chomp $entry;
163 :     ($id, $locus, $func) = split /\t/, $entry;
164 :    
165 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
166 :     && defined($contig) && $contig
167 :     && defined($beg) && $beg
168 :     && defined($end) && $end
169 :     && defined($len) && $len
170 :     && defined($strand) && $strand
171 :     )
172 :     {
173 :     $key = "$contig\t$strand$end";
174 :    
175 :     $type = q(peg); #...Until such time as RNAs are handled properly...
176 :     # if (($type eq 'peg') || ($type eq 'orf')) {
177 :     ++$num_pegs;
178 :     # }
179 :     # else {
180 :     # warn "Unknown feature type: $entry\n";
181 :     # }
182 :    
183 :     if (not defined($tbl->{$key})) {
184 :     $tbl->{$key} = [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $entry, ($func || q()), ];
185 :     }
186 :     else {
187 :     warn "Skipping same-STOP TBL entry for $file, $key:\n"
188 :     , "$tbl->{$key}->[ENTRY]\n$entry\n\n";
189 :     }
190 :     }
191 :     else {
192 :     warn "INVALID ENTRY:\t$entry\n";
193 :     }
194 :     }
195 :    
196 :     return ($tbl, $num_pegs);
197 :     }
198 :    
199 :     sub from_locus {
200 :     my ($locus) = @_;
201 :    
202 :     if ($locus) {
203 :     my ($contig, $left, $right, $dir) = SeedUtils::boundaries_of($locus);
204 :     my ($beg, $end) = ($dir eq q(+)) ? ($left, $right) : ($right, $left);
205 :    
206 :     return ($contig, $beg, $end, (1 + abs($right - $left)), $dir);
207 :     }
208 :     else {
209 :     die "Invalid locus $locus";
210 :     }
211 :    
212 :     return ();
213 :     }
214 :    
215 :     sub by_locus {
216 :     my ($x, $a, $b) = @_;
217 :    
218 :     my (undef, undef, $A_contig, $A_beg, $A_end, $A_len, $A_strand) = $x->[$a];
219 :     my (undef, undef, $B_contig, $B_beg, $B_end, $B_len, $B_strand) = $x->[$b];
220 :    
221 :     return ( ($A_contig cmp $B_contig)
222 :     || (&FIG::min($A_beg, $A_end) <=> &FIG::min($B_beg, $B_end))
223 :     || ($B_len <=> $A_len)
224 :     || ($A_strand cmp $B_strand)
225 :     );
226 :     }
227 :     sub by_key {
228 :     my ($x, $y) = @_;
229 :    
230 :     my ($X_contig, $X_strand, $X_end) = ($x =~ m/^([^\t]+)\t([+-])(\d+)/o);
231 :     my ($Y_contig, $Y_strand, $Y_end) = ($y =~ m/^([^\t]+)\t([+-])(\d+)/o);
232 :    
233 :     return ( ($X_contig cmp $Y_contig)
234 :     || ($X_end <=> $Y_end)
235 :     || ($X_strand cmp $Y_strand)
236 :     );
237 :     }
238 :    
239 :    
240 :     sub write_summary {
241 :     my ($old_pegs, $new_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost) = @_;
242 :    
243 :     print STDERR "old_num = $old_pegs PEGs\n";
244 :     print STDERR "new_num = $new_pegs PEGs\n\n";
245 :    
246 :     print STDERR ' Num %_Old %_New', qq(\n);
247 :     printf STDERR "same_stop = %4u %5.2f%% %5.2f%%\n"
248 :     , $same_stop, 100*$same_stop/$old_pegs, 100*$same_stop/$new_pegs;
249 :    
250 :    
251 :     printf STDERR "added = %4u %5.2f%% %5.2f%%\n"
252 :     , $added, 100*$added/$old_pegs, 100*$added/$new_pegs;
253 :    
254 :     printf STDERR "lost = %4u %5.2f%% %5.2f%%\n\n"
255 :     , $lost, 100*$lost/$old_pegs, 100*$lost/$new_pegs;
256 :    
257 :    
258 :     print STDERR ' Num %_Old %_New %_Common', qq(\n);
259 :     printf STDERR "identical = %4u %5.2f%% %5.2f%% %5.2f%%\n"
260 :     , $identical, 100*$identical/$old_pegs, 100*$identical/$new_pegs, 100*$identical/$same_stop;;
261 :    
262 :     printf STDERR "differ = %4u %5.2f%% %5.2f%% %5.2f%% %s\n"
263 :     , $differ, 100*$differ/$old_pegs, 100*$differ/$new_pegs, 100*$differ/$same_stop, q((Includes Lost and Added));
264 :    
265 :     printf STDERR "short = %4u %5.2f%% %5.2f%% %5.2f%%\n"
266 :     , $short, 100*$short/$old_pegs, 100*$short/$new_pegs, 100*$short/$same_stop;
267 :    
268 :     printf STDERR "long = %4u %5.2f%% %5.2f%% %5.2f%%\n\n"
269 :     , $long, 100*$long/$old_pegs, 100*$long/$new_pegs, 100*$long/$same_stop;
270 :    
271 :    
272 :     return 1;
273 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3