BioPerlTutorial - a tutorial for bioperl
1.2.3
Written by Peter Schattner <schattner@alum.mit.edu> Copyright Peter Schattner Contributions, additions and corrections have been made to this document by the following individuals: Jason Stajich Heikki Lehvaslaiho Brian Osborne Hilmar Lapp Chris Dagdigian Elia Stupka Ewan Birney
This tutorial includes "snippets" of code and text from various Bioperl documents including module documentation, example scripts and "t" test scripts. You may distribute this tutorial under the same terms as perl itself. This document is written in Perl POD (plain old documentation) format. You can run this file through your favorite pod translator (pod2html, pod2man, pod2text, etc.) if you would like a more convenient formatting. Table of Contents I. Introduction I.1 Overview I.2 Quick getting started scripts I.3 Software requirements I.3.1 For minimal bioperl installation I.3.2 For complete installation I.4 Installation procedures I.5 Additional comments for non-unix users I.6 Places to look for additional documentation II. Brief overview to bioperl's objects II.1 Sequence objects (Seq, PrimarySeq, LocatableSeq, RelSegment, LiveSeq, LargeSeq, RichSeq, SeqWithQuality, SeqI) II.2 Location objects (Simple, Split, Fuzzy) II.3 Interface objects and implementation objects III. Using bioperl III.1 Accessing sequence data from local and remote databases III.1.1 Accessing remote databases (Bio::DB::GenBank, etc) III.1.2 Indexing and accessing local databases (Bio::Index::*, bp_index.pl, bp_fetch.pl) III.2 Transforming formats of database/ file records III.2.1 Transforming sequence files (SeqIO) III.2.2 Transforming alignment files (AlignIO) III.3 Manipulating sequences III.3.1 Manipulating sequence data with Seq methods (Seq) III.3.2 Obtaining basic sequence statistics (SeqStats,SeqWord) III.3.3 Identifying restriction enzyme sites (Bio::Restriction) III.3.4 Identifying amino acid cleavage sites (Sigcleave) III.3.5 Miscellaneous sequence utilities: OddCodes, SeqPattern III.3.6 Converting coordinate systems (Coordinate::Pair, RelSegment) III.4 Searching for similar sequences III.4.1 Running BLAST remotely (using RemoteBlast.pm) III.4.2 Parsing BLAST and FASTA reports with Search and SearchIO III.4.3 Parsing BLAST reports with BPlite, BPpsilite, and BPbl2seq III.4.4 Parsing HMM reports (HMMER::Results, SearchIO) III.4.5 Running BLAST locally (StandAloneBlast) III.5 Manipulating sequence alignments (SimpleAlign) III.6 Searching for genes and other structures on genomic DNA (Genscan, Sim4, ESTScan, MZEF, Grail, Genemark, EPCR) III.7 Developing machine readable sequence annotations III.7.1 Representing sequence annotations (SeqFeature,RichSeq,Location) III.7.2 Representing sequence annotations (Annotation::Collection) III.7.3 Representing large sequences (LargeSeq) III.7.4 Representing changing sequences (LiveSeq) III.7.5 Representing related sequences - mutations, polymorphisms (Allele, SeqDiff) III.7.6 Incorporating quality data in sequence annotation (SeqWithQuality) III.7.7 Sequence XML representations - generation and parsing (SeqIO::game) III.7.8 Representing Sequence Features using GFF (Bio:Tools:GFF) III.8 Manipulating clusters of sequences (Cluster, ClusterIO) III.9 Representing non-sequence data in Bioperl: structures, trees, maps, graphics and bibliographic text III.9.1 Using 3D structure objects and reading PDB files (StructureI, Structure::IO) III.9.2 Tree objects and phylogenetic trees (Tree::Tree, TreeIO, PAML.pm ) III.9.3 Map objects for manipulating genetic maps (Map::MapI, MapIO) III.9.4 Bibliographic objects for querying bibliographic databases (Biblio) III.9.5 Graphics objects for representing sequence objects as images (Graphics) III.10 Bioperl alphabets III.10.1 Extended DNA / RNA alphabet III.10.2 Amino Acid alphabet IV. Auxiliary Bioperl Libraries (Bioperl-run, Bioperl-db, etc.) IV.1 Using the Bioperl Auxiliary Libraries IV.2 Running programs (Bioperl-run and Bioperl-ext) IV.2.1 Sequence manipulation using the Bioperl EMBOSS and PISE interfaces IV.2.2 Aligning 2 sequences with Blast using bl2seq and AlignIO IV.2.3 Aligning multiple sequences (Clustalw.pm, TCoffee.pm) IV.2.4 Aligning 2 sequences with Smith-Waterman (pSW) IV.3 Bioperl-db and BioSQL IV.4 Other Bioperl auxiliary libraries V. Appendices V.1 Finding out which methods are used by which Bioperl Objects V.2 Tutorial Demo Scripts
Bioperl is a collection of perl modules that facilitate the development of perl scripts for bioinformatics applications. As such, it does not include ready to use programs in the sense that many commercial packages and free web-based interfaces do (e.g. Entrez, SRS). On the other hand, bioperl does provide reusable perl modules that facilitate writing perl scripts for sequence manipulation, accessing of databases using a range of data formats and execution and parsing of the results of various molecular biology programs including Blast, clustalw, TCoffee, genscan, ESTscan and HMMER. Consequently, bioperl enables developing scripts that can analyze large quantities of sequence data in ways that are typically difficult or impossible with web based systems.
In order to take advantage of bioperl, the user needs a basic understanding of the perl programming language including an understanding of how to use perl references, modules, objects and methods. If these concepts are unfamiliar the user is referred to any of the various introductory or intermediate books on perl. We've liked S. Holzmer's Perl Core Language, Coriolis Technology Press, for example. This tutorial is not intended to teach the fundamentals of perl to those with little or no experience in the perl language. On the other hand, advanced knowledge of perl - such as how to write a object-oriented perl module - is not required for successfully using bioperl.
Bioperl is open source software that is still under active development. The advantages of open source software are well known. They include the ability to freely examine and modify source code and exemption from software licensing fees. However, since open source software is typically developed by a large number of volunteer programmers, the resulting code is often not as clearly organized and its user interface not as standardized as in a mature commercial product. In addition, in any project under active development, documentation may not keep up with the development of new features. Consequently the learning curve for actively developed, open source source software is sometimes steep.
This tutorial is intended to ease the learning curve for new users of bioperl. To that end the tutorial includes:
Descriptions of what bioinformatics tasks can be handled with bioperl
Directions on where to find the methods to accomplish these tasks within the bioperl package
Recommendations on where to go for additional information.
A runnable script, bptutorial.pl, which demonstrates many of the capabilities of Bioperl. Runnable example code can also be found in the scripts/ and examples/ directories. Summary descriptions of all of these scripts can be found in the file bioscripts.pod (or http://bioperl.org/Core/Latest/bioscripts.html). In addition, the POD documentation for many Bioperl modules should contain runnable code in the SYNOPSIS section which is meant to illustrate the use of a module and its methods. You will also find some interesting bits of code in the FAQ (http://bioperl.org/Core/Latest/faq.html).
Running the bptutorial.pl script while going through this tutorial - or better yet, stepping through it with an interactive debugger - is a good way of learning bioperl. The tutorial script is also a good place from which to cut-and-paste code for your scripts (rather than using the code snippets in this tutorial). Most of the scripts in the tutorial script should work on your machine - and if they don't it would probably be a good idea to find out why, before getting too involved with bioperl! Some of the demos require optional modules from the bioperl auxiliary libraries and/or external programs. These demos should be skipped if the demos are run and the required auxiliary programs are not found.
For newcomers and people who want to quickly evaluate whether this package is worth using in the first place, we have a very simple module which allows easy access to a small number of Bioperl's functionality in an easy to use manner. The Bio::Perl module provides some simple access functions, for example, this script will retrieve a swissprot sequence and write it out in fasta format
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);
Another example is the ability to blast a sequence using the facilities as NCBI. Please be careful not to abuse the compute that NCBI provides and so use this only for individual searches. If you want to do a large number of BLAST searches, please download the blast package locally.
use Bio::Perl; # this script will only work with an internet connection # on the computer it is run on $seq = get_sequence('swissprot',"ROA1_HUMAN"); # uses the default database - nr in this case $blast_result = blast_sequence($seq); write_blast(">roa1.blast",$blast_result);
Bio::Perl has a number of other easy-to-use functions, including
get_sequence - gets a sequence from standard, internet accessible databases read_sequence - reads a sequence from a file read_all_sequences - reads all sequences from a file new_sequence - makes a bioperl sequence just from a string write_sequence - writes a single or an array of sequence to a file translate - provides a translation of a sequence translate_as_string - provides a translation of a sequence, returning back just the sequence as a string blast_sequence - BLASTs a sequence against standard databases at NCBI write_blast - writes a blast report out to a file
Using the Bio::Perl.pm module, it is possible to manipulate sequence data in Bioperl without explicitly creating Seq or SeqIO objects described later in this tutorial. However, only limited data manipulation is supported in this mode.
Look at the documentation in Bio::Perl by going 'perldoc Bio::Perl' to learn more about these functions. In all these cases, Bio::Perl accesses a subset of the underlying Bioperl functions (for example, translation in Bioperl can handle many different translation tables and provides different options for stop codon processing) - in most cases, most users will migrate to using the underlying bioperl objects as their sophistication level increases, but Bio::Perl provides an easy on-ramp for newcomers and lazy programmers. Also see examples/bioperl.pl for more examples of usage of this module.
What's required to run bioperl.
For a minimal installation of bioperl, you will need to have perl itself installed as well as the bioperl "core modules". Bioperl has been tested primarily using perl 5.005, 5.6, and 5.8. The minimal bioperl installation should still work under perl 5.004. However, as increasing numbers of bioperl objects are using modules from CPAN (see below), problems have been observed for bioperl running under perl 5.004. So if you are having trouble running bioperl under perl 5.004, you should probably upgrade your version of perl.
In addition to a current version of perl, the new user of bioperl is encouraged to have access to, and familiarity with, an interactive perl debugger. Bioperl is a large collection of complex interacting software objects. Stepping through a script with an interactive debugger is a very helpful way of seeing what is happening in such a complex software system - especially when the software is not behaving in the way that you expect. The free graphical debugger ptkdb is highly recommended - it's available as Devel::ptkdb from CPAN. The standard perl distribution also contains a powerful interactive debugger with a command-line interface (use it like "perl -d <script>").
The Perl tool Data::Dumper used with the syntax:
use Data::Dumper; print Dumper($seqobj);
can also be helpful for obtaining debugging information on Bioperl objects.
Some of the capabilities of bioperl require software beyond that of the minimal installation. This additional software includes perl modules from CPAN, package-libraries from bioperl's auxiliary code-repositories, a bioperl xs-extension, and several standard compiled bioinformatics programs.
Perl - extensions
For a complete listing of external Perl modules required by bioperl please see the INSTALL file in the Bioperl package.
Bioperl auxiliary repositories
Some features of bioperl that require modules from bioperl's auxiliary code repositories. See section IV and references therein for further installation instructions for these modules.
Bioperl C extensions & external bioinformatics programs
Bioperl also uses several C programs for sequence alignment and local blast searching. To use these features of bioperl you will need an ANSI C or Gnu C compiler as well as the actual program available from sources such as:
for Smith-Waterman alignments- bioperl-ext-0.6 from http://bioperl.org/Core/external.shtml
for clustalw alignments- ftp://ftp.ebi.ac.uk/pub/software/unix/clustalw/ ftp://ftp-igbmc.u-strasbg.fr/pub/ClustalW/
for tcoffee alignments- http://igs-server.cnrs-mrs.fr/~cnotred/Projects_home_page/t_coffee_home_page.html
for local blast searching- ftp://ftp.ncbi.nih.gov/blast/executables/release/
for EMBOSS applications - http://www.emboss.org
The actual installation of the various system components is accomplished in the standard manner:
Locate the package on the network
Download
Decompress (with gunzip or a similiar utility)
Extract the file archive (e.g. with tar -xvf)
Create a Makefile with "perl Makefile.PL"
Run "make", "make test" and "make install". This procedure must be repeated for every CPAN module, bioperl-extension and external module to be installed. A helper module CPAN.pm is available from CPAN which automates the process for installing the perl modules.
The CPAN module can also be used to install all of the modules listed above in a single step as a "bundle" of modules, Bundle::BioPerl, eg
$>perl -MCPAN -e shell cpan>install Bundle::BioPerl <installation details....> cpan>install B/BI/BIRNEY/bioperl-1.2.2.tar.gz <installation details....> cpan>quit
Be advised that version numbers change regularly, so the number used above may not apply. A disadvantage of the "bundle" approach is that if there's a problem installing any individual module it may be a bit more difficult to isolate.
See bioperl's INSTALL file for more details.
For the external programs (clustal, Tcoffee, ncbi-blast), there is an extra step:
Set the relevant environmental variable (CLUSTALDIR, TCOFFEEDIR or BLASTDIR) to the directory holding the executable in your startup file - e.g. in .bashrc or .tcshrc. For running local blasts, it is also necessary that the name of local-blast database directory is known to bioperl. This will typically happen automatically, but in case of difficulty, refer to the documentation in Bio::Tools::Run::StandAloneBlast.
The only likely complication (at least on unix systems) that may occur is if you are unable to obtain system level writing privileges. For instructions on modifying the installation in this case and for more details on the overall installation procedure, see the INSTALL file in the bioperl distribution as well as the README files in the external programs you want to use (e.g. bioperl-ext, clustalw, TCoffee, NCBI-blast).
Bioperl has mainly been developed and tested under various unix environments, including Linux and MacOS X. In addition, this tutorial has been written largely from a Unix perspective.
Mac users may find Steve Cannon's installation notes and suggestions for Bioperl on OS X at http://www.tc.umn.edu/~cann0010/Bioperl_OSX_install.html helpful. Also Todd Richmond has written of his experiences with BioPerl on MacOS 9 (http://bioperl.org/Core/mac-bioperl.html).
The bioperl core has also been tested and should work under most versions of Microsoft Windows. For many windows users the perl and bioperl distributions from Active State, at http://www.activestate.com has been quite helpful. Other windows users have had success running bioperl under Cygwin (http://www.cygwin.com). See the package's INSTALL.WIN file for more details.
Many bioperl features require the use of CPAN modules, compiled extensions or external programs. These features probably will not work under some or all of these other operating systems. If a script attempts to access these features from a non-unix OS, bioperl is designed to simply report that the desired capability is not available. However, since the testing of bioperl in these environments has been limited, the script may well crash in a less graceful manner.
This tutorial does not intend to be a comprehensive description of all the objects and methods available in bioperl. For that the reader is directed to the documentation included with each of the modules. A very useful interface for finding one's way within all the module documentation can be found at http://doc.bioperl.org/bioperl-live/. This interface lists all bioperl modules and descriptions of all of their methods. In addition, beginner questions can often be answered by looking at the FAQ, INSTALL and README files (http://bioperl.org/Core/Latest/faq.html, http://bioperl.org/Core/Latest/INSTALL, http://bioperl.org/Core/Latest/README )in the top-level directory of the bioperl distribution.
One potential problem in locating the correct documentation is that multiple methods in different modules may all share the same name. Moreover, because of perl's complex method of inheritance it is not often clear which of the identically named methods is being called by a given object. One way to resolve this question is by using the software described in Appendix "V.1".
For those who prefer more visual descriptions, http://bioperl.org/Core/Latest/modules.html also offers links to PDF files which contain class diagrams that describe how many of the bioperl objects related to one another (Version 1.0 Class Diagrams).
In addition, a bioperl online course is available on the web at http://www.pasteur.fr/recherche/unites/sis/formation/bioperl. The user is also referred to numerous bioperl scripts in the scripts/ and examples/ directories (see bioscripts.pod for a description of these scripts, or http://bioperl.org/Core/Latest/bioscripts.html).
Another source of focussed documentation is the HOWTO files, found either in the bioperl doc/howto directory or at http://bioperl.org/HOWTOs/. Current topics include OBDA Access, SeqIO, SearchIO, and BioGraphics.
The purpose of this tutorial is to get you using bioperl to solve real-life bioinformatics problems as quickly as possible. The aim is not to explain the structure of bioperl objects or perl object-oriented programming in general. Indeed, the relationships among the bioperl objects is not simple; however, understanding them in detail is fortunately not necessary for successfully using the package.
Nevertheless, a little familiarity with the bioperl object bestiary can be very helpful even to the casual user of bioperl. For example there are (at least) eight different "sequence objects" - Seq, PrimarySeq, LocatableSeq, RelSegment, LiveSeq, LargeSeq, SeqI, and SeqWithQuality. Understanding the relationships among these objects - and why there are so many of them - will help you select the appropriate one to use in your script.
This section describes various Bioperl sequence objects. Many people using Bioperl will never know, or need to know, what kind of sequence object they are using. This is because the SeqIO module, section section "III.2.1", creates exactly the right type of object when given a file or a filehandle or a string. But if you're curious, or if you need to create a sequence object manually for some reason, then read on.
Seq is the central sequence object in bioperl. When in doubt this is probably the object that you want to use to describe a DNA, RNA or protein sequence in bioperl. Most common sequence manipulations can be performed with Seq. These capabilities are described in sections "III.3.1" and "III.7.1", or in Bio::Seq.
Seq objects may be created for you automatically when you read in a file containing sequence data using the SeqIO object. This procedure is described in section "III.2.1". In addition to storing its identification labels and the sequence itself, a Seq object can store multiple annotations and associated "sequence features", such as those contained in most Genbank and EMBL sequence files. This capability can be very useful - especially in development of automated genome annotation systems, see section "III.7.1".
On the other hand, if you need a script capable of simultaneously handling hundreds or thousands sequences at a time, then the overhead of adding annotations to each sequence can be significant. For such applications, you will want to use the PrimarySeq object. PrimarySeq is basically a stripped-down version of Seq. It contains just the sequence data itself and a few identifying labels (id, accession number, alphabet = dna, rna, or protein), and no features. For applications with hundreds or thousands or sequences, using PrimarySeq objects can significantly speed up program execution and decrease the amount of RAM the program requires. See Bio::PrimarySeq for more details.
RichSeq objects store additional annotations beyond those used by standard Seq objects. If you are using sources with very rich sequence annotation, you may want to consider using these objects which are described in section "III.7.1". RichSeq objects are created automatically when Genbank, EMBL, or Swissprot format files are read by SeqIO.
SeqWithQuality objects areu sed to manipulate sequences with quality data, like those produced by phred. These objects are described in section "III.7.6", Bio::Seq::RichSeqI, and in Bio::Seq::SeqWithQuality.
What is called a LocatableSeq object for historical reasons might be more appropriately called an "AlignedSeq" object. It is a Seq object which is part of a multiple sequence alignment. It has start and end positions indicating from where in a larger sequence it may have been extracted. It also may have gap symbols corresponding to the alignment to which it belongs. It is used by the alignment object SimpleAlign and other modules that use SimpleAlign objects (e.g. AlignIO.pm, pSW.pm).
In general you don't have to worry about creating LocatableSeq objects because they will be made for you automatically when you create an alignment (using pSW, Clustalw, Tcoffee, Lagan, or bl2seq) or when you input an alignment data file using AlignIO. However if you need to input a sequence alignment by hand (e.g. to build a SimpleAlign object), you will need to input the sequences as LocatableSeqs. Other sources of information include Bio::LocatableSeq, Bio::SimpleAlign, Bio::AlignIO, and Bio::Tools::pSW.
The RelSegment object is also a type of bioperl Seq object. RelSegment objects are useful when you want to be able to manipulate the origin of the genomic coordinate system. This situation may occur when looking at a sub-sequence (e.g. an exon) which is located on a longer underlying underlying sequence such as a chromosome or a contig. Such manipulations may be important, for example when designing a graphical genome browser. If your code may need such a capability, look at the documentation Bio::DB::GFF::RelSegment which describes this feature in detail.
A LargeSeq object is a special type of Seq object used for handling very long sequences (e.g. > 100 MB). If you need to manipulate such long sequences see section "III.7.3" which describes LargeSeq objects, or Bio::Seq::LargeSeq.
A LiveSeq object is another specialized object for storing sequence data. LiveSeq addresses the problem of features whose location on a sequence changes over time. This can happen, for example, when sequence feature objects are used to store gene locations on newly sequenced genomes - locations which can change as higher quality sequencing data becomes available. Although a LiveSeq object is not implemented in the same way as a Seq object, LiveSeq does implement the SeqI interface (see below). Consequently, most methods available for Seq objects will work fine with LiveSeq objects. Section "III.7.4" and Bio::LiveSeq contain further discussion of LiveSeq objects.
SeqI objects are Seq "interface objects" (see section "II.4" and Bio::SeqI). They are used to ensure bioperl's compatibility with other software packages. SeqI and other interface objects are not likely to be relevant to the casual bioperl user.
A Location object is designed to be associated with a Sequence Feature object in order to show where the feature is on a longer sequence. Location objects can also be standalone objects used to described positions. The reason why these simple concepts have evolved into a collection of rather complicated objects is that:
1) Some objects have multiple locations or sub-locations (e.g. a gene's exons may have multiple start and stop locations) 2) In unfinished genomes, the precise locations of features is not known with certainty.
Bioperl's various Location objects address these complications. In addition there are CoordinatePolicy objects that allow the user to specify how to measure the length of a feature if its precise start and end coordinates are not known. In most cases, you will not need to worry about these complications if you are using bioperl to handle simple features with well-defined start and stop locations. However, if you are using bioperl to annotate partially or unfinished genomes or to read annotations of such genomes with bioperl, understanding the various Location objects will be important. See the documentation of the various modules in the Bio::Locations directory or Bio::Location::CoordinatePolicyI or section "III.7.1" for more information.
One goal of the design of Bioperl is to separate interface and implementation objects. An interface is solely the definition of what methods one can call on an object, without any knowledge of how it is implemented. An implementation is an actual, working implementation of an object. In languages like Java, interface definition is part of the language. In Perl, you have to roll your own.
In bioperl, the interface objects usually have names like Bio::MyObjectI, with the trailing I indicating it is an interface object. The interface objects mainly provide documentation on what the interface is, and how to use it, without any implementations (though there are some exceptions). Although interface objects are not of much direct utility to the casual bioperl user, being aware of their existence is useful since they are the basis to understanding how bioperl programs can communicate with other bioinformatics projects and computer languages such as Ensembl and biopython and biojava.
For more discussion of design and development issues please see the biodesign.pod file in the package or biodesign.html (http://bioperl.org/Core/Latest/biodesign.html).
Bioperl provides software modules for many of the typical tasks of bioinformatics programming. These include:
Accessing sequence data from local and remote databases
Transforming formats of database/ file records
Manipulating individual sequences
Searching for similar sequences
Creating and manipulating sequence alignments
Searching for genes and other structures on genomic DNA
Developing machine readable sequence annotations
The following sections describe how bioperl can help perform all of these tasks.
Much of bioperl is focused on sequence manipulation. However, before bioperl can manipulate sequences, it needs to have access to sequence data. Now one can directly enter data sequence data into a bioperl Seq object, eg:
$seq = Bio::Seq->new(-seq => 'actgtggcgtcaact', -desc => 'Sample Bio::Seq object', -display_id => 'something', -accession_number => 'accnum', -alphabet => 'dna' );
However, in most cases, it is preferable to access sequence data from some online data file or database. Note that in common with conventional bioinformatics usage we will sometimes call a "database" what might be more appropriately referred to as an "indexed flat file".
Bioperl supports accessing remote databases as well as creating indices for accessing local databases. There are two general approaches to accomplishing this. If you know what kind of database the sequences are stored in (i.e. flat file, local relational database or a database accessed remotely over the internet), you can write a script that specifically accesses data from that kind of database. This approach is described in sections III.1.1 and III.1.2 for access from remote databases and local indexed flat files respectively. To explicitly access sequence data from a local relational database requires installing and setting up the modules in the bioperl-db library and the BioSQL schema, see "IV.3" for more information.
The other approach is to use the recently developed OBDA (Open Bioinformatics Data Access) Registry system. Using OBDA it is possible to import sequence data from a database without your needing to know whether the required database is flat-file or relational or even whether it is local or accessible only over the net. Descriptions of how to set up the necessary registry configuration file and access sequence data with the registry in described in BIODATABASE_ACCESS in the doc/howto subdirectory and won't be repeated here.
Accessing sequence data from the principal molecular biology databases is straightforward in bioperl. Data can be accessed by means of the sequence's accession number or id. Batch mode access is also supported to facilitate the efficient retrieval of multiple sequences. For retrieving data from genbank, for example, the code could be as follows:
$gb = new Bio::DB::GenBank(); # this returns a Seq object : $seq1 = $gb->get_Seq_by_id('MUSIGHBA1'); # this returns a Seq object : $seq2 = $gb->get_Seq_by_acc('AF303112'); # this returns a SeqIO object : $seqio = $gb->get_Stream_by_id(["J00522","AF303112","2981014"]);
See section "III.2.1" for information on using this SeqIO object.
Bioperl currently supports sequence data retrieval from the genbank, genpept, RefSeq, swissprot, and EMBL databases. See Bio::DB::GenBank, Bio::DB::GenPept, Bio::DB::SwissProt, Bio::DB::RefSeq and Bio::DB::EMBL for more information. A user can also specify a different database mirror for a database - this is especially relevent for the SwissProt resource where there are many ExPaSy mirrors. There are also configuration options for specifying local proxy servers for those behind firewalls.
The retrieval of NCBI RefSeqs sequences is supported through a special module called Bio::DB::RefSeq which actually queries an EBI server. Please see Bio::DB::RefSeq before using it as there are some caveats with RefSeq retrieval. RefSeq ids in Genbank begin with "NT_", "NC_", "NG_", "NM_", "NP_", "XM_", "XR_", or "XP_" (for more information see http://www.ncbi.nlm.nih.gov/LocusLink/refseq.html). Bio::DB::GenBank can be used to retrieve entries corresponding to these ids but bear in mind that these are not Genbank entries, strictly speaking. See Bio::DB::GenBank for special details on retrieving entries beginning with "NT_", these are specially formatted "CONTIG" entries.
Bioperl also supports retrieval from a remote Ace database. This capability requires the presence of the external AcePerl module. You need to download and install the aceperl module from http://stein.cshl.org/AcePerl/.
An additional module is available for accessing remote databases, BioFetch, which queries the dbfetch script at EBI. The available databases are EMBL, GenBank, or SWALL, and the entries can be retrieved in different formats as objects or streams (SeqIO objects), or as "tempfiles". See Bio::DB::BioFetch for the details.
Alternately, bioperl permits indexing local sequence data files by means of the Bio::Index or Bio::DB::Fasta objects. The following sequence data formats are supported by Bio::Index: genbank, swissprot, pfam, embl and fasta. Once the set of sequences have been indexed using Bio::Index, individual sequences can be accessed using syntax very similar to that described above for accessing remote databases. For example, if one wants to set up an indexed flat-file database of fasta files, and later wants then to retrieve one file, one could write scripts like:
# script 1: create the index use Bio::Index::Fasta; # using fasta file format use strict; # some users have reported that this is necessary my $Index_File_Name = shift; my $inx = Bio::Index::Fasta->new( -filename => $Index_File_Name, -write_flag => 1); $inx->make_index(@ARGV); # script 2: retrieve some files use Bio::Index::Fasta; use strict; # some users have reported that this is necessary my $Index_File_Name = shift; my $inx = Bio::Index::Fasta->new($Index_File_Name); foreach my $id (@ARGV) { my $seq = $inx->fetch($id); # Returns Bio::Seq object # do something with the sequence }
To facilitate the creation and use of more complex or flexible indexing systems, the bioperl distribution includes two sample scripts in the scripts/index directory, bp_index.PLS and bp_fetch.PLS. These scripts can be used as templates to develop customized local data-file indexing systems.
Bioperl also supplies Bio::DB::Fasta as a means to index and query Fasta format files. It's similar in spirit to Bio::Index::Fasta but offers more methods, e.g.
use Bio::DB::Fasta; use strict; my $db = Bio::DB::Fasta->new($file); # one file or many files my $seqstring = $db->seq($id); # get a sequence as string my $seqobj = $db->get_Seq_by_id($id); # get a PrimarySeq obj my $desc = $db->header($id); # get the header, or description line
See Bio::DB::Fasta for more information on this fully-featured module.
Both modules also offer the user the ability to designate a specific string within the fasta header as the desired id, such as the gi number within the string "gi|4556644|gb|X45555". Consider the following fasta-formatted sequence, in "test.fa":
>gi|523232|emb|AAC12345|sp|D12567 titin fragment MHRHHRTGYSAAYGPLKJHGYVHFIMCVVVSWWASDVVTYIPLLLNNSSAGWKRWWWIIFGGE GHGHHRTYSALWWPPLKJHGSKHFILCVKVSWLAKKERTYIPKKILLMMGGWWAAWWWI
By default Bio::Index::Fasta and Bio::DB::Fasta will use the first "word" they encounter in the fasta header as the retrieval key, in this case "gi|523232|emb|AAC12345|sp|D12567". What would be more useful as a key would be a single id. The code below will index the "test.fa" file and create an index file called "test.fa.idx" where the keys are the Swissprot, or "sp", identifiers.
$ENV{BIOPERL_INDEX_TYPE} = "SDBM_File"; # look for the index in the current directory $ENV{BIOPERL_INDEX} = "."; my $file_name = "test.fa"; my $inx = Bio::Index::Fasta->new( -filename => $file_name . ".idx", -write_flag => 1 ); # pass a reference to the critical function to the Bio::Index object $inx->id_parser(\&get_id); # make the index $inx->make_index($file_name); # here is where the retrieval key is specified sub get_id { my $header = shift; $header =~ /^>.*\bsp\|([A-Z]\d{5}\b)/; $1; }
Here is how you would retrieve the sequence, as a Bio::Seq object:
my $seq = $inx->fetch("D12567"); print $seq->seq;
What if you wanted to retrieve a sequence using either a Swissprot id or a gi number and the fasta header was actually a concatenation of headers with multiple gi's and Swissprots?
>gi|523232|emb|AAC12345|sp|D12567|gi|7744242|sp|V11223 titin fragment
Modify the function that's passed to the id_parser method:
sub get_id { my $header = shift; my (@sps) = $header =~ /^>.*\bsp\|([A-Z]\d{5})\b/g; my (@gis) = $header =~ /gi\|(\d+)\b/g; return (@sps,@gis); }
The Bio::DB::Fasta module uses the same principle, but the syntax is slightly different, for example:
my $db = Bio::DB::Fasta->new('test.fa', -makeid=>\&make_my_id); my $seqobj = $db->get_Seq_by_id($id); sub make_my_id { my $description_line = shift; $description_line =~ /gi\|(\d+)\|emb\|(\w+)/; ($1,$2); }
The core bioperl installation does not support accessing sequences and data stored in relational databases. However, this capability is available with the auxiliary bioperl-db library. See section "IV.3" for more information.
A common - and tedious - bioinformatics task is that of converting sequence data among the many widely used data formats. Bioperl's SeqIO object, however, makes this chore a breeze. SeqIO can read a stream of sequences - located in a single or in multiple files - in a number of formats: Fasta, EMBL, GenBank, Swissprot, PIR, GCG, SCF, phd/phred, Ace, fastq, exp, chado, or raw (plain sequence). SeqIO can also parse tracefiles in alf, ztr, abi, ctf, and ctr format Once the sequence data has been read in with SeqIO, it is available to bioperl in the form of Seq, PrimarySeq, or RichSeq objects, depending on what the sequence source is. Moreover, the sequence objects can then be written to another file (again using SeqIO) in any of the supported data formats making data converters simple to implement, for example:
use Bio::SeqIO; $in = Bio::SeqIO->new(-file => "inputfilename", -format => 'Fasta'); $out = Bio::SeqIO->new(-file => ">outputfilename", -format => 'EMBL'); while ( my $seq = $in->next_seq() ) {$out->write_seq($seq); }
In addition, the perl "tied filehandle" syntax is available to SeqIO, allowing you to use the standard <> and print operations to read and write sequence objects, eg:
$in = Bio::SeqIO->newFh(-file => "inputfilename" , -format => 'fasta'); $out = Bio::SeqIO->newFh(-format => 'embl'); print $out $_ while <$in>;
If the "-format" argument isn't used then Bioperl will try to determine the format based on the file's suffix, in a case-insensitive manner. If there's no suffix available then SeqIO will attempt to guess the format based on actual content. Here is the current set of suffixes:
Format Suffixes Comment fasta fasta|fast|seq|fa|fsa|nt|aa Fasta genbank gb|gbank|genbank|gbs|gbk Genbank scf scf SCF tracefile pir pir PIR embl embl|ebl|emb|dat EMBL raw txt plain gcg gcg GCG ace ace ACeDB bsml bsm|bsml BSML XML game GAME XML swiss swiss|sp SwissProt phd phd|phred Phred fastq fastq Fastq Locuslink LL_tmpl format qual Phred quality file chado Chado XML exp exp Staden experiment file abi* abi ABI tracefile alf* alf ALF tracefile ctf* ctf CTF tracefile ztr* ztr ZTR tracefile pln* pln Staden plain tracefile
* These formats require the bioperl-ext package and the io_lib library from the Staden package
For more information see Bio::SeqIO or the SeqIO HOWTO (http://bioperl.org/HOWTOs/html/SeqIO.html).
Data files storing multiple sequence alignments also appear in varied formats. AlignIO is the bioperl object for conversion of alignment files. AlignIO is patterned on the SeqIO object and its commands have many of the same names as the commands in SeqIO. Just as in SeqIO the AlignIO object can be created with "-file" and "-format" options:
use Bio::AlignIO; my $io = Bio::AlignIO->new(-file => "receptors.aln", -format => "clustalw" );
If the "-format" argument isn't used then Bioperl will try and determine the format based on the file's suffix, in a case-insensitive manner. Here is the current set of suffixes:
Format Suffixes Comment bl2seq clustalw aln emboss* water|needle fasta fasta|fast|seq|fa|fsa|nt|aa mase Seaview mega meg|mega meme meme metafasta msf msf|pileup|gcg GCG nexus nexus|nex pfam pfam|pfm phylip phylip|phlp|phyl|phy|phy|ph interleaved prodom psi psi PSI-BLAST selex selex|slx|selx|slex|sx HMMER stockholm
*water, needle, matcher, stretcher, merger, and supermatcher See "IV.2.1" on EMBOSS for more information
Unlike SeqIO AlignIO cannot create output files in every format. AlignIO currently supports output in these 6 formats: fasta, mase, selex, clustalw, msf/gcg, and phylip (interleaved).
Another significant difference between AlignIO and SeqIO is that AlignIO handles IO for only a single alignment at a time but SeqIO.pm handles IO for multiple sequences in a single stream. Syntax for AlignIO is almost identical to that of SeqIO:
use Bio::AlignIO; $in = Bio::AlignIO->new(-file => "inputfilename" , -format => 'fasta'); $out = Bio::AlignIO->new(-file => ">outputfilename", -format => 'pfam'); while ( my $aln = $in->next_aln() ) { $out->write_aln($aln); }
The only difference is that the returned object reference, $aln, is to a SimpleAlign object rather than to a Seq object.
AlignIO also supports the tied filehandle syntax described above for SeqIO. See Bio::AlignIO, Bio::SimpleAlign, and section "III.5" on SimpleAlign for more information.
Bioperl contains many modules with functions for sequence analysis. And if you cannot find the function you want in bioperl you may be able to find it in EMBOSS or PISE , which are accessible through the bioperl-run auxiliary library (see "IV.2.1").
OK, so we know how to retrieve sequences and access them as sequence objects. Let's see how we can use sequence objects to manipulate our sequence data and retrieve information. Seq provides multiple methods for performing many common (and some not-so-common) tasks of sequence manipulation and data retrieval. Here are some of the most useful:
These methods return strings or may be used to set values:
$seqobj->display_id(); # the human read-able id of the sequence $seqobj->seq(); # string of sequence $seqobj->subseq(5,10); # part of the sequence as a string $seqobj->accession_number(); # when there, the accession number $seqobj->alphabet(); # one of 'dna','rna','protein' $seqobj->primary_id(); # a unique id for this sequence irregardless # of its display_id or accession number $seqobj->desc(); # a description of the sequence
It is worth mentioning that some of these values correspond to specific fields of given formats. For example, the display_id method returns the LOCUS name of a Genbank entry, the (\S+) following the > character in a Fasta file, the ID from a SwissProt file, and so on. The desc() method will return the DEFINITION line of a Genbank file, the line following the display_id in a Fasta file, and the DE field in a SwissProt file.
The following methods return an array of Bio::SeqFeature objects:
$seqobj->get_SeqFeatures; # The 'top level' sequence features $seqobj->get_all_SeqFeatures; # All sequence features, including sub- # seq features
For a comment annotation, you can use:
use Bio::Annotation::Comment; $seq->annotation->add_Annotation('comment', Bio::Annotation::Comment->new(-text => 'some description');
For a reference annotation, you can use:
use Bio::Annotation::Reference; $seq->annotation->add_Annotation('reference', Bio::Annotation::Reference->new(-authors => 'author1,author2', -title => 'title line', -location => 'location line', -medline => 998122 );
Sequence features will be discussed further in section "III.7" on machine-readable sequence annotation. A general description of the object can be found in Bio::SeqFeature::Generic, and a description of related, top-level annotation is found in Bio::Annotation::Collection.
Additional sample code for obtaining sequence features can be found in the script gb2features.pl in the subdirectory examples/DB. Finally, there's a HOWTO on features and annotations (http://bioperl.org/HOWTOs/html/Feature-Annotation.html) and there's a section on features in the FAQ (http://bioperl.org/Core/Latest/faq.html#5).
The following methods returns new sequence objects, but do not transfer the features from the starting object to the resulting feature:
$seqobj->trunc(5,10); # truncation from 5 to 10 as new object $seqobj->revcom; # reverse complements sequence $seqobj->translate; # translation of the sequence
Note that some methods return strings, some return arrays and some return objects. See Bio::Seq for more information.
Many of these methods are self-explanatory. However, bioperl's flexible translation methods warrant further comment. Translation in bioinformatics can mean two slightly different things:
The bioperl implementation of sequence-translation does the first of these tasks easily. Any sequence object which is not of alphabet 'protein' can be translated by simply calling the method which returns a protein sequence object:
$translation1 = $my_seq_object->translate;
However, the translate method can also be passed several optional parameters to modify its behavior. For example, the first two arguments to translate() can be used to modify the characters used to represent stop (default '*') and unknown amino acid ('X'). (These are normally best left untouched.) The third argument determines the frame of the translation. The default frame is "0". To get translations in the other two forward frames, we would write:
$translation2 = $my_seq_object->translate(undef,undef,1); $translation3 = $my_seq_object->translate(undef,undef,2);
The fourth argument to translate() makes it possible to use alternative genetic codes. There are currently 16 codon tables defined, including tables for 'Vertebrate Mitochondrial', 'Bacterial', 'Alternative Yeast Nuclear' and 'Ciliate, Dasycladacean and Hexamita Nuclear' translation. These tables are located in the object Bio::Tools::CodonTable which is used by the translate method. For example, for mitochondrial translation:
$human_mitochondrial_translation = $seq_obj->translate(undef,undef,undef,2);
If we want to translate full coding regions (CDS) the way major nucleotide databanks EMBL, GenBank and DDBJ do it, the translate method has to perform more tricks. Specifically, 'translate' needs to confirm that the sequence has appropriate start and terminator codons at the beginning and the end of the sequence and that there are no terminator codons present within the sequence. In addition, if the genetic code being used has an atypical (non-ATG) start codon, the translate method needs to convert the initial amino acid to methionine. These checks and conversions are triggered by setting the fifth argument of the translate method to evaluate to "true".
If argument 5 is set to true and the criteria for a proper CDS are not met, the method, by default, issues a warning. By setting the sixth argument to evaluate to "true", one can instead instruct the program to die if an improper CDS is found, e.g.
$protein_object = $cds->translate(undef,undef,undef,undef,1,'die_if_errors');
See Bio::Tools::CodonTable for related details.
In addition to the methods directly available in the Seq object, bioperl provides various helper objects to determine additional information about a sequence. For example, SeqStats object provides methods for obtaining the molecular weight of the sequence as well the number of occurrences of each of the component residues (bases for a nucleic acid or amino acids for a protein.) For nucleic acids, SeqStats also returns counts of the number of codons used. For example:
use SeqStats; $seq_stats = Bio::Tools::SeqStats->new($seqobj); $weight = $seq_stats->get_mol_wt(); $monomer_ref = $seq_stats->count_monomers(); $codon_ref = $seq_stats->count_codons(); # for nucleic acid sequence
Note: sometimes sequences will contain ambiguous codes. For this reason, get_mol_wt() returns a reference to a two element array containing a greatest lower bound and a least upper bound of the molecular weight.
The SeqWords object is similar to SeqStats and provides methods for calculating frequencies of "words" (e.g. tetramers or hexamers) within the sequence. See Bio::Tools::SeqStats and Bio::Tools::SeqWords for more information.
Another common sequence manipulation task for nucleic acid sequences is locating restriction enzyme cutting sites. Bioperl provides the Bio::Restriction::Enzyme, Bio::Restriction::EnzymeCollection, and Bio::Restriction::Analysis objects for this purpose. These modules replace the older module Bio::Tools::RestrictionEnzyme. A new collection of enzyme objects would be defined like this:
use Bio::Restriction::EnzymeCollection; my $all_collection = Bio::Restriction::EnzymeCollection;
Bioperl's default Restriction::EnzymeCollection object comes with data for more than 500 different Type II restriction enzymes. A list of the available enzyme names can be accessed using the available_list() method, but these are just the names, not the functional objects. You also have access to enzyme subsets. For example to select all available Enzyme objects with recognition sites that are six bases long one could write:
my $six_cutter_collection = $all_collection->cutters(6); foreach my $enz ($six_cutter_collection){ print $enz->name,"\t",$enz->site,"\t",$enz->overhang_seq,"\n"; # prints name, recognition site, overhang }
There are other methods that can be used to select sets of enzyme objects, such as unique_cutters() and blunt_enzymes(). You can also select a Enzyme object by name, like so:
my $ecori_enzyme = $all_collection->get_enzyme('EcoRI');
Once an appropriate enzyme has been selected, the sites for that enzyme on a given nucleic acid sequence can be obtained using the fragments() method. The syntax for performing this task is:
use Bio::Restriction::Analysis; my $analysis = Bio::Restriction::Analysis->new(-seq => $seq); # where $seq is the Bio::Seq object for the DNA to be cut @fragments = $analysis->fragments($enzyme); # and @fragments will be an array of strings
To get information on isoschizomers, methylation sites, microbe source, vendor or availability you will need to create your EnzymeCollection directly from a REBASE file, like this:
use Bio::Restriction::IO; my $re_io = Bio::Restriction::IO->new(-file=>$file,-format=>'withrefm'); my $rebase_collection = $re_io->read;
A REBASE file in the correct format can be found at ftp://ftp.neb.com/pub/rebase, it will have a name like "withrefm.308". If need be you can also create new enzymes, like this:
my $re = new Bio::Restriction::Enzyme(-enzyme=>'BioRI',-seq=>'GG^AATTCC');
For more informatation see Bio::Restriction::Enzyme, Bio::Restriction::EnzymeCollection, Bio::Restriction::Analysis, and Bio::Restriction::IO.
For amino acid sequences we may be interested to know whether the amino acid sequence contains a cleavable signal sequence for directing the transport of the protein within the cell. SigCleave is a program (originally part of the EGCG molecular biology package) to predict signal sequences, and to identify the cleavage site based on the von Heijne algorithm.
The threshold setting controls the score reporting. If no value for threshold is passed in by the user, the code defaults to a reporting value of 3.5. SigCleave will only return score/position pairs which meet the threshold limit.
There are 2 accessor methods for this object. signals() will return a perl hash containing the sigcleave scores keyed by amino acid position. pretty_print() returns a formatted string similar to the output of the original sigcleave utility.
The syntax for using Sigcleave is as follows:
# create a Seq object, for example: $seqobj = Bio::Seq->new(-seq => "AALLHHHHHHGGGGPPRTTTTTVVVVVVVVVVVVVVV"); use Bio::Tools::Sigcleave; $sigcleave_object = new Bio::Tools::Sigcleave ( -seq => $seqobj, -threshold => 3.5, -desc => 'test sigcleave protein seq', -type => 'AMINO' ); %raw_results = $sigcleave_object->signals; $formatted_output = $sigcleave_object->pretty_print;
Note that the "type" in the Sigcleave object is "amino" whereas in a Seq object it would be called "protein". Please see Bio::Tools::Sigcleave for details.
OddCodes:
For some purposes it's useful to have a listing of an amino acid sequence showing where the hydrophobic amino acids are located or where the positively charged ones are. Bioperl provides this capability via the module Bio::Tools::OddCodes.
For example, to quickly see where the charged amino acids are located along the sequence we perform:
use Bio::Tools::OddCodes; $oddcode_obj = Bio::Tools::OddCodes->new($amino_obj); $output = $oddcode_obj->charge();
The sequence will be transformed into a three-letter sequence (A,C,N) for negative (acidic), positive (basic), and neutral amino acids. For example the ACDEFGH would become NNAANNC.
For a more complete chemical description of the sequence one can call the chemical() method which turns sequence into one with an 8-letter chemical alphabet { A (acidic), L (aliphatic), M (amide), R (aromatic), C (basic), H (hydroxyl), I (imino), S (sulfur) }:
$output = $oddcode_obj->chemical();
In this case the sample sequence ACDEFGH would become LSAARAC.
OddCodes also offers translation into alphabets showing alternate characteristics of the amino acid sequence such as hydrophobicity, "functionality" or grouping using Dayhoff's definitions. See the documentation in Bio::Tools::OddCodes for further details.
SeqPattern:
The SeqPattern object is used to manipulate sequences using perl regular expressions. A key motivation for SeqPattern is to have a way of generating a reverse complement of a nucleic acid sequence pattern that includes ambiguous bases and/or regular expressions. This capability leads to significant performance gains when pattern matching on both the sense and anti-sense strands of a query sequence are required. Typical syntax for using SeqPattern is shown below. For more information, there are several interesting examples in the script seq_pattern.pl in the examples/tools directory.
use Bio::Tools::SeqPattern; $pattern = '(CCCCT)N{1,200}(agggg)N{1,200}(agggg)'; $pattern_obj = new Bio::Tools::SeqPattern(-SEQ => $pattern, -TYPE => 'dna'); $pattern_obj2 = $pattern_obj->revcom(); $pattern_obj->revcom(1); # returns expanded rev complement pattern.
More detail can be found in Bio::Tools::SeqPattern.
Coordinate system conversion is a common requirement, for example, when one wants to look at the relative positions of sequence features to one another and convert those relative positions to absolute coordinates along a chromosome or contig. Although coordinate conversion sounds pretty trivial it can get fairly tricky when one includes the possibilities of switching to coordinates on negative (i.e. Crick) strands and/or having a coordinate system terminate because you have reached the end of a clone or contig. Bioperl has two different approaches to coordinate-system conversion (based on the modules Bio::Coordinate::Pair and Bio::DB::GFF::RelSegment, respectively).
The Coordinate::Pair approach is somewhat more "low level". With it, you define an input coordinate system and an output coordinate system, where in each case a coordinate system is a triple of a start position, end position and strand. The end position is especially important when dealing with unfinished assemblies where the coordinate system ends when one reaches the end of the sequence of a clone or contig. Once one has defined the two coordinate systems, one defines a Coordinate::Pair to map between them. Then one can map positions between the coordinates systems with code such as this:
$input_coordinates = Bio::Location::Simple->new (-seq_id => 'propeptide', -start => 1000, -end => 2000, -strand=>1 ); $output_coordinates = Bio::Location::Simple->new (-seq_id => 'peptide', -start => 1100, -end => 2100, -strand=>1 ); $pair = Bio::Coordinate::Pair->new (-in => $input_coordinates , -out => $output_coordinates ); $pos = Bio::Location::Simple->new (-start => 500, -end => 500 ); $res = $pair->map($pos); $converted_start = $res->start;
In this example $res is also a Bio::Location object, as you'd expect. See the documentation for Bio::Coordinate::Pair and Bio::Coordinate::GeneMapper for more details.
The Bio::DB::GFF::RelSegment approach is designed more for handling coordinate transformations of sequence features rather than for transforming arbitrary coordinate systems. With Bio::DB::GFF::RelSegment you define a coordinate system relative to a specific feature (called the "refseq"). You also have access to the absolute coordinate system (typically of the entire chromosome.) You can determine the position of a feature relative to some other feature simply by redefining the relevant reference feature (i.e. the "refseq") with code like this:
$db = Bio::DB::GFF->new(-dsn => 'dbi:mysql:elegans', -adaptor => 'dbi:mysqlopt'); $segment = $db->segment('ZK909'); $relative_start = $segment->start; # $relative_start = 1; # Now retrieve the start position of ZK909 relative to feature ZK337 $segment->refseq('ZK337'); $relative_start = $segment->start; # Now retrieve the start position of ZK909 relative to the entire chromosome $absolute_start = $segment->abs_start;
This approach is convenient because you don't have to keep track of coordinates directly, you just keep track of the name of a feature which in turn marks the coordinate-system origin. However, this approach does require that you have stored all the sequence features in GFF format. Moreover, Bio::DB::GFF::RelSegment has been principally developed and tested for applications where all the sequence features are stored in a Bioperl-db relational database. However, if one wants to use the Bio:DB::GFF machinery (including its coordinate transformation capabilities) without building a local relational database, this is possible by defining the 'database' as having an adaptor called 'memory', e.g.
$db = Bio::DB::GFF->new( '-adaptor' => 'memory' );
For more details on coordinate transformations and other GFF-related capabilities in Bioperl see Bio::DB::GFF::RelSegment, Bio::DB::GFF, and the test file t/BioDBGFF.t.
One of the basic tasks in molecular biology is identifying sequences that are, in some way, similar to a sequence of interest. The Blast programs, originally developed at the NCBI, are widely used for identifying such sequences. The bioperl and bioperl-run packages offer a number of modules to facilitate running Blast as well as to parse the often voluminous reports produced by Blast.
Bioperl supports remote execution of blasts at NCBI by means of the RemoteBlast object.
A skeleton script to run a remote blast might look as follows:
$remote_blast = Bio::Tools::Run::RemoteBlast->new ( -prog => 'blastp',-data => 'ecoli',-expect => '1e-10' ); $r = $remote_blast->submit_blast("t/data/ecolitst.fa"); while (@rids = $remote_blast->each_rid ) { foreach $rid ( @rids ) {$rc = $remote_blast->retrieve_blast($rid);} }
You may want to change some parameter of the remote job and this example shows how to change the matrix:
$Bio::Tools::Run::RemoteBlast::HEADER{'MATRIX_NAME'} = 'BLOSUM25';
For a description of the many CGI parameters see:
http://www.ncbi.nlm.nih.gov/BLAST/Doc/urlapi.html
Note that the script has to be broken into two parts. The actual Blast submission and the subsequent retrieval of the results. At times when the NCBI Blast is being heavily used, the interval between when a Blast submission is made and when the results are available can be substantial.
The object $rc would contain the blast report that could then be parsed with Bio::Tools::BPlite or Bio::SearchIO. The default object returned is SearchIO after version 1.0. The object type can be changed using the -readmethod parameter but bear in mind that the favored Blast parser is Bio::SearchIO, others won't be supported in later versions.
Note that to make this script actually useful, one should add details such as checking return codes from the Blast to see if it succeeded and a "sleep" loop to wait between consecutive requests to the NCBI server. See example 22 in the demonstration script in the appendix to see some working code you could use, or Bio::Tools::Run::RemoteBlast for details.
It should also be noted that the syntax for creating a remote blast factory is slightly different from that used in creating StandAloneBlast, Clustalw, and T-Coffee factories. Specifically RemoteBlast requires parameters to be passed with a leading hyphen, as in '-prog' => 'blastp', while the other programs do not pass parameters with a leading hyphen.
No matter how Blast searches are run (locally or remotely, with or without a perl interface), they return large quantities of data that are tedious to sift through. Bioperl offers several different objects - Search.pm/SearchIO.pm, and BPlite.pm (along with its minor modifications, BPpsilite and BPbl2seq) for parsing Blast reports. Search and SearchIO which are the principal Bioperl interfaces for Blast and FASTA report parsing, are described in this section. The older BPlite is described in section "III.4.3". We recommend you use SearchIO, it's certain to be supported in future releases.
The Search and SearchIO modules provide a uniform interface for parsing sequence-similarity-search reports generated by BLAST (in standard and BLAST XML formats), PSI-BLAST, RPS-BLAST, bl2seq and FASTA. The SearchIO modules also provide a parser for HMMER reports and in the future, it is envisioned that the Search/SearchIO syntax will be extended to provide a uniform interface to an even wider range of report parsers including parsers for Genscan.
Parsing sequence-similarity reports with Search and SearchIO is straightforward. Initially a SearchIO object specifies a file containing the report(s). The method next_result reads the next report into a Search object in just the same way that the next_seq method of SeqIO reads in the next sequence in a file into a Seq object.
Once a report (i.e. a SearchIO object) has been read in and is available to the script, the report's overall attributes (e.g. the query) can be determined and its individual hits can be accessed with the next_hit method. Individual high-scoring segment pairs for each hit can then be accessed with the next_hsp method. Except for the additional syntax required to enable the reading of multiple reports in a single file, the remainder of the Search/SearchIO parsing syntax is very similar to that of the BPlite object it is intended to replace. Sample code to read a BLAST report might look like this:
# Get the report $searchio = new Bio::SearchIO (-format => 'blast', -file => $blast_report); $result = $searchio->next_result; # Get info about the entire report $result->database_name; $algorithm_type = $result->algorithm; # get info about the first hit $hit = $result->next_hit; $hit_name = $hit->name ; # get info about the first hsp of the first hit $hsp = $hit->next_hsp; $hsp_start = $hsp->query->start;
For more details there is a good description of how to use SearchIO at http://www.bioperl.org/HOWTOs/html/SearchIO.html or in the docs/howto subdirectory of the distribution. Additional documentation can be found in Bio::SearchIO::blast, Bio::SearchIO::psiblast, Bio::SearchIO::blastxml, Bio::SearchIO::fasta, and Bio::SearchIO. There is also sample code in the examples/searchio directory which illustrates how to use SearchIO. And finally, there's a section with SearchIO questions in the FAQ (http://bioperl.org/Core/Latest/faq.html#3).
Bioperl's older BLAST report parsers - BPlite, BPpsilite, BPbl2seq and Blast.pm - are no longer supported but since legacy Bioperl scripts have been written which use these objects, they are likely to remain within Bioperl for some time.
Much of the user interface of BPlite is very similar to that of Search. However accessing the next hit or HSP uses methods called next_Sbjct and next_HSP, respectively - in contrast to Search's next_hit and next_hsp.
BPlite
The syntax for using BPlite is as follows where the method for retrieving hits is now called nextSbjct() (for "subject"), while the method for retrieving high-scoring-pairs is called nextHSP():
use Bio::Tools::BPlite; $report = new Bio::Tools::BPlite(-fh=>\*STDIN); $report->query; while(my $sbjct = $report->nextSbjct) { $sbjct->name; while (my $hsp = $sbjct->nextHSP) { print $hsp->score,"\n"; } }
A complete description of the module can be found in Bio::Tools::BPlite.
BPpsilite
BPpsilite and BPbl2seq are objects for parsing (multiple iteration) PSIBLAST reports and Blast bl2seq reports, respectively. They are both minor variations on the BPlite object. See Bio::Tools::BPbl2seq and Bio::Tools::BPpsilite for details.
The syntax for parsing a multiple iteration PSIBLAST report is as shown below. The only significant additions to BPlite are methods to determine the number of iterated blasts and to access the results from each iteration. The results from each iteration are parsed in the same manner as a (complete) BPlite object.
use Bio::Tools::BPpsilite; $report = new Bio::Tools::BPpsilite(-fh=>\*STDIN); $total_iterations = $report->number_of_iterations; $last_iteration = $report->round($total_iterations) while(my $sbjct = $last_iteration ->nextSbjct) { $sbjct->name; while (my $hsp = $sbjct->nextHSP) {$hsp->score; } }
See Bio::Tools::BPpsilite for details.
BPbl2seq
BLAST bl2seq is a program for comparing and aligning two sequences using BLAST. Although the report format is similar to that of a conventional BLAST, there are a few differences. Consequently, the standard bioperl parser BPlite ia unable to read bl2seq reports directly. From the user's perspective, one difference between bl2seq and other blast reports is that the bl2seq report does not print out the name of the first of the two aligned sequences. Consequently, BPbl2seq has no way of identifying the name of one of the initial sequence unless it is explicitly passed to constructor as a second argument as in:
use Bio::Tools::BPbl2seq; $report = Bio::Tools::BPbl2seq->new(-file => "t/data/dblseq.out", -queryname => "ALEU_HORVU"); $hsp = $report->next_feature; $answer=$hsp->score;
In addition, since there will only be (at most) one subject (hit) in a bl2seq report one should use the method $report->next_feature, rather than $report->nextSbjct->nextHSP to obtain the next high scoring pair. See Bio::Tools::BPbl2seq for more details.
Blast.pm
The Bio::Tools::Blast parser has been removed from Bioperl as of version 1.1. Consequently, the BPlite parser (described in the section "III.4.3") or the Search/SearchIO parsers (section "III.4.2") should be used for BLAST parsing within bioperl. SearchIO is the preferred approach and will be formally supported in future releases.
Blast is not the only sequence-similarity-searching program supported by bioperl. HMMER is a Hidden Markov Model (HMM) program that (among other capabilities) enables sequence similarity searching, from http://hmmer.wustl.edu. Bioperl does not currently provide a perl interface for running HMMER. However, bioperl does provide 2 HMMER report parsers, the recommended SearchIO HMMER parser and an older parser called HMMER::Results.
SearchIO can parse reports generated both by the HMMER program hmmsearch - which searches a sequence database for sequences similar to those generated by a given HMM - and the program hmmpfam - which searches a HMM database for HMMs which match domains of a given sequence. Sample usage for parsing a hmmsearch report might be:
use Bio::SearchIO; $in = new Bio::SearchIO(-format => 'hmmer',-file => '123.hmmsearch'); while ( $res = $in->next_result ){ # get a Bio::Search::Result::HMMERResult object print $res->query_name, " for HMM ", $res->hmm_name, "\n"; while ( $hit = $res->next_hit ){ print $hit->name, "\n"; while ( $hsp = $hit->next_hsp ){ print "length is ", $hsp->length, "\n"; } } }
Purists may insist that the term "hsp" is not applicable to hmmsearch or hmmpfam results and they may be correct - this is an unintended consequence of using the flexible and extensible SearchIO approach. See Bio::Search::Result::HMMERResult for more information.
For documentation on the older, unsupported HMMER parser, look at Bio::Tools::HMMER::Results.
There are several reasons why one might want to run the Blast programs locally - speed, data security, immunity to network problems, being able to run large batch runs, wanting to use custom or proprietary databases, etc. The NCBI provides a downloadable version of blast in a stand-alone format, and running blast locally without any use of perl or bioperl is completely straightforward. However, there are situations where having a perl interface for running the blast programs locally is convenient.
The module Bio::Tools::Run::StandAloneBlast offers the ability to wrap local calls to blast from within perl. All of the currently available options of NCBI Blast (e.g. PSIBLAST, PHIBLAST, bl2seq) are available from within the bioperl StandAloneBlast interface. Of course, to use StandAloneBlast, one needs to have installed BLAST from NCBI locally as well as one or more blast-readable databases.
Basic usage of the StandAloneBlast.pm module is simple. Initially, a local blast factory object is created.
@params = (program => 'blastn', database => 'ecoli.nt'); $factory = Bio::Tools::Run::StandAloneBlast->new(@params);
Any parameters not explicitly set will remain as the BLAST defaults. Once the factory has been created and the appropriate parameters set, one can call one of the supported blast executables. The input sequence(s) to these executables may be fasta file(s), a Seq object or an array of Seq objects, eg
$input = Bio::Seq->new(-id =>"test query", -seq =>"ACTAAGTGGGGG"); $blast_report = $factory->blastall($input);
The returned blast report will be in the form of a bioperl parsed-blast object. The report object may be either a SearchIO, BPlite, BPpsilite, BPbl2seq or Blast object depending on the type of blast search - the SearchIO object is returned by default. The raw blast report is also available.
The syntax for running PHIBLAST, PSIBLAST and bl2seq searches via StandAloneBlast is also straightforward. See Bio::Tools::Run::StandAloneBlast documentation for details. In addition, the script standaloneblast.pl in the examples/tools directory contains descriptions of various possible applications of the StandAloneBlast object. This script shows how the blast report object can access the SearchIO blast parser directly, e.g.
while (my $hit = $blast_report->next_hit){ while (my $hsp = $sbjct->next_hsp){ print $hsp->score," ",$hit->name,"\n"; } }
See the sections "III.4.2" and "III.4.3" for more details on parsing BLAST reports.
Once one has identified a set of similar sequences, one often needs to create an alignment of those sequences. Bioperl offers several perl objects to facilitate sequence alignment: pSW, Clustalw.pm, TCoffee.pm and the bl2seq option of StandAloneBlast. As of release 1.2 of bioperl, using these modules (except bl2seq) requires a bioperl auxiliary library (bioperl-ext for pSW, bioperl-run for the others) and are therefore described in section IV. Here we describe only the module within the bioperl core package for manipulating previously created alignments, namely the SimpleAlign module.
The script aligntutorial.pl in the examples/align/ subdirectory is another good source of information of ways to create and manipulate sequence alignments within bioperl.
SimpleAlign objects are produced by bioperl-run alignment creation objects (e.g. Clustalw.pm, BLAST's bl2seq, TCoffee.pm, Lagan.pm, or pSW from the bioperl-ext package) or they can be read in from files of multiple-sequence alignments in various formats using AlignIO.
Some of the manipulations possible with SimpleAlign include:
slice(): Obtaining an alignment "slice", that is, a subalignment inclusive of specified start and end columns. Sequences with no residues in the slice are excluded from the new alignment and a warning is printed.
column_from_residue_number(): Finding column in an alignment where a specified residue of a specified sequence is located.
consensus_string(): Making a consensus string. This method includes an optional threshold parameter, so that positions in the alignment with lower percent-identity than the threshold are marked by "?"'s in the consensus
percentage_identity(): A fast method for calculating the average percentage identity of the alignment
consensus_iupac(): Making a consensus using IUPAC ambiguity codes from DNA and RNA.
Skeleton code for using some of these features is shown below. More detailed, working code is in bptutorial.pl example 13 and in align_on_codons.pl in the examples/align directory. Additional documentation on methods can be found in Bio::SimpleAlign and Bio::LocatableSeq.
use Bio::SimpleAlign; $aln = Bio::SimpleAlign->new('t/data/testaln.dna'); $threshold_percent = 60; $consensus_with_threshold = $aln->consensus_string($threshold_percent); $iupac_consensus = $aln->consensus_iupac(); # dna/rna alignments only $percent_ident = $aln->percentage_identity; $seqname = '1433_LYCES'; $pos = $aln->column_from_residue_number($seqname, 14);
Automated searching for putative genes, coding sequences, sequence-tagged-sites (STS's) and other functional units in genomic and expressed sequence tag (EST) data has become very important as the available quantity of sequence data has rapidly increased. Many feature searching programs currently exist. Each produces reports containing predictions that must be read manually or parsed by automated report readers.
Parsers for six widely used gene prediction programs - Genscan, Sim4, Genemark, Grail, ESTScan and MZEF - are available in bioperl. The interfaces for these parsers are all similar. We illustrate the usage for Genscan and Sim4 here. The syntax is relatively self-explanatory; see Bio::Tools::Genscan, Bio::Tools::Genemark, Bio::Tools::Grail, Bio::Tools::ESTScan, Bio::Tools::MZEF, and Bio::Tools::Sim4::Results for further details.
use Bio::Tools::Genscan; $genscan = Bio::Tools::Genscan->new(-file => 'result.genscan'); # $gene is an instance of Bio::Tools::Prediction::Gene # $gene->exons() returns an array of Bio::Tools::Prediction::Exon objects while($gene = $genscan->next_prediction()) { @exon_arr = $gene->exons(); } $genscan->close();
See Bio::Tools::Prediction::Gene and Bio::Tools::Prediction::Exon for more details.
use Bio::Tools::Sim4::Results; $sim4 = new Bio::Tools::Sim4::Results(-file => 't/data/sim4.rev', -estisfirst => 0); # $exonset is-a Bio::SeqFeature::Generic with Bio::Tools::Sim4::Exons # as sub features $exonset = $sim4->next_exonset; @exons = $exonset->sub_SeqFeature(); # $exon is-a Bio::SeqFeature::FeaturePair $exon = 1; $exonstart = $exons[$exon]->start(); $estname = $exons[$exon]->est_hit()->seqname(); $sim4->close();
See Bio::SeqFeature::Generic and Bio::Tools::Sim4::Exons for more information.
A parser for the ePCR program is also available. The ePCR program identifies potential PCR-based sequence tagged sites (STSs) For more details see the documentation in Bio::Tools::EPCR. A sample skeleton script for parsing an ePCR report and using the data to annotate a genomic sequence might look like this:
use Bio::Tools::EPCR; use Bio::SeqIO; $parser = new Bio::Tools::EPCR(-file => 'seq1.epcr'); $seqio = new Bio::SeqIO(-format => 'fasta', -file => 'seq1.fa'); $seq = $seqio->next_seq; while( $feat = $parser->next_feature ) { # add EPCR annotation to a sequence $seq->add_SeqFeature($feat); }
Historically, annotations for sequence data have been entered and read manually in flat-file or relational databases with relatively little concern for machine readability. More recent projects - such as EBI's ENSEMBL project and the efforts to develop an XML molecular biology data specification - have begun to address this limitation. Because of its strengths in text processing and regular-expression handling, perl is a natural choice for the computer language to be used for this task. And bioperl offers numerous tools to facilitate this process - several of which are described in the following sub-sections.
In Bioperl, most sequence annotations are stored in sequence-feature (SeqFeature) objects, where the SeqFeature object is associated with a parent Seq object. A SeqFeature object generally has a description (e.g. "exon", "promoter"), a location specifying its start and end positions on the parent sequence, and a reference to its parent sequence. In addition, a Seq object can also have an Annotation object associated with it, which could be used to store database links, literature references and comments. Creating a new SeqFeature and Annotation and associating it with a Seq is accomplished with syntax like:
$feat = new Bio::SeqFeature::Generic(-start => 40, -end => 80, -strand => 1, -primary => 'exon', -source => 'internal' ); $seqobj->add_SeqFeature($feat); # Add the SeqFeature to the Seq object
Once the features and annotations have been associated with the Seq, they can be with retrieved, eg:
@topfeatures = $seqobj->top_SeqFeatures(); # just top level, or @allfeatures = $seqobj->all_SeqFeatures(); # descend into sub features $disease_annotation = $annotations->get_Annotations('disease');
The individual components of a SeqFeature can also be set or retrieved with methods including:
# methods which return numbers $feat->start # start position $feat->end # end position $feat->strand # 1 means forward, -1 reverse, 0 not relevant # methods which return strings $feat->primary_tag # the main 'name' of the sequence feature, # eg, 'exon' $feat->source_tag # where the feature comes from, e.g. 'BLAST' # methods which return Bio::PrimarySeq objects $feat->seq # the sequence between start and end $feat->entire_seq # the entire sequence $feat->spliced_seq # the "joined" sequence, when there are # multiple sub-locations # other useful methods include $feat->overlap($other) # do SeqFeature $feat and SeqFeature $other overlap? $feat->contains($other) # is $other completely within $feat? $feat->equals($other) # do $feat and $other completely agree? $feat->sub_SeqFeatures # create/access an array of subsequence features
It is worth mentioning that one can also retrieve the start and end positions of a feature using a Bio::LocationI object:
$location = $feat->location # $location is a Bio::LocationI object $location->start; # start position $location->end; # end position
This is useful because one can use a Bio::Location::SplitLocationI object in order to retrieve the split coordinates inside the Genbank or EMBL join() statements (e.g. "CDS join(51..142,273..495,1346..1474)"):
if ( $feat->location->isa('Bio::Location::SplitLocationI') && $feat->primary_tag eq 'CDS' ) { foreach $loc ( $feat->location->sub_Location ) { print $loc->start,"..",$loc->end,"\n"; } }
See Bio::LocationI and Bio::Location::SplitLocationI for more information.
If more detailed information is required than is currently available in Seq objects the RichSeq object may be used. It is applicable in particular to database sequences (EMBL, GenBank and Swissprot) with detailed annotations. Sample usage might be:
@secondary = $richseq->get_secondary_accessions; $division = $richseq->division; @dates = $richseq->get_dates; $seq_version = $richseq->seq_version;
See Bio::Seq::RichSeqI for more details.
Much of the interesting description of a sequence can be associated with sequence features but in sequence objects derived from Genbank or EMBL entries there can be useful information in other "annotation" sections, such as the COMMENTS section of a Genbank entry. In order to access this information you'll need to create an Annotation::Collection object. For example:
$db = new Bio::DB::GenBank; $seqobj = $db->get_Seq_by_acc("NM_125788"); $ann_coll = $seqobj->annotation;
This Collection object is just a container for other specialized objects, and its methods are described in Bio::Annotation::Collection. You can find the desired object within the Collection object by examining the "tagnames":
foreach $ann ($ann_coll->get_Annotations) { print "Comment: ",$ann->as_text if ($ann->tagname eq "comment"); }
Other possible tagnames include "date_changed", "keyword", and "reference". Objects with the "reference" tagname are Bio::Annotation::Reference objects and represent scientific articles. See Bio::Annotation::Reference for descriptions of the methods used to access the data in Reference objects. There is also a HOWTO on features and annotation (http://bioperl.org/HOWTOs/html/Feature-Annotation.html).
Very large sequences present special problems to automated sequence-annotation storage and retrieval projects. Bioperl's LargeSeq object addresses this situation.
A LargeSeq object is a SeqI compliant object that stores a sequence as a series of files in a temporary directory (see sect "II.1" or Bio::SeqI for a definition of SeqI objects). The aim is to enable storing very large sequences (e.g. > 100 MBases) without running out of memory and, at the same time, preserving the familiar bioperl Seq object interface. As a result, from the user's perspective, using a LargeSeq object is almost identical to using a Seq object. The principal difference is in the format used in the SeqIO calls. Another difference is that the user must remember to only read in small chunks of the sequence at one time. These differences are illustrated in the following code:
$seqio = new Bio::SeqIO(-format => 'largefasta', -file => 't/data/genomic-seq.fasta'); $pseq = $seqio->next_seq(); $plength = $pseq->length(); $last_4 = $pseq->subseq($plength-3,$plength); # this is OK # On the other hand, the next statement would # probably cause the machine to run out of memory # $lots_of_data = $pseq->seq(); # NOT OK for a large LargeSeq object
Data files with sequences that are frequently being updated present special problems to automated sequence-annotation storage and retrieval projects. Bioperl's LiveSeq object is designed to address this situation.
The LiveSeq object addresses the need for a sequence object capable of handling sequence data that may be changing over time. In such a sequence, the precise locations of features along the sequence may change. LiveSeq deals with this issue by re-implementing the sequence object internally as a "double linked chain." Each element of the chain is connected to other two elements (the PREVious and the NEXT one). There is no absolute position like in an array, hence if positions are important, they need to be computed (methods are provided). Otherwise it's easy to keep track of the elements with their "LABELs". There is one LABEL (think of it as a pointer) to each ELEMENT. The labels won't change after insertions or deletions of the chain. So it's always possible to retrieve an element even if the chain has been modified by successive insertions or deletions.
Although the implementation of the LiveSeq object is novel, its bioperl user interface is unchanged since LiveSeq implements a PrimarySeqI interface (recall PrimarySeq is the subset of Seq without annotations or SeqFeatures - see section "II.1" or Bio::PrimarySeq). Consequently syntax for using LiveSeq objects is familiar although a modified version of SeqIO called Bio::LiveSeq::IO::Bioperl needs to be used to actually load the data, e.g.:
$loader = Bio::LiveSeq::IO::BioPerl->load(-db => "EMBL", -file => "t/data/factor7.embl"); $gene = $loader->gene2liveseq(-gene_name => "factor7"); $id = $gene->get_DNA->display_id ; $maxstart = $gene->maxtranscript->start;
See Bio::LiveSeq::IO::BioPerl for more details.
A Mutation object allows for a basic description of a sequence change in the DNA sequence of a gene. The Mutator object takes in mutations, applies them to a LiveSeq gene and returns a set of Bio::Variation objects describing the net effect of the mutation on the gene at the DNA, RNA and protein level.
The objects in Bio::Variation and Bio::LiveSeq directory were originally designed for the "Computational Mutation Expression Toolkit" project at European Bioinformatics Institute (EBI). The result of using them to mutate a gene is a holder object, 'SeqDiff', that can be printed out or queried for specific information. For example, to find out if restriction enzyme changes caused by a mutation are exactly the same in DNA and RNA sequences, we can write:
use Bio::LiveSeq::IO::BioPerl; use Bio::LiveSeq::Mutator; use Bio::LiveSeq::Mutation; $loader = Bio::LiveSeq::IO::BioPerl->load(-file => "$filename"); $gene = $loader->gene2liveseq(-gene_name => $gene_name); $mutation = new Bio::LiveSeq::Mutation (-seq =>'G', -pos => 100 ); $mutate = Bio::LiveSeq::Mutator->new(-gene => $gene, -numbering => "coding" ); $mutate->add_Mutation($mutation); $seqdiff = $mutate->change_gene(); $DNA_re_changes = $seqdiff->DNAMutation->restriction_changes; $RNA_re_changes = $seqdiff->RNAChange->restriction_changes; $DNA_re_changes eq $RNA_re_changes or print "Different!\n";
For a complete working script, see the change_gene.pl script in the examples/liveseq directory. For more details on the use of these objects see Bio::LiveSeq::Mutator and Bio::LiveSeq::Mutation as well as the original documentation for the "Computational Mutation Expression Toolkit" project at http://www.ebi.ac.uk/mutations/toolkit/.
SeqWithQuality objects are used to describe sequences with very specific annotations - that is, data quality annotations. Data quality information is important for documenting the reliability of base calls in newly sequenced or otherwise questionable sequence data. The quality data is contained within a Bio::Seq::PrimaryQual object. Syntax for using SeqWithQuality objects is as follows:
# first, make a PrimarySeq object $seqobj = Bio::PrimarySeq->new( -seq => 'atcgatcg', -id => 'GeneFragment-12', -accession_number => 'X78121', -alphabet => 'dna'); # now make a PrimaryQual object $qualobj = Bio::Seq::PrimaryQual->new(-qual => '10 20 30 40 50 50 20 10', -id => 'GeneFragment-12', -accession_number => 'X78121', -alphabet => 'dna'); # now make the SeqWithQuality object $swqobj = Bio::Seq::SeqQithQuality->new(-seq => $seqobj, -qual => $qualobj); # Now we access the sequence with quality object $swqobj->id(); # the id of the SeqWithQuality object may not match the # id of the sequence or of the quality $swqobj->seq(); # the sequence of the SeqWithQuality object $swqobj->qual(); # the quality of the SeqWithQuality object
A SeqWithQuality object is created automatically when phred output, a *phd file, is read by SeqIO, e.g.
$seqio = Bio::SeqIO->new(-file=>"my.phd",-format=>"phd"); # or just 'Bio::SeqIO->new(-file=>"my.phd")' $seqWithQualObj = $seqio->next_seq;
See Bio::Seq::SeqWithQuality for a detailed description of the methods, Bio::Seq::PrimaryQual, and Bio::SeqIO::phd.
The previous subsections have described tools for automated sequence annotation by the creation of an object layer on top of a traditional database structure. XML takes a somewhat different approach. In XML, the data structure is unmodified, but machine readability is facilitated by using a data-record syntax with special flags and controlled vocabulary.
In order to transfer data with XML in biology, one needs an agreed upon a vocabulary of biological terms. Several of these have been proposed and bioperl has at least some support for three: GAME, BSML and AGAVE.
Once a vocabulary is agreed upon, it becomes possible to convert sequence XML sequence features can be turned into bioperl Annotation and SeqFeature objects. Conversely Seq object features and annotations can be converted to XML so that they become available to any other systems. Typical usage with GAME or BSML are shown below. No special syntax is required by the user. Note that some Seq annotation will be lost when using XML in this manner since generally XML does not support all the annotation information available in Seq objects.
$str = Bio::SeqIO->new(-file => 't/data/test.game', -format => 'game'); $seq = $str->next_primary_seq(); $id = $seq->id; @feats = $seq->all_SeqFeatures(); $first_primary_tag = $feats[0]->primary_tag; $str = Bio::SeqIO->new(-file => 'bsmlfile.xml', -format => 'bsml'); $seq = $str->next_primary_seq(); $id = $seq->id; @feats = $seq->all_SeqFeatures(); $first_primary_tag = $feats[0]->primary_tag;
Another format for transmitting machine-readable sequence-feature data is the Genome Feature Format (GFF). This file type is well suited to sequence annotation because it allows the ability to describe entries in terms of parent-child relationships (see http://www.sanger.ac.uk/software/GFF for details). Bioperl includes a parser for converting between GFF files and SeqFeature objects. Typical syntax looks like:
$gffio = Bio::Tools::GFF->new(-fh => \*STDIN, -gff_version => 2); # loop over the input stream while ($feature = $gffio->next_feature()) { # do something with feature } $gffio->close();
Further information can be found at Bio::Tools::GFF. Also see examples/tools/gff2ps.pl, examples/tools/gb_to_gff.pl, and the scripts in scripts/Bio-DB-GFF. Note: this module shouldn't be confused with the module Bio::DB::GFF which is for implementing relational databases when using bioperl-db.
Sequence alignments are not the only examples in which one might want to manipulate a group of sequences together. Such groups of related sequences are generally referred to as clusters. Examples include Unigene clusters and gene clusters resulting from clustering algorithms being applied to microarray data.
The bioperl Cluster and ClusterIO modules are available for handling sequence clusters. Currently, cluster input/output modules are available only for Unigene clusters. To read in a Unigene cluster (in the NCBI XML format) and then extract individual sequences for the cluster for manipulation might look like this:
my $stream = Bio::ClusterIO->new(-file => "Hs.data", -format => "unigene"); while ( my $in = $stream->next_cluster ) { print $in->unigene_id() . "\n"; while ( my $sequence = $in->next_seq ) { print $sequence->accession_number . "\n"; } }
See Bio::Cluster::UniGene for more details.
Though bioperl has its roots in describing and searching nucleotide and protein sequences it has also branched out into related fields of study, such as protein structure, phylogenetic trees and genetic maps.
A StructureIO object can be created from one or more 3D structures represented in Protein Data Bank, or pdb, format (see http://www.rcsb.org/pdb for details).
StructureIO objects allow access to a variety of related Bio:Structure objects. An Entry object consist of one or more Model objects, which in turn consist of one or more Chain objects. A Chain is composed of Residue objects, which in turn consist of Atom objects. There's a wealth of methods, here are just a few:
$structio = Bio::Structure::IO->new( -file => "1XYZ.pdb"); $struc = $structio->next_structure; # returns an Entry object $pseq = $struc->seqres; # returns a PrimarySeq object, thus $pseq->subseq(1,20); # returns a sequence string @atoms = $struc->get_atoms($res); # Atom objects, given a Residue @xyz = $atom->xyz; # the 3D coordinates of the atom
These lines show how one has access to a number of related objects and methods. For examples of typical usage of these modules, see the scripts in the examples/structure subdirectory. Also see Bio::Structure::IO, Bio::Structure::Entry, Bio::Structure::Model, Bio::Structure::Chain, Bio::Structure::Residue, and Bio::Structure::Atom for more information.
Bioperl Tree objects can store data for all kinds of computer trees and are intended especially for phylogenetic trees. Nodes and branches of trees can be individually manipulated. The TreeIO object is used for stream I/O of tree objects. Currently only phylip/newick tree format is supported. Sample code might be:
$treeio = new Bio::TreeIO( -format => 'newick', -file => $treefile); $tree = $treeio->next_tree; # get the tree @nodes = $tree->get_nodes; # get all the nodes $tree->get_root_node->each_Descendent; # get descendents of root node
See Bio::TreeIO and Bio::Tree::Tree for details.
Using the Bio::Tools::Phylo::PAML module one can also parse the results of the PAML tree-building programs codeml, baseml, basemlg, codemlsites and yn00. See Bio::Tools::Phylo::PAML or the PAML HOWTO (http://bioperl.org/HOWTOs/html/PAML.html) for more information.
Bioperl map objects can be used to describe any type of biological map data including genetic maps, STS maps etc. Map I/O is performed with the MapIO object which works in a similar manner to the SeqIO, SearchIO and similar I/O objects described previously. In principle, Map I/O with various map data formats can be performed. However currently only mapmaker format is supported. Manipulation of genetic map data with Bioperl Map objects might look like this:
$mapio = new Bio::MapIO(-format => 'mapmaker', -file => $mapfile); $map = $mapio->next_map; # get a map $maptype = $map->type ; foreach $marker ( $map->each_element ) { $marker_name = $marker->name ; # get the name of each map marker }
See Bio::MapIO and Bio::Map::SimpleMap for more information.
Bio::Biblio objects are used to query bibliographic databases, such as MEDLINE. The associated modules are built to work with OpenBQS-compatible databases (see http://industry.ebi.ac.uk/openBQS). A Bio::Biblio object can execute a query like:
my $collection = $biblio->find ('brazma', 'authors'); while ( $collection->has_next ) { print $collection->get_next; }
See Bio::Biblio, the scripts/biblio/biblio.PLS script, or the examples/biblio/biblio_examples.pl script for more information.
A user may want to represent sequence objects and their SeqFeatures graphically. The Bio::Graphics::* modules use Perl's GD.pm module to create a PNG or GIF image given the SeqFeatures (Section "III.7.1") contained within a Seq object.
These modules contain numerous methods to dictate the sizes, colors, labels, and line formats within the image. For information see the excellent Graphics-HOWTO (http://bioperl.org/HOWTOs/html/Graphics-HOWTO.html) or in the docs/howto subdirectory. Additional documentation can be found in Bio::Graphics, Bio::Graphics::Panel, and in the scripts in the examples/biographics/ and scripts/graphics directories in the Bioperl package.
Bioperl modules use the standard extended single-letter genetic alphabets to represent nucleotide and amino acid sequences.
In addition to the standard alphabet, the following symbols are also acceptable in a biosequence:
? (a missing nucleotide or amino acid) - (gap in sequence)
(includes symbols for nucleotide ambiguity) ------------------------------------------ Symbol Meaning Nucleic Acid ------------------------------------------ A A Adenine C C Cytosine G G Guanine T T Thymine U U Uracil M A or C R A or G W A or T S C or G Y C or T K G or T V A or C or G H A or C or T D A or G or T B C or G or T X G or A or T or C N G or A or T or C IUPAC-IUB SYMBOLS FOR NUCLEOTIDE NOMENCLATURE: Cornish-Bowden (1985) Nucl. Acids Res. 13: 3021-3030.
------------------------------------------ Symbol Meaning ------------------------------------------ A Alanine B Aspartic Acid, Asparagine C Cystine D Aspartic Acid E Glutamic Acid F Phenylalanine G Glycine H Histidine I Isoleucine K Lysine L Leucine M Methionine N Asparagine P Proline Q Glutamine R Arginine S Serine T Threonine V Valine W Tryptophan X Unknown Y Tyrosine Z Glutamic Acid, Glutamine * Terminator IUPAC-IUP AMINO ACID SYMBOLS: Biochem J. 1984 Apr 15; 219(2): 345-373 Eur J Biochem. 1993 Apr 1; 213(1): 2
Beyond the bioperl "core" distribution which you get with the "minimal" installation, bioperl contains numerous other modules in so-called auxiliary libraries. These auxiliary libraries include bioperl-run, bioperl-db, bioperl-pipeline, bioperl-microarray and bioperl-ext among others. Generally, modules are placed in an auxiliary library if either:
The module requires the installation of additional non-standard external programs or modules, or
The module is perceived to be of interest to only a small percentage of the bioinformatics community
However there are exceptions and it is not always obvious whether a given module will be found in the "core" or in an auxiliary library.
At present, modules in the auxiliary packages can be obtained only by means of the CVS system. To browse through the auxiliary libraries and to obtain the download files, go to:
http://cvs.bioperl.org/cgi-bin/viewcvs/viewcvs.cgi/?cvsroot=bioperl
Install much like Bioperl:
>perl Makefile.PL
and
>make
Then:
>make test
Once the auxiliary library has been installed in this manner, the modules can be used in exactly the same manner as if they were in the bioperl core.
It possible to run various external (to Bioperl) sequence alignment and sequence manipulation programs via a perl interface using bioperl. However in most cases this requires having the bioperl-run auxiliary library (some cases may require bioperl-ext). Prior to bioperl release 1.2, many of these features were available within the bioperl "core" release. However currently some of the required modules have been transferred out of the core library. Some of the more commonly used of these modules are described in this section.
EMBOSS (European Molecular Biology Open Source Software) is an extensive collection of sequence analysis programs written in the C programming language, from http://www.uk.embnet.org/Software/EMBOSS. There are a number of algorithms in EMBOSS that are not found in "Bioperl proper" (e.g. calculating DNA melting temperature, finding repeats, identifying prospective antigenic sites) so if you cannot find the function you want in bioperl you might be able to find it in EMBOSS.
EMBOSS programs are usually called from the command line but the bioperl-run auxiliary library provides a Perl wrapper for EMBOSS function calls so that they can be executed from within a Perl script. Of course, the EMBOSS package as well as the bioperl-run must be installed in order for the Bioperl wrapper to function.
In the future, it is planned that Bioperl EMBOSS objects will return appropriate Bioperl objects to the calling script in addition to generating standard EMBOSS reports. This functionality is being initially implemented with the EMBOSS sequence alignment programs, so that they will return SimpleAlign objects in a manner similar to the way the Bioperl-run modules TCoffee.pm and Clustalw.pm work (see section "III.5" for a discussion of SimpleAlign).
An example of the Bioperl EMBOSS wrapper where a file is returned would be:
$factory = new Bio::Factory::EMBOSS; $compseqapp = $factory->program('compseq'); %input = ( -word => 4, -sequence => $seqObj, -outfile => $compseqoutfile ); $compseqapp->run(\%input); $seqio = Bio::SeqIO->new( -file => $compseqoutfile ); # etc...
Note that a Seq object was used as input. The EMBOSS object can also accept a file name as input, eg
-sequence => "inputfasta.fa"
Some EMBOSS programs will return strings, others will create files that can be read directly using Bio::SeqIO (section "III.2.1"), as in the example above. It's worth mentioning that another way to align sequences in bioperl is to run a program from the EMBOSS suite, such as 'matcher'. This can produce an output file that bioperl can read in using AlignIO:
my $factory = new Bio::Factory::EMBOSS; my $prog = $factory->program('matcher'); $prog->run({ -sequencea => Bio::Seq->new(-id => "seq1", -seq => $seqstr1), -sequenceb => Bio::Seq->new(-id => "seq2", -seq => $seqstr2), -aformat => "pair", -alternatives => 2, -outfile' => $outfile}); my $alignio_fmt = "emboss"; my $align_io = new Bio::AlignIO(-format => $alignio_fmt, -file => $outfile);
The Pise interface is another way of extending Bioperl's sequence analysis capabilities. To use EMBOSS programs within Bioperl you need to have EMBOSS locally installed, as well as the bioperl-run library. In contrast, with Pise you only need to install bioperl-run, since the actual analysis programs reside at the Pise site. Advantages of Pise include not having to load additional programs locally and having access to an extraordinary variety of programs, including EMBOSS. However Pise has the disadvantages of lower performance and decreased security since the data is transmitted over the net. Consider this example:
my $genscan = Pise::genscan->new( "http://bioweb.pasteur.fr/cgi-bin/seqanal/genscan.pl",'letondal@pasteur.fr', seq => $seq, parameter_file => "Arabidopsis.smat"); my $job = $genscan->run; my $parser = Bio::Tools::Genscan->new(-fh => $job->fh('genscan.out') ); while(my $gene = $parser->next_prediction()) { my $prot = $gene->predicted_protein; print $prot->seq, "\n"; }
Extremely simple! For more information on the Bioperl Pise interface see http://www-alt.pasteur.fr/~letondal/Pise/ or the documentation in the bioperl-run package.
As an alternative to Smith-Waterman, two sequences can also be aligned in Bioperl using the bl2seq option of Blast within the StandAloneBlast object. To get an alignment - in the form of a SimpleAlign object - using bl2seq, you need to parse the bl2seq report with the Bio::AlignIO file format reader as follows:
$factory = Bio::Tools::Run::StandAloneBlast->new('outfile' => 'bl2seq.out'); $bl2seq_report = $factory->bl2seq($seq1, $seq2); # Use AlignIO.pm to create a SimpleAlign object from the bl2seq report $str = Bio::AlignIO->new(-file => 'bl2seq.out', -format => 'bl2seq'); $aln = $str->next_aln();
For aligning multiple sequences (i.e. two or more), bioperl offers a perl interface to the bioinformatics-standard clustalw and tcoffee programs. Clustalw has been a leading program in global multiple sequence alignment (MSA) for several years. TCoffee is a relatively recent program - derived from clustalw - which has been shown to produce better results for local MSA.
To use these capabilities, the clustalw and/or tcoffee programs themselves need to be installed on the host system. In addition, the environmental variables CLUSTALDIR and TCOFFEEDIR need to be set to the directories containg the executables. See section "I.4" and the Bio::Tools::Run::Alignment::Clustalw and Bio::Tools::Run::Alignment::TCoffee for information on downloading and installing these programs.
From the user's perspective, the bioperl syntax for calling Clustalw.pm or TCoffee.pm is almost identical. The only differences are the names of the modules themselves appearing in the initial "use" and constructor statements and the names of the some of the individual program options and parameters.
In either case, initially, a factory object must be created. The factory may be passed most of the parameters or switches of the relevant program. In addition, alignment parameters can be changed and/or examined after the factory has been created. Any parameters not explicitly set will remain as the underlying program's defaults. Clustalw.pm/TCoffee.pm output is returned in the form of a SimpleAlign object. It should be noted that some Clustalw and TCoffee parameters and features (such as those corresponding to tree production) have not been implemented yet in the Perl interface.
Once the factory has been created and the appropriate parameters set, one can call the method align() to align a set of unaligned sequences, or profile_align() to add one or more sequences or a second alignment to an initial alignment. Input to align() consists of a set of unaligned sequences in the form of the name of file containing the sequences or a reference to an array of Seq objects. Typical syntax is shown below. (We illustrate with Clustalw.pm, but the same syntax - except for the module name - would work for TCoffee.pm)
use Bio::Tools::Run::Alignment::Clustalw; @params = ('ktuple' => 2, 'matrix' => 'BLOSUM'); $factory = Bio::Tools::Run::Alignment::Clustalw->new(@params); $ktuple = 3; $factory->ktuple($ktuple); # change the parameter before executing $seq_array_ref = \@seq_array; # where @seq_array is an array of Bio::Seq objects $aln = $factory->align($seq_array_ref);
Clustalw.pm/TCoffee.pm can also align two (sub)alignments to each other or add a sequence to a previously created alignment by using the profile_align method. For further details on the required syntax and options for the profile_align method, the user is referred to Bio::Tools::Run::Alignment::Clustalw and Bio::Tools::Run::Alignment::TCoffee. The user is also encouraged to examine the script clustalw.pl in the examples/align directory.
The Smith-Waterman (SW) algorithm is the standard method for producing an optimal local alignment of two sequences. Bioperl supports the computation of SW alignments via the pSW object with the auxiliary bioperl-ext library. Note that pSW only supports the alignment of protein sequences, not nucleotide.
The SW algorithm itself is implemented in C and incorporated into bioperl using an XS extension. This has significant efficiency advantages but means that pSW will not work unless you have compiled the bioperl-ext auxiliary library. If you have compiled the bioperl-ext package, usage is simple, where the method align_and_show displays the alignment while pairwise_alignment produces a (reference to) a SimpleAlign object.
use Bio::Tools::pSW; $factory = new Bio::Tools::pSW( '-matrix' => 'blosum62.bla', '-gap' => 12, '-ext' => 2, ); $factory->align_and_show($seq1, $seq2, STDOUT); $aln = $factory->pairwise_alignment($seq1, $seq2);
SW matrix, gap and extension parameters can be adjusted as shown. Bioperl comes standard with blosum62 and gonnet250 matrices. Others can be added by the user. For additional information on accessing the SW algorithm via pSW see the script psw.pl in the examples/tools directory and the documentation in Bio::Tools::pSW.
The bioperl-db package is intended to enable the easy access and manipulation of biology relational databases via a perl interface. Obviously it requires having administrative access to a relational database. Currently the bioperl-db interface is implemented to support databases in the Mysql, Postgres and Oracle formats. More details on bioperl-db can be found in the bioperl-db CVS directory at http://cvs.bioperl.org/cgi-bin/viewcvs/viewcvs.cgi/bioperl-db/?cvsroot=bioperl.
The database schema itself is not specified in the bioperl-db package but in the BioSQL package, available at http://obda.open-bio.org/. It is worth mentioning that most of the bioperl objects mentioned above map directly to tables in the Biosql schema. Therefore object data such as sequences, their features, and annotations can be easily loaded into the databases, as in
$loader->store($newid,$seqobj)
Similarly one can query the database in a variety of ways and retrieve arrays of Seq objects. See biodatabases.pod, Bio::DB::SQL::SeqAdaptor, Bio::DB::SQL::QueryConstraint, and Bio::DB::SQL::BioQuery for examples. The README file in the bioperl-db package has a helpful overview of the approach used in bioperl-db.
There a several other auxiliary libraries in the bioperl CVS repository including bioperl-microarray, bioperl-pedigree, bioperl-gui, bioperl-pipeline, bioperl-das-client and bioperl-corba-client. They are typically for specialized uses and/or require multiple external programs to run and/or are still pretty new and undeveloped. But if you have a need for any of these capabailities, it is easy to take a look at them at: http://cvs.bioperl.org/cgi-bin/viewcvs/viewcvs.cgi/?cvsroot=bioperl and see if they might be of use to you.
At numerous places in the tutorial, the reader is directed to the "documentation included with each of the modules." As was mentioned in the introduction, it is sometimes not easy in perl to determine the appropriate documentation because objects inherit methods from other objects (and the relevant documentation will be stored in the object from which the method was inherited.)
For example, say you wanted to find documentation on the parse() method of the module Genscan.pm. You would not find this documentation in the code for Genscan.pm, but rather in the code for AnalysisResult.pm from which Genscan.pm inherits the parse method!
So how would you know to look in AnalysisResult.pm for this documentation? The easy way is to use the special function "option 100" in the bptutorial script. Specifically if you run:
> perl -w bptutorial.pl 100 Bio::Tools::Genscan
you will receive the following output:
***Methods for Object Bio::Tools::Genscan ******** Methods taken from package Bio::Root::IO catfile close gensym new qualify qualify_to_ref rmtree tempdir tempfile ungensym Methods taken from package Bio::Root::RootI DESTROY stack_trace stack_trace_dump throw verbose warn Methods taken from package Bio::SeqAnalysisParserI carp confess croak next_feature Methods taken from package Bio::Tools::AnalysisResult analysis_date analysis_method analysis_method_version analysis_query analysis_subject parse Methods taken from package Bio::Tools::Genscan next_prediction
From this output, it is clear exactly from which object each method of Genscan.pm is taken, and, in particular that parse() is taken from the package Bio::Tools::AnalysisResult.
With this approach you can easily determine the source of any method in any bioperl object.
The following scripts demonstrate many of the features of bioperl. To run all the core demos, run:
> perl -w bptutorial.pl 0
To run a subset of the scripts do
> perl -w bptutorial.pl
and use the displayed help screen.
It may be best to start by just running one or two demos at a time. For example, to run the basic sequence manipulation demo, do:
> perl -w bptutorial.pl 1
Some of the later demos require that you have an internet connection and/or that you have an auxiliary bioperl library and/or external cpan module and/or external program installed. They may also fail if you are not running under Linux or Unix. In all of these cases, the script should fail gracefully simply saying the demo is being skipped. However if the script crashes, simply run the other demos individually (and perhaps send an email to bioperl-l@bioperl.org detailing the problem :-).
To install Bio::Seq, copy and paste the appropriate command in to your terminal.
cpanm
cpanm Bio::Seq
CPAN shell
perl -MCPAN -e shell install Bio::Seq
For more information on module installation, please visit the detailed CPAN module installation guide.