utils/load_affymetrix_probe_annotations.pl

   1 #!/usr/bin/perl
   2 # load_affymetrix_probe_annotations.pl loads affymetrix probe annotations
   3 # and is released under the terms of the GNU GPL version 3, or any
   4 # later version, at your option. See the file README and COPYING for
   5 # more information.
   6 # Copyright 2013 by Don Armstrong <don@donarmstrong.com>.
   7
   8
   9 use warnings;
  10 use strict;
  11
  12 use Getopt::Long;
  13 use Pod::Usage;
  14
  15 =head1 NAME
  16
  17 load_affymetrix_probe_annotations.pl - loads affymetrix probe annotations
  18
  19 =head1 SYNOPSIS
  20
  21 load_affymetrix_probe_annotations.pl [options] [annotation files]
  22
  23  Options:
  24   --service, -s pgsql service
  25   --progress, -p show progress bar
  26   --debug, -d debugging level (Default 0)
  27   --help, -h display this help
  28   --man, -m display manual
  29
  30 =head1 OPTIONS
  31
  32 =over
  33
  34 =item B<--service,-s>
  35
  36 Postgresql service
  37
  38 =item B<--progress,-p>
  39
  40 Show a progress bar
  41
  42 =item B<--debug, -d>
  43
  44 Debug verbosity. (Default 0)
  45
  46 =item B<--help, -h>
  47
  48 Display brief usage information.
  49
  50 =item B<--man, -m>
  51
  52 Display this manual.
  53
  54 =back
  55
  56 =head1 EXAMPLES
  57
  58 load_affymetrix_probe_annotations.pl
  59
  60 =cut
  61
  62
  63 use vars qw($DEBUG);
  64 use DBI;
  65 use Term::ProgressBar;
  66 use Fcntl qw(:seek);
  67 use Text::CSV;
  68
  69
  70 my %options = (debug           => 0,
  71                help            => 0,
  72                man             => 0,
  73                service         => 'snp',
  74                progress        => 1,
  75               );
  76
  77 GetOptions(\%options,
  78            'service|s=s',
  79            'progress|p!',
  80            'debug|d+','help|h|?','man|m');
  81
  82 pod2usage() if $options{help};
  83 pod2usage({verbose=>2}) if $options{man};
  84
  85 $DEBUG = $options{debug};
  86
  87 my @USAGE_ERRORS;
  88 # if (1) {
  89 #      push @USAGE_ERRORS,"You must pass something";
  90 # }
  91
  92 pod2usage(join("\n",@USAGE_ERRORS)) if @USAGE_ERRORS;
  93
  94 my $dbh = DBI->connect("dbi:Pg:service=$options{service}",
  95                        '','',{AutoCommit => 0}) or
  96     die "Unable to connect to database: ".$DBI::errstr;
  97
  98 my %sth;
  99 $sth{insert_annotation} = $dbh->prepare(<<'END') // die "Unable to prepare insert annotation statement: ".$dbh->errstr;
 100 INSERT INTO affy_annotation
 101 (probe,gene_symbol,
 102 gene_name,species,
 103 array_name,entrez_id,
 104 refseq_prot,refseq_transcript) VALUES ($1,$2,$3,$4,$5,$6,$7,$8);
 105 END
 106
 107 $sth{delete_annotation_id} = $dbh->prepare(<<'END') // die "Unable to prepare delete annotation id statement: ".$dbh->errstr;
 108 DELETE FROM affy_annotation aa WHERE aa.probe = $1;
 109 END
 110
 111 $sth{select_annotation_id} = $dbh->prepare(<<'END') // die "Unable to prepare select annotation id statement: ".$dbh->errstr;
 112 SELECT aa.id FROM affy_annotation aa WHERE aa.probe = $1;
 113 END
 114
 115
 116 $sth{select_affy_probe_id} = $dbh->prepare(<<'END') // die "Unable to prepare select annotation id statement: ".$dbh->errstr;
 117 SELECT id FROM affy_probe ap WHERE ap.probe = $1;
 118 END
 119
 120 $sth{insert_affy_probe_id} = $dbh->prepare(<<'END') // die "Unable to prepare insert affy probe id statement: ".$dbh->errstr;
 121 INSERT INTO affy_probe (probe) VALUES ($1);
 122 END
 123
 124
 125
 126 my @ifh;
 127 for my $ifn (@ARGV) {
 128     my $ifh = IO::File->new($ifn,'r') or
 129         die "Unable to open $ifn for reading: $!";
 130     push @ifh,$ifh;
 131 }
 132
 133 if (not @ARGV) {
 134     push @ifh,\*STDIN;
 135 }
 136
 137 my %header_regex =
 138     (probe => qr/(?i)Probe\s*Set\s*ID/,
 139      gene_symbol => qr/(?i)Gene\s*Symbol/,
 140      gene_name => qr/(?i)Gene\s*Title/,
 141      species => qr/(?i)Species\s*Scientific\s*Name/,
 142      array_name => qr/(?i)GeneChip\s*Array/,
 143      entrez_id => qr/(?i)Entrez\s*Gene/,
 144      refseq_prot => qr/(?i)RefSeq\s*Protein\s*ID/,
 145      refseq_transcript => qr/(?i)RefSeq\s*Transcript\s*ID/,
 146     );
 147
 148 for my $ifh (@ifh) {
 149     my $p;
 150     if ($options{progress}) {
 151         if ($ifh->seek(0,SEEK_END)) {
 152             $p = Term::ProgressBar->new({count => $ifh->tell,
 153                                          remove => 1,
 154                                          ETA=>'linear'});
 155             $ifh->seek(0,SEEK_SET);
 156         }
 157     }
 158     my @header;
 159     my $csv = Text::CSV->new({sep_char=>','});
 160     my %headers;
 161     my %important_headers;
 162     while (<$ifh>) {
 163         chomp;
 164         next if /^#/;
 165         if (not $csv->parse($_)) {
 166             die "Unable to parse line $. of file: ".$csv->error_diag();
 167         }
 168         my @row = $csv->fields();
 169         if (not @header) {
 170             @header = @row;
 171             @headers{@header} = 0..$#row;
 172             for my $header (keys %header_regex) {
 173                 my @match =
 174                     grep { $_ =~ $header_regex{$header}
 175                        } keys %headers;
 176                 $important_headers{$header} = $headers{$match[0]} if @match;
 177                 if (not @match) {
 178                     use Data::Printer;
 179                     p %important_headers;
 180                     p @header;
 181                     p %headers;
 182                     die "unable to find header match for $header";
 183                 }
 184             }
 185             next;
 186         }
 187         insert_annotation($dbh,\%sth,
 188                           {fixup_row(\%important_headers,\@row)
 189                           },
 190                          );
 191         if (defined $p) {
 192             $p->update($ifh->tell);
 193         }
 194     }
 195     $dbh->commit();
 196 }
 197
 198 sub insert_annotation {
 199     my ($dbh,$sth,$annot) = @_;
 200     # find the probe id
 201     $annot->{probe_id} = select_affy_probe_id(@_);
 202     # see if this annotation already exists
 203     return unless defined $annot->{probe_id};
 204     my $annot_id = select_annotation_id(@_);
 205     # if not, insert it
 206     if (not defined $annot_id) {
 207         my $rv = $sth->{insert_annotation}->execute(@{$annot}{(qw(probe_id gene_symbol gene_name species array_name entrez_id),
 208                                                                qw(refseq_prot refseq_transcript),
 209                                                               )})
 210             // die "Unable to execute statement properly: ".$dbh->errstr;
 211         $sth->{insert_annotation}->finish;
 212     } else {
 213         print STDERR "probe: $annot->{probe} is already annotated ($annot_id)\n" if $DEBUG;
 214     }
 215 }
 216
 217 sub select_annotation_id {
 218     my ($dbh,$sth,$annot) = @_;
 219     if (not defined $annot->{probe_id}) {
 220         $annot->{probe_id} = select_affy_probe_id(@_);
 221         return unless defined $annot->{probe_id};
 222     }
 223     my $rv = $sth->{select_annotation_id}->execute($annot->{probe_id}) //
 224         die "Unable to execute statement properly: ".$dbh->errstr;
 225     my ($sample_id) = map {ref $_ ?@{$_}:()}
 226         map {ref $_ ?@{$_}:()} $sth->{select_annotation_id}->fetchall_arrayref([0]);
 227     $sth->{select_annotation_id}->finish;
 228     return ($sample_id);
 229 }
 230
 231 sub select_affy_probe_id {
 232     my ($dbh,$sth,$annot) = @_;
 233     my $rv = $sth->{select_affy_probe_id}->execute($annot->{probe}) //
 234         die "Unable to execute statement properly: ".$dbh->errstr;
 235     my ($probe_id) = map {ref $_ ?@{$_}:()}
 236         map {ref $_ ?@{$_}:()} $sth->{select_affy_probe_id}->fetchall_arrayref([0]);
 237     $sth->{select_affy_probe_id}->finish;
 238     return ($probe_id);
 239 }
 240
 241 sub fixup_row {
 242     my ($ih,$r) = @_;
 243     my %r; # return
 244     for my $h (keys %{$ih}) {
 245         $r{$h} = $r->[$ih->{$h}];
 246     }
 247     return %r;
 248 }
 249
 250
 251 __END__