Phat (a Pretty Handy Annotation Tool)
Anthony Wirth, Simon Cawley and Terry Speed
(Questions and comments are welcome, email simon_cawley@affymetrix.com)
INTRODUCTION
Phat is a HMM-based genefinder, originally developed for genefinding in
Plasmodium falciparum. Care has been taken to make it easy to supply different
parameter files to run it on different organisms. This release comes with
parameter sets ready-made for:
Homo sapiens
Plasmodium falciparum
Plasmodium vivax
Phat has been tested on sequences up to just over 1M in length, and in principle
there is no limit to the size of the input sequence (though memory requirements
and execution time increase linearly with sequence length).
Phat comes with code to re-estimate parameters for other organisms.
(The code is in the subdirectory "Retrain" and some documentation may be
found in the file Retrain/DOC).
Phat comes in two flavours, one which operates on each strand independently
and one which processes both strands at the same time. These are called
halfphat and fullphat respectively. In the case of Plasmodium falciparum
fullphat performs significantly better. Coding regions on one strand look
similar to coding sequence on the reverse strand and so halfphat often predicts
a gene on each strand for each real gene present. Reducing these overlapping
predictions to a non-overlapping set is a difficult problem. Fullphat gets
around this by considering genes on both strands in its model. It takes
more time and memory to run, but comes up with better predictions.
UPDATES
release 20001212
- Changed format of the posterior coding probability files. The old
format gave a probability for every base being coding, which could be
a big file if the input sequence is big. The new format only records
the points at which the posterior probability changes. Since it doesn't
change much, this is a much more compact file. For backwards compatability,
the old format can still be produced by using the -oldprb option.
EXTRACTING AND COMPILING
After downloading the code, extract it by typing
gzcat phat.tar.gz | tar xvf -
cd into the new directory Phat and compile the code by typing
make
The code has been successfully compiled on a number of systems and
hopefully will compile with no trouble. If you experience problems
compiling please email me.
RUNNING
PART 1: running phat directly from the command line.
Note that there is a Perl script phat.pl which extends the functionality
of the program to include more input formats and multiple sequences per
file. See part 2 for details (requires both Perl and Bioperl)
Phat takes its parameters from a bunch of files all located in the same
directory. One such directory is supplied with the code, it is called
Pars.P.falciparum. (These parameters were estimated from the annotation
of chromosomes 2 and 3 of Plasmodium falciparum, as found in Genbank with
accessions NC_000910 and NT_002521)
There are two ways to tell the program which set of parameters to use.
One is to specify the directory on the command line:
fullphat -p Pars.P.falciparum myseq.fasta
halfphat -p Pars.P.falciparum myseq.gbk
Another way is to set the environment variable PHATPARS
setenv PHATPARS Pars.P.falciparum
fullphat myseq.fasta
The -p option will override the environment variable if it is already set.
As the examples imply, input sequences can be in fasta or genbank format.
If a multiple fasta file is passed in only the first sequence in the
file will be analyzed.
Code to estimate Phat parameters from annotated genbank files is currently
being developed, and will shortly be part of the release. The aim is to
make it easy for users to retrain Phat for their own particular uses.
AN EXAMPLE
There should be a directory called Sequence in the Phat directory, it
contains a genbank P. falciparum sequence of about 50kb in length
(AF030694.gbk), and a GFF file specifying the correct exon locations
in the file (AF030694.true.gff).
After compiling the code, cd to the Sequence directory and run fullphat
on the sequence:
../fullphat -p ../Pars.P.falciparum AF030694.gbk -noextras
This may take up to a minute or two, depending on the speed of the machine
you're using. The -noextras option limits the amount of output data
fullphat produces.
There should now be a file called AF030694.gff containing the fullphat
predictions. If you have a copy of gff2ps handy, you can make a postscript
plot of the predicted and true exons with the following command:
gff2ps -v AF030694.true.gff AF030694.gff > AF030694.ps
(download gff2ps from http://www1.imim.es/~jabril/GFFTOOLS/GFF2PS.html)
The plot is written to AF030694.ps and can be viewed with a postscript
viewer such as gv or ghostview.
gv AF030694.ps
PHAT OUTPUT
Phat writes a number of output files and prints a summary of the predicted
exons to the screen. The base name of the output files will be the same as
the base name of the input sequence. The output files are:
.rna A multiple fasta format file of predicted spliced mRNAs (without UTRs).
.pep A multiple fasta format file of predicted protein sequences.
.gff A GFF (General Features Format) file of the exon predictions.
(for a definition of GFF, see www.sanger.ac.uk/Software/formats/GFF/)
The GFF format is very useful, a bunch of tools for processing it
have been developed. In particular, there's an excellent free program
gff2ps for converting .gff formatted predictions to postscript plots.
(see http://www1.imim.es/~jabril/GFFTOOLS/GFF2PS.html)
.prb/.prf Probabilities of being coding for each base. Phat computes the
probability of each base being some kind of coding exon. In the future
a plot of these probabilities will be produced rather than a text file
of the values. These values offer another method identifying coding
regions, though they won't give a parse. They could be useful in
identifying candidate regions in which to design PCR primers. For
fullphat just one file (.prb) is produced. For halfphat, .prf and .prb
will contain the probabilities for the forward and backward strands
respectively.
Finally, phat will print output to the screen, something like the following:
Gene Exon Dir Type Start Stop Length Fr Ph Pr
1 1 + init 1849 3279 1431 0 0 0.83
1 2 + term 3417 3539 123 2 0 0.84
2 1 + init 3716 3832 117 1 0 0.54
2 2 + term 4103 4909 807 1 0 1.00
3 1 - sing 5852 8053 2202 1 0 1.00
4 1 - sing 8139 8597 459 2 0 0.79
5 1 - term 9741 10469 729 2 0 0.99
5 2 - init 10596 10625 30 2 0 0.96
6 1 - sing 11754 12377 624 2 0 0.93
7 1 + sing 16000 18597 2598 0 0 0.76
There is one line for each predicted exon. The columns are defined as follows:
Gene: Index of the gene to which the exon belongs.
Exon: Index of the exon in its gene.
Dir: Direction, + for the direct strand, - for the reverse.
Type: One of sing/init/intl/term for a single/initial/internal/terminal exon.
Start: The start coordinate of the exon.
Stop: The stop coordinate of the exon. (Start is always less that Stop).
Length: Exon length = Stop-Start+1
Fr: Frame: if the "rightmost" base of one of the exon's codons ends at
position k, the frame is (k mod 3).
Phase: The number of bases to skip in the exon (going 5' -> 3') before
reaching the first base of the next codon.
Pr: Exon probability: the probability of the prediction being a real
exon, taking into account the entire sequence. This is a good
indicator of the exon being correct.
PHAT COMMAND-LINE OPTIONS
[ -2 ] (Only works with halfphat) Predict on both strands. The default
is to predict in only the direct strand).
[ -acc accfile ] Use the acceptor splice site model in accfile, which should
be in the parameter directory being used. The default is "acceptor.dat"
Allows the use of alternative acceptor site models.
[ -debug ] Debug mode - prints *lots* of output to the screen. This option
was mainly used for debugging the code. You probably don't want to use it,
unless you're tinkering around with the code.
[ -don donfile ] Use the donor splice site model in donfile, which should
be in the parameter directory being used. The default is "donor.dat"
Allows the use of alternative donor site models.
[ -f ] Reverse complement ("flip") the sequence before predicting.
[ -geo ] Use geometric length distributions for exons, rather than the
general probability distributions specified in the parameter files
exon.len.*.dat located in the parameter directory. The mean of each
geometric distribution used is chosen to be the same as the mean of
the distribution supplied in the exon.len.*.dat file. This doesn't
make things run any faster or slower, its just available to make it
possible to compare the performance of HMM's with Generalized HMMs.
[ -seq seqname ] Use seqname (instead of the base name for output files)
as the name of the sequence in the .gff file. Can be useful (along
with the -o option) if you're processing a large sequence in separate
batches and then concatenating the resulting .gff files.
[ -noext ] Don't add the ".number" extension to the names given the objects
(proteins, mRNAs, etc) in the output files. The default is that
they will be given a numeric extension.
[ -norna ] Don't print the .rna file containing the predicted spliced
mRNA (without untranslated regions). The default is to print
it (in multiple fasta format).
[ -noextras ] Don't print "extra" output files, just write the .gff file
and print the predictions to the screen. If processing a large
sequence the output files can be very big, especially the exon
probability file.
[ -nopep ] Don't print the .pep file containing the predicted peptide
sequences. The default is to print it (in multiple fasta format).
[ -noprob ] Don't print the .prb file (in the case of fullphat) or the
.prf and .prb files (in the case of halfphat). These files contain
probabilities of each base in the sequence being some kind of
coding exon, and since there is one value per base they can be
quite large. The default is to print it or them.
[ -o filename ] Any output files will be written with the base name filename.
The default is to use the base name of the sequence file being
processed.
[ -offset offset ] Adds offset to all coordinates in the output written to
stdout and to the .gff file. Can be useful for changing coordinate
systems. Reported frame information will be relative to the new
coordinate system.
[ -okstops ] (fullphat only) If invoked, in-frame stop codons are given
very low probabilities, but they're not completely forbidden. You
might want to use this if the sequence being analyzed is prone to
errors, it could stop an erroneous stop codon breaking up an
otherwise likely exon. Without this option in-frame stops are
completely forbidden in fullphat. halfphat will be modified to
operate the same way in the future, currently it will allow
in-frame stops. For both phats, a stop codon spanning an intron
can't be ruled out. This is an unfortunate inherent feature of
all current HMM genefinders.
[ -p dir ] Specifies the directory where the parameter files are located.
dir can be an absolute or relative pathname and it doesn't matter
whether or not it is terminated with a "/". In addition to using
organism-specific parameters, you might want to use different
parameter sets for different regions.
[ -r ] (Only works with halfphat) Reverse complement the sequence before
predicting. All output coordinates will be in terms of the reverse
complemented sequence.
[ -v ] Run in verbose mode. A fair bit of output will be printed to the
screen. Most users can ignore this.
PART 2: running phat.pl
phat.pl is a perl script which uses bioperl modules to make it easier to
run fullphat/halfphat, and to extend their functionality. Advantages of
using the phat.pl script are:
1) Handles any input format that Bioperl can handle (Fasta,GenBank,GCG,etc)
2) Handles multiple sequence entries per input file.
3) Supports a simple type of parallel processing.
You will need to have recent versions of Perl and Bioperl installed to run
phat.pl. Tests to see the absolute minimum required versions of Perl and
Bioperl have not been conducted, but Perl 5.005_03 and Bioperl 0.6.1 or
later should do. They can be downloaded (for free):
Perl http://www.perl.com
Bioperl http://bio.perl.org
The phat.pl script will probably need to be edited a little to reflect where
you choose to store the fullphat/halfphat binaries on your system. It also
uses some rules determining which format is specified by which file
extension, these are easy to customise in the script. The default number
of CPUs is just 1, if you want to process sequences in parallel change
the value either in the script or when calling phat.pl with the "-cpu"
option.
Output naming conventions: Phat output files will be written in the current
directory. If the input file contains one sequence only the phat output
files will use the same base name. If the input file (called "inFile")
contains multiple sequences (named "seq1", "seq2" and so on) then phat
output files will have the base name inFile.seq1, inFile.seq2 and so on.
If there are any whitespace characters in the name, the name will be
shortened by dropping everything after and including the first space.
A word about parallel processing: it would be nicest to implement this
using Perl threads, but this is still an experimental feature in Perl
and is often not compiled in to Perl when it is installed. Instead we
choose to batch up to nCPUs jobs at a time and wait for them all to
finish. True parallel processing wouldn't wait around for the other
jobs to finish. Nevertheless, what is implemented here should speed
things up. In the future this will be improved.
EXAMPLES
phat.pl Sequence/AF030694.gbk -ph "-noprob -norna"
Runs fullphat on AF030694.gbk, suppresses generation
of .prb and .rna files.
phat.pl Sequence/AE00*.gbk -ph "-noextras" -cpu 4
Runs fullphat on AE00*.gbk sequences, suppresses generation of
everything other than the .gff files, processes up to 4 sequences
in parallel.
phat.pl Sequence/AF030694.gbk -ph "-noprob -nopep" -type half
Runs halfphat on AF030694.gbk, suppresses generation
of .prb and .pep files.
PHAT.PL COMMAND-LINE OPTIONS (call with no options to see help info)
[ -v ] Verbose mode
[ -h ] Help - prints out some usage info
[ -type half/full ] Tells phat.pl whether to use halfphat or fullphat
[ -ph "-opt1 -opt2 ... " ] Passes on the specified options list to the
call to halfphat or fullphat.
[ -cpu nCPUs ] The max number of sequence that will be analyzed in parallel.