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