scripts/DB-HIV/hivq.PLS

   1 #!perl
   2 #$Id: hivq.PL 362 2009-01-23 22:45:51Z maj $#
   3
   4 # command-line script for HIV sequence queries
   5 # using HIV.pm and HIVQuery.pm
   6
   7 =head1 NAME
   8
   9 hivq.PL - an interactive command-line interface to L<Bio::DB::HIV> and L<Bio::DB::Query::HIVQuery>
  10
  11 =head1 SYNOPSIS
  12
  13  $ perl hivq.PL
  14  hivq> query C[subtype] SI[phenotype]
  15  hivq> prerun
  16  80 sequences returned
  17  Query: C[subtype] SI[phenotype]
  18  hivq> outfile csi.fas
  19  hivq> run
  20  Download complete.
  21  hivq> outfile dsi.fas
  22  hivq> run D[subtype] SI[phenotype]
  23  Download complete.
  24  hivq> count
  25  25 sequences returned
  26  Query: D[subtype] SI[phenotype]
  27  hivq> exit
  28  $
  29
  30 =head1 DESCRIPTION
  31
  32 The BioPerl modules L<Bio::DB::HIV> and L<Bio::DB::Query::HIVQuery>
  33 together allow batch queries against the Los Alamos National
  34 Laboratories' HIV Sequence Database using a simple query
  35 language. C<hivq.PL> provides both an example of the use of these
  36 modules, and a standalone interactive command-line interface to the
  37 LANL HIV DB. Simple commands allow the user to retrieve HIV sequences
  38 and annotations using the query language implemented in
  39 L<Bio::DB::Query::HIVQuery>. Visit the man pages for those modules for
  40 more details.
  41
  42 =head1 USAGE
  43
  44 Run the script using C<perl hivq.PL> or, in Unix, C<./hivq.PL>. You will see
  45 the
  46
  47  hivq>
  48
  49 prompt. Type commands with queries to retrieve sequence and annotation
  50 data.  See the SYNOPSIS for a sample session. Available commands
  51 are described below.
  52
  53 =head2 TIPS
  54
  55 The LANL database is pretty complex and extensive. Use the C<find> facility to
  56 explore the available database tables and fields. To identify aliases for a particular field, use C<find alias [fieldname]>. For example, to find a short alias to the
  57 weirdly named field C<seq_sample.ssam_second_receptor>, do
  58
  59  hivq> find alias seq_sample.ssam_second_receptor
  60
  61 which returns
  62
  63  coreceptor             second_receptor
  64
  65 Now, instead of the following query
  66
  67  hivq> run C[subtype] CCR5[seq_sample.ssam_second_receptor]
  68
  69 you know you can do
  70
  71  hivq> run C[subtype] CCR5[coreceptor]
  72
  73 Use the C<outfile> command to set the file that receives the retrieved
  74 sequences. You can change the current output file simply by issuing a
  75 new C<outfile> command during the session. The output file defaults to
  76 standard output.
  77
  78 Use the C<query> command to validate a query without hitting the
  79 database. Use the C<prerun> or C<count> commands to get a count of
  80 sequence hits for a query without retrieving the data. Use C<run> or
  81 C<do> to perform a complete query, retrieving sequence data into the
  82 currently set output files.
  83
  84 To process C<hivq.PL> commands in batch, create a text file
  85 (C<hivq.cmd>, for example) containing desired commands one per
  86 line. Then execute the following from the shell:
  87
  88  $ cat hivq.cmd | perl hivq.PL
  89
  90 =head1 COMMANDS
  91
  92 Here is a complete list of commands. Options in single brackets (C<[req_option]>) are
  93 required; options in double brackets (C<[[opt_option]]>) are optional.
  94
  95  confirm            : Toggle interactive confirmation before
  96                       executing queries
  97  exit               : Quit script
  98  find               : Explore database schema
  99   find tables                 Display all database tables
 100   find fields                 Display all database fields (columns)
 101   find fields [table]         Display all fields in [table]
 102   find alias [field]          Display valid aliases for [field]
 103  help [[command]]   : Show command help
 104                       if [[command]] not specified, list all
 105                       available commands
 106  id                 : Display current session id
 107  outfile [filename] : Set file for collecting retrieved data
 108  ping               : Check if LANL DB is available
 109  prerun [[query]]   : Execute query but retreive hit count only
 110                       if [[query]] not specified, use current query
 111  query [query]      : Validate and set current query
 112  run [[query]]      : Execute query and retrieve data
 113                       if [[query]] not specified, use current query
 114  state              : Display current state of the script
 115
 116  bye                : Alias for 'exit'
 117  config             : Alias for 'state'
 118  count              : Alias for 'prerun'
 119  do                 : Alias for 'run'
 120  out                : Alias for 'outfile'
 121  quit               : Alias for 'exit'
 122
 123 =head1 OPTIONS
 124
 125  -v : verbose; turns on the internal debug() function
 126
 127 =head1 FEEDBACK
 128
 129 =head2 Mailing Lists
 130
 131 User feedback is an integral part of the evolution of this and other
 132 Bioperl modules. Send your comments and suggestions preferably to
 133 the Bioperl mailing list.  Your participation is much appreciated.
 134
 135   bioperl-l@bioperl.org                  - General discussion
 136   http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
 137
 138 =head2 Reporting Bugs
 139
 140 Report bugs to the Bioperl bug tracking system to help us keep track
 141 of the bugs and their resolution. Bug reports can be submitted via the
 142 web:
 143
 144   https://redmine.open-bio.org/projects/bioperl/
 145
 146 =head1 AUTHOR - Mark A. Jensen
 147
 148 Mark A. Jensen  E<lt>maj@fortinbras.usE<gt>
 149
 150 =cut
 151
 152 $|=1;
 153 use strict;
 154 use Term::ReadLine;
 155 use Bio::DB::HIV;
 156 use Bio::DB::Query::HIVQuery;
 157 use Bio::SeqIO;
 158 use Error qw(:try);
 159
 160 use Getopt::Std;
 161 our ($opt_v);
 162 getopts('v');
 163 my $db = new Bio::DB::HIV;
 164
 165 my $q = new Bio::DB::Query::HIVQuery();
 166 my $schema = $q->_schema;
 167 my $seqio;
 168 my @cmds = qw(
 169 bye
 170 config
 171 confirm
 172 count
 173 do
 174 exit
 175 find
 176 help
 177 id
 178 out
 179 outfile
 180 ping
 181 prerun
 182 query
 183 quit
 184 run
 185 state
 186 );
 187
 188 my %state = (
 189     'confirm'    => 0,
 190     'outfile'    => {'text'=>'[stdout]','handle'=>\*STDOUT},
 191     'query'      => '',
 192     'session_id' => '',
 193     'timeout'    => '',
 194     'curct'      => 0
 195 );
 196
 197 my %help = (
 198     'bye' => 'bye: End script',
 199     'config' => 'config: Print current configuration state of hivq',
 200     'confirm' => 'confirm: Toggle user confirmation of query execution',
 201     'count' => 'count [[query string]]: Get number of sequences available for query',
 202     'do' => 'do [[query string]]: Run query to completion',
 203     'exit' => 'exit: End script',
 204     'find' => join("\n","find tables: Print all LANL table names",
 205                "find fields: Print all field names",
 206                "find fields [tablename]: Print all fields in specified table",
 207                "find aliases [fieldname]: Print all valid aliases for specified field",
 208                    "find options [fieldname/alias]: Print valid match options for specified field"),
 209     'help' => 'help [[command]]: Get description of command',
 210     'id' => 'id: Print current session id',
 211     'out' => 'out [output filename]: Set filename to hold returned sequences',
 212     'outfile' => 'outfile [output filename]: Set filename to hold returned sequences',
 213     'ping' => 'ping: Ping LANL database',
 214     'prerun' => 'prerun [[query string]]: Get number of sequences available',
 215     'query' => 'query [[query string]]: Validate and set current query string',
 216     'quit' => 'quit: End script',
 217     'run' => 'run [[query string]]: Run query to completion',
 218     'state' => 'state: Print current configuration state of hivq',
 219     );
 220 my $done = 0;
 221 debug("We have arrived.");
 222 #### main loop
 223 my ($cmd, @tok);
 224 my $prompt = "hivq> ";
 225 my $term = new Term::ReadLine 'hivq';
 226 CMD:
 227 while ( !$done ) {
 228     @tok=();$cmd='';
 229     for ( $cmd = $term->readline($prompt) ) {
 230         chomp;
 231         s/^\s+//; #trim whsp
 232         $term->addhistory($_);
 233         @tok = split(/\s+/);
 234         $tok[0] = lc $tok[0]; # normalize
 235         for ($tok[0]) {
 236             (!$_ || /^\#/) && do {
 237                 next CMD;
 238             };
 239             !(grep(/^$tok[0]$/, @cmds)) && do {
 240                 error("Unrecognized command \'$tok[0]\'");
 241                 next CMD;
 242             };
 243             (m/^quit$/ || m/^bye$/ || m/^exit$/) && do {
 244                 my $fh = $state{outfile}->{handle};
 245                 unless ($fh == \*STDOUT) {
 246                     close $fh;
 247                 }
 248                 $done=1;
 249                 next CMD;
 250             };
 251             m/^confirm$/ && do {
 252                 $state{confirm} = !$state{confirm};
 253                 print "User confirmation before running query: ". ($state{confirm} ? "yes" : "no")."\n";
 254                 next CMD;
 255             };
 256             (m/^outfile$/ || m/^out$/) && do {
 257                 unless ($tok[1]) {
 258                     error('Output filename not specified');
 259                     next CMD;
 260                 }
 261                 my $fh;
 262                 eval {
 263                     open $fh, ">$tok[1]" or die $!;
 264                 };
 265                 if ($@) {
 266                     error($@);
 267                 }
 268                 else {
 269                     my $oldfh = $state{outfile}->{handle};
 270                     if ($state{outfile}->{text} ne '[stdout]') {
 271                         close($oldfh);
 272                     }
 273                     $state{outfile}->{handle} = $fh;
 274                     $state{outfile}->{text} = $tok[1];
 275                 }
 276                 next CMD;
 277             };
 278             m/^query$/ && do {
 279                 my $query = join(' ', @tok[1..$#tok]);
 280                 if ($query) {
 281                     $q->_reset;
 282                     $q->query($query);
 283                 }
 284                 else {
 285                     if ($state{query}) {
 286                         print "Current query: $state{query}\n";
 287                     }
 288                     else {
 289                         print "No query set\n";
 290                     }
 291                     next CMD;
 292                 }
 293                 try { #catch errors
 294                     $q->_do_query(0);
 295                 }
 296                 catch Error with {
 297                     my $E = shift;
 298                     error($E->text);
 299                     next CMD;
 300                 };
 301                 $state{session_id} = $q->_session_id;
 302                 $state{query} = $query;
 303                 next CMD;
 304             };
 305             (m/^prerun$/ || m/^count$/) && do {
 306                 my $query = join(' ', @tok[1..$#tok]);
 307                 if (!$query) {
 308                     if (!$q->query()) {
 309                         error('No query currently set');
 310                         next CMD;
 311                     }
 312                     elsif ($q->{'_RUN_LEVEL'} >= 1 &&
 313                         ($q->query() eq $state{query})) {
 314                         outputPrerun();
 315                         next CMD;
 316                     }
 317                     # else, fallthrough to prerun current query
 318                 }
 319                 else { # new query
 320                     $q->_reset;
 321                     $q->query($query);
 322                     try {
 323                         $q->_do_query(0);
 324                     }
 325                     catch Error with {
 326                         my $E = shift;
 327                         error($E->text);
 328                         next CMD;
 329                     };
 330                 }
 331                 # on success:
 332                 $state{confirm} && do { next CMD unless getConfirm(); };
 333                 $state{session_id} = $q->_session_id;
 334                 $state{query} = $query;
 335                 try {
 336                     $q->_do_query(1);
 337                 }
 338                 catch Error with {
 339                     my $E = shift;
 340                     error($E->text);
 341                     next CMD;
 342                 };
 343                 # on success:
 344                 $state{curct} = $q->count;
 345                 $state{query} = $q->query;
 346                 outputPrerun();
 347                 next CMD;
 348
 349             };
 350             (m/^run$/ || m/^do$/) && do {
 351                 my $query = join(' ', @tok[1..$#tok]);
 352                 if (!$query) {
 353                     unless (defined $q->{'_RUN_LEVEL'} && $q->{'_RUN_LEVEL'} < 2) {
 354                         error('No query currently set');
 355                         next CMD;
 356                     }
 357                 }
 358                 else { # new query
 359                     $q->_reset;
 360                     $q->query($query);
 361                     try {
 362                         $q->_do_query(0);
 363                     }
 364                     catch Error with {
 365                         my $E = shift;
 366                         error($E->text);
 367                         next CMD;
 368                     };
 369                     $state{query} = $query;
 370                     $state{curct} = 0;
 371                     $state{session_id} = $q->_session_id;
 372                 }
 373                 # on success:
 374                 $state{confirm} && do { next CMD unless getConfirm(); };
 375                 try {
 376                     $q->_do_query(2);
 377                     $seqio = $db->get_Stream_by_query($q);
 378                 }
 379                 catch Error with {
 380                     my $E = shift;
 381                     error($E->text);
 382                     next CMD;
 383                 };
 384                 # on success:
 385                 $state{curct} = $q->count;
 386                 # output the seqs
 387                 if ($q->count) {
 388                     outputSeqs($seqio);
 389                 }
 390                 else {
 391                     error('No sequences returned');
 392                 }
 393                 next CMD;
 394             };
 395             m/^find$/ && do {
 396                 outputFind(@tok);
 397                 next CMD;
 398             };
 399             m/^help$/ && do {
 400                 unless ($tok[1]) {
 401                     outputUsage();
 402                     next CMD;
 403                 }
 404                 if (grep(/$tok[1]/, @cmds)) {
 405                     print $help{$tok[1]}, "\n";
 406                 }
 407                 else {
 408                     error("Command \'$tok[1]\' unrecognized\n");
 409                 }
 410                 next CMD;
 411             };
 412             (m/^state$/ || m/^config$/) && do {
 413                 foreach my $k (sort keys %state) {
 414                     if (ref($state{$k}) eq 'HASH') {
 415                         print "$k: ".$state{$k}->{text}."\n";
 416                     }
 417                     else {
 418                         print "$k: ";
 419                         if ($k eq 'confirm') {
 420                             print $state{$k} ? "yes" : "no";
 421                             print "\n";
 422                         }
 423                         else {
 424                             print "$state{$k}\n";
 425                         }
 426                     }
 427                 }
 428                 next CMD;
 429             };
 430             do {
 431                 error("Command currently unimplemented\n");
 432                 next CMD;
 433             };
 434         }
 435     }
 436 }
 437
 438 ### end main
 439
 440 sub error {
 441     my $msg = shift;
 442     if ($msg =~ /MSG:/) {
 443         ($msg) = grep (/^MSG:/, split(/\n|\r/,$msg));
 444         $msg =~ s/^MSG: *//;
 445         $msg =~ s/\sat\s.*$//;
 446     }
 447     print STDERR "hivq: $msg\n";
 448     return 1;
 449 }
 450
 451 sub outputPrerun {
 452     print (($state{curct} ? $state{curct} : "No")
 453         . " sequence"
 454         . ($state{curct}>1 ? "s" : "")
 455         . " returned\n");
 456     print "Query: ".$state{query}."\n";
 457     return 1;
 458 }
 459
 460 sub outputSeqs {
 461     my ($seqio) = @_;
 462     my @ret;
 463     while (my $seq = $seqio->next_seq) {
 464         my $nameline = ">";
 465         $nameline .= $seq->annotation->get_value('Special', 'accession').
 466             '('.$seq->id.') ';
 467         # loop through categories:
 468         foreach my $cat ($seq->annotation->get_keys()) {
 469             foreach my $an ($seq->annotation->get_keys($cat)) {
 470                 next if ($an eq 'accession');
 471                 my $value = $seq->annotation->get_value($cat, $an);
 472                 # next line: kludge to skip if there's an annotation
 473                 # object instead of a value (I believe this is a bug)
 474                 next if ref($value);
 475                 $nameline .= "\t".join('=', "'$an'", "'".$value."'");
 476             }
 477         }
 478         push @ret, $nameline."\n";
 479         push @ret, $seq->seq()."\n";
 480     }
 481     doOutputSeqs(@ret);
 482 }
 483
 484 sub doOutputSeqs {
 485     my @lines = @_;
 486     if (!@lines) {
 487         error("No sequences to output");
 488     }
 489     my $fh = $state{outfile}->{handle};
 490     print $fh @lines;
 491     print "Download complete\n";
 492     return 1;
 493 }
 494
 495 sub getConfirm {
 496     print "Run query? [y/n]";
 497     my $ans = <STDIN>;
 498     ($ans =~ /^[yY]/) && do {return 1;};
 499     return 0;
 500 }
 501
 502 sub outputUsage {
 503     print "Available commands:\n";
 504     outputInColumns(\@cmds);
 505     return 1;
 506 }
 507
 508 sub outputFind {
 509     my @tok = @_;
 510     for my $arg ($tok[1]) {
 511         !$arg && do {
 512             print $help{find},"\n";
 513             return;
 514         };
 515         ($arg =~ m/^tables$/) && do {
 516             outputInColumns([$schema->tables]);
 517             return;
 518         };
 519         ($arg =~ m/^fields$/) && do {
 520             my $tbl = $tok[2];
 521             if (!$tbl) {
 522                 outputInColumns([$schema->fields]);
 523             }
 524             else {
 525                 unless (grep /^$tbl$/, $schema->tables) {
 526                     error("Table \'$tbl\' not valid");
 527                     return;
 528                 }
 529                 outputInColumns([grep /^$tbl/, $schema->fields()]);
 530             }
 531             return;
 532         };
 533         ($arg =~ /^options$/) && do {
 534             my $fld = $tok[2];
 535             my @aliases = $schema->aliases($fld);
 536             if (!@aliases) {
 537                 unless (grep /^$fld$/, $schema->aliases) {
 538                     error("Field \'$tok[2]\' not valid");
 539                     return;
 540                 }
 541                 foreach ($schema->fields) {
 542                     if ( grep /$fld/, $schema->aliases($_) ) {
 543                         $fld = $_;
 544                         last;
 545                     }
 546                 }
 547                 # on success: $fld is set to valid field
 548             }
 549             my @options = sort {$a cmp $b}  $schema->options($fld);
 550             @options = (@options ? @options : ('Free text'));
 551             outputInColumns(\@options);
 552             return;
 553         };
 554         do {
 555             ($arg =~ /^alias/) && ($arg = $tok[2]);
 556             if (grep /$arg/, $schema->fields) {
 557                 my @aliases = sort $schema->aliases($arg);
 558                 if (@aliases) {
 559                     outputInColumns(\@aliases);
 560                 }
 561                 else {
 562                     error("No aliases to field \'$arg\'");
 563                 }
 564             }
 565             else {
 566                 error("Field \'$arg\' not valid");
 567             }
 568             return;
 569         };
 570     }
 571 }
 572
 573
 574
 575
 576 sub outputInColumns {
 577     my ($items, $n, $w) = @_;
 578     my $i = 0;
 579     my @items = @$items;
 580     my $maxl = length([sort {length($b)<=>length($a)} @items]->[0]);
 581     $n ||= int(60/($maxl+3)) || 1;
 582     $w ||= $maxl+3;
 583     my $coll = int(@items/$n);
 584     $coll == @items/$n ? $coll : ++$coll;
 585     my @t;
 586     for my $j (0..$n-1) {
 587         if ($j<$n-1) {
 588             $t[$j] = [@items[$j*$coll..$j*$coll+$coll-1]];
 589         }
 590         else {
 591             $t[$j] = [@items[$j*$coll..$#items]];
 592         }
 593     }
 594     @items = map { my $j = $_; map { $t[$_]->[$j] || () } (0..$#t) } (0..scalar(@{$t[0]})-1);
 595     foreach (@items) {
 596         $_ .= (' ' x ($w-length($_)));
 597     }
 598     while ($i < @items) {
 599         print join('', map { $_ || ()  } @items[$i..$i+$n-1] ),"\n";
 600         $i+=$n;
 601     }
 602     return 1;
 603 }
 604
 605 sub debug {
 606     print STDERR shift()."\n" if $opt_v;
 607     return;
 608 }