[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.2 - (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 : gdpusch 1.2 use Data::Dumper;
55 :     #use Carp;
56 : gdpusch 1.1
57 : gdpusch 1.2 use GenomeTypeObject;
58 : gdpusch 1.1 use SeedUtils;
59 :     use YAML::Any;
60 :    
61 : gdpusch 1.2 #use Bio::KBase::GenomeAnnotation::Client;
62 :     #use JSON::XS;
63 :    
64 : gdpusch 1.1
65 :     my $help;
66 :     my $first_file;
67 :     my $second_file;
68 :     my $output_file;
69 :     my $summary_yaml;
70 :    
71 :     my $trouble;
72 :     use Getopt::Long;
73 :     my $rc = GetOptions('help' => \$help,
74 :     'file1=s' => \$first_file,
75 :     'file2=s' => \$second_file,
76 :     'output=s' => \$output_file,
77 :     'summary=s' => \$summary_yaml,
78 :     );
79 :    
80 :     if (@ARGV >= 2) {
81 :     $first_file ||= shift @ARGV;
82 :     if (!-s $first_file) {
83 :     $trouble = 1;
84 :     warn "ERROR: file1=\'$first_file\' does not exist or is empty\n";
85 :     }
86 :    
87 :     $second_file ||= shift @ARGV;
88 :     if (!-s $second_file) {
89 :     $trouble = 1;
90 :     warn "ERROR: file2=\'$second_file\' does not exist or is empty\n";
91 :     }
92 :     }
93 :    
94 :    
95 :     if (@ARGV == 1) {
96 :     $summary_yaml = shift @ARGV;
97 :     }
98 :    
99 :     if (!$rc || $help || $trouble || @ARGV != 0) {
100 :     seek(DATA, 0, 0);
101 :     while (<DATA>) {
102 :     last if /^=head1 COMMAND-LINE /;
103 :     }
104 :     while (<DATA>) {
105 :     last if (/^=/);
106 :     print $_;
107 :     }
108 :     exit($help ? 0 : 1);
109 :     }
110 :    
111 :    
112 :     my $out_fh;
113 :     if ($output_file) {
114 :     open($out_fh, ">", $output_file) or die "Cannot open $output_file: $!";
115 :     } else { $out_fh = \*STDOUT; }
116 :    
117 :    
118 : gdpusch 1.2 my ($old_tbl, $old_num_pegs, $old_metrics) = &load_gto($first_file);
119 :     my ($new_tbl, $new_num_pegs, $new_metrics) = &load_gto($second_file);
120 :     #die (Dump($old_metrics), Dump($new_metrics));
121 :     #die Dumper($old_num_pegs, $new_num_pegs, $old_tbl, $new_tbl);
122 : gdpusch 1.1
123 :     use constant FID => 0;
124 :     use constant LOCUS => 1;
125 :     use constant CONTIG => 2;
126 :     use constant STRAND => 3;
127 :     use constant START => 4;
128 :     use constant STOP => 5;
129 :     use constant LEFT => 6;
130 :     use constant RIGHT => 7;
131 :     use constant LEN => 8;
132 :     use constant TYPE => 9;
133 :     use constant FUNC => 10;
134 :     use constant ALT_IDS => 11;
135 :     use constant ENTRY => 12;
136 :    
137 :    
138 :     my $identical = 0;
139 :     my $same_stop = 0;
140 :     my $differ = 0;
141 :     my $short = 0;
142 :     my $long = 0;
143 :     my $added = 0;
144 :     my $lost = 0;
145 :    
146 :    
147 :     my (%keys, @keys);
148 :     foreach my $key ((keys %$old_tbl), (keys %$new_tbl)) {
149 :     $keys{$key} = 1;
150 :     }
151 :     @keys = sort { &by_key($a,$b) } (keys %keys);
152 : gdpusch 1.2 my $num_union = @keys;
153 : gdpusch 1.1
154 :     print STDERR (q(Num keys = ), (scalar @keys), qq(\n\n)) if $ENV{VERBOSE};
155 :    
156 :     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));
157 :    
158 :     foreach my $key (sort { &by_key($a,$b) } @keys) {
159 :     my $case = q();
160 :    
161 :     my $old_fid = q();
162 :     my $old_func = q();
163 :     my $old_loc = q();
164 :     my $old_len = 0;
165 :     my $old_alt = q();
166 :     if (defined($old_tbl->{$key})) {
167 :     $old_fid = $old_tbl->{$key}->[FID];
168 :     $old_func = $old_tbl->{$key}->[FUNC];
169 :     $old_loc = $old_tbl->{$key}->[LOCUS];
170 :     $old_len = $old_tbl->{$key}->[LEN];
171 :     $old_alt = $old_tbl->{$key}->[ALT_IDS];
172 :     }
173 :    
174 :     my $new_fid = q();
175 :     my $new_func = q();
176 :     my $new_loc = q();
177 :     my $new_len = 0;
178 :     my $new_alt = q();
179 :     if (defined($new_tbl->{$key})) {
180 :     $new_fid = $new_tbl->{$key}->[FID];
181 :     $new_func = $new_tbl->{$key}->[FUNC];
182 :     $new_loc = $new_tbl->{$key}->[LOCUS];
183 :     $new_len = $new_tbl->{$key}->[LEN];
184 :     $new_alt = $new_tbl->{$key}->[ALT_IDS];
185 :     }
186 :    
187 :     if (defined($old_tbl->{$key})) {
188 :     if (not defined($new_tbl->{$key})) {
189 :     $case = q(lost);
190 :    
191 :     ++$lost;
192 :     ++$differ;
193 :     die Dumper($old_tbl->{$key}) unless $old_len;
194 :     }
195 :     else {
196 :     ++$same_stop;
197 :     if ($old_len == $new_len) {
198 :     $case = q(ident);
199 :    
200 :     ++$identical;
201 :     }
202 :     elsif ($old_len > $new_len) {
203 :     $case = q(short);
204 :    
205 :     ++$short;
206 :     ++$differ;
207 :     }
208 :     elsif ($old_len < $new_len) {
209 :     $case = q(long);
210 :    
211 :     ++$long;
212 :     ++$differ;
213 :     }
214 :     else {
215 :     die "Could not handle $key";
216 :     }
217 :     }
218 :     }
219 :     else {
220 :     $case = q(added);
221 :    
222 :     ++$added;
223 :     ++$differ;
224 :     }
225 :     my $diff = $new_len - $old_len;
226 :    
227 :     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));
228 :     }
229 :    
230 :     if (defined($summary_yaml))
231 :     {
232 :     if (open(my $fh, ">", $summary_yaml))
233 :     {
234 : gdpusch 1.2 &write_summary_yaml($fh, $old_metrics, $new_metrics,
235 :     $old_num_pegs, $new_num_pegs, $num_union,
236 :     $identical, $same_stop, $differ,
237 :     $short, $long, $added, $lost);
238 : gdpusch 1.1 }
239 :     else
240 :     {
241 :     die "Error opening $summary_yaml for writing: $!";
242 :     }
243 :     }
244 : gdpusch 1.2
245 :     &write_summary($old_metrics, $new_metrics,
246 :     $old_num_pegs, $new_num_pegs, $num_union,
247 :     $identical, $same_stop, $differ,
248 :     $short, $long, $added, $lost);
249 : gdpusch 1.1
250 :     exit(0);
251 :    
252 :    
253 :    
254 :     sub load_gto {
255 :     my ($filename) = @_;
256 :     my ($tbl, $num_pegs);
257 :     my ($key, $id, $locus, $func, $contig, $beg, $end, $left, $right, $len, $strand, $type, $alt_ids);
258 :    
259 : gdpusch 1.2 #my $fh;
260 :     #open($fh, "<", $filename) or die "Cannot open $filename: $!";
261 :     #my $json = JSON::XS->new;
262 :    
263 :     my $gto = GenomeTypeObject->new({ file => $filename });
264 :     my $metrics = $gto->metrics();
265 :     $metrics->{gc_content} = $gto->compute_contigs_gc();
266 : gdpusch 1.1
267 : gdpusch 1.2 my $total_coding = 0;
268 :     foreach my $feature (@ { $gto->{features} })
269 : gdpusch 1.1 {
270 : gdpusch 1.2 ($id, $contig, $strand,
271 :     $left, $right,
272 :     $beg, $end,
273 :     $len, $locus) = &feature_bounds($feature);
274 :    
275 :     $type = $feature->{type};
276 :     next unless ($type =~ m/^peg|CDS$/i);
277 :    
278 :     $total_coding += $len;
279 :    
280 :     $func = $feature->{function} || q();
281 :    
282 :     $alt_ids = join(',', (defined($feature->{aliases}) ? @ { $feature->{aliases} } : q()));
283 :    
284 :     if (defined($contig) && $contig
285 :     && defined($beg) && $beg
286 :     && defined($end) && $end
287 :     && defined($len) && $len
288 :     && defined($strand) && $strand
289 :     )
290 : gdpusch 1.1 {
291 :    
292 : gdpusch 1.2 $key = join("\t", ($contig, $strand.$end));
293 : gdpusch 1.1
294 : gdpusch 1.2 if (not defined($tbl->{$key})) {
295 :     ++$num_pegs;
296 :     $tbl->{$key} = [ $id, $locus, $contig, $strand,
297 :     $beg, $end, $left, $right, $len,
298 :     $type, $func, $alt_ids, $feature ];
299 : gdpusch 1.1 }
300 :     else {
301 :     warn ("Skipping same-STOP TBL entry for $filename, $key:\n"
302 :     , Dumper($feature),
303 :     , "\n\n");
304 :     }
305 : gdpusch 1.2 }
306 :     else {
307 :     warn ("INVALID ENTRY:\n", Dumper($feature), "\n\n");
308 : gdpusch 1.1 }
309 :     }
310 : gdpusch 1.2 $metrics->{coding_fraction} = 100.0 * $total_coding / $metrics->{totlen};
311 : gdpusch 1.1
312 : gdpusch 1.2 return ($tbl, $num_pegs, $metrics);
313 : gdpusch 1.1 }
314 :    
315 :     sub feature_bounds {
316 :     my ($feature) = @_;
317 :    
318 :     my $location = $feature->{location};
319 :    
320 :     my ($feature_contig,
321 :     $feature_strand,
322 :     $feature_left,
323 :     $feature_right,
324 :     ) = &parse_exon($location->[0]);
325 :    
326 :     foreach my $exon (@$location) {
327 :     my ($contig, $strand, $left, $right) = &parse_exon($exon);
328 :    
329 :     if ($feature_contig) {
330 :     if ($contig ne $feature_contig) {
331 :     warn ("Malformed feature --- \'$contig\' ne \'$feature_contig\':\n", Dumper($feature));
332 :     return ();
333 :     }
334 :     }
335 :    
336 :     if ($feature_strand) {
337 :     if ($strand ne $feature_strand) {
338 :     warn ("Malformed feature --- \'$strand\' ne \'$feature_strand\':\n", Dumper($feature));
339 :     return ();
340 :     }
341 :     }
342 :    
343 :     $feature_left = &SeedUtils::min($left, $feature_left);
344 :     $feature_right = &SeedUtils::max($right, $feature_right);
345 :     }
346 :    
347 :     my ($feature_beg, $feature_end) = ($feature_strand eq q(+))
348 :     ? ($feature_left, $feature_right)
349 :     : ($feature_right, $feature_left);
350 :    
351 :     my $feature_length = $feature_right - $feature_left + 1;
352 :    
353 :     my $feature_locus = join(',', map { join('', ($_->[0], q(_), $_->[1], $_->[2], $_->[3])) } @$location);
354 :    
355 :     return ($feature->{id}, $feature_contig, $feature_strand,
356 :     $feature_left, $feature_right,
357 :     $feature_beg, $feature_end,
358 :     $feature_length, $feature_locus);
359 :     }
360 :    
361 :     sub parse_exon {
362 :     my ($exon) = @_;
363 :     my ($contig, $beg, $strand, $len) = @$exon;
364 :     my $end = ($strand eq '+') ? $beg + ($len-1) : $beg - ($len-1);
365 :    
366 :     return ($contig, $strand, &SeedUtils::min($beg,$end), &SeedUtils::max($beg,$end));
367 :     }
368 :    
369 :    
370 :    
371 :     sub load_tbl
372 :     {
373 :     my ($file) = @_;
374 :     my ($tbl, $num_pegs);
375 :     my ($key, $entry, $id, $locus, $func, $rest, $contig, $beg, $end, $len, $strand, $type);
376 :    
377 :     open(TBL, "<$file") || die "Could not read-open \'$file\'";
378 :     while (defined($entry = <TBL>))
379 :     {
380 :     next if ($entry =~ m/^\#/);
381 :    
382 :     chomp $entry;
383 :     my @fields = split /\t/, $entry, -1;
384 :     $id = shift @fields;
385 :     $func = pop @fields;
386 :     $locus = pop @fields;
387 :    
388 :     $rest = join(q(,), (grep { $_ } @fields) );
389 :    
390 :     if ((($contig, $beg, $end, $len, $strand) = &from_locus($locus))
391 :     && defined($contig) && $contig
392 :     && defined($beg) && $beg
393 :     && defined($end) && $end
394 :     && defined($len) && $len
395 :     && defined($strand) && $strand
396 :     )
397 :     {
398 :     $key = "$contig\t$strand$end";
399 :    
400 :     $type = q(peg); #...Until such time as RNAs are handled properly...
401 :     # if (($type eq 'peg') || ($type eq 'orf')) {
402 :     ++$num_pegs;
403 :     # }
404 :     # else {
405 :     # warn "Unknown feature type: $entry\n";
406 :     # }
407 :    
408 :     if (not defined($tbl->{$key})) {
409 :     $tbl->{$key} = [ $id, $locus, $contig, $beg, $end, $len, $strand, $type, $entry, ($func || q()), $rest ];
410 :     }
411 :     else {
412 :     warn "Skipping same-STOP TBL entry for $file, $key:\n"
413 :     , "$tbl->{$key}->[ENTRY]\n$entry\n\n";
414 :     }
415 :     }
416 :     else {
417 :     warn "INVALID ENTRY:\t$entry\n";
418 :     }
419 :     }
420 :    
421 :     return ($tbl, $num_pegs);
422 :     }
423 :    
424 :     sub from_locus {
425 :     my ($locus) = @_;
426 :    
427 :     if ($locus) {
428 :     my ($contig, $left, $right, $dir) = SeedUtils::boundaries_of($locus);
429 :     my ($beg, $end) = ($dir eq q(+)) ? ($left, $right) : ($right, $left);
430 :    
431 :     if ($contig && $left && $right && $dir) {
432 :     return ($contig, $beg, $end, (1 + abs($right - $left)), $dir);
433 :     }
434 :     else {
435 :     die "Invalid locus $locus";
436 :     }
437 :     }
438 :     else {
439 :     die "Missing locus";
440 :     }
441 :    
442 :     return ();
443 :     }
444 :    
445 :     sub by_locus {
446 :     my ($x, $a, $b) = @_;
447 :    
448 :     my (undef, undef, $A_contig, $A_beg, $A_end, $A_len, $A_strand) = $x->[$a];
449 :     my (undef, undef, $B_contig, $B_beg, $B_end, $B_len, $B_strand) = $x->[$b];
450 :    
451 :     return ( ($A_contig cmp $B_contig)
452 :     || (&FIG::min($A_beg, $A_end) <=> &FIG::min($B_beg, $B_end))
453 :     || ($B_len <=> $A_len)
454 :     || ($A_strand cmp $B_strand)
455 :     );
456 :     }
457 :     sub by_key {
458 :     my ($x, $y) = @_;
459 :    
460 :     my ($X_contig, $X_strand, $X_end) = ($x =~ m/^([^\t]+)\t([+-])(\d+)/o);
461 :     my ($Y_contig, $Y_strand, $Y_end) = ($y =~ m/^([^\t]+)\t([+-])(\d+)/o);
462 :    
463 :     return ( ($X_contig cmp $Y_contig)
464 :     || ($X_end <=> $Y_end)
465 :     || ($X_strand cmp $Y_strand)
466 :     );
467 :     }
468 :    
469 :    
470 :     sub write_summary {
471 : gdpusch 1.2 my ($old_metrics, $new_metrics,
472 :     $old_pegs, $new_pegs, $num_union,
473 :     $identical, $same_stop, $differ,
474 :     $short, $long, $added, $lost
475 :     ) = @_;
476 :     my $num_common = $same_stop;
477 :    
478 :     my $pct_same_old = 100.0 * $same_stop / $old_pegs;
479 :     my $pct_same_new = 100.0 * $same_stop / $new_pegs;
480 :     my $pct_same_union = 100.0 * $same_stop / $num_union;
481 :    
482 :     my $pct_added_old = 100.0 * $added / $old_pegs;
483 :     my $pct_added_new = 100.0 * $added / $new_pegs;
484 :     my $pct_added_union = 100.0 * $added / $num_union;
485 :    
486 :     my $pct_lost_old = 100.0 * $lost / $old_pegs;
487 :     my $pct_lost_new = 100.0 * $lost / $new_pegs;
488 :     my $pct_lost_union = 100.0 * $lost / $num_union;
489 :    
490 :     my $pct_identical_old = 100.0 * $identical / $old_pegs;
491 :     my $pct_identical_new = 100.0 * $identical / $new_pegs;
492 :     my $pct_identical_common = 100.0 * $identical / $num_common;
493 :     my $pct_identical_union = 100.0 * $identical / $num_union;
494 :    
495 :     my $pct_differ_old = 100.0 * $differ / $old_pegs;
496 :     my $pct_differ_new = 100.0 * $differ / $new_pegs;
497 :     my $pct_differ_common = 100.0 * $differ / $num_common;
498 :     my $pct_differ_union = 100.0 * $differ / $num_union;
499 :    
500 :     my $pct_short_old = 100.0 * $short / $old_pegs;
501 :     my $pct_short_new = 100.0 * $short / $new_pegs;
502 :     my $pct_short_common = 100.0 * $short / $num_common;
503 :     my $pct_short_union = 100.0 * $short / $num_union;
504 :    
505 :     my $pct_long_old = 100.0 * $long / $old_pegs;
506 :     my $pct_long_new = 100.0 * $long / $new_pegs;
507 :     my $pct_long_common = 100.0 * $long / $num_common;
508 :     my $pct_long_union = 100.0 * $long / $num_union;
509 : gdpusch 1.1
510 : gdpusch 1.2 {
511 :     format STDERR =
512 :     genome_size = @####### bp
513 :     $old_metrics->{totlen}
514 :     GC_Content = @##.##%
515 :     $old_metrics->{gc_content}
516 :    
517 :     Num_PEGs Coding_Fraction
518 :     old_num = @##### @##.##%
519 :     $old_pegs, $old_metrics->{coding_fraction}
520 :     new_num = @##### @##.##%
521 :     $new_pegs, $new_metrics->{coding_fraction}
522 :     num_union = @#####
523 :     $num_union
524 :    
525 :     Num %_Old %_New %_Union
526 :     same_stop = @#### @##.##% @##.##% @##.##%
527 :     $same_stop, $pct_same_old, $pct_same_new, $pct_same_union
528 :     added = @#### @##.##% @##.##% @##.##%
529 :     $added, $pct_added_old, $pct_added_new, $pct_added_union
530 :     lost = @#### @##.##% @##.##% @##.##%
531 :     $lost, $pct_lost_old, $pct_lost_new, $pct_lost_union
532 :    
533 :     Num %_Old %_New %_Common %_Union
534 :     identical = @#### @##.##% @##.##% @##.##% @##.##%
535 :     $identical, $pct_identical_old, $pct_identical_new, $pct_identical_common, $pct_identical_union
536 :     differ = @#### @##.##% @##.##% @##.##% @##.##% (Includes Lost and Added)
537 :     $differ, $pct_differ_old, $pct_differ_new, $pct_differ_common, $pct_differ_union
538 :     short = @#### @##.##% @##.##% @##.##% @##.##%
539 :     $short, $pct_short_old, $pct_short_new, $pct_short_common, $pct_short_union
540 :     long = @#### @##.##% @##.##% @##.##% @##.##%
541 :     $long, $pct_long_old, $pct_long_new, $pct_long_common, $pct_long_union
542 :    
543 :     .
544 : gdpusch 1.1
545 : gdpusch 1.2 write STDERR;
546 :     }
547 : gdpusch 1.1
548 :     return 1;
549 :     }
550 :    
551 :    
552 :     sub write_summary_yaml {
553 : gdpusch 1.2 my ($fh, $old_metrics, $new_metrics,
554 :     $old_pegs, $new_pegs, $num_union,
555 :     $identical, $same_stop, $differ,
556 :     $short, $long, $added, $lost
557 :     ) = @_;
558 : gdpusch 1.1
559 :     my $dat = {
560 : gdpusch 1.2 genome_size => $old_metrics->{totlen},
561 :     gc_content => $old_metrics->{gc_content},
562 : gdpusch 1.1 old_num => $old_pegs,
563 :     new_num => $new_pegs,
564 : gdpusch 1.2 num_common => $same_stop,
565 :     num_union => $num_union,
566 :     old_coding_fraction => $old_metrics->{coding_fraction},
567 :     new_coding_fraction => $new_metrics->{coding_fraction},
568 : gdpusch 1.1 };
569 :    
570 : gdpusch 1.2 foreach my $what (qw(same_stop added lost identical differ short long))
571 : gdpusch 1.1 {
572 :     my $val = eval "\$$what";
573 :     $dat->{$what} = $val;
574 : gdpusch 1.2 $dat->{"${what}_pct_old"} = 100 * $val / $old_pegs;
575 :     $dat->{"${what}_pct_new"} = 100 * $val / $new_pegs;
576 :     $dat->{"${what}_pct_union"} = 100 * $val / $num_union;
577 : gdpusch 1.1 }
578 : gdpusch 1.2
579 :     foreach my $what (qw(identical differ short long))
580 : gdpusch 1.1 {
581 :     my $val = eval "\$$what";
582 :     $dat->{"${what}_pct_common"} = 100 * $val / $same_stop;
583 :     }
584 : gdpusch 1.2
585 : gdpusch 1.1 print $fh Dump($dat);
586 : gdpusch 1.2 #print STDERR Dumper($dat);
587 :    
588 : gdpusch 1.1 return 1;
589 :     }
590 :    
591 :     __DATA__

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3