[ Program Manual | User's Guide | Data Files | Databases ]
TFastA does a Pearson and Lipman search for similarity between a query peptide sequence and any group of nucleotide sequences. TFastA translates the nucleotide sequences in all six reading frames before performing the comparison. It is designed to answer the question, "What implied peptide sequences in a nucleotide sequence database are similar to my peptide sequence?"
TFastA uses the method of Pearson and Lipman (Proc. Natl. Acad. Sci. USA 85; 2444-2448 (1988)) to search for similarities between a query protein sequence and any group of nucleotide sequences. TFastA translates the nucleotide sequences in all six frames before performing the comparison. Each translated reading frame is treated as a separate sequence to be searched. In the first step of this search, the comparison can be viewed as a set of dot plots, with the query as the vertical sequence and the group of sequences to which the query is being compared as the different horizontal sequences. This first step finds the registers of comparison (diagonals) having the largest number of short perfect matches (words) for each comparison. In the second step, these "best" regions are rescored using a scoring matrix that allows conservative replacements, ambiguity symbols, and runs of identities shorter than the size of a word. In the third step, the program checks to see if some of these initial highest-scoring diagonals can be joined together. Finally, the search set sequences with the highest scores are aligned to the query sequence for display.
A word is any short sequence (n-mer or k-tuple) where you have set n to some small integer less than or equal to six. The word GGATGG is one of the 4,096 possible words of length six that can be created from an alphabet consisting of the four letters G, A, T, and C. The word QL is one of the 400 possible words of length two that you can make with the 20 letters of the amino acid alphabet.
Here is a session using TFastA to identify sequences in the GenEMBL nucleotide sequence database that contain translated regions similar to a human globin protein:
% tfasta
TFASTA with what query sequence ? ggamma.pep
Begin (* 1 *) ?
End (* 148 *) ?
Search for query in what sequence(s) (* GenEMBL:* *) ?
What word size (* 2 *) ?
Don't show scores whose E() value exceeds: (* 10.0 *):
What should I call the output file (* ggamma.tfasta *) ?
Don't show scores whose E() value exceeds: (* 10.0 *):
1 Sequences 1,349 aa searched GB_BA:A16STM214
101 Sequences 153,771 aa searched GB_BA:ABCPRECA
/////////////////////////////////////////////////////////
CPU time used:
Database scan: 0:12:04.0
Post-scan processing: 0:00: 2.2
Total CPU time: 0:12:06.8
Output file: ggamma.tfasta
%
The output from FastA is a list file, and is suitable for input to any GCG program that allows indirect file specifications. (For information about indirect file specification, see Chapter 2, Using Sequence Files and Databases of the User's Guide.)
Here is some of the output file:
!!SEQUENCE_LIST 1.0
(Peptide) TFASTA of: ggamma.pep from: 1 to: 148 September 18, 1996 10:47
TRANSLATE of: gamma.seq check: 6474 from: 2179 to: 2270
and of: gamma.seq check: 6474 from: 2393 to: 2615
and of: gamma.seq check: 6474 from: 3502 to: 3630
generated symbols 1 to: 148.
Human fetal beta globins G and A gamma
from Shen, Slightom and Smithies, Cell 26; 191-203. . . .
TO: GenEMBL:* Sequences: 260,418 Symbols: 350,478,837 Word Size: 2
Sequences too short to analyze: 20 (84 symbols)
Databases searched:
GenBank, Release 95.0, Released on 19Jun96, Formatted on 19Jul1996
GenBank, Release 95.0, Released on 21Jun96, Formatted on 21Jul1996
EMBL, Release 47.0, Released on 22Jun96, Formatted on 22Jul1996
EMBL, Release 47.0, Released on 23Jun96, Formatted on 23Jul1996
Searching all six frames.
Scoring matrix: GenRunData:blosum50.cmp
Variable pamfactor used
Gap creation penalty: 16 Gap extension penalty: 4
Histogram Key:
Each histogram symbol represents 2705 search set sequences
Each inset symbol represents 25 search set sequences
z-scores computed from opt scores
z-score obs exp
(=) (*)
< 20 69060 0 :*=========================
22 11 0 :*
24 31 1 :*
26 55 31 :*
28 150 334 :*
30 542 2032 :*
32 1989 7857 := *
34 6700 21307 :=== *
36 20529 43759 :======== *
38 48452 72318 :================== *
40 84574 100877 :================================ *
42 120981 123309 :=============================================*
44 146732 136022 :==================================================*====
46 162247 138542 :===================================================*========
48 161214 132638 :=================================================*==========
50 149022 121032 :============================================*===========
52 125208 106408 :=======================================*=======
54 100266 90891 :=================================*====
56 78114 75922 :============================*
58 58357 62330 :====================== *
60 47900 50491 :==================*
62 34216 40479 :============= *
64 26386 32193 :========== *
66 20852 25444 :======== *
68 16592 20014 :=======*
70 12461 15684 :=====*
72 9490 12256 :====*
74 7429 9555 :===*
76 6279 7437 :==*
78 5161 5781 :==*
80 4310 4489 :=*
82 3703 3434 :=*
84 2883 2720 :=*
86 2460 2105 :*
88 2167 1628 :*
90 1530 1260 :*
92 1230 975 :* :======================================*=
94 1078 754 :* :==============================*=========
96 742 584 :* :=======================*======
98 676 452 :* :==================*=========
100 512 349 :* :=============*=======
102 430 270 :* :==========*=======
104 294 209 :* :========*===
106 236 162 :* :======*===
108 191 125 :* :====*===
110 133 97 :* :===*==
112 110 75 :* :==*==
114 99 58 :* :==*=
116 66 45 :* :=*=
118 67 35 :* :=*=
>120 1129 27 :* :=*======================================
Results sorted and z-values calculated from opt score
1955 scores saved that exceeded 106
111754 optimizations performed
Joining threshold: 36, optimization threshold: 36, opt. width: 16
The best scores are: frame init1 initn opt z-sc E(...)..
GB_PR:HUMHBGG Begin: 18 End: 458
! M15386 Human glycine-gamma-globin, ...(3) 971 971 971 1707.5 0
GB_PR:HSGGGPHG Begin: 2 End: 382
! X55656 H.sapiens mRNA for gamma-G g...(2) 843 843 843 1483.4 0
GB_RO:MUSGLOBEP Begin: 18 End: 458
! M26897 Mouse epsilon-globin mRNA, c...(3) 760 760 764 1344.0 0
/////////////////////////////////////////////////////////////////////////////
\\End of List
ggamma.pep
GB_PR:HUMHBGG
M15386 Human glycine-gamma-globin, 3' end. 11/94
LOCUS HUMHBGG 545 bp mRNA PRI 08-NOV-1994
DEFINITION Human glycine-gamma-globin, 3' end.
ACCESSION M15386
NID g183884
KEYWORDS gamma-globin; globin. . . .
SCORES (3) Init1: 971 Initn: 971 Opt: 971 z-score: 1707.5 E(): 0
100.0% identity in 147 aa overlap
10 20 30 40 50
ggamma.pep MGHFTEEDKATITSLWGKVNVEDAGGETLGRLLVVYPWTQRFFDSFGNLSSASAI
|||||||||||||||||||||||||||||||||||||||||||||||||||||||
HUMHBGG PSPDAMGHFTEEDKATITSLWGKVNVEDAGGETLGRLLVVYPWTQRFFDSFGNLSSASAI
10 20 30 40 50 60
60 70 80 90 100 110
ggamma.pep MGNPKVKAHGKKVLTSLGDAIKHLDDLKGTFAQLSELHCDKLHVDPENFKLLGNVLVTVL
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
HUMHBGG MGNPKVKAHGKKVLTSLGDAIKHLDDLKGTFAQLSELHCDKLHVDPENFKLLGNVLVTVL
70 80 90 100 110 120
120 130 140
ggamma.pep AIHFGKEFTPEVQASWQKMVTGVASALSSRYHX
|||||||||||||||||||||||||||||||||
HUMHBGG AIHFGKEFTPEVQASWQKMVTGVASALSSRYHXARCPXCRASRIGFILQAIQIINLFCXE
130 140 150 160 170 180
///////////////////////////////////////////////////////////////////
! CPU time used:
! Database scan: 0:12:04.5
! Post-scan processing: 0:00:02.2
! Total CPU time: 0:12:06.8
! Output File: ggamma.tfasta
The first part of the output file contains a histogram showing the distribution of the z-scores between the query and search set sequences. (See the ALGORITHM topic for an explanation of z-score.) The histogram is composed of bins of size 2 that are labeled according to the higher score for that bin (the leftmost column of the histogram). For example, the bin labeled 24 stores the number of sequence pairs that had scores of 23 or 24.
The next two columns of the histogram list the number of z-scores that fell within each bin. The second column lists the number of z-scores observed in the search and the third column lists the number of z-scores that were expected.
The body of the histogram displays a graphical representation of the score distributions. Equal signs (=) indicate the number of scores of that magnitude that were observed during the search, while asterisks (*) plot the number of scores of that magnitude that were expected.
At the bottom of the histogram is a list of some of the parameters pertaining to the search. These are displayed even if the histogram itself has been suppressed by -NOHIStogram.
Below the histogram, TFastA displays a listing of the best scores. This listing includes the reading frame in the original nucleotide sequence from which the reported translated sequence is derived.
Following the list of best scores, TFastA displays the alignments of the regions of best overlap between the query and search sequences. In these alignments, stop codons are represented by the letter X.
This program displays only the region of overlap between the two aligned sequences (plus some residues on either side of the region to provide context for the alignment) unless you put -SHOWall on the command line. The display of identities and conservative replacements between the aligned sequences depends on the value of the -MARKx command-line parameter. By default ( -MARKx=3), the pipe character (|) is used to denote identities and the colon (:) to denote conservative replacements.
TFastA accepts a single protein sequence as the query sequence. The search set is either a single nucleic acid sequence or multiple nucleic acid sequences. You can specify multiple sequences in a number of ways: by using a list file, for example @project.list; by using an MSF or RSF file, for example project.msf{*}; or by using a sequence specification with an asterisk (*) wildcard, for example GenEMBL:*. If TFastA rejects your protein sequence, turn to Appendix VI to see how to change or set the type of a sequence.
WordSearch identifies sequences in the database that share large numbers of common words in the same register of comparison with your query sequence. The output of WordSearch can be displayed with Segments.
Segments aligns and displays the segments of similarity found by WordSearch.
If you run Compare with the command-line parameter -WORd, the program calculates the points for a dot plot that shows where common words between two sequences occur.
ProfileSearch uses a profile (representing a group of aligned sequences) as a query to search the database for new sequences with similarity to the group. The profile is created with the program ProfileMake.
BLAST searches for sequences similar to a query sequence. The query and the database searched can be either peptide or nucleic acid in any combination. BLAST can search databases on your own computer or databases maintained at the National Center for Biotechnology Information (NCBI) in Bethesda, Maryland, USA.
FastA does a Pearson and Lipman search for similarity between a query sequence and a group of sequences of the same type (nucleic acid or protein). For nucleotide searches, FastA may be more sensitive than BLAST.
FrameSearch searches a group of protein sequences for similarity to one or more nucleotide query sequences, or searches a group of nucleotide sequences for similarity to one or more protein query sequences. For each sequence comparison, the program finds an optimal alignment between the protein sequence and all possible codons on each strand of the nucleotide sequence. Optimal alignments may include reading frame shifts.
FindPatterns, StringSearch, LookUp, and Names are other sequence-identification programs.
PepData translates DNA sequence(s) in all six frames.
The query sequence may not be longer than 32,000 symbols. You cannot select a list size of more than 1,000 best scores nor view more than 1,000 alignments. The word size must be either 1 or 2.
For the estimates of statistical significance to be valid, the search set must contain a large sample of unrelated sequences. The statistical estimates will not be calculated at all if there are fewer than 60 scores saved (10 sequences in the search set when all six frames are searched, or 20 sequences when three frames are searched).
If -NOOPTall is specified on the command line, the estimates of statistical significance will not be accurate.
For a description of the algorithm, see the FastA program documentation.
The E() values are affected by the similarities in sequence composition between the query sequence and the search set sequence. Unrelated sequences may have "significant" scores because of composition bias.
TFastA treats each reading frame as a different sequence. If a nucleotide sequence contains a gene coding for a protein similar to your query, but with an intervening sequence that changes the reading frame, the output displays different alignments, one for each reading frame in which there is overlap.
If there is a database sequence with several overlaps in each reading frame, only the best overlap in each reading frame appears in the alignment display.
TFastA translates stop codons in search set sequences to the sequence symbol X.
There are two ways to control the size of the list of best scores. By default TFastA will list scores until a specific E() value is reached. You may set the value by typing it in at the prompt or by using the -EXPect parameter; otherwise the program uses 10.0 for protein searches, 2.0 for nucleic acid searches. (If you are running the program interactively, it will show no more than 40 scores initially, and ask if you want to see more scores if there are any more that are less than the E() value.)
If -NOOPTall is on the command line or if the list size is specified on the command line (for example, -LIStsize=40), the E() value is ignored, and the program will list either the number of scores you requested or 40 scores if -NOOPTall is specified alone. If you are running the program interactively, it will then ask if you want to see more scores, up to the maximum of 1000 scores.
You can control the number of alignments using the -NOALIgn and -ALIgn= command-line parameters. The program behaves differently depending on whether it is being run noninteractively (in batch or with -Default on the command line) or interactively. In the noninteractive case, the program displays the number of alignments set by the -ALIgn= parameter. (If this is not present, it shows 40 alignments or the number of scores that were listed, whichever is smaller.) If you run the program interactively, it displays the list of best scores, then asks you how many alignments you want to see. This allows you to override the -ALIgn= command-line parameter if you see that you need more (or fewer) alignments than you had originally thought. (This prompt does not appear if -NOALIgn is on the command line.)
By default, TFastA uses a word size of 2. Using a word size of 1 increases the sensitivity at the expense of dramatically increasing the amount of CPU time required to run the program. A word size of 1 should be used if the query sequence is a short peptide.
TFastA chooses default gap creation and extension penalties that are appropriate for the scoring matrix it reads. If you select a different scoring matrix with the -MATRix command-line parameter, the program will adjust the default gap penalties accordingly. (See Appendix VII for information about how to set the default gap penalties for any scoring matrix.) The histogram display gives a qualitative view of the quality of fit between the actual distribution of scores and the expected distribution of scores. This information may indicate whether or not suitable gap creation and extension penalties were used for the search. When the histogram shows poor agreement between the actual distribution and the theoretical distribution, you might consider using -GAPweight and -LENgthweight to specify higher gap creation and extension penalties, respectively. For example, you might increase the gap creation penalty from 16 to 20 and the gap extension penalty from 4 to 6.
If you want to search a single database division instead of an entire database, see the "Using Database Sequences" topic of Chapter 2, Using Sequence Files and Databases of the User's Guide for a list of the logical names used for the databases and the divisions of each database. The search set can also consist of a group of sequence files that are not in a database. Use a multiple sequence specification to name these. For information about naming groups of sequences for the search set, see the topics "Specifying Files" and "Using Wildcards" in Chapter 1, Getting Started, and "Using Database Sequences," "Using Multiple Sequence Format (MSF) Files", "Using Rich Sequence Format (RSF) Files", and "Using List Files" in Chapter 2, Using Sequence Files and Databases of the User's Guide.
TFastA is one of the few programs in the Wisconsin Package(TM) that can take more than a few minutes to run. Most comparisons should probably be run in the batch queue. You can specify that this program run at a later time in the batch queue by using the command-line parameter -BATch. Run this way, the program prompts you for all the required parameters and then automatically submits itself to the batch or at queue. For more information, see "Using the Batch Queue" in Chapter 3, Using Programs in the User's Guide. Very large comparisons may exceed the CPU limit set by some systems.
You can type <Ctrl>C to interrupt a search and see the results from the part of the search that has already been completed.
FastA and TFastA were written by Professor William Pearson of the University of Virginia Department of Biochemistry (Pearson and Lipman, Proc. Natl. Acad. Sci., USA 85; 2444-2448 (1988)). In collaboration with Professor Pearson, they were modified and documented for distribution with GCG Version 6.1 by Mary Schultz and Irv Edelman, and for Versions 8 and 9 by Sue Olson.
All parameters for this program may be put on the command line. Use the parameter -CHEck to see the summary below and to have a chance to add things to the command line before the program executes. In the summary below, the capitalized letters in the parameter names are the letters that you must type in order to use the parameter. Square brackets ([ and ]) enclose parameter values that are optional. For more information, see "Using Program Parameters" in Chapter 3, Using Programs in the User's Guide.
Minimal Syntax: % tfasta [-INfile1=]ggamma.pep -Default
Prompted Parameters:
[-INfile2=]GenEMBL:* search set (all of GenEMBL)
[-OUTfile=]ggamma.tfasta output file name
-BEGin=1 -END=148 range of interest
-WORdsize=2 word size
-EXPect=2.0 lists scores until E() value reaches 2.0
Local Data Files:
-MATRix=blosum50.cmp scoring matrix for peptides
Optional Parameters:
-GAPweight=16 gap creation penalty
-LENgthweight=4 gap extension penalty
-SINce=6.90 limits search to sequences dated on or after June 1990
-THREEFrames translates and searches only the three forward reading
frames
-FRAme=1 translates and searches only the frame specified.
-NOPAMfactor uses a constant factor to calculate initial diagonal scores
-LIStsize=40 shows the best 40 scores (overrides EXPect)
-NOATTRibutes suppresses writing the Begin, End, and Strand
list attributes to the list of best scores
-ALIgn=20 shows the best 20 alignments
-NOALIgn suppresses sequence alignments
-OPTall=20 immediately computes opt score when the initn score is 20
or higher; sorts on opt score
-NOOPTall doesn't compute opt score during search; sorts on initn
-SWalign does final alignment as Smith-Waterman
-SHOWall shows complete sequences in alignment, not just overlaps
-MARKx=3 determines the alignment display mode
-NOHIStogram suppresses printing the histogram
-LINEsize=60 number of sequence symbols per line of the alignment
-NODOCLines suppresses sequence documentation in the alignment
-NOMONitor suppresses the screen trace for each search set sequence
-BATch submits the program to run in the batch queue
-MINLength=1000 searches only sequences of 1000 or more residues
-MAXLength=5000 searches only sequences of 5000 or fewer residues
The files described below supply auxiliary data to this program. The program automatically reads them from a public data directory unless you either 1) have a data file with exactly the same name in your current working directory; or 2) name a file on the command line with an expression like -DATa1=myfile.dat. For more information see Chapter 4, Using Data Files in the User's Guide.
This program reads one or more scoring matrices for the comparison of sequence characters. The program automatically reads the program default scoring matrix file in a public data directory unless you either 1) have a data file with exactly the same name as the program default scoring matrix in your current working directory; or 2) have a data file with exactly the same name as the program default scoring matrix in the directory with the logical name MyData; or 3) name a file on the command line with an expression like -MATRix=mymatrix.cmp. If you don't include a directory specification when you name a file on the command line with -MATRix, the program searches for the file first in your local directory, then in the directory with the logical name MyData, then in the public data directory with the logical name GenMoreData, and finally in the public data directory with the logical name GenRunData. For more information see "Using a Special Kind of Data File: A Scoring Matrix" in Chapter 4, Using Data Files in the User's Guide.
TFastA reads a scoring matrix containing the values for every possible match from your working directory or the public database. The default matrix is blosum50.cmp, which is a BLOSUM50 matrix. You can use the Fetch program to obtain a copy of this file if you need to modify it for your own needs.
The parameters listed below can be set from the command line. For more information, see "Using Program Parameters" in Chapter 3, Using Programs in the User's Guide.
allows you to specify a scoring matrix file name other than the program default. If you don't include a directory specification when you name a file on the command line with -MATRix, the program searches for the file first in your local directory, then in the directory with the logical name MyData, then in the public data directory with the logical name GenMoreData, and finally in the public data directory with the logical name GenRunData. For more information see the Local Scoring Matrices topic above.
limits the search to sequences that have been entered into the database or modified since June 1990. As this is being written, only the EMBL, GenBank, and SWISS-PROT databases support this parameter.
translates and searches only the three forward reading frames.
limits the translation and search to the first reading frame. The default is to translate and search all six frames.
uses a constant factor for the calculation of initial diagonal scores, instead of using the identical match scores from the scoring matrix.
specifies the gap creation penalty that is subtracted from the alignment score whenever a gap is created.
specifies the gap extension penalty that is subtracted from the alignment score for each residue added to an existing gap.
shows all scores whose E() value is less than 2.0. Ignored if -LIStsize or -NOOPTall is on the command line.
shows the best 40 scores. Overrides -EXPect.
suppresses writing to the list of best scores the Begin, End, and Strand attributes that indicate the region of the search set sequence that was aligned with the query sequence.
limits the number of alignments to display in the output file to the 10 best-scoring regions in the list.
suppresses the sequence alignments in the output file. The resulting output file can be used as a list file for input to other Wisconsin Package programs.
immediately performs an alignment and calculates the opt score when the initn score is greater than the specified threshold score. This parameter allows you to override the default threshold calculated by the program. Scores are sorted and saved by opt score during the search.
doesn't compute the opt score until the search is complete. Scores are sorted and saved by initn score instead of by opt score.
does an unlimited Smith-Waterman alignment as the final alignment for nucleotide searches and TFastA searches, instead of the "alignment in a band" version of Smith-Waterman. (Note: this can be very slow.)
shows entire sequences in the alignment display, instead of just the best region of overlap and its surroundings.
determines the alignment display mode -- especially the symbols that identify matches and mismatches. The default value, 3, uses a pipe character (|) to show identities and a colon (:) to show conservative replacements. -MARKx=0 uses a colon to show identities and a period (.) to show conservative replacements. -MARKx=1 will not mark identities; instead, conservative replacements are connected with a lowercase x, and non-conservative substitutions are connected with an uppercase X. If -MARKx=2, the residues in the second sequence are shown only if they differ from the first sequence.
Use -MARKx=10 to get aligned sequences in the FastA "parsable" output format. A document describing this format appears after FastA in the Program Manual.
suppresses printing the histogram.
lets you set the number of sequence symbols in each line of the alignment to any number between 60 and 200.
suppresses the documentation from the search set sequence accompanying the alignment in the output file. Use -DOCLines=5 to copy only five non-blank lines of documentation.
restricts the search to search set sequences that are equal to or longer than 1000 residues.
restricts the search to search set sequences that are equal to or shorter than 5000 residues.
submits the program to the batch queue for processing after prompting you for all required user inputs. Any information that would normally appear on the screen while the program is running is written into a log file. Whether that log file is deleted, printed, or saved to your current directory depends on how your system manager has set up the command that submits this program to the batch queue. All output files are written to your current directory, unless you direct the output to another directory when you specify the output file.
monitors this program's progress on your screen. Use this parameter to see this same monitor in the log file for a batch process. If the monitor is slowing down the program because your terminal is connected to a slow modem, suppress it with -NOMONitor.
The monitor is updated every time the program processes 100 sequences or files. You can use a value after the parameter to set this monitoring interval to some other number.
[ Program Manual | User's Guide | Data Files | Databases ]
Documentation Comments: doc-comments@gcg.com
Technical Support: help@gcg.com
Copyright (c) 1982, 1983, 1985, 1986, 1987, 1989, 1991, 1994, 1995, 1996, 1997 Genetics Computer Group, Inc. a wholly owned subsidiary of Oxford Molecular Group, Inc. All rights reserved.
Licenses and Trademarks Wisconsin Package is a trademark of Genetics Computer Group, Inc. GCG and the GCG logo are registered trademarks of Genetics Computer Group, Inc.
All other product names mentioned in this documentation may be trademarks, and if so, are trademarks or registered trademarks of their respective holders and are used in this documentation for identification purposes only.