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 $max_weight = max(values %auto_weight);
 280 for my $keyword (keys %auto_weight) {
 281      $auto_weight{$keyword} = $auto_weight{$keyword}/$max_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__