skip ENSG results, and use average weight
[function2gene.git] / bin / combine_results
1 #! /usr/bin/perl
2
3 # combine_results, is part of the gene search suite, and is released
4 # under the terms of the GPL version 2, or any later version, at your
5 # option. See the file README and COPYING for more information.
6
7 # Copyright 2006,2007 by Don Armstrong <don@donarmstrong.com>.
8
9
10 use warnings;
11 use strict;
12
13
14 use Getopt::Long;
15 use Pod::Usage;
16
17 =head1 NAME
18
19   combine_results -- combines parsed result files; outputs to stdout.
20
21 =head1 SYNOPSIS
22
23  combine_results parsed_results_1.txt [parsedresultfiles ...]
24
25  Options:
26   --keywords newline delineated list of keywords to search for
27   --results file to store results in
28   --results-table optional file to write result summary table to
29   --debug, -d debugging level [default 0]
30   --help, -h display this help
31   --man, -m display manual
32
33 =head1 OPTIONS
34
35 =over
36
37 =item B<--keywords>
38
39 A file which contains a newline delinated list of keywords to search
40 for. Can be specified multiple times. Lines starting with # or ; are
41 ignored. An optional weight can be specified after the keyword, which
42 is separated from the keyword by a tab. (If not specified, 1 is
43 assumed.)
44
45 =item B<--results>
46
47 A file in which to store the combined results; defaults to STDOUT
48
49 =item B<--results-table>
50
51 A file in which to store a summary table of results. If not provided,
52 no summary table is written.
53
54 =item B<--debug, -d>
55
56 Debug verbosity. (Default 0)
57
58 =item B<--help, -h>
59
60 Display brief useage information.
61
62 =item B<--man, -m>
63
64 Display this manual.
65
66 =back
67
68 =head1 EXAMPLES
69
70   combine_results foo_1.txt
71
72 Will pretty much do what you want
73
74 =cut
75
76
77
78 use vars qw($DEBUG);
79
80 BEGIN{
81      $DEBUG = 0 unless defined $DEBUG;
82 }
83
84 use XML::Parser::Expat;
85 use IO::File;
86 use List::Util qw(sum max);
87
88 # XXX parse config file
89
90 my %options = (debug    => 0,
91                help     => 0,
92                man      => 0,
93                keywords => [],
94                results  => undef,
95                results_table => undef,
96               );
97
98 GetOptions(\%options,'keywords|k=s@',
99            'results=s','results_table|results-table=s',
100            'debug|d+','help|h|?','man|m');
101
102
103 pod2usage() if $options{help};
104 pod2usage({verbose=>2}) if $options{man};
105
106 $DEBUG = $options{debug};
107
108 my $results_fh = \*STDOUT;
109 my $results_table_fh = undef;
110
111 if ($options{results}) {
112      $results_fh = IO::File->new($options{results},'w') or
113          die "Unable to open results file $options{results} for writing";
114 }
115 if ($options{results_table}) {
116      $results_table_fh = IO::File->new($options{results_table},'w') or
117          die "Unable to open results table file $options{results_table} for writing";
118 }
119
120 # CSV columns
121 use constant {NAME        => 0,
122               REFSEQ      => 1,
123               LOCATION    => 2,
124               ALIAS       => 3,
125               FUNCTION    => 4,
126               DESCRIPTION => 5,
127               KEYWORD     => 6,
128               DBNAME      => 7,
129               FILENAME    => 8,
130              };
131
132 my @csv_fields = qw(name hits rzscore weightedscore autoscore refseq location alias database terms description function);
133
134 my %genes;
135
136 my %keyword_weight;
137
138 if (@{$options{keywords}}) {
139      for my $keyword_file (@{$options{keywords}}) {
140           my $keyword_fh = IO::File->new($keyword_file,'r') or
141                die "Unable to open $keyword_file for reading: $!";
142           while (<$keyword_fh>) {
143                next if /^\s*[#;]/;
144                next unless /\w+/;
145                chomp;
146                my ($keyword,$weight) = split /\t/, $_;
147                $weight = 1 if not defined $weight;
148                $keyword_weight{$keyword} = $weight;
149           }
150      }
151 }
152
153
154
155 my %alias_reverse;
156
157 for my $file_name (@ARGV) {
158      my $file = new IO::File $file_name, 'r' or die "Unable to open file $file_name $!";
159      while (<$file>) {
160           next if /^"Name"/;
161           my @gene = map {s/^\"//; s/\"$//; $_;} split /(?<=\")\,(?=\")/, $_;
162           # check to see if there's a different name we should be using
163           if (not exists $genes{$gene[NAME]}) {
164                # if the gene has a valid name, we do at least one test.
165                my $num_tested = $gene[NAME] ne 'NO NAME' ? 1 : 0;
166                my %candidates;
167                if ($gene[NAME] ne 'NO NAME' and exists $alias_reverse{$gene[NAME]} and $alias_reverse{$gene[NAME]} ne '') {
168                     $candidates{$alias_reverse{$gene[NAME]}}++;
169                }
170                else {
171                     for my $alias (grep {$_ !~ /^NO (ALIASES|NAME)$/} split(/; /, $gene[ALIAS])) {
172                          if (exists $alias_reverse{$alias} and $alias_reverse{$alias} ne '') {
173                               $candidates{$alias_reverse{$alias}}++;
174                          }
175                          $num_tested++;
176                     }
177                }
178                #print STDERR "Choosing $alias_reverse{$gene[NAME]} for $gene[NAME]\n";
179                for my $candidate (keys %candidates) {
180                     if ($candidates{$candidate} > ($num_tested/2)) {
181                          print STDERR "Choosing $candidate for '$gene[NAME]', as it matched $candidates{$candidate} of $num_tested tests\n";
182                          $gene[NAME] = $candidate;
183                     }
184                }
185           }
186           $genes{$gene[NAME]}{name} = $gene[NAME];
187           $genes{$gene[NAME]}{database}{$gene[DBNAME]}++;
188           $genes{$gene[NAME]}{hits}++;
189           $genes{$gene[NAME]}{terms}{$gene[KEYWORD]}++;
190           $genes{$gene[NAME]}{terms}{$gene[KEYWORD].'['.$gene[DBNAME].']'}++;
191           add_unique_parts($genes{$gene[NAME]},'refseq',$gene[REFSEQ]);
192           add_if_better($genes{$gene[NAME]},'description',$gene[DESCRIPTION]);
193           add_if_better($genes{$gene[NAME]},'location',$gene[LOCATION]);
194           add_unique_parts($genes{$gene[NAME]},'function',split(/; /, $gene[FUNCTION]));
195           my @aliases = grep {$_ ne 'NO ALIASES'} split(/; /, $gene[ALIAS]);
196           add_unique_parts($genes{$gene[NAME]},'alias', @aliases);
197           if ($gene[NAME] ne 'NO NAME') {
198                for my $alias (@aliases) {
199                     if (not exists $alias_reverse{$alias}) {
200                          $alias_reverse{$alias} = $gene[NAME];
201                     }
202                     elsif ($alias_reverse{$alias} ne $gene[NAME]) {
203                          print STDERR "Alias $alias for $gene[NAME] also points at $alias_reverse{$alias} [".join(',',@aliases).".]\n";
204                          $alias_reverse{$alias} = '';
205                     }
206                }
207           }
208      }
209 }
210
211 my %databases;
212 my %terms;
213 my %auto_weight;
214 my %keyword_keyword;
215 for my $gene (keys %genes) {
216      my %term_temp;
217      my %db_temp;
218      my %gene_temp;
219      my %gene_temp2;
220      for my $term (keys %{$genes{$gene}{terms}}) {
221           if ($term =~ /\[/) {
222                my ($keyword,$database) = $term =~ /([^[]+)\[([^\]]+)\]/;
223                my $hits = $genes{$gene}{terms}{$term};
224                $keyword =~ s/[-+_]/ /g;
225                $keyword =~ s/\s*$//;
226                $keyword =~ s/[*]//;
227                $gene_temp{$keyword}{$database} = 1;
228                $gene_temp2{$database}{$keyword} = 1;
229                $databases{$database}{$keyword}{count}++;
230                $db_temp{$database}++;
231                $terms{$keyword}{$database}{count}++;
232           }
233           else {
234                my $keyword = $term;
235                my $hits = $genes{$gene}{terms}{$term};
236                $keyword =~ s/[-+_]/ /g;
237                $keyword =~ s/\s*$//;
238                $keyword =~ s/[*]//;
239                $terms{$keyword}{total}{count}++;
240           }
241      }
242      if (keys %gene_temp == 1) {
243           $terms{[keys %gene_temp]->[0]}{total}{unique}++;
244           if (keys %{$gene_temp{[keys %gene_temp]->[0]}} == 1) {
245                $databases{total}{total}{unique}++
246           }
247      }
248      if (keys %gene_temp2 == 1) {
249           $databases{[keys %gene_temp2]->[0]}{total}{unique}++;
250      }
251      for my $keyword (keys %gene_temp) {
252           if (keys %{$gene_temp{$keyword}} == 1) {
253                $terms{$keyword}{[keys %{$gene_temp{$keyword}}]->[0]}{unique}++;
254           }
255           for my $keyword2 (keys %gene_temp) {
256                $keyword_keyword{$keyword}{$keyword2}++
257           }
258      }
259      for my $database (keys %db_temp) {
260           $databases{$database}{total}{count}++;
261      }
262      $databases{total}{total}{count}++;
263
264 }
265
266 for my $keyword (keys %keyword_keyword) {
267      # the autoweight table is the diagonal over the sum of the column of the keyword/keyword table
268      # we use max here to avoid 0/0 problems.
269      my $results_by_this_keyword = max(1,$keyword_keyword{$keyword}{$keyword});
270      my $results_combined = max(1,grep {defined $_}
271                                 sum(map {$keyword_keyword{$keyword}{$_}}
272                                     grep {$_ ne $keyword}
273                                     keys %{$keyword_keyword{$keyword}}
274                                    )
275                                );
276      $auto_weight{$keyword} = $results_by_this_keyword/$results_combined;
277 }
278
279 my $avg_weight = sum(values %auto_weight) / scalar keys %auto_weight;
280 for my $keyword (keys %auto_weight) {
281      $auto_weight{$keyword} = $auto_weight{$keyword}/$avg_weight;
282 }
283
284 print {$results_fh} join(',',map {qq("$_")} @csv_fields),qq(\n);
285 for my $gene (keys %genes) {
286      $genes{$gene}{rzscore} = scalar grep {$_ !~ /\[/} keys %{$genes{$gene}{terms}};
287      $genes{$gene}{weightedscore}= sum(0,
288                                        map {defined $keyword_weight{$_}?$keyword_weight{$_}:1}
289                                        grep {$_ !~ /\[/} keys %{$genes{$gene}{terms}}
290                                       );
291      $genes{$gene}{autoscore}= sum(0,
292                                    map {defined $auto_weight{$_}?$auto_weight{$_}:1}
293                                    grep {$_ !~ /\[/} keys %{$genes{$gene}{terms}}
294                                   );
295 }
296
297 my $sort = 'autoscore';
298 if (scalar grep {$_ != 1 } values %keyword_weight) {
299      $sort='weightedscore';
300 }
301 for my $gene (sort {$genes{$b}{$sort} <=> $genes{$a}{$sort}} keys %genes) {
302      print {$results_fh} join (',',
303                         map {s/"//g; qq("$_")}
304                         map {
305                              my $value = $_;
306                              if (ref $value eq 'HASH') {
307                                   join('; ',map {qq($_:$$value{$_})} keys %$value);
308                              }
309                              elsif (ref $value eq 'ARRAY') {
310                                   join('; ', @$value);
311                              }
312                              else {
313                                   $value;
314                              }
315                         } @{$genes{$gene}}{@csv_fields}
316                        ), qq(\n);
317 }
318
319
320 sub add_unique_parts{
321      my ($hr,$key,@values) = @_;
322      if (not defined $$hr{key}) {
323           $$hr{$key} = [@values];
324      }
325      else {
326           return unless @values;
327           my %temp_hash;
328           @temp_hash{@{$$hr{$key}}} = (1) x scalar @{$$hr{$key}};
329           $temp_hash{@values} = (1) x scalar @values;
330           $$hr{$key} = [keys %temp_hash];
331      }
332 }
333
334 sub add_if_better{
335      my ($hash, $key, $value) = @_;
336      if (not defined $$hash{$key}) {
337           $$hash{$key} = $value;
338      }
339      elsif (length $$hash{$key} < length $value and $value !~ /^NO\s+(LOCATION|DESCRIPTION)+$/) {
340           $$hash{$key} = $value;
341      }
342 }
343
344 sub space_fill{
345      my ($value,$length,$right) = @_;
346      $right ||= 0;
347      if (length($value) > $length) {
348           $value =~ m/(.{$length})/;
349           return $1;
350      }
351      if (length($value) == $length) {
352           return $value
353      }
354      if ($right) {
355           return join('',
356                       ' ' x ($length - length($value)),
357                       $value,
358                      );
359      }
360      else {
361           return join('',
362                       $value,
363                       ' ' x ($length - length($value)),
364                      );
365      }
366 }
367
368 sub results_table_line {
369      my ($keyword,@fields) = @_;
370      return join( ' & ',
371                   space_fill($keyword,23),
372                   map {space_fill($_,11,1)} @fields
373                 )."\n";
374 }
375
376 my @database_order = grep {lc($_) ne 'total'} keys %databases;
377 if (defined $results_table_fh) {
378      my $keyword;
379      print {$results_table_fh} results_table_line('Keyword','Weight','Autoweight',
380                                                   map {ucfirst($_)} @database_order,
381                                                   'Total',
382                                                  );
383
384      for $keyword (sort keys %terms) {
385           my @fields =
386                map {
387                     if (not defined $_) {
388                          '$-$';
389                     }
390                     else {
391                          $_->{unique} ||= 0;
392                          "$_->{count} ($_->{unique})";
393                     }
394                } @{$terms{$keyword}}{@database_order,'total'};
395           unshift @fields, $auto_weight{$keyword};
396           unshift @fields, $keyword_weight{$keyword} || 1;
397           print {$results_table_fh} results_table_line($keyword,
398                                                        @fields
399                                                       );
400
401      }
402      $keyword = 'Total';
403      my @fields = ('','');
404      push @fields,
405           map {
406                if (not defined $_) {
407                     '$-$';
408                }
409                else {
410                     $_->{unique} ||= 0;
411                     "$_->{count} ($_->{unique})";
412                }
413           } map {$_->{total}} @databases{@database_order,'total'};
414      print {$results_table_fh} results_table_line($keyword,
415                                                   @fields
416                                                  );
417 }
418
419 __END__