9.3 Testing Bioperl
To check that things were working okay with my new Bioperl installation, I first wrote a little test to see if a Perl program could find the Bio::Perl module:
use Bio::Perl;
I ran it by putting it in a file bp0.pl and giving it to the Perl interpreter:
[tisdall]# perl bp0.pl [tisdall]#
As you see, it didn't complain, which means Perl found the Bio::Perl module. If it can't find it, it will complain. When I copied the bp0.pl program to a file called bp0.pl.broken and changed the module call of Bio::Perl to a call for the nonexistent module Bio::Perrl, I got the following (slightly truncated) error output:
[tisdall@coltrane development]$ perl bp0.pl.broken Can't locate Bio/Perrl.pm in @INC BEGIN failed--compilation aborted at bp0.pl.broken line 1.
9.3.1 Second Test
Now I knew that Bio::Perl could be found and loaded. I next tried a couple of test programs that are given in the bptutorial.pl document. I found a link to that tutorial document on the Web from the http://bioperl.org page. Alternately, I could have opened a window and typed at the command prompt:
perldoc bptutorial.pl
I went to the section near the beginning of the document called "I.2 Quick getting started scripts" and created a file tut1.pl on my computer by pasting in the text of the first tutorial script:
use Bio::Perl;
# this script will only work with an internet connection
# on the computer it is run on
$seq_object = get_sequence('swissprot',"ROA1_HUMAN");
write_sequence(">roa1.fasta",'fasta',$seq_object);
I then ran the program and looked at the output file:
[tisdall]$ perl tut1.pl [tisdall]$ cat roa1.fasta >ROA1_HUMAN Heterogeneous nuclear ribonucleoprotein A1 (Helix-destabilizing protein) (Single-strand binding protein) (hnRNP core protein A1). SKSESPKEPEQLRKLFIGGLSFETTDESLRSHFEQWGTLTDCVVMRDPNTKRSRGFGFVT YATVEEVDAAMNARPHKVDGRVVEPKRAVSREDSQRPGAHLTVKKIFVGGIKEDTEEHHL RDYFEQYGKIEVIEIMTDRGSGKKRGFAFVTFDDHDSVDKIVIQKYHTVNGHNCEVRKAL SKQEMASASSSQRGRSGSGNFGGGRGGGFGGNDNFGRGGNFSGRGGFGGSRGGGGYGGSG DGYNGFGNDGGYGGGGPGYSGGSRGYGSGGQGYGNQGSGYGGSGSYDSYNNGGGRGFGGG SGSNFGGGGSYNDFGNYNNQSSNFGPMKGGNFGGRSSGPYGGGGQYFAKPRNQGGYGGSS SSSSYGSGRRF [tisdall]$
That seemed to work perfectly.
9.3.2 Third Test
I tried the next short script from the same section of the tutorial, pasting it into a file called tut2.pl:
[tisdall]$ cat tut2.pl
use Bio::Perl;
# this script will only work with an internet connection
# on the computer it is run on
$seq_object = get_sequence('swissprot',"ROA1_HUMAN");
# uses the default database - nr in this case
$blast_result = blast_sequence($seq);
write_blast(">roa1.blast",$blast_report);
[tisdall]$ perl tut2.pl
-------------------- WARNING ---------------------
MSG: req was POST http://www.ncbi.nlm.nih.gov/blast/Blast.cgi
User-Agent: libwww-perl/5.69
Content-Length: 178
Content-Type: application/x-www-form-urlencoded
...
---------------------------------------------------
Submitted Blast for [blast-sequence-temp-id]
[tisdall]$ cat roa1.blast
[tisdall]$ ls -l roa1.blast
-rw-rw-r-- 1 tisdall tisdall 0 Apr 30 11:28 roa1.blast
[tisdall]$
Here, I experienced a problem. Running the program created a page of error output, mostly formatted in HTML, which I've truncated in the output. When I tried to cat the output file, there was nothing in it, which I verified (on Linux) by examining the file with ls -l, which showed that it had 0 bytes in it.
This wasn't very encouraging; the second of the two short example scripts was failing.
Looking at the two programs tut1.pl and tut2.pl, you can see that the second is an extension of the first; after getting the sequence, the second program also tries to blast it, and this is where it is failing. It's running (printing those ... dots took a while because the program was waiting for a reply from NCBI over the Internet), but it's not printing out the BLAST report to the file roa1.blast. What's going wrong?
I started my investigation by looking at the documentation for the Bio::Perl module by typing:
perldoc Bio::Perl
and searching for the function that is failing, blast_sequence. Here's what I found:
blast_sequence
Title : blast_sequence
Usage : $blast_result = blast_sequence($seq)
$blast_result = blast_sequence('MFVEGGTFASEDDDSASAEDE');
Function: If the computer has Internet accessibility, blasts
the sequence using the NCBI BLAST server against nrdb.
It choose the flavour of BLAST on the basis of the sequence.
This function uses Bio::Tools::Run::RemoteBlast, which itself
use Bio::SearchIO - as soon as you want to more, check out
these modules
Returns : Bio::Search::Result::GenericResult.pm
Args : Either a string of protein letters or nucleotides, or a
Bio::Seq object
That last sentence about the Args gave me a clue. I went back and looked at the failing program, which I placed in a file called tut2.pl, and, sure enough, the argument for the blast_sequence function, $seq, is neither a string of sequence nor a Bio::Seq object. In fact, it's not even defined! Clearly, what the tutorial writer meant instead of $seq, was $seq_object, which is defined and populated by the earlier get_sequence method call.
I was visited by a brainwave, and I decided to put a use strict and a use warnings into the program and to declare variables with my to see what ensued:
[tisdall]$ cat tut2.pl
use Bio::Perl;
use strict;
use warnings;
# this script will only work with an internet connection
# on the computer it is run on
my $seq_object = get_sequence('swissprot',"ROA1_HUMAN");
# uses the default database - nr in this case
my $blast_result = blast_sequence($seq);
write_blast(">roa1.blast",$blast_report);
[tisdall]$ perl tut2.pl
Global symbol "$seq" requires explicit package name at tut2.pl line 11.
Global symbol "$blast_report" requires explicit package name at tut2.pl line 13.
Execution of tut2.pl aborted due to compilation errors.
[tisdall]$
Well, that's pretty clear. The variables $seq and $blast_report are wrong; apparently, the author intended to reuse the variables $seq_object and $blast_result instead. So, I edited the file and ran it again:
[tisdall]$ cat tut2.pl
use Bio::Perl;
use strict;
use warnings;
# this script will only work with an internet connection
# on the computer it is run on
my $seq_object = get_sequence('swissprot',"ROA1_HUMAN");
# uses the default database - nr in this case
my $blast_result = blast_sequence($seq_object);
write_blast(">roa1.blast",$blast_result);
[tisdall]$ perl tut2.pl
[tisdall]$ perl tut2.pl
Submitted Blast for [ROA1_HUMAN] ...................
[tisdall]$ ls -l roa1.blast
-rw-rw-r-- 1 tisdall tisdall 56888 May 5 15:05 roa1.blast
[tisdall]$
Examining the roa1.blast file convinced me that the program had successfully called blast. I decided that was a success, albeit a qualified one: this was a spot in the documentation that could use a little attention, clearly. By the time you're reading this, it may well have been fixed.
9.3.3 Fourth Test
Now, let me show you one more test of my new Bioperl installation.
Looking at the Bio::Perl documentation, I found the following example and discussion at the very beginning.
Bio::Perl(3) User Contributed Perl Documentation Bio::Perl(3)
NAME
Bio::Perl - Functional access to BioPerl for people who
don't know objects
SYNOPSIS
use Bio::Perl;
# will guess file format from extension
$seq_object = read_sequence($filename);
# forces genbank format
$seq_object = read_sequence($filename,'genbank');
# reads an array of sequences
@seq_object_array = read_all_sequences($filename,'fasta');
# sequences are Bio::Seq objects, so the following methods work
# for more info see L<Bio::Seq>, or do 'perldoc Bio/Seq.pm'
print "Sequence name is ",$seq_object->display_id,"\n";
print "Sequence acc is ",$seq_object->accession_number,"\n";
print "First 5 bases is ",$seq_object->subseq(1,5),"\n";
# get the whole sequence as a single string
$sequence_as_a_string = $seq_object->seq( );
# writing sequences
write_sequence(">$filename",'genbank',$seq_object);
write_sequence(">$filename",'genbank',@seq_object_array);
# making a new sequence from just strings you have
# from something else
$seq_object = new_sequence("ATTGGTTTGGGGACCCAATTTGTGTGTTATATGTA",
"myname","AL12232");
# getting a sequence from a database (assumes internet connection)
$seq_object = get_sequence('swissprot',"ROA1_HUMAN");
$seq_object = get_sequence('embl',"AI129902");
$seq_object = get_sequence('genbank',"AI129902");
# BLAST a sequence (assummes an internet connection)
$blast_report = blast_sequence($seq_object);
write_blast(">blast.out",$blast_report);
DESCRIPTION
Easy first time access to BioPerl via functions
FEEDBACK
Mailing Lists
User feedback is an integral part of the evolution of this
and other Bioperl modules. Send your comments and sugges-
tions preferably to one of the Bioperl mailing lists.
Your participation is much appreciated.
bioperl-l@bio.perl.org
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://bugzilla.bioperl.org/
AUTHOR - Ewan Birney
Email bioperl-l@bio.perl.org
Describe contact details here
APPENDIX
The rest of the documentation details each of the object
methods. Internal methods are usually preceded with a _
...
I took the example code from the SYNOPSIS section and put it in a file called bp1.pl. I noticed that the code was not meant to be a running program. The variable $filename refers to some sequence file in the first line of code without being initialized, and then the same variable name is used in several other lines for different purposes. This is not an uncommon situation in such SYNOPSIS sections; the point is to show how individual calls can be made to methods in the class, not to demonstrate a complete working program (although sometimes the code can be run exactly as is).
I had to make some changes to make this code a working, runnable program. I declared the strict and warnings pragmas and declared each variable with my. I created variables for the different sequence filenames I needed as both input and output, and I made sure that the input sequence files were on disk and of the correct variety (for example, I created the array.fasta file with three FASTA headers and sequences one after the other). In the end I had this code:
use Bio::Perl;
use strict;
use warnings;
my $gbfilename = 'AI129902.genbank';
# will guess file format from extension
my $seq_object0 = read_sequence($gbfilename);
# forces genbank format
my $seq_object1 = read_sequence($gbfilename,'genbank');
my $fastafilename = 'array.fasta';
# reads an array of sequences
my @seq_object_array = read_all_sequences($fastafilename,'fasta');
# sequences are Bio::Seq objects, so the following methods work
# for more info see L<Bio::Seq>, or do 'perldoc Bio/Seq.pm'
print "Sequence name is ",$seq_object1->display_id,"\n";
print "Sequence acc is ",$seq_object1->accession_number,"\n";
print "First 5 bases is ",$seq_object1->subseq(1,5),"\n";
# get the whole sequence as a single string
my $sequence_as_a_string = $seq_object1->seq( );
# writing sequences
my $gbfilenameout = 'bpout.genbank';
write_sequence(">$gbfilenameout",'genbank',$seq_object1);
write_sequence(">$gbfilenameout",'genbank',@seq_object_array);
# making a new sequence from just strings you have
# from something else
my $seq_object2 = new_sequence("ATTGGTTTGGGGACCCAATTTGTGTGTTATATGTA",
"myname","AL12232");
# getting a sequence from a database (assumes internet connection)
my $seq_object3 = get_sequence('swissprot',"ROA1_HUMAN");
my $seq_object4 = get_sequence('embl',"AI129902");
my $seq_object5 = get_sequence('genbank',"AI129902");
# BLAST a sequence (assummes an internet connection)
my $blast_report = blast_sequence($seq_object3);
write_blast(">blast.out",$blast_report);
As you see, I didn't check on the output of all the method calls, but the real purpose of this test of my newly installed Bioperl modules was just to see if all the modules and methods could be found and would run when called.
As mentioned previously, I also added use strict; and use warnings; and declared all the variables with my. Often, you don't see that usage in this kind of documentation; it's not required in Perl, and it may distract some readers of the documentation from the main point of showing how the methods are called. This is not an unusual omission to find in Perl class documentation, but of course, it's recommended and often very helpful, as you saw in the last section.
So, I ran my slightly edited version of the example code from the beginning of the Bio::Perl manpage, with the following results:
[tisdall]$ perl bp1.pl Sequence name is AI129902 Sequence acc is AI129902 First 5 bases is CTCCG Submitted Blast for [ROA1_HUMAN] ......... [tisdall]$
I looked at each of the input and output files and verified the expected contents. The following output from a listing of the files will give you an idea of what to expect (although you may get different results by running the program with different input files or database lookups):
-rw-rw-r-- 1 tisdall tisdall 2391 May 5 10:37 AI129902.genbank -rw-rw-r-- 1 tisdall tisdall 2485 May 5 13:00 array.fasta -rw-rw-r-- 1 tisdall tisdall 1513 May 5 13:02 bp1.pl -rw-rw-r-- 1 tisdall tisdall 3653 May 5 13:02 bpout.genbank -rw-rw-r-- 1 tisdall tisdall 56888 May 5 13:04 blast.out
Notice that at the end of the program I call the blast_sequence method on the protein sequence object $seq_object3. I discovered that trying to run the blast_sequence method on a nucleotide sequence object (such as $seq_object5) failed. Although the documentation for the method said that the sequence type would be examined and the appropriate BLAST program called (for example, blastp for protein sequence and blastx for nucleotide sequence, against the nr nonredundant protein database), it always seemed to call blastp no matter what the input sequence, and therefore it failed unless called with protein sequence as I had stored in $seq_object3. Perhaps this bug, a disconnect between the code and the documentation, has been fixed by the time you read this.
