A local alignment searches for regions of local similarity between two sequences and need not include the entire length of the sequences. Local alignment methods are very useful for scanning databases or other circumsatnces when you wish to find matches between small regions of sequences, for example between protein domains.
Dynamic programming methods ensure the optimal local alignment by exploring all possible alignments and choosing the best. It does this by reading in a scoring matrix that contains values for every possible residue or nucleotide match. water finds an alignment with the maximum possible score where the score of an alignment is equal to the sum of the matches taken from the scoring matrix.
An important problem is the treatment of gaps, i.e., spaces inserted to optimise the alignment score. A penalty is subtracted from the score for each gap opened (the 'gap open' penalty) and a penalty is subtracted from the score for the total number of gap spaces multiplied by a cost (the 'gap extension' penalty).
Typically, the cost of extending a gap is set to be 5-10 times lower than the cost for opening a gap.
There are two ways to compute a penalty for a gap of n positions :
gap opening penalty + (n - 1) * gap extension penalty gap penalty + n * gap length penalty
The first way is used by EMBOSS and WU-BLAST
The second way is used by NCBI-BLAST, GCG, Staden and CLUSTAL. Fasta used it for a long time the first way, but Prof. Pearson decided recently to shift to the second.
The two methods are basically equivalent.
The Smith-Waterman algorithm contains no negative scores in the path matrix it creates. The algorithm starts the alignment at the highest path matrix score and works backwards until a cell contains zero.
See the Reference Smith et al. for details.
% water tsw:hba_human tsw:hbb_human Smith-Waterman local alignment. Gap opening penalty [10.0]: Gap extension penalty [0.5]: Output alignment [hba_human.water]:
Go to the input files for this example
Go to the output files for this example
Standard (Mandatory) qualifiers: [-asequence] sequence Sequence USA [-bsequence] seqall Sequence database USA -gapopen float The gap open penalty is the score taken away when a gap is created. The best value depends on the choice of comparison matrix. The default value assumes you are using the EBLOSUM62 matrix for protein sequences, and the EDNAFULL matrix for nucleotide sequences. -gapextend float The gap extension penalty is added to the standard gap penalty for each base or residue in the gap. This is how long gaps are penalized. Usually you will expect a few long gaps rather than many short gaps, so the gap extension penalty should be lower than the gap penalty. An exception is where one or both sequences are single reads with possible sequencing errors in which case you would expect many single base gaps. You can get this result by setting the gap open penalty to zero (or very low) and using the gap extension penalty to control gap scoring. [-outfile] align Output alignment file name Additional (Optional) qualifiers: -datafile matrixf This is the scoring matrix file used when comparing sequences. By default it is the file 'EBLOSUM62' (for proteins) or the file 'EDNAFULL' (for nucleic sequences). These files are found in the 'data' directory of the EMBOSS installation. Advanced (Unprompted) qualifiers: -[no]brief boolean Brief identity and similarity Associated qualifiers: "-asequence" associated qualifiers -sbegin1 integer Start of the sequence to be used -send1 integer End of the sequence to be used -sreverse1 boolean Reverse (if DNA) -sask1 boolean Ask for begin/end/reverse -snucleotide1 boolean Sequence is nucleotide -sprotein1 boolean Sequence is protein -slower1 boolean Make lower case -supper1 boolean Make upper case -sformat1 string Input sequence format -sdbname1 string Database name -sid1 string Entryname -ufo1 string UFO features -fformat1 string Features format -fopenfile1 string Features file name "-bsequence" associated qualifiers -sbegin2 integer Start of each sequence to be used -send2 integer End of each sequence to be used -sreverse2 boolean Reverse (if DNA) -sask2 boolean Ask for begin/end/reverse -snucleotide2 boolean Sequence is nucleotide -sprotein2 boolean Sequence is protein -slower2 boolean Make lower case -supper2 boolean Make upper case -sformat2 string Input sequence format -sdbname2 string Database name -sid2 string Entryname -ufo2 string UFO features -fformat2 string Features format -fopenfile2 string Features file name "-outfile" associated qualifiers -aformat3 string Alignment format -aextension3 string File name extension -adirectory3 string Output directory -aname3 string Base file name -awidth3 integer Alignment width -aaccshow3 boolean Show accession number in the header -adesshow3 boolean Show description in the header -ausashow3 boolean Show the full USA in the alignment -aglobal3 boolean Show the full sequence in alignment General qualifiers: -auto boolean Turn off prompts -stdout boolean Write standard output -filter boolean Read standard input, write standard output -options boolean Prompt for standard and additional values -debug boolean Write debug output to program.dbg -verbose boolean Report some/full command line options -help boolean Report command line options. More information on associated and general qualifiers can be found with -help -verbose -warning boolean Report warnings -error boolean Report errors -fatal boolean Report fatal errors -die boolean Report deaths
|Standard (Mandatory) qualifiers||Allowed values||Default|
|Sequence USA||Readable sequence||Required|
|Sequence database USA||Readable sequence(s)||Required|
|-gapopen||The gap open penalty is the score taken away when a gap is created. The best value depends on the choice of comparison matrix. The default value assumes you are using the EBLOSUM62 matrix for protein sequences, and the EDNAFULL matrix for nucleotide sequences.||Number from 0.000 to 100.000||10.0 for any sequence|
|-gapextend||The gap extension penalty is added to the standard gap penalty for each base or residue in the gap. This is how long gaps are penalized. Usually you will expect a few long gaps rather than many short gaps, so the gap extension penalty should be lower than the gap penalty. An exception is where one or both sequences are single reads with possible sequencing errors in which case you would expect many single base gaps. You can get this result by setting the gap open penalty to zero (or very low) and using the gap extension penalty to control gap scoring.||Number from 0.000 to 10.000||0.5 for any sequence|
|Output alignment file name||Alignment output file|
|Additional (Optional) qualifiers||Allowed values||Default|
|-datafile||This is the scoring matrix file used when comparing sequences. By default it is the file 'EBLOSUM62' (for proteins) or the file 'EDNAFULL' (for nucleic sequences). These files are found in the 'data' directory of the EMBOSS installation.||Comparison matrix file in EMBOSS data path||EBLOSUM62 for protein
EDNAFULL for DNA
|Advanced (Unprompted) qualifiers||Allowed values||Default|
|-[no]brief||Brief identity and similarity||Boolean value Yes/No||Yes|
ID HBA_HUMAN STANDARD; PRT; 141 AA. AC P01922; DT 21-JUL-1986 (Rel. 01, Created) DT 21-JUL-1986 (Rel. 01, Last sequence update) DT 15-JUL-1999 (Rel. 38, Last annotation update) DE HEMOGLOBIN ALPHA CHAIN. GN HBA1 AND HBA2. OS Homo sapiens (Human), Pan troglodytes (Chimpanzee), and OS Pan paniscus (Pygmy chimpanzee) (Bonobo). OC Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Mammalia; OC Eutheria; Primates; Catarrhini; Hominidae; Homo. RN  RP SEQUENCE FROM N.A. (ALPHA-1). RX MEDLINE; 81088339. RA MICHELSON A.M., ORKIN S.H.; RT "The 3' untranslated regions of the duplicated human alpha-globin RT genes are unexpectedly divergent."; RL Cell 22:371-377(1980). RN  RP SEQUENCE FROM N.A. (ALPHA-2). RX MEDLINE; 81175088. RA LIEBHABER S.A., GOOSSENS M.J., KAN Y.W.; RT "Cloning and complete nucleotide sequence of human 5'-alpha-globin RT gene."; RL Proc. Natl. Acad. Sci. U.S.A. 77:7054-7058(1980). RN  RP SEQUENCE FROM N.A. (ALPHA-2). RX MEDLINE; 80137531. RA WILSON J.T., WILSON L.B., REDDY V.B., CAVALLESCO C., GHOSH P.K., RA DERIEL J.K., FORGET B.G., WEISSMAN S.M.; RT "Nucleotide sequence of the coding portion of human alpha globin RT messenger RNA."; RL J. Biol. Chem. 255:2807-2815(1980). RN  RP SEQUENCE FROM N.A. (ALPHA-1 AND ALPHA-2). RA FLINT J., HIGGS D.R.; RL Submitted (JAN-1997) to the EMBL/GenBank/DDBJ databases. RN  RP SEQUENCE. RA BRAUNITZER G., GEHRING-MULLER R., HILSCHMANN N., HILSE K., HOBOM G., RA RUDLOFF V., WITTMANN-LIEBOLD B.; RT "The constitution of normal adult human haemoglobin."; RL Hoppe-Seyler's Z. Physiol. Chem. 325:283-286(1961). RN  RP SEQUENCE. RA HILL R.J., KONIGSBERG W.; RT "The structure of human hemoglobin: IV. The chymotryptic digestion of RT the alpha chain of human hemoglobin."; RL J. Biol. Chem. 237:3151-3156(1962). RN  [Part of this file has been deleted for brevity] FT /FTId=VAR_002841. FT VARIANT 130 130 A -> D (IN YUDA; O2 AFFINITY DOWN). FT /FTId=VAR_002842. FT VARIANT 131 131 S -> P (IN QUESTEMBERT; HIGHLY UNSTABLE; FT CAUSES ALPHA-THALASSEMIA). FT /FTId=VAR_002843. FT VARIANT 133 133 S -> R (IN VAL DE MARNE; O2 AFFINITY UP). FT /FTId=VAR_002844. FT VARIANT 135 135 V -> E (IN PAVIE). FT /FTId=VAR_002845. FT VARIANT 136 136 L -> M (IN CHICAGO). FT /FTId=VAR_002846. FT VARIANT 136 136 L -> P (IN BIBBA; UNSTABLE; FT CAUSES ALPHA-THALASSEMIA). FT /FTId=VAR_002847. FT VARIANT 138 138 S -> P (IN ATTLEBORO; O2 AFFINITY UP). FT /FTId=VAR_002848. FT VARIANT 139 139 K -> E (IN HANAKAMI; O2 AFFINITY UP). FT /FTId=VAR_002849. FT VARIANT 139 139 K -> T (IN TOKONAME; O2 AFFINITY UP). FT /FTId=VAR_002850. FT VARIANT 140 140 Y -> H (IN ROUEN; O2 AFFINITY UP). FT /FTId=VAR_002851. FT VARIANT 141 141 R -> C (IN NUNOBIKI; O2 AFFINITY UP). FT /FTId=VAR_002852. FT VARIANT 141 141 R -> L (IN LEGNANO; O2 AFFINITY UP). FT /FTId=VAR_002853. FT VARIANT 141 141 R -> H (IN SURESNES; O2 AFFINITY UP). FT /FTId=VAR_002854. FT VARIANT 141 141 R -> P (IN SINGAPORE). FT /FTId=VAR_002855. FT HELIX 4 35 FT HELIX 37 42 FT TURN 44 45 FT TURN 50 51 FT HELIX 53 71 FT TURN 72 74 FT HELIX 76 79 FT TURN 80 80 FT HELIX 81 89 FT TURN 90 91 FT TURN 95 95 FT HELIX 96 112 FT TURN 114 116 FT HELIX 119 136 FT TURN 137 139 SQ SEQUENCE 141 AA; 15126 MW; 5EC7DB1E CRC32; VLSPADKTNV KAAWGKVGAH AGEYGAEALE RMFLSFPTTK TYFPHFDLSH GSAQVKGHGK KVADALTNAV AHVDDMPNAL SALSDLHAHK LRVDPVNFKL LSHCLLVTLA AHLPAEFTPA VHASLDKFLA SVSTVLTSKY R //
ID HBB_HUMAN STANDARD; PRT; 146 AA. AC P02023; DT 21-JUL-1986 (Rel. 01, Created) DT 21-JUL-1986 (Rel. 01, Last sequence update) DT 15-JUL-1999 (Rel. 38, Last annotation update) DE HEMOGLOBIN BETA CHAIN. GN HBB. OS Homo sapiens (Human), Pan troglodytes (Chimpanzee), and OS Pan paniscus (Pygmy chimpanzee) (Bonobo). OC Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Mammalia; OC Eutheria; Primates; Catarrhini; Hominidae; Homo. RN  RP SEQUENCE. RC SPECIES=HUMAN; RA BRAUNITZER G., GEHRING-MULLER R., HILSCHMANN N., HILSE K., HOBOM G., RA RUDLOFF V., WITTMANN-LIEBOLD B.; RT "The constitution of normal adult human haemoglobin."; RL Hoppe-Seyler's Z. Physiol. Chem. 325:283-286(1961). RN  RP SEQUENCE FROM N.A. RC SPECIES=HUMAN; RX MEDLINE; 81064667. RA LAWN R.M., EFSTRATIADIS A., O'CONNELL C., MANIATIS T.; RT "The nucleotide sequence of the human beta-globin gene."; RL Cell 21:647-651(1980). RN  RP SEQUENCE OF 121-146 FROM N.A. RC SPECIES=HUMAN; RX MEDLINE; 85205333. RA LANG K.M., SPRITZ R.A.; RT "Cloning specific complete polyadenylylated 3'-terminal cDNA RT segments."; RL Gene 33:191-196(1985). RN  RP X-RAY CRYSTALLOGRAPHY (2.5 ANGSTROMS) OF DEOXYHEMOGLOBIN. RC SPECIES=HUMAN; RX MEDLINE; 76027820. RA FERMI G.; RT "Three-dimensional fourier synthesis of human deoxyhaemoglobin at RT 2.5-A resolution: refinement of the atomic model."; RL J. Mol. Biol. 97:237-256(1975). RN  RP SEQUENCE. RC SPECIES=P.TROGLODYTES; RX MEDLINE; 66071496. RA RIFKIN D.B., KONIGSBERG W.; RT "The characterization of the tryptic peptides from the hemoglobin of RT the chimpanzee (Pan troglodytes)."; RL Biochim. Biophys. Acta 104:457-461(1965). RN  [Part of this file has been deleted for brevity] FT VARIANT 140 140 A -> T (IN ST JACQUES: O2 AFFINITY UP). FT /FTId=VAR_003081. FT VARIANT 140 140 A -> V (IN PUTTELANGE; POLYCYTHEMIA; FT O2 AFFINITY UP). FT /FTId=VAR_003082. FT VARIANT 141 141 L -> R (IN OLMSTED; UNSTABLE). FT /FTId=VAR_003083. FT VARIANT 142 142 A -> D (IN OHIO; O2 AFFINITY UP). FT /FTId=VAR_003084. FT VARIANT 143 143 H -> D (IN RANCHO MIRAGE). FT /FTId=VAR_003085. FT VARIANT 143 143 H -> Q (IN LITTLE ROCK; O2 AFFINITY UP). FT /FTId=VAR_003086. FT VARIANT 143 143 H -> P (IN SYRACUSE; O2 AFFINITY UP). FT /FTId=VAR_003087. FT VARIANT 143 143 H -> R (IN ABRUZZO; O2 AFFINITY UP). FT /FTId=VAR_003088. FT VARIANT 144 144 K -> E (IN MITO; O2 AFFINITY UP). FT /FTId=VAR_003089. FT VARIANT 145 145 Y -> C (IN RAINIER; O2 AFFINITY UP). FT /FTId=VAR_003090. FT VARIANT 145 145 Y -> H (IN BETHESDA; O2 AFFINITY UP). FT /FTId=VAR_003091. FT VARIANT 146 146 H -> D (IN HIROSHIMA; O2 AFFINITY UP). FT /FTId=VAR_003092. FT VARIANT 146 146 H -> L (IN COWTOWN; O2 AFFINITY UP). FT /FTId=VAR_003093. FT VARIANT 146 146 H -> P (IN YORK; O2 AFFINITY UP). FT /FTId=VAR_003094. FT VARIANT 146 146 H -> Q (IN KODAIRA; O2 AFFINITY UP). FT /FTId=VAR_003095. FT HELIX 5 15 FT TURN 16 17 FT HELIX 20 34 FT HELIX 36 41 FT HELIX 43 45 FT HELIX 51 55 FT TURN 56 56 FT HELIX 58 75 FT TURN 76 77 FT HELIX 78 94 FT TURN 95 96 FT TURN 100 100 FT HELIX 101 121 FT HELIX 124 142 FT TURN 143 144 SQ SEQUENCE 146 AA; 15867 MW; EC9744C9 CRC32; VHLTPEEKSA VTALWGKVNV DEVGGEALGR LLVVYPWTQR FFESFGDLST PDAVMGNPKV KAHGKKVLGA FSDGLAHLDN LKGTFATLSE LHCDKLHVDP ENFRLLGNVL VCVLAHHFGK EFTPPVQAAY QKVVAGVANA LAHKYH //
The output is a standard EMBOSS alignment file.
The results can be output in one of several styles by using the command-line qualifier -aformat xxx, where 'xxx' is replaced by the name of the required format. Some of the alignment formats can cope with an unlimited number of sequences, while others are only for pairs of sequences.
The available multiple alignment format names are: unknown, multiple, simple, fasta, msf, trace, srs
The available pairwise alignment format names are: pair, markx0, markx1, markx2, markx3, markx10, srspair, score
See: http://emboss.sf.net/docs/themes/AlignFormats.html for further information on alignment formats.
######################################## # Program: water # Rundate: Fri Jul 15 2005 12:00:00 # Align_format: srspair # Report_file: hba_human.water ######################################## #======================================= # # Aligned_sequences: 2 # 1: HBA_HUMAN # 2: HBB_HUMAN # Matrix: EBLOSUM62 # Gap_penalty: 10.0 # Extend_penalty: 0.5 # # Length: 145 # Identity: 63/145 (43.4%) # Similarity: 88/145 (60.7%) # Gaps: 8/145 ( 5.5%) # Score: 293.5 # # #======================================= HBA_HUMAN 2 LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHF-DLS- 49 |:|.:|:.|.|.|||| :..|.|.|||.|:.:.:|.|:.:|..| ||| HBB_HUMAN 3 LTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRFFESFGDLST 50 HBA_HUMAN 50 ----HGSAQVKGHGKKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDP 95 .|:.:||.|||||..|.::.:||:|::....:.||:||..||.||| HBB_HUMAN 51 PDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDP 100 HBA_HUMAN 96 VNFKLLSHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKY 140 .||:||.:.|:..||.|...||||.|.|:..|.:|.|:..|..|| HBB_HUMAN 101 ENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVANALAHKY 145 #--------------------------------------- #---------------------------------------
The Identity: is the percentage of identical matches between the two sequences over the reported aligned region (including any gaps in the length).
The Similarity: is the percentage of matches between the two sequences over the reported aligned region (including any gaps in the length).
EMBOSS data files are distributed with the application and stored in the standard EMBOSS data directory, which is defined by the EMBOSS environment variable EMBOSS_DATA.
To see the available EMBOSS data files, run:
% embossdata -showall
To fetch one of the data files (for example 'Exxx.dat') into your current directory for you to inspect or modify, run:
% embossdata -fetch -file Exxx.dat
Users can provide their own data files in their own directories. Project specific files can be put in the current directory, or for tidier directory listings in a subdirectory called ".embossdata". Files for all EMBOSS runs can be put in the user's home directory, or again in a subdirectory called ".embossdata".
The directories are searched in the following order:
water is for aligning the best matching subsequences of two sequences. It does not necessarily align whole sequences against each other; you should use needle if you wish to align closely related sequences along their whole lengths.
A true Smith Waterman implementation like water needs memory proportional to the product of the sequence lengths. For two sequences of length 10,000,000 and 1,000 it therefore needs memory proportional to 10,000,000,000 characters. Two arrays of this size are produced, one of ints and one of floats so multiply that figure by 8 to get the memory usage in bytes. That doesn't include other overheads. Therefore only use water and needle for accurate alignment of reasonably short sequences.
If you run out of memory, try using supermatcher or matcher.
Uncaught exception Assertion failed raised at ajmem.c:xxx
Probably means you have run out of memory. Try using supermatcher or matcher if this happens.
|matcher||Finds the best local alignments between two sequences|
|seqmatchall||All-against-all comparison of a set of sequences|
|supermatcher||Match large sequences against one or more other sequences|
|wordmatch||Finds all exact matches of a given size between 2 sequences|
matcher is a local alignment program that gives as good an alignment as
waterbut it uses far less memory. However,
waterruns twice as fast as matcher.
supermatcher is designed for local alignments of very large sequences. It gives good results as long as there is not significant amounts of insertion or deletion in the alignment.
supermatcher Finds a match of a large sequence against one or more sequences matcher Finds the best local alignments between two sequences
Modified 27th July 1999 - tweaking scoring.
Modified 22 Oct 2000 - added ID and Similarity scores.