[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.4 - (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 : olson 1.2 $results->{num_calls_lost} = 0;
104 :     $results->{num_calls_gained} = 0;
105 :     $results->{num_calls_same} = 0;
106 :     $results->{num_calls_diff} = 0;
107 :    
108 : olson 1.1 #
109 :     # Read first run.
110 :     #
111 :     open(F, "<", $f1) or die "Cannot open $f1: $!";
112 :    
113 :     my %r1;
114 :     my $num_assigned_1 = 0;
115 :     while (<F>)
116 :     {
117 :     chomp;
118 :     my($id, $func) = split(/\t/);
119 :     $func = SeedUtils::strip_func($func);
120 :     $r1{$id} = $func;
121 :     $num_assigned_1++;
122 :     }
123 :    
124 :     $results->{prev_num_assigned} = $num_assigned_1;
125 : olson 1.4
126 : olson 1.1 #
127 :     # Read second run.
128 :     #
129 :     open(F, "<", $f2) or die "Cannot open $f2: $!";
130 :     my $num_assigned_2 = 0;
131 :     while (<F>)
132 :     {
133 :     chomp;
134 :     my($id, $func) = split(/\t/);
135 :     $func = SeedUtils::strip_func($func);
136 :    
137 : olson 1.4 my $pfunc = delete $r1{$id};
138 : olson 1.1
139 :     if ($func eq '' && $pfunc eq '')
140 :     {
141 :     # neither assigned, don't care.
142 :     }
143 :     elsif ($func eq '' && $pfunc ne '')
144 :     {
145 :     # last assigned, we did not
146 :     $results->{num_calls_lost}++;
147 :     push(@{$results->{calls_lost}}, [$id, $pfunc, $func]);
148 :     }
149 :     elsif ($func ne '' && $pfunc eq '')
150 :     {
151 :     # we assigned, last did not
152 :     $results->{num_calls_gained}++;
153 :     push(@{$results->{calls_gained}}, [$id, $pfunc, $func]);
154 :     }
155 :     else
156 :     {
157 :     # we both assigned, are same?
158 :     if ($func eq $pfunc)
159 :     {
160 :     $results->{num_calls_same}++;
161 :     }
162 :     else
163 :     {
164 :     $results->{num_calls_diff}++;
165 :     push(@{$results->{calls_diff}}, [$id, $pfunc, $func]);
166 :     }
167 :     }
168 :    
169 :     $num_assigned_2++;
170 :     }
171 :    
172 : olson 1.4 for my $id (keys %r1)
173 :     {
174 :     # last assigned, we did not
175 :     $results->{num_calls_lost}++;
176 :     my $pfunc = $r1{$id};
177 :     push(@{$results->{calls_lost}}, [$id, $pfunc, undef]);
178 :     }
179 :    
180 : olson 1.1 $results->{this_num_assigned} = $num_assigned_2;
181 :    
182 :     close(F);
183 :    
184 : olson 1.2 #
185 :     # Compute differences in metabolic reconstructions.
186 :     # They are in
187 :     #
188 :     my $meta_old = $f1;
189 :     $meta_old =~ s/\.out$/.reconstruction/;
190 :     my $meta_new = $f2;
191 :     $meta_new =~ s/\.out$/.reconstruction/;
192 :    
193 :     my ($in_rec_old, $ss_old) = &count_lines($meta_old);
194 :     my ($in_rec_new, $ss_new) = &count_lines($meta_new);
195 :    
196 :     $results->{prev_ss_in_reconstruction} = $ss_old;
197 :     $results->{this_ss_in_reconstruction} = $ss_new;
198 :     $results->{ss_in_reconstruction_gained} = int($ss_new - $ss_old);
199 :    
200 :     $results->{prev_pegs_in_reconstruction} = $in_rec_old;
201 :     $results->{this_pegs_in_reconstruction} = $in_rec_new;
202 :     $results->{pegs_in_reconstruction_gained} = int($in_rec_new - $in_rec_old);
203 :    
204 : olson 1.1 open(OUT, ">", "$out_dir/$k.$genome") or die "Cannot write $out_dir/$k.$genome: $!";
205 :     print OUT Dump($results);
206 :     close(OUT);
207 :     }
208 :    
209 : olson 1.2 sub count_lines
210 :     {
211 :     my($f) = @_;
212 :     if (!open(F, "<", $f))
213 :     {
214 :     return (0, 0);
215 :     }
216 :     my $n = 0;
217 :     my %ss;
218 :     while (<F>)
219 :     {
220 :     chomp;
221 :     my($id, $fn, $subsys) = split(/\t/);
222 :     $ss{$subsys}++;
223 :     $n++;
224 :     }
225 :     close(F);
226 :     return $n, int(keys %ss);
227 :     }
228 : olson 1.1
229 :     sub validate
230 :     {
231 :     my($dir) = @_;
232 :    
233 :     -d $dir or die "$dir not a directory\n";
234 :     -d "$dir/tests" or die "$dir has no test data available\n";
235 :    
236 :     my $base = basename($dir);
237 : olson 1.3 if ($base =~ /^Release(\S+)$/)
238 : olson 1.1 {
239 :     return $1;
240 :     }
241 :    
242 :     die "$dir doesn't have the right directory name (should be ReleaseXX)\n";
243 :     }

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3