COPYRIGHT NOTICE
Copyright 1988, 1991, 1992, 1994, 1995, 1996 by William R.
Pearson and the University of Virginia. All rights reserved. The
FASTA program and documentation may not be sold or incorporated
into a commercial product, in whole or in part, without written
consent of William R. Pearson and the University of Virginia.
For further information regarding permission for use or
reproduction, please contact: David Hudson, Assistant Provost for
Research, University of Virginia, P.O. Box 9025, Charlottesville,
VA 22906-9025, (804) 924-6853
The FASTA program package
Introduction
This documentation describes the version 2.0x of the FASTA
program package (see W. R. Pearson and D. J. Lipman (1988),
"Improved Tools for Biological Sequence Analysis", PNAS 85:2444-
2448, and W. R. Pearson (1990) "Rapid and Sensitive Sequence
Comparison with FASTP and FASTA" Methods in Enzymology 183:63-
98). Version 2.0 modifies version 1.8 to include explicit
statistical estimates for similarity scores based on the extreme
value distribution. In addition, FASTA protein alignments now
use the Smith-Waterman algorithm with no limitation on gap size.
FASTA and SSEARCH now use the BLOSUM50 matrix by default, with
options to change gap penalties on the command line. Version 1.7
replaces rdf2 and rss with prdf and prss, which use the extreme-
value distribution to calculate accurate probability estimates.
Although there are a large number of programs in this package,
they belong to four groups:
Library search programs: FASTA, FASTX, TFASTA, TFASTX, SSEARCH
Local homology programs: LFASTA, PLFASTA, LALIGN, PLALIGN, FLALIGN
Statistical significance: PRDF, RELATE, PRSS, RANDSEQ
Global alignment: ALIGN
In addition, I have included several programs for protein
sequence analysis, including a Kyte-Doolittle hydropathicity
plotting program (GREASE, TGREASE), and a secondary structure
prediction package (GARNIER).
The FASTA sequence comparison programs on this disk are
improved versions of the FASTP program, originally described in
Science (Lipman and Pearson, (1985) Science 227:1435-1441). We
have made several improvements. First, the library search
programs use a more sensitive method for the initial comparison
of two sequences which allows the scores of several similar
regions to be combined. As a result, the results of a library
search are now given with three scores, initn (the new initial
score which may include several similar regions), init1 (the old
fastp initial score from the best initial region), and opt (the
old fastp optimized score allowing gaps in a 32 residue wide
band).
These programs have also been modified to become "universal"
(hence FAST-A, for FASTA-All, as opposed to FAST-P (protein) or
FAST-N (nucleotides)); by changing the environment variable
SMATRIX, the programs can be used to search protein sequences,
DNA sequences, or whatever you like. By default, FASTA, LFASTA,
and the PRDF programs automatically recognize protein and DNA
sequences. Sequences are first read as amino acids, and then
converted to nucleotides if the sequence is greater than 85%
A,C,G,T (the '-n' option can be used to indicate DNA sequences).
TFASTA compares protein sequences to a translated DNA sequence.
Alternative scoring matrices can also be used. In addition to
the BLOSUM50 matrix for proteins, the PAM250 matrix or matrices
based on simple identities or the genetic code can also be used
for sequence comparisons or evaluation of significance. Several
different protein sequence matrices have been included;
instructions for constructing your own scoring matrix are
included in the file FORMAT.DOC.
The remainder of this document is divided into three sections:
(1) a brief history of the changes to the FASTA package; (2) A
guide to installing the programs and databases; (3) A guide to
using the FASTA programs. The programs are very easy to use, so
if you are using them on a machine that is administered by
someone else, you may want to skip to section (3) to learn how to
use the programs, and then read section (1) to look at some of
the more recent changes. If you are installing the programs on
your own machine, you will need to read section (2) carefully.
1. Revision History
1.1. Changes with version 2.0u
Version 2.0u provides several major improvements over
previous versions of FASTA (and SSEARCH). The most important is
the incorporation of explicit statistical estimates and
appropriate normalization of similarity scores. This improvement
is discussed in more detail below in the section entitled
Statistical Significance. In addition, all of the protein
comparison programs now use the BLOSUM50 matrix, with gap
penalties of -12, -2, by default. BLOSUM50 performs
significantly better than the older PAM250 matrix. PAM250 can
still be used with the command line option: -s 250. (DNA
sequence comparisons use a more stringent gap penalty of -16, -4,
which produces excellent statistical estimates when optimized
scores are used. TFASTA uses -16, -4 as well.)
The quality of the fit of the extreme value distribution to
the actual distribution of similarity scores is summarized with
the Kolmogorov-Smirnov statistic. The acceptance limits for this
statistic can be found in many statistics books. In general,
values <0.10 (N=30) indicate excellent agreement between the
actual and theoretical distributions. If this statistic is >
0.2, consider using a higher (more stringent) gap penalty, e.g.
-16, -4 rather than -12, -2. The default scoring matrix for DNA
has been changed to score +5 for an identity and -4 for a
mismatch. These are the same scores used by BLASTN.
With explicit expectation calculations, the program now
shows all scores and alignments with expectations less than 10.0
(with optimized scores, 2.0 without optimization) when the "-Q"
(quiet) mode is used. The expectation threshold can be changed
with the "-E" option.
Finally, the algorithm used to produce the final alignments
of protein sequences is now a full Smith-Waterman, with unlimited
gaps. (The older band-limited alignments are used for DNA
sequences and TFASTA by default, because Smith-Waterman
alignments are very slow for long sequences.) Both the optimized
and Smith-Waterman scores are reported; if the Smith-Waterman
score is higher, then additional gaps allowed a better alignment
and similarity score to be calculated.
FASTA searches now optimize similarity scores by default
(this slows searches about 2-fold (worst case) for ktup=2). Thus,
the meaning of the "-o" option has been reversed; "-o" now turns
off optimization and reports results sorted by "initn" scores.
Optimization significantly improves the sensitivity of FASTA, so
that it almost matches Smith-Waterman. With version 2.0, the
default band width used for optimized calculations can be varied
with the "-y" option. For proteins with ktup=2, a width of 16
(-y 16) is used; 16 is also used for DNA sequences. For proteins
and ktup=1, a width of 32 is used. Searches that disable
optimization with the "-o" option will work fine for sequences
that share 25% or more identity in general, but to detect
evolutionary relationships with 20% - 25% identity, the more
sensitive default optimization is often required. Optimization
is required for accurate statistical estimates with either
protein or DNA sequences.
The FASTA package now includes FASTX, a program that
compares a DNA sequence to a protein sequence database by
translating the DNA sequence in three frames (the reverse frames
are selected with the -i option) and aligning the three-frame
translation with the sequences in the protein database.
Alignment scores allow frameshifts so that a cDNA or EST sequence
with insertion/deletion errors can be aligned with its homologues
from beginning to end.
With release 20u6, there is also a TFASTX program, which is
a replacement for TFASTA. TFASTA treats each of the six reading
frames of a DNA library sequence as a different sequence; TFASTX
compares a protein sequence against only two sequences from each
DNA sequence - the forward and reverse orientation. For a given
orientation, TFASTX calculates a similarity score for alignments
that allow frameshifts, thus considering all possible reading
frames.
Another new program is included - randseq - which will
produce a randomly shuffled (uniform or local shuffle) from an
input sequence. This randomly shuffled sequence can be used to
evaluate the statistical estimates produced by FASTA, SSEARCH, or
BLAST.
1.2. Changes with version 1.7
Version 1.7 has been released to provide the PRDF and PRSS
programs for shuffling sequences and estimating accurately the
probabilities of the unshuffled-sequence scores.
PRDF a version of RDF2 that uses calculates the probability
of a similarity score more accurately by using a fit to
an extreme value distribution. Code to fit the extreme
value distribution parameters and the impetus to update
RDF2 was provided by Phil Green, U. of Washington.
PRSS a version of PRDF that uses a rigorous Smith-Waterman
calculation to score similarities
1.3. Changes with version 1.6
FASTA version 1.6 uses a new method for calculating optimal
scores in a band (the optimization or last step in the FASTA
algorithm). In addition, it uses a linear-space method for
calculating the actual alignments. FASTA v1.6 package includes
several new programs:
SSEARCH a program to search a sequence database using the
rigorous Smith-Waterman algorithm (this program is
about 100-fold slower than FASTA with ktup=2 (for
proteins).
LALIGN A rigorous local sequence alignment program that will
display the N-best local alignments (N=10 by default).
PLALIGN a version of lalign that plots the local alignments to
a tektronix display.
FLALIGN a version of lalign that plots the local alignments to
a GCG Figure file.
The LALIGN/PLALIGN/FLALIGN programs incorporate the "sim"
algorithm described by Huang and Miller (1991) Adv. Appl. Math.
12:337-357. The SSEARCH and PRSS programs incorporate algorithms
described by Huang, Hardison, and Miller (1990) CABIOS 6:373-381.
LFASTA and PLFASTA now calculate a different number of local
similarities; they now behave more like LALIGN/PLALIGN. Since
local alignments of identical sequences produce "mirror-image"
alignments, lalign and lfasta consider only one-half of the
potential alignments between sequences from identical file names.
Thus
lfasta mchu.aa mchu.aa
Displays only two alignments, with earlier versions of the
program, it would have displayed five, including the identity
alignment. PLFASTA does display five alignments; when two
identical filenames are given, it draws the identity alignment,
calculates the two unique local alignments, draws them, and draws
their mirror images. LFASTA/PLFASTA and LALIGN/PLALIGN use the
filenames, rather than the actual sequences, to determine whether
sequences are identical; you can "trick" the programs into
behaving the old way by putting the same sequence in two
different files.
1.4. Changes with version 1.5
FASTA version 1.5 includes a number of substantial revisions
to improve the performance and sensitivity of the program. It is
now possible to tell the program to optimize all of the initn
scores greater than a threshold. The threshold is set at the
same value as the old FASTA cutoff score. Alternatively, you can
tell FASTA to sort the results by the init1, rather than the
initn, score by using the -1 option. FASTA -1 ... will report
the results the way the older FASTP program did.
A new method has been provided for selecting libraries. In
the past, one could enter the name of a sequence file to be
searched or a single letter that would specify a library from the
list included in the $FASTLIBS file. Now, you can specify a set
of library files with a string of letters preceded by a '%'.
Thus, if the FASTLIBS file has the lines:
Genbank 70 primates$1P/seqlib/gbpri.seq 1
Genbank 70 rodents$1R/seqlib/gbrod.seq 1
Genbank 70 other mammals$1M/seqlib/gbmam.seq 1
Genbank 70 vertebrates $1B/seqlib/gbvrt.seq 1
Then the string: "%PRMB" would tell FASTA to search the four
libraries listed above. The %PRMB string can be entered either
on the command line or when the program asks for a filename or
library letter.
FASTA1.5 also provides additional flexibility for specifying
the number of results and alignments to be displayed with the -Q
(quiet) option. The -b number option allows you to specify the
number of sequence scores to show when the search is finished.
Thus
FASTA -b 100 ...
tells the program to display the top 100 sequence scores. In the
past, if you displayed 100 scores (in -Q mode), you would also
have store 100 alignments. The -d option allows you to limit the
number of alignments shown. FASTA -b 100 -d 20 would show 100
scores and 20 alignments.
Finally, FASTA can provide a complete list of all of the
sequences and scores calculated to a file with the -r (results)
option. FASTA -r results.out ... creates a file with a list of
scores for every sequence in the library. The list is not
sorted, and only includes those scores calculated during the
initial scan of the library.
2. Installing the FASTA package
2.1. Installing the programs
2.1.1. Unix version
The FASTA distribution comes with several makefile's that
can be used to compile the FASTA programs. Over the years, as
ATT Unix System 5 and BSD unix have converged, these files have
become very similar. To begin with, I recommend using the
standard Makefile. There are two values in the makefile that
should be checked against the values used on your system: the HZ
value, which is the frequency in ticks per second used by the
times() system call, this value can usually be found by running:
grep HZ /usr/include/sys/*
and the functions available to return random numbers. If you
have a rand48() function that returns a 32-bit random number, use
it and use the lines:
NRAND=nrand48
RANFLG= -DRAND32
If not, you will need to use the rand() function call and
determine whether it returns a 16-bit or a 32-bit value. These
functions are used by PRDF and PRSS. If you have problems
compiling the programs, you may want to examine the makefile.unx
and makefile.sun files, to look for differences. I have tried to
use very standard unix functions in these programs, and they have
been successfully compiled, with very small changes to the
Makefile, on Sun's (Sun OS 4.1), IBM RS/6000's (AIX), and MIPS
machines (under the BSD environment).
2.1.2. IBM-PC/DOS version
For the IBM-PC/DOS version, the FASTA source code disk
contains the complete source code to all of the programs on the
other disks. The programs were compiled with Borland's Turbo
'C++', using Borland's MAKE utility. The graphics programs
(PLFASTA, TGREASE) use the graphics device drivers supplied with
the Turbo 'C' V2.0 package. Also included are the documentation
files PROGRAMS.DOC and FORMAT.DOC. You do not need any of the
files the source code disk to run the programs. The files on
this disk are identical to the UNIX and VMS versions that run on
larger machines. Also included is the code to compile
ALIGN0.EXE. ALIGN0 is the same as ALIGN, but does not penalize
for end-gaps.
If you have the DOS or Macintosh version of the FASTA
package, to install the programs you should:
(1) Make a new directory (folder) for the FASTA programs.
This need not be the same as the directory for your
sequence databases.
(2) Copy the files from the FASTA source disk to the new
directory.
(3) (DOS only) Edit your AUTOEXEC.BAT file to (a) modify your
PATH command to include the FASTA directory and (b) add
the line:
set FASTLIBS=c:\yourfastadirectory\fastgbs
On the Macintosh, you may need to edit the "environment"
file and change the line that reads:
FASTLIBS=fastgbs
to indicate the full directory path for the fastgbs file,
for example:
FASTLIBS=Q105:FASTA:fastgbs
(4) Finally, you will need to edit the fastgbs file. This is
usually the most confusing part of the installation. An
example of this file is shown below; to customize this
file for your machine, you will need to change the file
names from those provided in the fastgbs file to ones that
reflect the directory names and file names you use on your
machine. This is explained in more detail below. In
addition, some entries in the fastgbs file refer to other
files of file names. These files of file names (as
opposed to actual database files) may also need to be
edited.
2.2. Installing the libraries
2.2.1. The NBRF protein sequence library
The FASTA program package does not include any protein or
DNA sequence libraries. You can obtain the PIR protein sequence
database from:
National Biomedical Research Foundation
Georgetown University Medical Center
3900 Reservoir Rd, N.W.
Washington, D.C. 20007
In addition, this database is available via anonymous ftp from
the host "ftp.bchs.uh.edu". It is available in two formats, VMS
and CODATA format. The "VMS" format (library type 5 below) can
be searched much faster, can be easily reformatted for use by the
"BLAST" rapid searching program, and is compatible with the
Genetics Computer Group package of programs. The CODATA format
is used by the EUGENE/MBIR computing package from Baylor (library
type 2).
2.2.2. The GENBANK DNA sequence library
FASTA, and TFASTA search sequences from the GENBANK
"flatfile" (not ASN.1) DNA sequence library in the flat-file
format distributed by the National Center for Biotechnology
Information and the PIR format used by EBI/EMBL. CD-ROMs can be
obtained from:
Genbank
National Center for Biotechnology Information
National Library of Medicine
National Institutes of Health
8600 Rockville Pike
Bethesda, MD 20894
The GenBank DNA sequence library is also available via
anonymous FTP from ncbi.nlm.nih.gov.
2.2.3. The EBI/EMBL CD-ROM libraries
The European Bioinformatics Institute (EBI) is now
distributing the EMBL CD-ROM that contains both the complete EMBL
DNA sequence database (which should be essentially identical to
the GenBank DNA sequence database) and the SWISS-PROT protein
sequence database. SWISS-PROT is derived from the NBRF Protein
sequence database with additions from the EBI/EMBL DNA sequence
database. This CD-ROM is a "best-buy," since it provides both
DNA and protein sequence libraries. It is available from:
European Bioinformatics Institute
Hinxton Genome Campus, Hinxton Hall
Hinxton, Cambridge CB10 1RQ,
United Kingdom
Tel: +44 1223 4944
Fax: +44 1223 494468
Email: DATALIB@ebi.ac.uk
In addition, the SWISS-PROT protein sequence database is
available via anonymous FTP from ncbi.nlm.nih.gov.
2.3. Finding the libraries: FASTLIBS
FASTA and TFASTA use the environment variable FASTLIBS to
find the protein and DNA sequence libraries. The FASTLIBS
variable contains the name of a file that has the actual
filenames of the libraries. The FASTGBS file on is an example of
a file that can be referred to by FASTLIBS. To use the FASTGBS
file, type:
setenv FASTLIBS /usr/lib/fasta/fastgbs (BSD UNIX/csh)
or
export FASTLIBS=/usr/lib/fasta/fastgbs (SysV UNIX/ksh)
Then edit the FASTGBS file to indicate where the protein and DNA
sequence libraries can be found. If you have a hard disk and
your protein sequence library is kept in the file
/usr/lib/aabank.lib and your Genbank DNA sequence library is kept
in the directory: /usr/lib/genbank, then fastgbs might contain:
NBRF Protein$0P/usr/lib/seq/aabank.lib 0
SWISS PROT 10$0S/usr/lib/vmspir/swiss.seq 5
GB Primate$1P@/usr/lib/genbank/gpri.nam
GB Rodent$1R@/usr/lib/genbank/grod.nam
GB Mammal$1M@/usr/lib/genbank/gmammal.nam
^ 1 ^^^^ 4 ^ ^
23 (5)
The first line of this file says that there is a copy of the NBRF
protein sequence database (which is a protein database) that can
be selected by typing "P" on the command line or when the
database menu is presented in the file /usr/lib/seq/aabank.lib.
Note that there are 4 or 5 fields in the lines in fastgbs.
The first field is the description of the library which will be
displayed by FASTA; it ends with a '$'. The second field (1
character), is a 0 if the library is a protein library and 1 if
it is a DNA library. The third field (1 character) is the
character to be typed to select the library.
The fourth field is the name of the library file. In the
example above, the /usr/lib/seq/aabank.lib file contains the
entire protein sequence library. However the DNA library file
names are preceded by a '@', because these files (gpri.nam,
grod.nam, gmammal.nam) do not contain the sequences; instead they
contain the names of the files which contain the sequences. This
is done because the GENBANK DNA database is broken down in to a
large number of smaller files. In order to search the entire
primate database, you must search more than a dozen files.
In addition, an optional fifth field can be used to specify
the format of the library file. Alternatively, you can specify
the library format in a file of file names (a file preceded by an
'@'). This field must be separated from the file name by a space
character (' ') from the filename. In the example above, the
aabank.lib file is in Pearson/FASTA format, while the swiss.seq
file is in PIR/VMS format (from the EMBL CD-ROM). Currently,
FASTA can read the following formats:
0 Pearson/FASTA (>SEQID - comment/sequence)
1 Uncompressed Genbank (LOCUS/DEFINITION/ORIGIN)
2 NBRF CODATA (ENTRY/SEQUENCE)
3 EMBL/SWISS-PROT (ID/DE/SQ)
4 Intelligenetics (;comment/SEQID/sequence)
5 NBRF/PIR VMS (>P1;SEQID/comment/sequence)
6 GCG (version 8.0) Unix Protein and DNA (compressed)
11 NCBI Blast1.3.2 format (unix only)
In particular, this version will work with the EMBL and PIR VMS
formats that are distributed on the EMBL CD-ROM. The latter
format (PIR VMS) is much faster to search than EMBL format. This
release also works with the protein and DNA database formats
created for the BLASTP and BLASTN programs by SETDB and PRESSDB
and with the new NCBI search format. If a library format is not
specified, for example, because you are just comparing two
sequences, Pearson/FASTA (format 0) is used by default. To
change this default, you may set the LIBTYPE environment variable
to a number. For example,
setenv LIBTYPE 1
would cause the program to use the GenBank LOCUS format by
default for libraries (or the second sequence file), but the
Pearson/FASTA format would still be used for the query sequence.
You can specify a group of library files by putting a '@'
symbol before a file that contains a list of file names to be
searched. For example, if @gmam.nam is in the fastgbs file, the
file "gmam.nam" might contain the lines:
' in the first
column. (3) distributed sequence libraries (this is a broad class
that includes the NBRF/PIR VMS and blocked ascii formats, Genbank
flat-file format, EMBL flat-file format, and Intelligenetics
format. All of the files that you create should be of type (1)
or (2). Type (2) files (ones with a be used as query or library
sequence files by all of the programs.
I have included several sample test files, *.AA. The first
line may begin with a '>' or ';' followed by a comment. The
text after ';' in other lines will be ignored. Spaces and
tabs (and anything else that is not an amino-acid code) are
ignored.
Library files should have the form:
>Sequence name and identifier
A F A S Y T .... actual sequence.
F S S .... second line of sequence.
>Next sequence name and identifier
This is often referred to as "FASTA" or "Pearson" format. You
can build your own library by concatenating several sequence
files. Just be sure that each sequence is preceded by a line
beginning with a '>' with a sequence name.
The test file should not have lines longer than 120
characters, and sequences entered with word processors should use
a document mode, with normal carriage returns at the end of
lines.
Program Summary
3.3. Sequence search programs
FASTA universal sequence comparison. Defaults to comparing
protein sequences; if the sequences are > 85% A+C+G+T
or the -n option is used, a DNA sequence is assumed.
FASTX Search a protein sequence library using amino acid
sequence comparison to the forward three frames of a
translated DNA query sequence. (The reverse frames are
specified with the -i option.) Alignment scores allow
frameshifts; the final alignment uses a Smith-Waterman
type alignment routine (no limit on gaps) that allows
frameshifts.
TFASTA Search DNA library for a protein sequence by
translating the DNA sequence to protein in all six
frames (three forward frames with the -3 command line
option). TFASTA with ktup=2 is about as fast as a DNA
FASTA with ktup=4, and is substantially more sensitive.
(also reads the GENBANK library)
TFASTX Search DNA library for a protein sequence by
translating the DNA sequence to protein in all six
frames (three forward frames with the -3 command line
option) calculating similarity scores that allow
frameshifts. TFASTX produces an optimal Smith-Waterman
alignment of the query and translated-library sequence.
SSEARCH Universal sequence comparison using the Smith-Waterman
algorithm ( T. F. Smith and M. S. Waterman (1981) J.
Mol. Biol. 147:195-197). This program uses code
developed by Huang and Miller (X. Huang, R. C.
Hardison, W. Miller (1990) CABIOS 6:373-381) for
calculating the local similarity score and code from
the ALIGN program (see below) for calculating the local
alignment. SSEARCH is about 50-times slower than FASTA
with ktup=2 (for proteins).
ALIGN optimal global alignment of two sequences with no
short-cuts. This program is a slightly modified
version of one taken from E. Myers and W. Miller. The
algorithm is described in E. Myers and W. Miller,
"Optimal Alignments in Linear Space" (CABIOS (1988)
4:11-17).
3.4. Local similarity programs
LFASTA local similarity searches showing local alignments.
The algorithm used to calculate the local alignment in
a band has been improved (Chao, Pearson, and Miller,
submitted).
PLFASTA local similarity searches with plot output (on the IBM,
this program requires that the environment variable
BGIDIR be set).
PCLFASTA (unix only) local similarity searches with plot output
using pic commands.
LALIGN Calculates the N-best local alignments using a rigorous
algorithm. (N=10 by default.) The algorithm was
developed by Huang and Miller (X. Huang and W. Miller
(1991) Adv. Appl. Math. 12:337-357), which is a
linear-space version of an algorithm described by M. S.
Waterman and M. Eggert (J. Mol. Biol. 197:723-728).
Like SSEARCH, LALIGN is rigorous, but also very slow.
PLALIGN A version of LALIGN that plots its output to a screen
or to a Tektronix terminal emulator.
3.5. Statistical Significance
With version 2.0 of the FASTA program distribution, FASTA,
TFASTA, and SSEARCH now provide estimates of statistical
significance for library searches. Work by Altschul, Arratia,
Karlin, Mott, Waterman, and others (see Altschul et al. (1994)
Nature Genetics 6:119 for an excellent review) suggests that
local sequence similarity scores follow the extreme value
distribution, so that P(s > x) = 1 - exp(-exp(-lambda(x-u)) where
u = ln(Kmn)/lambda and m,m are the lengths of the query and
library sequence. This formula can be rewritten as: 1 - exp(-Kmn
exp(-lambda x), which shows that the average score for an
unrelated library sequence increases with the logarithm of the
length of the library sequence. FASTA and SSEARCH use simple
linear regression against the the log of the library sequence
length to calculate a normalized "z-score" with mean 50,
regardless of library sequence length, and variance 10. These
z-scores can then be used with the extreme value distribution and
the poisson distribution (to account for the fact that each
library sequence comparison is an independent test) to calculate
the number of library sequences to obtain a score greater than or
equal to the score obtained in the search. The original idea and
routines to do the linear regression on library sequence length
were provided Phil Green, U. Washington. This version of FASTA
and SSEARCH uses a slightly different strategy for fitting the
data than those originally provided by Dr. Green.
The expected number of sequences is plotted in the histogram
using an "*". Since the parameters for the extreme value
distribution are not calculated directly from the distribution of
similarity scores, the pattern of "*'s" in the histogram gives a
qualitative view of how well the statistical theory fits the
similarity scores calculated by FASTA and SSEARCH. For FASTA, if
optimized scores are calculated for each sequence in the database
(the default), the agreement between the actual distribution of
"z-scores" and the expected distribution based on the length
dependence of the score and the extreme value distribution is
usually very good. Likewise, the distribution of SSEARCH Smith-
Waterman scores typically agrees closely with the actual
distribution of "z-scores." The agreement with unoptimized
scores, ktup=2, is often not very good, with too many high
scoring sequences and too few low scoring sequences compared with
the predicted relationship between sequence length and similarity
score. In those cases, the expectation values may be
overestimates.
The statistical routines assume that the library contains a
large sample of unrelated sequences. If this is not the case,
then the expectation values are meaningless. Likewise, if there
are fewer than 20 sequences in the library, the statistical
calculations are not done.
For protein searches, library sequences with E() values <
0.01 for searches of a 10,000 entry protein database are almost
always homologous. Frequently sequences with E()-values from 1 -
10 are related as well. Remember, however, that these E() values
also reflect differences between the amino acid composition of
the query sequence and that of the "average" library sequence.
Thus, when searches are done with query sequences with "biased"
amino-acid composition, unrelated sequences may have
"significant" scores because of sequence bias. The programs
below, PRDF and PRSS, can address this problem by calculating
similarity scores for random sequences with the same length and
amino acid composition.
If optimization is not used ("-o"), E-values for DNA
sequences overestimate the significance of the scores that are
obtained and unrelated sequences frequently have E()-values <
0.0005. With optimization, the agreement between E()-value
compares favorably with protein sequence comparison. This is in
part due to the use of more stringent gap penalties for DNA
sequence comparison, -16, -4 rather than -12, -2. With the
latter penalties, many unrelated sequences appear to have
significant similarity. Nevertheless, since protein sequence
comparison is much more sensitive, DNA sequence comparison should
not be used to identify sequences that encode protein. Even with
ktup=6, optimization rarely increases run-times more than 50%
with mRNA-size query sequences. Optimization should be used
whenever possible.
Similar comments apply to TFASTA, where higher gap
penalties (-16,-4) are required for accurate statistical
estimates. Because TFASTA produces so many artificial "coding"
sequences with atypical amino acid compositions, the statistical
estimates with TFASTA are often over estimates. With optimized
scores, ktup=1, and gap penalties of -16, -4, unrelated sequences
will sometimes have E() values of 0.1. If initn scores are used,
unrelated sequences may have have E() values < 0.01.
PRDF improved version of RDF program that includes accurate
probability estimates for all three scoring methods
(includes local or window shuffle routine)
PRSS A version of PRDF that uses the rigorous Smith-Waterman
calculation used by SSEARCH.
RANDSEQ produces a randomly shuffled sequence from a query
sequence.
RELATE significance program described by Dayhoff (Atlas of
Protein Sequence and Structure, Vol. 5, Supplement 3).
Each chunk of 25 residues in one sequence is compared
to every 25 residue fragment of the second sequence.
Sequences which are genuinely related will have a large
number of scores greater than 3 standard deviations
above the mean score of all of the comparisons.
3.6. Other analysis programs
AACOMP calculate the amino acid composition and molecular
weight of a sequence.
BESTSCOR calculate the best self-comparison score.
GREASE Kyte-Doolittle hydropathicity profile
TGREASE graphic plot of Kyte-Doolittle profile
FROMGB convert from GenBank LOCUS format (also used by the
IBI-Pustell programs) to Pearson/FASTA format.
GARNIER A secondary structure prediction program using the
method of Garnier, Osgusthorpe, and Robson, J. Mol.
Biol., (1978) 120:97-120.
3.7. Options
These programs have a number of output options, which are
invoked by the environment variables LINLEN, SHOWALL, and MARKX.
Alternatively, these values can be controlled by command line
options. The number of sequence residues per output line is now
adjustable by setting the environment variable LINLEN, or the
command line option -w. LINLEN is normally 60, to change it set
LINLEN=80 before running the program or add -w 80 to the command
line. LINLEN can be set up to 200. SHOWALL (-a) determines
whether all, or just a portion, of the aligned sequences are
displayed. Previously, FASTP would show the entire length of
both sequences in an alignment while FASTN would only show the
portions of the two sequences that overlapped. Now the default is
to show only the overlap between the two sequences, to show
complete sequences, set SHOWALL=1, or use the -a option on the
command line.
The differences between the two aligned sequences can be
highlighted in three different ways by changing the environment
variable MARKX or the -m option. Normally (MARKX=0) the program
uses ':' do denote identities and '.' to denote conservative
replacements. If MARKX=1, the program will not mark identities;
instead conservative replacements are denoted by a 'x' and non-
conservative substitutions by a 'X'. If MARKX=2, the residues in
the second sequence are only shown if they are different from the
first. MARKX=3 displays the aligned library sequences without the
query sequence; these can be used to build a primitive multiple
alignment. MARKX=4 provides a graphical display of the
boundaries of the alignments. Thus the five options are:
MARKX=0 MARKX=1 MARKX=2 MARKX=3 MARKX=4
MWRTCGPPYT MWRTCGPPYT MWRTCGPPYT MWRTCGPPYT
::..:: ::: xx X ..KS..Y... MWKSCGYPYT ----------
MWKSCGYPYT MWKSCGYPYT
(fasta20u4, Feb. 1996) In addition MARKX=10 is a new, parseable
format for use with other programs. See the file"readme.v20u4"
for a more complete description.
3.8. Command line options
It is now possible to specify several options on the
command line, instead of using environment variables. The
command line options are preceded by a dash; the following
options are available:
-a same as showall=1
-A force Smith-Waterman alignments for DNA sequences and
TFASA. By default, only FASTA protein sequence
comparisons use Smith-Waterman alignments.
-b # Number of sequence scores to be shown on output. In
the absence of this option, fasta (and tfasta and
ssearch) display all library sequences obtaining
similarity scores with expectations less than 10.0 if
optimized score are used, or 2.0 if they are not. The
-b option can limit the display further, but it will
not cause additional sequences to be displayed.
-c # Threshold score for optimization (OPTCUT). Set "-c 1"
to optimize every sequence in a database. (This slows
the program down about 5-fold).
-E # Limit the number of scores and alignments shown based
on the expected number of scores. Used to override the
expectation value of 10.0 used by default. When used
with -Q, -E 2.0 will show all library sequences with
scores with an expectation value <= 2.0.
-d # Number of alignments to be reported by default. (Used
in conjunction with -Q). No longer necessary, see "-b"
above.
-f Penalty for the first residue in a gap (-12 by default
for proteins, -16 for DNA or for TFASTA).
-g Penalty for additional residues in a gap (-2 by default
for proteins, -4 for DNA and TFASTA ).
-h Penalty for frameshift (FASTX, TFASTX only).
-H Omit histogram.
-i Invert (reverse complement) the query sequence if it is
DNA. For TFASTX, search the reverse complement of the
library sequence only.
-k # Threshold for joining init1 segments to build an initn
score (GAPCUT).
-l file Location of library menu file (FASTLIBS).
-L Display more information about the library sequence in
the alignment.
-m # MARKX = # (0, 1, 2, 3, 4, 10)
-n Force the query sequence to be treated as a DNA
sequence. This is particularly useful for query
sequences that contain a large number of ambiguous
residues, e.g. transcription factor binding sites.
-O Send copy of results to "filename." Helpful for
environments without STDOUT.
-o Turn off default optimization of all scores greater
than OPTCUT. Sort results by "initn" scores.
-Q,-q Quiet - does not prompt for any input. Writes scores
and alignments to the terminal or standard output file.
-r file Save a results summary line for every sequence in the
sequence library. The summary line includes the
sequence identifier, superfamily number (if available)
position in the library, and the similarity scores
calculated. This option can be used to evaluate the
sensitivity and selectivity of different search
strategies (see W. R. Pearson (1991) Genomics 11:635-
650.)
-s file SMATRIX is read from file. Several SMATRIX files are
provided with the standard distribution. For protein
sequences: codaa.mat - based on minimum mutation
matrix; idnaa.mat - identity matrix; pam250.mat - the
PAM250 matrix developed by Dayhoff et al (Atlas of
Protein Sequence and Structure, vol. 5, suppl. 3,
1978); pam120.mat - a PAM120 matrix. The default
scoring matrix is BLOSUM50, PAM250 is available with
"-s 250", BLOSUM62 ("-s BL62") is also available.
-v (LINEVAL) values used for line styles in plfasta
-w # Line length (width) = number (<200)
-x Specifies offsets for the beginning of the query and
library sequence. For example, if you are comparing
upstream regions for two genes, and the first sequence
contains 500 nt of upstream sequence while the second
contains 300 nt of upstream sequence, you might try:
fasta -x "-500 -300" seq1.nt seq2.nt
If the -x option is not used, FASTA assumes numbering
starts with 1. This option will not work properly with
the translated library sequence with tfasta. (You
should double check to be certain the negative
numbering works properly.)
-y Set the width of the band used for calculating
"optimized" scores. For proteins and ktup=2, the width
is 16. For proteins with ktup=1, the width is 32 by
default. For DNA the width is 16.
-z Turn off statistical calculations.
-1 sort output by init1 score (as FASTP used to do).
-3 (TFASTA, TFASTX only) translate only three forward
frames
For example:
fasta -w 80 -a seq1.aa seq.aa
would compare the sequence in seq1.aa to that in seq2.aa and
display the results with 80 residues on an output line, showing
all of the residues in both sequences. Be sure to enter the
options before entering the file names, or just enter the options
on the command line, and the program will prompt for the file
names.
Not all of these options are appropriate for all of the
programs. The options above are used by FASTA and TFASTA. RELATE
uses the -s option, ALIGN uses the -w, -m, and -s options, and
the PRDF program uses -c, -f, -k, and -s.
4. Environment variable summary
Environment variables allow you to set search parameters
that will be used frequently when you run a program; for example,
if you prefer to use the PAM250 scoring matrix, you might "set
SMATRIX=250." Command line parameters, if used, always override
environment variable settings. The following environment
variables are used by this program:
AABANK the file name of the default sequence library.
FASTLIBS the location of the file which contains the list of
library files to be searched.
GAPCUT threshold used for joining init1 regions in the second
step of FASTA. Normally set based on sequence length
and ktup.
LIBTYPE used to specify the format of the library sequence for
FASTA and TFASTA.
LINLEN output line length - can go up to 200
LINEVAL used by plfasta to determine the relationship between
line style and similarity score (-v). This should be a
string of three numbers, e.g. "200 100 50"
MARKX symbol for denoting matches, mismatches. Note that this
symbol is only used across the optimized local region;
sequences that are outside this region are not marked.
OPTCUT Set the threshold to be used for optimization in a band
around the best initial region. Normally the OPTCUT
value is calculated from the length of the sequence and
the ktup value (for a 200 residue sequence, it is about
28). If OPTCUT=1, every sequence in the database will
be optimized. This is the most sensitive option.
PAMFACT This version of fasta uses a more sensitive method for
identifying initial regions. Instead of using a
constant factor (fact) for each match in a ktup, it
uses the scoring matrix (PAM) scores. While this works
well for protein sequences, it has not been as
carefully tested for DNA sequences, so by default, this
modification is used for proteins but not for DNA.
Setting the PAMFACT environment variable to 1 forces
the option on; PAMFACT=0 turns it off.
SHOWALL on output, show the complete sequence instead of just
the overlap of the two aligned sequences.
SMATRIX alternative scoring matrix file.
TEKPLOT (IBM-PC only, Unix and VMS versions generate Tektronix
graphics by default) Generate Tektronix output.
Normally, PLFASTA and TGREASE plot graphs using the
Turbo C graphics library. Unfortunately, often these
plots cannot be printed out without special programs.
However, if you set TEKPLOT=1, tektronix graphics
commands will be used. Tektronix commands can be used
together with the PLOTDEV program, available from
Microplot Systems. They no lonter sell this program,
but it can be downloaded from
http://iquest.com/~microplt/index1.html. PLOTDEV also
allows you to print out graphics on the screen.
As always, please inform me of bugs as soon as possible.
William R. Pearson
Department of Biochemistry
Box 440, Jordan Hall
U. of Virginia
Charlottesville, VA
wrp@virginia.EDU