Bioperl: Extracting genes from GenBank flat file (fwd)

Ewan Birney birney@ebi.ac.uk
Wed, 14 Jun 2000 13:27:03 +0100 (GMT)


  This message is in MIME format.  The first part should be readable text,
  while the remaining parts are likely unreadable without MIME-aware tools.
  Send mail to mime@docserver.cac.washington.edu for more info.

--------------24D1E38F782C0746943A7695
Content-Type: TEXT/PLAIN; CHARSET=US-ASCII
Content-ID: <Pine.LNX.4.10.10006141326462.30975@riker.ebi.ac.uk>



Imre meant to send this to bioperl, but send it to me instead ;)


-----------------------------------------------------------------
Ewan Birney. Mobile: +44 (0)7970 151230, Work: +44 1223 494420
<birney@ebi.ac.uk>. 
-----------------------------------------------------------------

---------- Forwarded message ----------
Date: Wed, 14 Jun 2000 13:56:50 +0300
From: Imre Vastrik <imre.vastrik@helsinki.fi>
To: Ewan Birney <birney@ebi.ac.uk>
Subject: Re: Bioperl: Extracting genes from GenBank flat file

Sorry for joining this thread so late, but it was just now that I
noticed few things about using Bio::DB::GenBank module. In it's present
form (in bioperl-0.6.1 and live) this module fetches fasta-formatted
sequence. That's the reason why features/annotations are not there.
However, for getting the thing work properly one has to do just 4 small
changes:

on line 115 change 'dopt=f' to 'dopt=g'
on line 223 change 'FORMAT=1' to 'FORMAT=6'
on line 246 change 'Fasta' to 'genbank'
on line 277 change 'Fasta' to 'genbank'

