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

Annotation of /FigKernelScripts/compare_gto_called_CDSs.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : gdpusch 1.1 #!/usr/bin/env perl
2 :     # -*- perl -*-
3 :     # This is a SAS Component.
4 :     ########################################################################
5 :     # Copyright (c) 2003-2006 University of Chicago and Fellowship
6 :     # for Interpretations of Genomes. All Rights Reserved.
7 :     #
8 :     # This file is part of the SEED Toolkit.
9 :     #
10 :     # The SEED Toolkit is free software. You can redistribute
11 :     # it and/or modify it under the terms of the SEED Toolkit
12 :     # Public License.
13 :     #
14 :     # You should have received a copy of the SEED Toolkit Public License
15 :     # along with this program; if not write to the University of Chicago
16 :     # at info@ci.uchicago.edu or the Fellowship for Interpretation of
17 :     # Genomes at veronika@thefig.info or download a copy from
18 :     # http://www.theseed.org/LICENSE.TXT.
19 :     ########################################################################
20 :    
21 :     =head1 NAME
22 :    
23 :     compare-called-CDSs
24 :    
25 :     =head1 SYNOPSIS
26 :    
27 :     usage: compare-called-CDSs old.gto new.gto [summary.yaml] > comparison.tab 2> summary.txt
28 :    
29 :     =head1 DESCRIPTION
30 :    
31 :     detailed_description_of_purpose_of_script
32 :    
33 :     Example:
34 :    
35 :     example_of_use
36 :    
37 :     example_description
38 :    
39 :     =head1 COMMAND-LINE OPTIONS
40 :    
41 :     Usage: short_usage_msg
42 :    
43 :     -opt1
44 :    
45 :     =head1 AUTHORS
46 :    
47 :     L<The SEED Project|http://www.theseed.org>
48 :    
49 :     =cut
50 :    
51 :    
52 :     use strict;
53 :     use warnings;
54 :    
55 :     use SeedUtils;
56 :     use Data::Dumper;
57 :     use YAML::Any;
58 :     # use Carp;
59 :    
60 :     use Bio::KBase::GenomeAnnotation::Client;
61 :     use JSON::XS;
62 :    
63 :     my $help;
64 :     my $first_file;
65 :     my $second_file;
66 :     my $output_file;
67 :     my $summary_yaml;
68 :    
69 :     my $trouble;
70 :     use Getopt::Long;
71 :     my $rc = GetOptions('help' => \$help,
72 :     'file1=s' => \$first_file,
73 :     'file2=s' => \$second_file,
74 :     'output=s' => \$output_file,
75 :     'summary=s' => \$summary_yaml,
76 :     );
77 :    
78 :     if (@ARGV >= 2) {
79 :     $first_file ||= shift @ARGV;
80 :     if (!-s $first_file) {
81 :     $trouble = 1;
82 :     warn "ERROR: file1=\'$first_file\' does not exist or is empty\n";
83 :     }
84 :    
85 :     $second_file ||= shift @ARGV;
86 :     if (!-s $second_file) {
87 :     $trouble = 1;
88 :     warn "ERROR: file2=\'$second_file\' does not exist or is empty\n";
89 :     }
90 :     }
91 :    
92 :    
93 :     if (@ARGV == 1) {
94 :     $summary_yaml = shift @ARGV;
95 :     }
96 :    
97 :     if (!$rc || $help || $trouble || @ARGV != 0) {
98 :     seek(DATA, 0, 0);
99 :     while (<DATA>) {
100 :     last if /^=head1 COMMAND-LINE /;
101 :     }
102 :     while (<DATA>) {
103 :     last if (/^=/);
104 :     print $_;
105 :     }
106 :     exit($help ? 0 : 1);
107 :     }
108 :    
109 :    
110 :     my $out_fh;
111 :     if ($output_file) {
112 :     open($out_fh, ">", $output_file) or die "Cannot open $output_file: $!";
113 :     } else { $out_fh = \*STDOUT; }
114 :    
115 :    
116 :     my ($old_tbl, $old_num_pegs) = &load_gto($first_file);
117 :     my ($new_tbl, $new_num_pegs) = &load_gto($second_file);
118 :    
119 :    
120 :     use constant FID => 0;
121 :     use constant LOCUS => 1;
122 :     use constant CONTIG => 2;
123 :     use constant STRAND => 3;
124 :     use constant START => 4;
125 :     use constant STOP => 5;
126 :     use constant LEFT => 6;
127 :     use constant RIGHT => 7;
128 :     use constant LEN => 8;
129 :     use constant TYPE => 9;
130 :     use constant FUNC => 10;
131 :     use constant ALT_IDS => 11;
132 :     use constant ENTRY => 12;
133 :    
134 :    
135 :     my $identical = 0;
136 :     my $same_stop = 0;
137 :     my $differ = 0;
138 :     my $short = 0;
139 :     my $long = 0;
140 :     my $added = 0;
141 :     my $lost = 0;
142 :    
143 :    
144 :     my (%keys, @keys);
145 :     foreach my $key ((keys %$old_tbl), (keys %$new_tbl)) {
146 :     $keys{$key} = 1;
147 :     }
148 :     @keys = sort { &by_key($a,$b) } (keys %keys);
149 :    
150 :     print STDERR (q(Num keys = ), (scalar @keys), qq(\n\n)) if $ENV{VERBOSE};
151 :    
152 :     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));
153 :    
154 :     foreach my $key (sort { &by_key($a,$b) } @keys) {
155 :     my $case = q();
156 :    
157 :     my $old_fid = q();
158 :     my $old_func = q();
159 :     my $old_loc = q();
160 :     my $old_len = 0;
161 :     my $old_alt = q();
162 :     if (defined($old_tbl->{$key})) {
163 :     $old_fid = $old_tbl->{$key}->[FID];
164 :     $old_func = $old_tbl->{$key}->[FUNC];
165 :     $old_loc = $old_tbl->{$key}->[LOCUS];
166 :     $old_len = $old_tbl->{$key}->[LEN];
167 :     $old_alt = $old_tbl->{$key}->[ALT_IDS];
168 :     }
169 :    
170 :     my $new_fid = q();
171 :     my $new_func = q();
172 :     my $new_loc = q();
173 :     my $new_len = 0;
174 :     my $new_alt = q();
175 :     if (defined($new_tbl->{$key})) {
176 :     $new_fid = $new_tbl->{$key}->[FID];
177 :     $new_func = $new_tbl->{$key}->[FUNC];
178 :     $new_loc = $new_tbl->{$key}->[LOCUS];
179 :     $new_len = $new_tbl->{$key}->[LEN];
180 :     $new_alt = $new_tbl->{$key}->[ALT_IDS];
181 :     }
182 :    
183 :     if (defined($old_tbl->{$key})) {
184 :     if (not defined($new_tbl->{$key})) {
185 :     $case = q(lost);
186 :    
187 :     ++$lost;
188 :     ++$differ;
189 :     die Dumper($old_tbl->{$key}) unless $old_len;
190 :     }
191 :     else {
192 :     ++$same_stop;
193 :     if ($old_len == $new_len) {
194 :     $case = q(ident);
195 :    
196 :     ++$identical;
197 :     }
198 :     elsif ($old_len > $new_len) {
199 :     $case = q(short);
200 :    
201 :     ++$short;
202 :     ++$differ;
203 :     }
204 :     elsif ($old_len < $new_len) {
205 :     $case = q(long);
206 :    
207 :     ++$long;
208 :     ++$differ;
209 :     }
210 :     else {
211 :     die "Could not handle $key";
212 :     }
213 :     }
214 :     }
215 :     else {
216 :     $case = q(added);
217 :    
218 :     ++$added;
219 :     ++$differ;
220 :     }
221 :     my $diff = $new_len - $old_len;
222 :    
223 :     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));
224 :     }
225 :    
226 :     if (defined($summary_yaml))
227 :     {
228 :     if (open(my $fh, ">", $summary_yaml))
229 :     {
230 :     &write_summary_yaml($fh, $old_num_pegs, $new_num_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost);
231 :     }
232 :     else
233 :     {
234 :     die "Error opening $summary_yaml for writing: $!";
235 :     }
236 :     }
237 :     else
238 :     {
239 :     &write_summary($old_num_pegs, $new_num_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost);
240 :     }
241 :    
242 :     exit(0);
243 :    
244 :    
245 :    
246 :     sub load_gto {
247 :     my ($filename) = @_;
248 :     my ($tbl, $num_pegs);
249 :     my ($key, $id, $locus, $func, $contig, $beg, $end, $left, $right, $len, $strand, $type, $alt_ids);
250 :    
251 :     my $fh;
252 :     open($fh, "<", $filename) or die "Cannot open $filename: $!";
253 :    
254 :     my $json = JSON::XS->new;
255 :     my $gto;
256 :     {
257 :     local $/;
258 :     undef $/;
259 :     my $gto_txt = <$fh>;
260 :     $gto = $json->decode($gto_txt);
261 :     foreach my $feature (@ { $gto->{features} })
262 :     {
263 :     ($id, $contig, $strand,
264 :     $left, $right,
265 :     $beg, $end,
266 :     $len, $locus) = &feature_bounds($feature);
267 :    
268 :     $type = $feature->{type};
269 :     next unless ($type =~ m/^peg|CDS$/i);
270 :    
271 :     $func = $feature->{function} || q();
272 :    
273 :     $alt_ids = join(',', (defined($feature->{aliases}) ? @ { $feature->{aliases} } : q()));
274 :    
275 :     if (defined($contig) && $contig
276 :     && defined($beg) && $beg
277 :     && defined($end) && $end
278 :     && defined($len) && $len
279 :     && defined($strand) && $strand
280 :     )
281 :     {
282 :    
283 :     $key = join("\t", ($contig, $strand.$end));
284 :    
285 :     if (not defined($tbl->{$key})) {
286 :     ++$num_pegs;
287 :     $tbl->{$key} = [ $id, $locus, $contig, $strand, $beg, $end, $left, $right, $len, $type, $func, $alt_ids, $feature ];
288 :     }
289 :     else {
290 :     warn ("Skipping same-STOP TBL entry for $filename, $key:\n"
291 :     , Dumper($feature),
292 :     , "\n\n");
293 :     }
294 :     }
295 :     else {
296 :     warn ("INVALID ENTRY:\n", Dumper($feature), "\n\n");
297 :     }
298 :     }
299 :     }
300 :    
301 :     return ($tbl, $num_pegs);
302 :     }
303 :    
304 :     sub feature_bounds {
305 :     my ($feature) = @_;
306 :    
307 :     my $location = $feature->{location};
308 :    
309 :     my ($feature_contig,
310 :     $feature_strand,
311 :     $feature_left,
312 :     $feature_right,
313 :     ) = &parse_exon($location->[0]);
314 :    
315 :     foreach my $exon (@$location) {
316 :     my ($contig, $strand, $left, $right) = &parse_exon($exon);
317 :    
318 :     if ($feature_contig) {
319 :     if ($contig ne $feature_contig) {
320 :     warn ("Malformed feature --- \'$contig\' ne \'$feature_contig\':\n", Dumper($feature));
321 :     return ();
322 :     }
323 :     }
324 :    
325 :     if ($feature_strand) {
326 :     if ($strand ne $feature_strand) {
327 :     warn ("Malformed feature --- \'$strand\' ne \'$feature_strand\':\n", Dumper($feature));
328 :     return ();
329 :     }
330 :     }
331 :    
332 :     $feature_left = &SeedUtils::min($left, $feature_left);
333 :     $feature_right = &SeedUtils::max($right, $feature_right);
334 :     }
335 :    
336 :     my ($feature_beg, $feature_end) = ($feature_strand eq q(+))
337 :     ? ($feature_left, $feature_right)
338 :     : ($feature_right, $feature_left);
339 :    
340 :     my $feature_length = $feature_right - $feature_left + 1;
341 :    
342 :     my $feature_locus = join(',', map { join('', ($_->[0], q(_), $_->[1], $_->[2], $_->[3])) } @$location);
343 :    
344 :     return ($feature->{id}, $feature_contig, $feature_strand,
345 :     $feature_left, $feature_right,
346 :     $feature_beg, $feature_end,
347 :     $feature_length, $feature_locus);
348 :     }
349 :    
350 :     sub parse_exon {
351 :     my ($exon) = @_;
352 :     my ($contig, $beg, $strand, $len) = @$exon;
353 :     my $end = ($strand eq '+') ? $beg + ($len-1) : $beg - ($len-1);
354 :    
355 :     return ($contig, $strand, &SeedUtils::min($beg,$end), &SeedUtils::max($beg,$end));
356 :     }
357 :    
358 :    
359 :    
360 :     sub load_tbl
361 :     {
362 :     my ($file) = @_;
363 :     my ($tbl, $num_pegs);
364 :     my ($key, $entry, $id, $locus, $func, $rest, $contig, $beg, $end, $len, $strand, $type);
365 :    
366 :     open(TBL, "<$file") || die "Could not read-open \'$file\'";
367 :     while (defined($entry = <TBL>))
368 :     {
369 :     next if ($entry =~ m/^\#/);
370 :    
371 :     chomp $entry;
372 :     my @fields = split /\t/, $entry, -1;
373 :     $id = shift @fields;
374 :     $func = pop @fields;
375 :     $locus = pop @fields;
376 :    
377 :     $rest = join(q(,), (grep { $_ } @fields) );
378 :    
379 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
380 :     && defined($contig) && $contig
381 :     && defined($beg) && $beg
382 :     && defined($end) && $end
383 :     && defined($len) && $len
384 :     && defined($strand) && $strand
385 :     )
386 :     {
387 :     $key = "$contig\t$strand$end";
388 :    
389 :     $type = q(peg); #...Until such time as RNAs are handled properly...
390 :     # if (($type eq 'peg') || ($type eq 'orf')) {
391 :     ++$num_pegs;
392 :     # }
393 :     # else {
394 :     # warn "Unknown feature type: $entry\n";
395 :     # }
396 :    
397 :     if (not defined($tbl->{$key})) {
398 :     $tbl->{$key} = [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $entry, ($func || q()), $rest ];
399 :     }
400 :     else {
401 :     warn "Skipping same-STOP TBL entry for $file, $key:\n"
402 :     , "$tbl->{$key}->[ENTRY]\n$entry\n\n";
403 :     }
404 :     }
405 :     else {
406 :     warn "INVALID ENTRY:\t$entry\n";
407 :     }
408 :     }
409 :    
410 :     return ($tbl, $num_pegs);
411 :     }
412 :    
413 :     sub from_locus {
414 :     my ($locus) = @_;
415 :    
416 :     if ($locus) {
417 :     my ($contig, $left, $right, $dir) = SeedUtils::boundaries_of($locus);
418 :     my ($beg, $end) = ($dir eq q(+)) ? ($left, $right) : ($right, $left);
419 :    
420 :     if ($contig && $left && $right && $dir) {
421 :     return ($contig, $beg, $end, (1 + abs($right - $left)), $dir);
422 :     }
423 :     else {
424 :     die "Invalid locus $locus";
425 :     }
426 :     }
427 :     else {
428 :     die "Missing locus";
429 :     }
430 :    
431 :     return ();
432 :     }
433 :    
434 :     sub by_locus {
435 :     my ($x, $a, $b) = @_;
436 :    
437 :     my (undef, undef, $A_contig, $A_beg, $A_end, $A_len, $A_strand) = $x->[$a];
438 :     my (undef, undef, $B_contig, $B_beg, $B_end, $B_len, $B_strand) = $x->[$b];
439 :    
440 :     return ( ($A_contig cmp $B_contig)
441 :     || (&FIG::min($A_beg, $A_end) <=> &FIG::min($B_beg, $B_end))
442 :     || ($B_len <=> $A_len)
443 :     || ($A_strand cmp $B_strand)
444 :     );
445 :     }
446 :     sub by_key {
447 :     my ($x, $y) = @_;
448 :    
449 :     my ($X_contig, $X_strand, $X_end) = ($x =~ m/^([^\t]+)\t([+-])(\d+)/o);
450 :     my ($Y_contig, $Y_strand, $Y_end) = ($y =~ m/^([^\t]+)\t([+-])(\d+)/o);
451 :    
452 :     return ( ($X_contig cmp $Y_contig)
453 :     || ($X_end <=> $Y_end)
454 :     || ($X_strand cmp $Y_strand)
455 :     );
456 :     }
457 :    
458 :    
459 :     sub write_summary {
460 :     my ($old_pegs, $new_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost) = @_;
461 :    
462 :     print STDERR "old_num = $old_pegs PEGs\n";
463 :     print STDERR "new_num = $new_pegs PEGs\n\n";
464 :    
465 :     print STDERR ' Num %_Old %_New', qq(\n);
466 :     printf STDERR "same_stop = %4u %5.2f%% %5.2f%%\n"
467 :     , $same_stop, 100*$same_stop/$old_pegs, 100*$same_stop/$new_pegs;
468 :    
469 :    
470 :     printf STDERR "added = %4u %5.2f%% %5.2f%%\n"
471 :     , $added, 100*$added/$old_pegs, 100*$added/$new_pegs;
472 :    
473 :     printf STDERR "lost = %4u %5.2f%% %5.2f%%\n\n"
474 :     , $lost, 100*$lost/$old_pegs, 100*$lost/$new_pegs;
475 :    
476 :    
477 :     print STDERR ' Num %_Old %_New %_Common', qq(\n);
478 :     printf STDERR "identical = %4u %5.2f%% %5.2f%% %5.2f%%\n"
479 :     , $identical, 100*$identical/$old_pegs, 100*$identical/$new_pegs, 100*$identical/$same_stop;;
480 :    
481 :     printf STDERR "differ = %4u %5.2f%% %5.2f%% %5.2f%% %s\n"
482 :     , $differ, 100*$differ/$old_pegs, 100*$differ/$new_pegs, 100*$differ/$same_stop, q((Includes Lost and Added));
483 :    
484 :     printf STDERR "short = %4u %5.2f%% %5.2f%% %5.2f%%\n"
485 :     , $short, 100*$short/$old_pegs, 100*$short/$new_pegs, 100*$short/$same_stop;
486 :    
487 :     printf STDERR "long = %4u %5.2f%% %5.2f%% %5.2f%%\n\n"
488 :     , $long, 100*$long/$old_pegs, 100*$long/$new_pegs, 100*$long/$same_stop;
489 :    
490 :    
491 :     return 1;
492 :     }
493 :    
494 :    
495 :     sub write_summary_yaml {
496 :     my ($fh, $old_pegs, $new_pegs, $identical, $same_stop, $differ, $short, $long, $added, $lost) = @_;
497 :    
498 :     my $dat = {
499 :     old_num => $old_pegs,
500 :     new_num => $new_pegs,
501 :     };
502 :    
503 :     for my $what (qw(same_stop added lost identical differ short long))
504 :     {
505 :     my $val = eval "\$$what";
506 :     $dat->{$what} = $val;
507 :     $dat->{"${what}_pct_old"} = 100 * $val / $old_pegs;
508 :     $dat->{"${what}_pct_new"} = 100 * $val / $new_pegs;
509 :     }
510 :     for my $what (qw(identical differ short long))
511 :     {
512 :     my $val = eval "\$$what";
513 :     $dat->{"${what}_pct_common"} = 100 * $val / $same_stop;
514 :     }
515 :    
516 :     print $fh Dump($dat);
517 :     return 1;
518 :     }
519 :    
520 :     __DATA__

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3