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

Annotation of /FigKernelScripts/FFB2_compare_tests.pl

Parent Directory Parent Directory | Revision Log Revision Log


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

1 : olson 1.1 #
2 :     # Compare the results of two tested FIGfam dirs, writing the results into the
3 :     # newer one.
4 :     #
5 :    
6 :     use Data::Dumper;
7 :     use YAML::Any;
8 :     use strict;
9 :     use File::Basename;
10 :     use SeedUtils;
11 :    
12 :     @ARGV == 2 or die "Usage: FFB2_compare_tests FFdir1 FFdir2\n";
13 :    
14 :     my $ff_dir_1 = shift;
15 :     my $ff_dir_2 = shift;
16 :    
17 :     my $ff_rel_1 = validate($ff_dir_1);
18 :     my $ff_rel_2 = validate($ff_dir_2);
19 :    
20 :     print "Comparing release $ff_rel_1 ($ff_dir_1) to release $ff_rel_2 ($ff_dir_2)\n";
21 :    
22 :     my $out_dir = "$ff_dir_2/comparisons/$ff_rel_1";
23 :     &SeedUtils::verify_dir($out_dir);
24 :    
25 :     #
26 :     # Find the set of kmers & genomes in common between the two test runs.
27 :     #
28 :    
29 :     my %all_genomes;
30 :     my @to_compare;
31 :    
32 :     for my $f1kdir (<$ff_dir_1/tests/*>)
33 :     {
34 :     my $k = basename($f1kdir);
35 :     my $f2kdir = "$ff_dir_2/tests/$k";
36 :     if (! -d $f2kdir)
37 :     {
38 :     warn "Skipping $k, not present in release $ff_rel_2\n";
39 :     next;
40 :     }
41 :    
42 :     for my $outfile (<$f1kdir/*.out>)
43 :     {
44 :     my $outbase = basename($outfile);
45 :    
46 :     if (! -f "$f2kdir/$outbase")
47 :     {
48 :     warn "Skipping $outbase, not present in release $ff_rel_2\n";
49 :     next;
50 :     }
51 :    
52 :     my ($genome) = $outbase =~ /^(\d+\.\d+)/;
53 :     $all_genomes{$genome}++;
54 :     push(@to_compare, [$k, $genome, $outfile, "$f2kdir/$outbase"]);
55 :     }
56 :     }
57 :    
58 :     #
59 :     # Find genome peg counts (common data).
60 :     #
61 :    
62 :     my %genome_peg_count;
63 :     for my $genome (keys %all_genomes)
64 :     {
65 :     open(F, "<", "/vol/figfam-prod/test_genomes/$genome.peg.fasta") or die "cannot open /vol/figfam-prod/test_genomes/$genome.peg.fasta: $!";
66 :     my $c = 0;
67 :     while (<F>)
68 :     {
69 :     /^>/ and $c++;
70 :     }
71 :     close(F);
72 :     $genome_peg_count{$genome} = $c;
73 :     }
74 :    
75 :     for my $comp (sort { $a->[0] <=> $b->[0] or $a->[1] <=> $b->[1] } @to_compare)
76 :     {
77 :     my($k, $genome, $f1, $f2) = @$comp;
78 :    
79 :     print "Compare @$comp\n";
80 :    
81 :     &compare($k, $genome, $f1, $f2);
82 :     }
83 :    
84 :     #
85 :     # Compare two runs of the kmers.
86 :     #
87 :     # Compute a bunch of metrics, put in data structure, and
88 :     # write the YAML version into $ff_dir_2/comparisons/$ff_rel_1/$k.$genome.
89 :     #
90 :     sub compare
91 :     {
92 :     my($k, $genome, $f1, $f2) = @_;
93 :    
94 :     my $results = {};
95 :    
96 :     $results->{peg_count} = $genome_peg_count{$genome};
97 :     $results->{kmer} = $k;
98 :     $results->{genome} = $genome;
99 :     $results->{this_version} = $ff_rel_2;
100 :     $results->{prev_version} = $ff_rel_1;
101 :     $results->{this_data} = $f2;
102 :     $results->{prev_data} = $f1;
103 :    
104 :     #
105 :     # Read first run.
106 :     #
107 :     open(F, "<", $f1) or die "Cannot open $f1: $!";
108 :    
109 :     my %r1;
110 :     my $num_assigned_1 = 0;
111 :     while (<F>)
112 :     {
113 :     chomp;
114 :     my($id, $func) = split(/\t/);
115 :     $func = SeedUtils::strip_func($func);
116 :     $r1{$id} = $func;
117 :     $num_assigned_1++;
118 :     }
119 :    
120 :     $results->{prev_num_assigned} = $num_assigned_1;
121 :    
122 :     #
123 :     # Read second run.
124 :     #
125 :     open(F, "<", $f2) or die "Cannot open $f2: $!";
126 :     my $num_assigned_2 = 0;
127 :     while (<F>)
128 :     {
129 :     chomp;
130 :     my($id, $func) = split(/\t/);
131 :     $func = SeedUtils::strip_func($func);
132 :    
133 :     my $pfunc = $r1{$id};
134 :    
135 :     if ($func eq '' && $pfunc eq '')
136 :     {
137 :     # neither assigned, don't care.
138 :     }
139 :     elsif ($func eq '' && $pfunc ne '')
140 :     {
141 :     # last assigned, we did not
142 :     $results->{num_calls_lost}++;
143 :     push(@{$results->{calls_lost}}, [$id, $pfunc, $func]);
144 :     }
145 :     elsif ($func ne '' && $pfunc eq '')
146 :     {
147 :     # we assigned, last did not
148 :     $results->{num_calls_gained}++;
149 :     push(@{$results->{calls_gained}}, [$id, $pfunc, $func]);
150 :     }
151 :     else
152 :     {
153 :     # we both assigned, are same?
154 :     if ($func eq $pfunc)
155 :     {
156 :     $results->{num_calls_same}++;
157 :     }
158 :     else
159 :     {
160 :     $results->{num_calls_diff}++;
161 :     push(@{$results->{calls_diff}}, [$id, $pfunc, $func]);
162 :     }
163 :     }
164 :    
165 :     $num_assigned_2++;
166 :     }
167 :    
168 :     $results->{this_num_assigned} = $num_assigned_2;
169 :    
170 :     close(F);
171 :    
172 :     open(OUT, ">", "$out_dir/$k.$genome") or die "Cannot write $out_dir/$k.$genome: $!";
173 :     print OUT Dump($results);
174 :     close(OUT);
175 :     }
176 :    
177 :    
178 :     sub validate
179 :     {
180 :     my($dir) = @_;
181 :    
182 :     -d $dir or die "$dir not a directory\n";
183 :     -d "$dir/tests" or die "$dir has no test data available\n";
184 :    
185 :     my $base = basename($dir);
186 :     if ($base =~ /^Release(\d+)$/)
187 :     {
188 :     return $1;
189 :     }
190 :    
191 :     die "$dir doesn't have the right directory name (should be ReleaseXX)\n";
192 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3