bin/snp_study_table

   1 #! /usr/bin/perl
   2 # snp_study_table outputs a table of study results, and is released
   3 # under the terms of the GPL version 2, or any later version, at your
   4 # option. See the file README and COPYING for more information.
   5 # Copyright 2011 by Don Armstrong <don@donarmstrong.com>.
   6 # $Id: perl_script 1825 2011-01-02 01:53:43Z don $
   7
   8
   9 use warnings;
  10 use strict;
  11
  12 use Getopt::Long;
  13 use Pod::Usage;
  14
  15 =head1 NAME
  16
  17 snp_study_table - output an xls table of study results
  18
  19 =head1 SYNOPSIS
  20
  21 snp_study_table [options] [snps]
  22
  23  Options:
  24   --output output filename (Default STDOUT)
  25   --debug, -d debugging level (Default 0)
  26   --help, -h display this help
  27   --man, -m display manual
  28
  29 =head1 SNPS
  30
  31 SNPS can be specified as rs12345 or as chr-pos; in the later case, the
  32 ga_chr_pos table will be used, in the former, the ga_snp table will be
  33 used.
  34
  35 =head1 OPTIONS
  36
  37 =over
  38
  39 =item B<--output>
  40
  41 Output filename (defaults to STDOUT)
  42
  43 =item B<--debug, -d>
  44
  45 Debug verbosity. (Default 0)
  46
  47 =item B<--help, -h>
  48
  49 Display brief usage information.
  50
  51 =item B<--man, -m>
  52
  53 Display this manual.
  54
  55 =back
  56
  57 =head1 EXAMPLES
  58
  59
  60 =cut
  61
  62
  63 use vars qw($DEBUG);
  64
  65 use DBI;
  66 use Spreadsheet::WriteExcel;
  67
  68 my %options = (debug           => 0,
  69                help            => 0,
  70                man             => 0,
  71                );
  72
  73 GetOptions(\%options,
  74            'output|o=s',
  75            'debug|d+','help|h|?','man|m');
  76
  77 pod2usage() if $options{help};
  78 pod2usage({verbose=>2}) if $options{man};
  79
  80 $DEBUG = $options{debug};
  81
  82 my @USAGE_ERRORS;
  83
  84 pod2usage(join("\n",@USAGE_ERRORS)) if @USAGE_ERRORS;
  85
  86 my $out_fh = \*STDOUT;
  87 if (exists $options{output}) {
  88     $out_fh = IO::File->new($options{output},'w') or
  89         die "Unable to open $options{output} for writing: $!";
  90 }
  91 my $wb = Spreadsheet::WriteExcel->new($out_fh);
  92 my $ws = $wb->add_worksheet('snp_results');
  93
  94 my $dbh = DBI->connect("dbi:Pg:service=snp",'','',{AutoCommit => 0}) or
  95     die "Unable to connect to database: ".$DBI::errstr;
  96
  97 my $sth_ga_info = $dbh->prepare(<<'END') // die "Unable to prepare ga_snp info statement: ".$dbh->errstr;
  98 SELECT * FROM ga_snp WHERE snp_id=?;
  99 END
 100
 101 my $sth_ga_cp_info = $dbh->prepare(<<'END') // die "Unable to prepare ga_snp info statement: ".$dbh->errstr;
 102 SELECT * FROM ga_chr_pos WHERE chr=? AND pos=?;
 103 END
 104
 105
 106 if (not @ARGV) {
 107     while (<STDIN>) {
 108         next if /^#/;
 109         chomp;
 110         push @ARGV,$_;
 111     }
 112 }
 113
 114 my %studies;
 115 my %snps;
 116 for my $snp (@ARGV) {
 117     my $results;
 118     if ($snp =~ /^(?:rs)(\d+)$/) {
 119         $results = get_snp_results(dbh => $dbh,
 120                                    sth => $sth_ga_info,
 121                                    bind => [$1],
 122                                   );
 123     }
 124     elsif ($snp =~ /^(\d+)-(\d+)$/) {
 125         $results = get_snp_results(dbh => $dbh,
 126                                    sth => $sth_ga_cp_info,
 127                                    bind => [$1,$2],
 128                                   );
 129     }
 130     else {
 131         print STDERR "Invalid snp format '$snp'\n";
 132         next;
 133     }
 134     $snps{$snp} = $results;
 135     for my $result (@{$results}) {
 136         $studies{$result->{study_name}}{$result->{subpart_name}} = 1;
 137     }
 138 }
 139 $dbh->disconnect();
 140
 141 #column names
 142 my @cn = ('A'..'Z','AA'..'ZZ');
 143
 144 my $r = 1;
 145 my $c = 0;
 146 my $f_center = $wb->add_format(align => 'center');
 147 $ws->write($cn[$c++].($r+2),'SNP');
 148 for my $study (sort keys %studies) {
 149     my $n_substudy = keys %{$studies{$study}};
 150     $ws->write($cn[$c].$r,$study);
 151     $ws->merge_range($r-1,$c,$r-1,$c+$n_substudy*3-1,$study,$f_center);
 152     for my $substudy (sort keys %{$studies{$study}}) {
 153         $ws->write($cn[$c].($r+1),$substudy);
 154         $ws->merge_range($r,$c,$r,$c+2,$substudy,$f_center);
 155         $ws->write($cn[$c].($r+2),'P value');
 156         $ws->write($cn[$c+1].($r+2),'Q value');
 157         $ws->write($cn[$c+2].($r+2),'FDR');
 158         $c+=3;
 159     }
 160 }
 161 $r+=3;
 162 $c=0;
 163 for my $snp (sort keys %snps) {
 164     $ws->write($cn[$c++].$r,$snp);
 165     my %snp_studies;
 166     for my $row (@{$snps{$snp}}) {
 167         $snp_studies{$row->{study_name}}{$row->{subpart_name}} = [@{$row}{qw(pvalue qvalue fdr)}];
 168     }
 169     for my $study (sort keys %studies) {
 170         for my $substudy (sort keys %{$studies{$study}}) {
 171             if (not exists $snp_studies{$study} or
 172                 not exists $snp_studies{$study}{$substudy}
 173                ) {
 174                 $c+=3;
 175             }
 176             else {
 177                 for (0..2) {
 178                     $ws->write($cn[$c++].$r,defined $snp_studies{$study}{$substudy}[$_]?$snp_studies{$study}{$substudy}[$_]:'');
 179                 }
 180             }
 181         }
 182     }
 183     $r++;
 184     $c=0;
 185 }
 186
 187 sub get_snp_results{
 188     my %param = @_;
 189     my $rv = $param{sth}->execute(@{$param{bind}}) or
 190         die "Unable to execute statement properly: ".$param{dbh}->errstr;
 191     my $results = $param{sth}->fetchall_arrayref({});
 192     if ($param{sth}->err) {
 193         print STDERR $param{sth}->errstr;
 194         return {};
 195     }
 196     return $results;
 197 }
 198
 199
 200 __END__