I include the the modified version of GenBank.pm at as an attachment
(don't know if this is PC these days ;-)). Maybe someone with cvs
writing rights can put it to the appropriate location.

imre

--------------------------------------------
- Imre Vastrik, Ph.D., tel:+359-9-19126886 -
-    http://www.hi.helsinki.fi/bioinfo     -
--------------------------------------------

gert thijs wrote:
> 
> I am trying to develop an some perl script to extract the location of CDS or
> genes in a DNA sequence based on the features described in the genbank flat
> file.
> I have two question:
> - Has anyone ever tried such things before? I haven't found anything browsing
> through the mail archive.
> - I use Bio::DB::GenBank to download a sequence from genbank, now I was
> wondering if the features in the genbank description are also downloaded and
> stored somewhere. I thought I could use the annotation object to extract this
> information but this seems not to be the case.
> I use something like this:
> $seq = $gb->get_Seq_by_acc($Accn[$i]);
> $ann = $seq->annotation();
> print $ann->description();
> 
> But then $ann->description() seems to be an empty string.

--------------24D1E38F782C0746943A7695
Content-Type: TEXT/PLAIN; CHARSET=us-ascii; NAME="GenBank.pm"
Content-ID: <Pine.LNX.4.10.10006141326463.30975@riker.ebi.ac.uk>
Content-Description: 
Content-Disposition: INLINE; FILENAME="GenBank.pm"

# test
#
# BioPerl module for Bio::DB::GenBank
#
# Cared for by Aaron Mackey <amackey@virginia.edu>
#
# Copyright Aaron Mackey
#
# You may distribute this module under the same terms as perl itself

# POD documentation - main docs before the code

=head1 NAME

Bio::DB::GenBank - Database object interface to GenBank

=head1 SYNOPSIS

    $gb = new Bio::DB::GenBank;

    $seq = $gb->get_Seq_by_id('MUSIGHBA1'); # Unique ID

    # or ...

    $seq = $gb->get_Seq_by_acc('J00522'); # Accession Number

=head1 DESCRIPTION

Allows the dynamic retrieval of Sequence objects (Bio::Seq) from the GenBank
database at NCBI, via an Entrez query.

WARNING: Please do NOT spam the Entrez web server with multiple requests.
NCBI offers Batch Entrez for this purpose.  Batch Entrez support will likely
be supported in a future version of DB::GenBank.

=head1 FEEDBACK

=head2 Mailing Lists

User feedback is an integral part of the
evolution of this and other Bioperl modules. Send
your comments and suggestions preferably to one
of the Bioperl mailing lists. Your participation
is much appreciated.

  vsns-bcd-perl@lists.uni-bielefeld.de          - General discussion
  vsns-bcd-perl-guts@lists.uni-bielefeld.de     - Technically-oriented discussion
  http://bio.perl.org/MailList.html             - About the mailing lists

=head2 Reporting Bugs

Report bugs to the Bioperl bug tracking system to
help us keep track the bugs and their resolution.
Bug reports can be submitted via email or the
web:

  bioperl-bugs@bio.perl.org
  http://bio.perl.org/bioperl-bugs/

=head1 AUTHOR - Aaron Mackey

Email amackey@virginia.edu

=head1 APPENDIX

The rest of the documentation details each of the
object methods. Internal methods are usually
preceded with a _

=cut

# Let the code begin...

package Bio::DB::GenBank;
use strict;
use vars qw(@ISA);

# Object preamble - inherits from Bio::DB::RandomAccessI

use Bio::DB::RandomAccessI;
use Bio::SeqIO;
use IO::Socket;
use IO::File;

@ISA = qw(Bio::Root::Object Bio::DB::RandomAccessI);

# new() is inherited from Bio::Root::Object

# _initialize is where the heavy stuff will happen when new is called

sub _initialize {
  my($self,@args) = @_;

  my $make = $self->SUPER::_initialize;
  
# set stuff in self from @args  
  return $make; # success - we hope!
}

=head2 get_Seq_by_id

 Title   : get_Seq_by_id
 Usage   : $seq = $db->get_Seq_by_id($uid);
 Function: Gets a Bio::Seq object by its unique identifier/name
 Returns : a Bio::Seq object
 Args    : $uid : the id (as a string) of the desired sequence entry

=cut

sub get_Seq_by_id {

  my $self = shift;
  my $uid = shift or $self->throw("Must supply an identifier!\n");

  my $entrez = "db=n&form=6&dopt=g&html=no&title=no&uid=$uid";

  my $stream = $self->_get_stream($entrez);
  my $seq = $stream->next_seq();
  $self->throw("Unable to get seq for id $uid, is it really a genbank id?\n") 
      if ( !defined $seq );
  return $seq;
}

=head2 get_Seq_by_acc

  Title   : get_Seq_by_acc
  Usage   : $seq = $db->get_Seq_by_acc($acc);
  Function: Gets a Bio::Seq object by its accession number
  Returns : a Bio::Seq object
  Args    : $acc : the accession number of the desired sequence entry
  Note    : For GenBank, this just calls the same code for get_Seq_by_id()

=cut

sub get_Seq_by_acc {

  my $self = shift;
  my $acc = shift or $self->throw("Must supply an accesion number!\n");
  
  return $self->get_Seq_by_id($acc);
}

=head2 get_Stream_by_id

  Title   : get_Stream_by_id
  Usage   : $stream = $db->get_Stream_by_id( [$uid1, $uid2] );
  Function: Gets a series of Seq objects by unique identifiers
  Returns : a Bio::SeqIO stream object
  Args    : $ref : a reference to an array of unique identifiers for
                   the desired sequence entries


=cut

sub get_Stream_by_id {

  my $self = shift;
  my $id = shift or $self->throw("Must supply a unique identifier!\n");
  ref($id) eq "ARRAY" or $self->throw("Must supply an array ref!\n");

  my $uid = join(',', @{$id});
  my $entrez = "db=n&form=6&dopt=f&html=no&title=no&uid=$uid" ;

  return $self->_get_stream($entrez);

}

=head2 get_Stream_by_acc

  Title   : get_Stream_by_acc
  Usage   : $seq = $db->get_Seq_by_acc($acc);
  Function: Gets a series of Seq objects by accession numbers
  Returns : a Bio::SeqIO stream object
  Args    : $ref : a reference to an array of accession numbers for
                   the desired sequence entries
  Note    : For GenBank, this just calls the same code for get_Stream_by_id()

=cut

sub get_Stream_by_acc {

  my $self = shift;
  my $acc = shift or $self->throw("Must supply an accession number!\n");

  return $self->get_Seq_by_id($acc);
}

=head2 get_Stream_by_batch

  Title   : get_Stream_by_batch
  Usage   : $seq = $db->get_Stream_by_batch($ref);
  Function: Retrieves Seq objects from Entrez 'en masse', rather than one
            at a time.  For large numbers of sequences, this is far superior
            than get_Stream_by_[id/acc]().
  Example :
  Returns : a Bio::SeqIO stream object
  Args    : $ref : either an array reference, a filename, or a filehandle
            from which to get the list of unique id's/accession numbers.

=cut

sub get_Stream_by_batch {
   my $self = shift;
   my $ref = shift or $self->throw("Must supply an argument!\n");
   my $which = ref($ref);
   my $fh;
   my $filename;
   if ( $which eq 'ARRAY') { # $ref is an array reference
       $fh = new_tmpfile IO::File;
       for ( @{$ref} ) {
	   print $fh $_ . "\n";
       }
       seek $fh, 0, 0;
       $filename = "tempfile.txt";
   } elsif ( $which eq '') { # $ref is a filename
       $fh = new IO::File $ref, "r";
       $filename = $ref;
   } elsif ( $which eq 'GLOB' or $which eq 'IO::File') { # $ref is assumed to be a filehandle
       $fh = $ref;
       $filename = "tempfile.txt";
   }

   my $wwwbuf = "DB=n&REQUEST_TYPE=LIST_OF_GIS&FORMAT=6&HTML=FALSE&SAVETO=FALSE&NOHEADER=TRUE&UID=" . join(',', grep { chomp; } <$fh> );

   my $sock = $self->_get_sock();

   select $sock;
   print "POST /cgi-bin/Entrez/qserver.cgi HTTP/1.0\015\012";
   print "Host: www.ncbi.nlm.nih.gov\015\012";
   print "User-Agent: $0::Bio::DB::GenBank\015\012";
   print "Connection: Keep-Alive\015\012";
   print "Content-type: application/x-www-form-urlencoded\015\012";
   print "Content-length: " . length($wwwbuf) . "\015\012";
   print "\015\012";
   print $wwwbuf;

   while (<$sock>) {
       if ( m,^HTTP/\d+\.\d+\s+(\d+)[^\012]\012, ) {
	   my $code = $1;
	   return undef unless $code =~ /^2/;
       }
       $self->throw("Entrez Error - check query sequences!\n") if m/^ERROR/i;
       last if m/Batch Entrez results/;
   }

   return Bio::SeqIO->new('-fh' => $sock, '-format' => 'genbank');

}



sub _get_stream {

  my($self, $entrez) = @_;

# most of this socket stuff is borrowed heavily from LWP::Simple, by
# Gisle Aas and Martijn Koster.  They copyleft'ed it, but we should give
# them full credit for this little diddy.

  my $sock = $self->_get_sock();

  print $sock join("\015\012" =>
		   "GET /htbin-post/Entrez/query?$entrez HTTP/1.0",
		   "Host: www.ncbi.nlm.nih.gov",
		   "User-Agent: $0::Bio::DB::GenBank",
		   "", "");

  while(<$sock>) {
    if ( m,^HTTP/\d+\.\d+\s+(\d+)[^\012]\012, ) {
      my $code = $1;
      return undef unless $code =~ /^2/;
    }
    $self->throw("Entrez Error - check query sequences!\n") if m/^ERROR/i;
    last if m/^------/; # Kludgy, but it's how L. Stein does Boulder too
  }

  return Bio::SeqIO->new('-fh' => $sock, '-format' => 'genbank');

}

sub _get_sock {
    my $self = shift;
    my $sock = IO::Socket::INET->new(PeerAddr => 'www.ncbi.nlm.nih.gov',
				     PeerPort => 80,
				     Proto    => 'tcp',
				     Timeout  => 60
				     );
    unless ($sock) {
	$@ =~ s/^.*?: //;
	$self->throw("Can't connect to GenBank ($@)\n");
    }
    $sock->autoflush();		# just for safety's sake if they have old IO::Socket

    return $sock;
}


1;
__END__












--------------24D1E38F782C0746943A7695--
=========== Bioperl Project Mailing List Message Footer =======
Project URL: http://bio.perl.org/
For info about how to (un)subscribe, where messages are archived, etc:
http://www.techfak.uni-bielefeld.de/bcd/Perl/Bio/vsns-bcd-perl.html
====================================================================