Vol.8, no 4. 1992 Pages 3 7 7 - 3 8 3
CABIOS
Sequence pattern matching on a massively parallel computer Robert Jones Abstract
Introduction Consensus patterns serve as compact representations of sequence features and highlight the most important elements. Searching for matches to these patterns within a sequence is an important method in computational molecular biology. A newly sequenced protein of unknown function can be examined for the presence of a variety of ligand-binding sites, with matches serving as important clues to structure and function. Scanning a eukaryotic DNA sequence for matches to consensus donor and acceptor sites for RNA splicing can help uncover intron and exon structure. Comparing a single sequence against a database using dynamic programming is the most sensitive method available for finding weak relationships between sequences but is computationally very expensive (Smith and Waterman, 1981). An alternative is to produce short consensus sequence patterns that are diagnostic for particular protein families and then to use these to identify other members of those families. Several groups have published databases of these consensus patterns (Bairoch, 1990; Smith and Smith, 1990). In many ways sequence pattern matching resembles the problem of finding regular expressions in text, and many of those algorithms have been adapted to the domain of biological sequences. However, there are important features found in sequence patterns that do not occur in regular text. In many positions of a consensus pattern it is necessary to represent conservation of some physical or chemical property, rather than the identity of a specific amino acids. Here one must score a match if any one of a class of characters matches at the given position in the pattern. In addition two regions of a consensus
Thinking Machines Corporation, 245 First Street. Cambridge. MA 02142. USA
© Oxford University Press
Methods In this section a syntax is presented in which complex sequence patterns can be described, followed by the algorithms that enable them to be matched with a sequence database. At present there is no standard sequence pattern syntax in use and no one existing syntax could be found in which to express all the features described here. Accordingly a novel syntax is presented, although this could be readily modified to conform to any emerging standard. The pattern syntax used here can be viewed as having three levels, with the lowest level representing the basic elements. The simplest pattern consists of a contiguous string of character codes for explicit amino acids or nucleotides, 'DEGAG' for example. A match to this pattern is reported when all the individual positions match in the same spatial order to a contiguous region of a database sequence. In the current work, every matching position has a score of one, with mismatched positions scoring zero. Hence an exact match to a pattern of length N has a score of N. Ambiguity can be introduced by defining a threshold score for the pattern and reporting a match only if the accumulated score equals or exceeds the threshold. A threshold score is written after the pattern following a colon. For example the pattern 'DEGAG:4' will match to any region of five amino acids in which four of the five positions match exactly. This ability to express a degree of amgibuity 377
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
A method is described for finding all occurrences of a sequence pattern within a database of molecular sequences. Implementation of this on a massively parallel computer allows the user to perform very fast database searches using complex patterns. In particular, the software supports approximate pattern matching with score thresholds for either the entire pattern or specified elements thereof. Matches to individual elements can be linked by variable length gaps within userspecified limits.
pattern may exhibit variation in their spacing in the examples used to derive the pattern, usually constrained within certain bounds. The syntax of most text-searching methods is unable to represent this. There have been a number of tools developed for searching specifically with sequence patterns. In the area of protein structure prediction very flexible pattern specification languages and search tools have been described (Cohen et al, 1983, 1986; Taylor, 1988, 1989). Several tools have been reported for searching databases with consensus patterns characteristic of protein families, and searching with short sequence fragments forms the basis of a fast database comparison method (Altschul et al, 1990; Smith and Smith, 1990; G.Mehldau and E.Myers, submitted). The current work describes a novel approach to searching sequence databases that uses a massively parallel computer, the Connection Machine CM-2 ®, to achieve very high performance. A rich pattern syntax is provided in which to describe the diversity of observed sequence patterns.
RJones
378
pattern '$ATP-binding' expands to the explicit pattern 'G . . G . GK'. The algorithm described here exploits the architecture of a massively parallel computer, the Connection Machine CM-2, to achieve very high performance,. The CM-2 can be configured with up to 65,536 single-bit processors and up to 128 kbytes of memory per processor. Instructions are broadcast to the processors from a front-end serial computer in a single instruction stream and these are performed in parallel on the data local to each processor. All processors receive all the instructions but only those in which the 'context' flag is set actually operate on their data. Conditional programming constructs are implemented through manipulation of this flag. For many problems the number of items of data greatly exceeds the maximum number of physical processors and this is certainly the case in database searching. The software of the CM-2 can emulate very large machines by having each physical processor simulate a larger number of virtual processors (VPs). This abstraction is made transparent to the programmer and though the execution time increases as the number of virutal processors grows, the user's program remains unchanged. The physical processors are linked by a hypercube communications network. In software, the processors can be configured as Ndimensional arrays, referred to as geometries, and the same physical processors may participate in several different geometries during the course of a program. Communication between processors makes use of the hypercube network in two ways. If the pattern of communication is regular, for instance when all processor / send a value to all processor i — 1, very efficient communications primitives are available. If the pattern is irregular, with the possibility of several processors sending data to a single recipient, then a more general method must be used. The software described here was written in the language C*, a version of the C language that has been extended to exploit massively parallel architectures. In this application the machine is configured as a onedimensional array of virtual processors onto which the sequence database is laid out with one amino acid per processor. Adjacent sequences are delimited by a special character. Consider the task of finding positions in the database that are an exact match to a short simple pattern. Each processor has a score that is initially set to zero. The character at the end of the sequence pattern is broadcast to all processors and those that contain this character increment their score by one. Processors holding the delimiter character reset their score to zero. All processors pass their current score to their neighor below them in order of the array. A fast regular communication primitive can be used as all processors are passing data to the same relative address. The cycle is repeated by broadcasting the second to last character in the pattern to all processors. Those that match the character increment their score, delimiting positions reset their score to zero and all processors pass their score to their lower neighbor.
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
in a pattern match is found to be very useful in the domain of sequence analysis. A position in the pattern that is poorly conserved can be deemphasized by using a period in that position. Any character can match to a period and the match does not contribute to the scoring. Fixed length gaps between pattern elements can be presented by a series of periods, e.g. 'DE . . . GAG'. A character class, which allows a limited group of characters to match at a given pattern position, is represented by enclosing the acceptable characters in square brackets. For example, the pattern [ILV] matches to a single sequence position containing I, L or V. A character class that excludes a set of characters from a position is denoted by preceding the contents of the class with the character '!'. The parser expands this into a character class containing all amino acids not specifically excluded. The second level of syntax allows blocks of these basic pattern elements to be combined in a very flexible manner. This level of the syntax is not explicitly supported by most text-searching algorithms. The use of score thresholds permits ambiguity at the level of the entire pattern. For example, one might wish to combine an exact match to one region of a pattern with an approximate match to a different region. This is achieved by splitting the pattern into a number of 'blocks'. Enclosing part of the pattern in parentheses marks it as a block and it may then be assigned a 'local' score threshold. For example, the pattern '(DEGAG:4) (GT)' requires an exact match to 'GT' and an exact match in at least four positions to the five characters of 'DEGAG'. If the context in which a block is found is unambiguous, then the parsing function is able to recognize the implicit block structure and parentheses are not required. The concept of pattern blocks allows variable spacing between blocks to be defined. This is represented by angle brackets enclosing the inclusive bounds of an acceptable gap. For example the pattern 'DEG AG' will match to a sequence containing the pattern ' D E C separated by between two and five other amino acids from the pattern 'AG'. The third level of pattern syntax provides a number of convenience features. Character classes are used extensively in consensus patterns and in many cases the class of amino acids represents a particular conserved chemical characteristic, e.g. the class '[KR]' can represent conservation of a positive charge and the class '[AGS]' can represent conservation of a small side chain. Rather than having to specify these classes explicitly every time they are employed, the user can define named character classes and store these in a file. Named classes can then be incorporated into patterns by preceding the name with a ' # ' character. For example, the definition '# aliphatic' is expanded to the explicit class '[ILV]'. Entire patterns can be assigned names and then referred to by preceding the name with a '$' character. Predefined patterns are provided for a collection of protein ligand binding sites as well as 'fingerprint' patterns for a number of protein families. For example, the named
Pattern matching on a parallel computer
Pattern : Q S T Processor
10
11
12
13
14
Q
S
15
Sequence Score Step MalchT
2:
ShW score
3:
Match"?
4:
Shin score
5:
Match YJ
Processors with a score equal to the length of the pattern (3) indicate the start of exact matches to the pattern
Fig. 1. Matching the pattern 'QST' to a sequence.
The process continues until the entire pattern has been broadcast. Scores are not passed to neighboring processors in the final iteration. Figure 1 shows the steps involved for a simple example involving 16 processors and a pattern of three characters. For an exact match the threshold score is simply the length of the pattern. All processors with scores less than that are disabled, leaving only those that mark the start of exact matches to the pattern. Accordingly, for simple patterns like this of length n, the algorithm requires n broadcast instructions, n character matches and n — 1 local communications operations. The complexity of the algorithm is O(nP/N), where n is the pattern length, N is the number of physical processors and P is the length of the sequence database. Where an exact match is required, the threshold is simply the number of effective positions in the pattern. Ambiguity is expressed by explicitly specifying a lower threshold score. Processors with scores equal or greater than this threshold mark the start of regions that contain up to the specified number of mismatched positions. These mismatches can occur anywhere within a single pattern block. For example in Figure 1 processor 5 has a final score of 2, indicating that it marks the start of a region (QVT) that contains one mismatch to the supplied pattern (QST). Note that there is no additional cost whatsoever involved in identifying partial matches compared to exact matches. This feature is particularly important for techniques that involve approximate pattern matching, such as the BLAST sequence comparison algorithm (Altschul et ai, 1990). A character class in the pattern is handled by broadcasting and matching all members of the class in turn and only then passing on the current score to the neighboring processors. Accordingly for a class that contains C members there are C broadcast and matching operations but only one communica-
tion step. Positions in the pattern that can match to any character, indicated by a period in the syntax, do not contribute to the threshold score. When such a character is encountered no matching occurs and the current score is passed directly to the neighboring processors. These simple operations implement all the features in the first level of the pattern syntax. The next stage is to return the locations of matching positions to the front-end computer. This requires that the locations and scores of the N matches to the patterns be represented in the first N processors of the array, as opposed to being dispersed throughout the machine. The process is illustrated in Figure 2. The subset of active processors, at which matches start, is enumerated and each processor uses its position within this order as the address of a destination processor, to which it sends both its score and its absolute address in the processor array. For example in Figure 2, processors 1, 9 and 13 all mark the start of a match. Following the enumeration step they carry addresses of 0, 1 and 2 respectively. Scores and absolute coordinates are transferred, using a 'send' operation, to the computed destination processors such that in the figure processor 1 sends its data to processor 0, etc. After this operation the first N processors contain the scores and positions of the N matches to the pattern. These are then read from the machine to the frontend serial computer using an efficient array transfer operation. The titles and sequences of the matching database sequence are then referenced and the search results are displayed for the user. The observed complexity of this stage is approximately O{mlN), where m is the number of matches and N is the number of physical processors. The second level of the pattern syntax provided a way of representing independent pattern blocks that may carry local
379
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
1:
RJones
Processor
0
1
Sequence
D
Disable VPs with scores below the threshold
2".
Enumerate active VPs tor use as a send address
3
Q
0[3i
Rnal Score 1:
2
X
4
5
6
S
T
L
Q
0
0
0
2
X
X
X
7 V
0
X
8 T
9 R
10 Q
11
12
13
14
15
T
K
Q
S
T
0
0
X
X
S
0
0
0
X
X
X
0
X
0
X
3: Compute VP coordinates within the array
1
9
Only the shaded processors are active during these steps
13
The positions of N matches are now In the first N VPs Fig. 2 . Gathering the locations of matches.
Search lor pattern 'QST LIL EDLV
D O S T L Q V T R Q S T K Q S T N L 1: Match block 3 (EDLV)
•
•
2; Match block 2 (LIL)
•
• 1 •
3: Match block 1 (QST)
• «3 '•
IL A C H L
S E D L V F Q S R
1 1
•
1 •
2
4: Unk blocks 2 and 3
5: Link blocks 1 and 2+3
Match 1 Match 2
Q S T
Fig. 3 . Matching pattern blocks linked by variable length gaps.
threshold scores and that may be separated by variable-length spacers. These are first parsed into a series of independent simple patterns, with the relative order and spacing restrictions being recorded. Each basic pattern block is matched to the database independently as described above. Matching positions for each block are flagged and the scores stored in an array on each processor rather than being sent to the front-end. Figure 3 illustrates how these independent matching blocks become linked. First the spacing limits between the end of one block and the start of another are converted into spacing limits between the start positions of adjacent blocks. As with the positions of an individual block, the blocks of a complex pattern are handled
380
in reverse order so that the starts of overall matches are returned at.the end of the procedure. The scores and coordinates of processors that match to the rightmost block (i) are passed in a series of nearest-neighbor communication steps to those processors that fall within the spacing limits for the next lowest pattern block (/ - 1). Any processors within this region that are flagged as matching to block i — 1 are therefore legitmate matches to the pattern segment consisting of blocks i and i — 1. Their score is incremented by that of block /. These processors then continue the procedure by passing their scores and coordinates, as well as the relative coordinates of the associated matches to block j , to those processors that fall
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
4: Send coordinate to send address
Pattern matching on a parallel computer
A Simple Sequence Pattern:
F...D.D..G.X...EL t*cl»
TTHUCS TPTCCI NCUO (CM
CtlBOdull
HCDO HCOO HCUUA HCUUC
CalModuli - Dlctyoic*: IUB ditcoldau* C-lxxh.ll - Olotyoct*. lua discoldaua * P « 2* P ottln - Purpi • M grcblo tpn 2c P ocaln - Pur* • M archIn P«rv«lbuL n beta - Edit 1* froq Parr* 11MM. n beta - b i t l c CKMl P*r**lb-ja o - Tbornback ray C*lclu»-b
CUaoduJi
FVTO
FVCO PVITC RLCII KLM»
sr*
3 •
OMIUTKAA EtHIUTKAA
« 1f * f
uoiurirxA TDtttEIJUA UQIWKU OOSUCLItA
• 1 •
noiuriXA TDTTEE1UXA a WAIT KM
• ctowvrus 7 CTSAADAJOCV 1 ATSADELOL 7 MSOAELAII EHCAKCrifM 11 « XHCOKLDOU If 4 P«tt*rn.
) . c . .o.GJciAcorri |ACG*TI . . . . (ILVI
same protein
Thysldln* Adcoylit*
KIUTt RIHUA KltC A TVHU11A TVKVM TVBYSR
c h U co, • ) - Hi* - lacharlchU coll 1 polyprocvln - 1 — 11 • I roraat v l r p r o t a l n t (tr«n>poBMi Tn7)
bratn
Cilclwfa
PaMM
=j
UttMngSw^i mo* rmrDAOcccoisvuL nwimKicDirrm FSLTDXOCOCTITTIIL ntVTD KDCHCT I >AAI L r ALTD XDCOCTI r m L n v n ROCHC L i I A U L r s LTD KDCOCS i r m L nCVTD KDCMCTICAAIL FKCCCTOCOCK1TS[[L FUIDAOCOGXITPEE L niLDRDUCriEOOIL n i AD ED n c r I u c c L FWVU>GDQ*GY I E V t l L nKYDOOCHCYlDCJICL Ft LYDQOGOCYI DDCt L F. . . D . D . . C . I
CTVTWOCQT CTVHHKLGQT CTVHULCOtl ~ MVHTWLCM CTVmSLCQN MVHTWLCII _ OTVMMLCON ~ MVHTSLGf K _ RAATUIEID KAAntflEID CLTLQWRAJ lUlIATAAD KKTLKCrSDC QAUJCDLCtK DALLXDLCIK
Two domains In the same protein
• 1
21 1
I » t 9
$ 1 SO It
» n
PVHrpACILT cjcuumciir mil KTTttLV VtAKTtTKLV rnfLRTTurv LSIKQGTVIO lAPYlTTWa UMTAOSLL
E V L V V V V r I
OAFGVCVTT GCPCKICT aAKACJCOT CACCVCJtW OAMtVCUA GOGCVGMA GASaSGKST CVPGSOXSA CCSOSatTT
HLHI L RAV f ODLTIWFEF Q C M I VQI Y GYTHLrTCOt. 0 * 0 f 1 HIM Y OltQISTGDM rVDCrOPTIE LT1Q L ION rVOCTDPTIE LTIQ L I OK rVOCTOPTIE LTIQ L 1 0 3 IPDfOOVUD LTKL I OUT DLVTMKXEJI t i u L VTK PCVIYUULH ILHR I LAT
..CL
Fig. 4 . Partial output from a search with a simple sequence pattern.
within the next set of spacing limits. Processors within this region that indicate a match to block / — 2 then increment their score and record the relative coordinates of the previous matching blocks. At each stage only those processors that are legitimate starts for the current aggregate pattern are able to initiate the local searching for matches to the next pattern block. When all the blocks have been linked, the scores and sets of coordinates are returned to the front-end. The relative coordinates of the blocks for an individual match are used to reconstruct the actual spacing when the match is displayed. Execution time for this stage of the algorithm increases linearly with the sum of differences between the upper and lower bounds of the spacings between pattern elements. In the example shown in Figure 3 there is only one processor that can link pattern blocks 2 and 3. However, there are two instances of matches to block 1 that fall within the spacing limits between blocks 1 and 2. Both of these are returned as legitmate matches to the entire pattern. In this way multiple start positions within a sequence can share the same matching blocks with higher coordinates. The current method does not permit a single start position to show ambiguity in blocks with higher coordinates. In other words, each start position returned specifies a unique arrangement of matching blocks. This restriction could be lifted by linking the blocks from both directions, as opposed to the current top-down approach, and then removing duplicate matches. The convenience features that comprise the third level of the syntax described earlier are handled entirely within the parsing routines that reside on the serial front-end computer. The method described forms the core of an interactive search tool, called 'cmsearch'. The program is invoked with the name of a sequence database, which is read from disk and loaded
PMttn:
(C..G.CK.T)
i: utct I- Ontrt l! wncr*
(R
V.Q1
M pwraMM I m i Cl|r«lM
OUCSCItT ie*tcscoT JCLJCICBIT
uj«ucrriDgvLiDC«iAjiiioMurvt
ld«VCWLCi (MVUJMS11 rr«LTMvTfli unnmwnv wuiwivrfli i r t u r a r v
Fig. 5 . (a) A search with a pattern that contains character classes and predefined
names, (b) A search with a pattern that contains a variable
spacer between blocks.
into CM memory. The user can then repeatedly enter patterns at the prompt with the results being dispalyed within a few seconds. The display takes the form of a list of names of matching database entries followed by the matchingregionsof sequence. These regions are displayed in their context, usually with ten characters on either side, which is found to be very useful in the interactive refinement of a pattern. The location of the pattern start within the sequence is displayed as well as the score of the match. With patterns that contain variable spacing between blocks the display is padded to align matches to the individual blocks, using therelativeoffsets between blocks that were recorded during the search. Interpretation of matches to such patterns without this feature is extremely difficult. Results The results described here were obtained with an 8192 processor CM-2 and employed a version of the PIR database that contained 2.1 million amino acids (Protein Identification Resource, Release 25.0, National Biomedical Research Foundation, Washington, DC). In several of the figures the list of matching sequences has been edited for reasons for space.
381
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
7
~] Two domains in the
R.Jones
Pattern: IYGKGGIGKST 1: MICLT* 2: HIKBTP 3t HIAVF
Threshold - 11 (Exact Match)
Search Time versus Length of Pattern (Exact Match)
H l t r o g a n j s a Iron p r o t a l n - c l o s t r l d l u a pastaurljnui m t r o o a n a s a Iron p r o t a l n - K l a b s l a l U poauvonLja l*Uroaana»a Iron p r o t a l n - Axotob«ctar v l n a U n d l l
1.501: 2: 3|
11 11 11
< 7 I
HRQVA :IYGKGCIGXST: TTWLTSGLH TMBQCA :ITCKCCICKSTI TTOILVAALA KAMRQCA :IYCXCCICKIT: TT0NLVHA1A
Pattern: IYGKGGK3KST
1.25-
Threshotd-10, 9, or 8
Additional Matches: KIAir niitni HIXRFT HILVC
«I 5: (l 7t
10 10 10 10
m t r o g a n t s a Iron p r o t a t n - An«b*an« s p . Mltrooan«*a Iron p r o t a l n - Rhlxoblua M l l l o t l Mltrooarusa Iron protaln - fthlzoblua laoiMlnoaaru frxC protaln - Llv«n*ort (H«rch«ntlA polyaorpnt)
11 » • 3
HTSeMROIA HAALAOIA HAAIJtOIA H(IA
Pattern: IYGKGGIGKST
:FICXCCia(JT: iFTCKSCICMTl :FTCXGGICUT: :VTCXCCIGKST:
T3QWTLAAMA TSCWTLAALV T30HTLAALV TSCHllIALA
Threshold-7
Additional Matches : ll «l
10: 11: 12: 13:
I: »: 10: lit 12: 13:
BVAIA1 KIBETU TVDOR1 OOTWT 8VECUA ErBTlA
7 7 7 7
7 7
1• ( 9 29 1 37 (111 ?•
0.50-
0.25-
10
20
30
40
50
Pattern Length
aroM protaln - EMrlcall* nldul«n* Thy*ldlna kln«sa - Fallna harpacvlrus 1 Transforming. proLaln (rji) - Dlctyoatellu
Fig. 7. The effect of the length of a simple pattern on search performance.
uvrA p r o t a l n - E s c h a r l c h U c o l l E l o n o « d o n faccor l - « l p h * A - YaaaL
GPORCKAJIr KSOVHTVSIT KTKYKLV LOARAGDVIS LTLPVCLTTC IGMTTTGHL
(34 a j L C h a i L o t ' l
IICMRCACKST: ACHWVSKAUi 1DCATCICXSL: TAKYLVRAOE lVCCCCVCXSA: LTioLiomr IICStCSCXST: FLRCIHFLCK ITCVSCSdUT: LIHOTLTPIA ITXCSCIDIIRT: IEKFEXEAAE
- output a d l t a d f o r c l i r l t y )
Fig. 6. An example of the effect of pattern threshold scores.
Figure 4 shows an example of the output from the program with a simple pattern as input. The pattern 'F . . . D . D . . G . I . . . EL' is diagnostic of a subset of calcium-binding sites and requires an exact match to all seven specified amino acids. Although too restrictive to match all known calcium-binding proteins, this pattern matched 57 sequences in the database, all of which are known calciumbinding proteins. A number of these are shown in Figure 4 to illustrate the output of the program. Note that in the case of the calmodulin family the pattern has matched to two of the four binding sites present in each sequence. Figure 5(a) illustrates the use of character classes within a pattern, as well as the use of predefined names for both classes and entire pattern blocks. The original pattern consists of a named pattern block, for an ATP-binding site, and five named character classes. During parsing of the pattern these are expanded into an explicit sequence pattern, as shown, and then matched to the database. Matches to the four blocks of the pattern are shown separated by single space columns from the non-matching regions and the flanking sequence. This helps to highlight which positions in a disperse pattern have contributed to a match. In the first matching column of the figure there are instances of all members of the specificed character class
382
0.75-
(aliphatic). An example of specifying spacing limits between pattern blocks is given in Figure 5(b). Matches to blocks of the pattern are shown in aligned columns, but the spacer regions are simply left-justified. The use of threshold scores for matches to a sequence block is illustrated in Figure 6. The original pattern was an 11 amino acid fragment from a single sequence. Three sequences are an exact match to the pattern. Relaxing the stringency of the match, to allow one, two or three mismatches anywhere within the 11 character pattern, results in seven matches. When a threshold of seven is used, allowing up to four mismatches,the number of matches jumps to 34. Among the matching sequences are many examples of known ATP-binding sites from a variety of proteins families. Score thresholds provide a versatile way to find sequences with similar, but not identical, patterns and to help determine which elements of a pattern are the most conserved amongst a diverse set of proteins. The search algorithm is well suited to a massively parallel computer and achieves very high performance for complex patterns. Loading the database from aremotedisk into processor memory takes of the order of 5 s. The elapsed time between entering a pattern at the keyboard and completion of the screen display of the results typically takes between 2 and 10 s, depending on the number of matches that are found. This performance results in a fully interactive search tool, allowing the user to experiment with variations of an original pattern and to see the output immediately. Simple patterns, to which exact matches are required, are searched in time proportional to the pattern length. Figure 7 shows the time taken for the inner loop of the algorithm for a range of pattern lengths. Consensus patterns typically contain
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
4: Si (t 71
Pattern matching on a parallel computer
patterns from sets of related protein sequences. Proc. Nail. Acad. Sci. (USA), 87, 118-122. Smith.T.F. and Waterman,M.S. (1981) Identification of common molecular subsequences. J. Mol. Biol.. 147, 195-197. Taylor.W.R. (1988) Pattern matching methods in protein sequence comparison and structure prediction. Prot. Eng.. 2, 77-86. Taylor.W.R. (1989) A template based method of pattern matching in protein sequences. Prog. Biophys. Mol. Biol., 54, 159-252. Received on October 17, 1991; accepted on January 12, 1992
Circle No. 9 on Reader Enquiry Card Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
of the order of 10 informative characters and the figure shows that the preloaded database can be searched at a speed of - 6 Mbytes/s for patterns of this length, on an 8192 processor CM-2. Performance scales linearly with increasing numbers of processors such that a 64K processor CM-2 can be predicted to perform the same task at around 48 Mbytes/s for 10 character patterns. In order to search very large databases, the titles of matching database entries are retrieved from file, rather than being stored in the memory. Displaying the search results therefore involves disk access, usually across a network, and is often the limiting stage in the process. Given that the overall time for most searches is only a few seconds, this does not seem to be a serious problem. Discussion High-performance database search tools offer one solution to the problem of searching a body of data that is growing rapidly. The combination of the rich pattern syntax described here and the truly interactive response of the search tool facilitates experimentation with sequence patterns. It allows a user to explore the available sequence data in a manner that has not previously been available. Putative diagnostic patterns can be rapidly tested and the effects of variants examined. Ambiguity in a pattern can be introduced at specific positions through the use of character classes or through the use of threshold scores. The method could be extended to include the use of an amino acid similarity matrix in character matching and to permit the use of amino acid frequency matrices as elements of the syntax. These features would increase the ways in which ambiguity could be expressed in patterns. The performance of the current method is independent of whether or not mismtaches are allowed within a pattern block. This could be used to advantage by sequence comparison methods that employ word neighborhoods to express ambiguity in matching, such as the program BLAST (Altschul etai, 1990). The software described there is available from the author or by anonymous ftp from the Internet node 'think.com'. Connection Machine CM-2 is a trademark of Thinking Machines Corporation. References Altschul.S.F., Gish.W., Miller,W., Myers.E.W. and Lipman.DJ. (1990) Basic local alignment search tool. J. Mol. Biol., 215, 403-410. Bairoch.A. (1990) PROSITE: A Dictionary of Protein Sites and Patterns. University of Geneva. Cohen.F.E., Arbancl.R.M., Kuntz.I.D. and Flctterick.R.D. (1983) Secondary structure assignment for alpha/beta proteins by a combinatorial approach. Biochemistry. 22, 4894-4904. Cohen.F.E., Arbanel.R.M., Kuntz.I.D. and Fletterick.R.D. (1986) Turn prediction in proteins using a pattern-matching approach. Biochemistry, 25, 266-275. Smith,R.F. and Smith.T.F. (1990) Automatic generation of primary sequence
383
CABIOS
Sequence pattern matching on a massively parallel computer Robert Jones Abstract
Introduction Consensus patterns serve as compact representations of sequence features and highlight the most important elements. Searching for matches to these patterns within a sequence is an important method in computational molecular biology. A newly sequenced protein of unknown function can be examined for the presence of a variety of ligand-binding sites, with matches serving as important clues to structure and function. Scanning a eukaryotic DNA sequence for matches to consensus donor and acceptor sites for RNA splicing can help uncover intron and exon structure. Comparing a single sequence against a database using dynamic programming is the most sensitive method available for finding weak relationships between sequences but is computationally very expensive (Smith and Waterman, 1981). An alternative is to produce short consensus sequence patterns that are diagnostic for particular protein families and then to use these to identify other members of those families. Several groups have published databases of these consensus patterns (Bairoch, 1990; Smith and Smith, 1990). In many ways sequence pattern matching resembles the problem of finding regular expressions in text, and many of those algorithms have been adapted to the domain of biological sequences. However, there are important features found in sequence patterns that do not occur in regular text. In many positions of a consensus pattern it is necessary to represent conservation of some physical or chemical property, rather than the identity of a specific amino acids. Here one must score a match if any one of a class of characters matches at the given position in the pattern. In addition two regions of a consensus
Thinking Machines Corporation, 245 First Street. Cambridge. MA 02142. USA
© Oxford University Press
Methods In this section a syntax is presented in which complex sequence patterns can be described, followed by the algorithms that enable them to be matched with a sequence database. At present there is no standard sequence pattern syntax in use and no one existing syntax could be found in which to express all the features described here. Accordingly a novel syntax is presented, although this could be readily modified to conform to any emerging standard. The pattern syntax used here can be viewed as having three levels, with the lowest level representing the basic elements. The simplest pattern consists of a contiguous string of character codes for explicit amino acids or nucleotides, 'DEGAG' for example. A match to this pattern is reported when all the individual positions match in the same spatial order to a contiguous region of a database sequence. In the current work, every matching position has a score of one, with mismatched positions scoring zero. Hence an exact match to a pattern of length N has a score of N. Ambiguity can be introduced by defining a threshold score for the pattern and reporting a match only if the accumulated score equals or exceeds the threshold. A threshold score is written after the pattern following a colon. For example the pattern 'DEGAG:4' will match to any region of five amino acids in which four of the five positions match exactly. This ability to express a degree of amgibuity 377
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
A method is described for finding all occurrences of a sequence pattern within a database of molecular sequences. Implementation of this on a massively parallel computer allows the user to perform very fast database searches using complex patterns. In particular, the software supports approximate pattern matching with score thresholds for either the entire pattern or specified elements thereof. Matches to individual elements can be linked by variable length gaps within userspecified limits.
pattern may exhibit variation in their spacing in the examples used to derive the pattern, usually constrained within certain bounds. The syntax of most text-searching methods is unable to represent this. There have been a number of tools developed for searching specifically with sequence patterns. In the area of protein structure prediction very flexible pattern specification languages and search tools have been described (Cohen et al, 1983, 1986; Taylor, 1988, 1989). Several tools have been reported for searching databases with consensus patterns characteristic of protein families, and searching with short sequence fragments forms the basis of a fast database comparison method (Altschul et al, 1990; Smith and Smith, 1990; G.Mehldau and E.Myers, submitted). The current work describes a novel approach to searching sequence databases that uses a massively parallel computer, the Connection Machine CM-2 ®, to achieve very high performance. A rich pattern syntax is provided in which to describe the diversity of observed sequence patterns.
RJones
378
pattern '$ATP-binding' expands to the explicit pattern 'G . . G . GK'. The algorithm described here exploits the architecture of a massively parallel computer, the Connection Machine CM-2, to achieve very high performance,. The CM-2 can be configured with up to 65,536 single-bit processors and up to 128 kbytes of memory per processor. Instructions are broadcast to the processors from a front-end serial computer in a single instruction stream and these are performed in parallel on the data local to each processor. All processors receive all the instructions but only those in which the 'context' flag is set actually operate on their data. Conditional programming constructs are implemented through manipulation of this flag. For many problems the number of items of data greatly exceeds the maximum number of physical processors and this is certainly the case in database searching. The software of the CM-2 can emulate very large machines by having each physical processor simulate a larger number of virtual processors (VPs). This abstraction is made transparent to the programmer and though the execution time increases as the number of virutal processors grows, the user's program remains unchanged. The physical processors are linked by a hypercube communications network. In software, the processors can be configured as Ndimensional arrays, referred to as geometries, and the same physical processors may participate in several different geometries during the course of a program. Communication between processors makes use of the hypercube network in two ways. If the pattern of communication is regular, for instance when all processor / send a value to all processor i — 1, very efficient communications primitives are available. If the pattern is irregular, with the possibility of several processors sending data to a single recipient, then a more general method must be used. The software described here was written in the language C*, a version of the C language that has been extended to exploit massively parallel architectures. In this application the machine is configured as a onedimensional array of virtual processors onto which the sequence database is laid out with one amino acid per processor. Adjacent sequences are delimited by a special character. Consider the task of finding positions in the database that are an exact match to a short simple pattern. Each processor has a score that is initially set to zero. The character at the end of the sequence pattern is broadcast to all processors and those that contain this character increment their score by one. Processors holding the delimiter character reset their score to zero. All processors pass their current score to their neighor below them in order of the array. A fast regular communication primitive can be used as all processors are passing data to the same relative address. The cycle is repeated by broadcasting the second to last character in the pattern to all processors. Those that match the character increment their score, delimiting positions reset their score to zero and all processors pass their score to their lower neighbor.
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
in a pattern match is found to be very useful in the domain of sequence analysis. A position in the pattern that is poorly conserved can be deemphasized by using a period in that position. Any character can match to a period and the match does not contribute to the scoring. Fixed length gaps between pattern elements can be presented by a series of periods, e.g. 'DE . . . GAG'. A character class, which allows a limited group of characters to match at a given pattern position, is represented by enclosing the acceptable characters in square brackets. For example, the pattern [ILV] matches to a single sequence position containing I, L or V. A character class that excludes a set of characters from a position is denoted by preceding the contents of the class with the character '!'. The parser expands this into a character class containing all amino acids not specifically excluded. The second level of syntax allows blocks of these basic pattern elements to be combined in a very flexible manner. This level of the syntax is not explicitly supported by most text-searching algorithms. The use of score thresholds permits ambiguity at the level of the entire pattern. For example, one might wish to combine an exact match to one region of a pattern with an approximate match to a different region. This is achieved by splitting the pattern into a number of 'blocks'. Enclosing part of the pattern in parentheses marks it as a block and it may then be assigned a 'local' score threshold. For example, the pattern '(DEGAG:4) (GT)' requires an exact match to 'GT' and an exact match in at least four positions to the five characters of 'DEGAG'. If the context in which a block is found is unambiguous, then the parsing function is able to recognize the implicit block structure and parentheses are not required. The concept of pattern blocks allows variable spacing between blocks to be defined. This is represented by angle brackets enclosing the inclusive bounds of an acceptable gap. For example the pattern 'DEG AG' will match to a sequence containing the pattern ' D E C separated by between two and five other amino acids from the pattern 'AG'. The third level of pattern syntax provides a number of convenience features. Character classes are used extensively in consensus patterns and in many cases the class of amino acids represents a particular conserved chemical characteristic, e.g. the class '[KR]' can represent conservation of a positive charge and the class '[AGS]' can represent conservation of a small side chain. Rather than having to specify these classes explicitly every time they are employed, the user can define named character classes and store these in a file. Named classes can then be incorporated into patterns by preceding the name with a ' # ' character. For example, the definition '# aliphatic' is expanded to the explicit class '[ILV]'. Entire patterns can be assigned names and then referred to by preceding the name with a '$' character. Predefined patterns are provided for a collection of protein ligand binding sites as well as 'fingerprint' patterns for a number of protein families. For example, the named
Pattern matching on a parallel computer
Pattern : Q S T Processor
10
11
12
13
14
Q
S
15
Sequence Score Step MalchT
2:
ShW score
3:
Match"?
4:
Shin score
5:
Match YJ
Processors with a score equal to the length of the pattern (3) indicate the start of exact matches to the pattern
Fig. 1. Matching the pattern 'QST' to a sequence.
The process continues until the entire pattern has been broadcast. Scores are not passed to neighboring processors in the final iteration. Figure 1 shows the steps involved for a simple example involving 16 processors and a pattern of three characters. For an exact match the threshold score is simply the length of the pattern. All processors with scores less than that are disabled, leaving only those that mark the start of exact matches to the pattern. Accordingly, for simple patterns like this of length n, the algorithm requires n broadcast instructions, n character matches and n — 1 local communications operations. The complexity of the algorithm is O(nP/N), where n is the pattern length, N is the number of physical processors and P is the length of the sequence database. Where an exact match is required, the threshold is simply the number of effective positions in the pattern. Ambiguity is expressed by explicitly specifying a lower threshold score. Processors with scores equal or greater than this threshold mark the start of regions that contain up to the specified number of mismatched positions. These mismatches can occur anywhere within a single pattern block. For example in Figure 1 processor 5 has a final score of 2, indicating that it marks the start of a region (QVT) that contains one mismatch to the supplied pattern (QST). Note that there is no additional cost whatsoever involved in identifying partial matches compared to exact matches. This feature is particularly important for techniques that involve approximate pattern matching, such as the BLAST sequence comparison algorithm (Altschul et ai, 1990). A character class in the pattern is handled by broadcasting and matching all members of the class in turn and only then passing on the current score to the neighboring processors. Accordingly for a class that contains C members there are C broadcast and matching operations but only one communica-
tion step. Positions in the pattern that can match to any character, indicated by a period in the syntax, do not contribute to the threshold score. When such a character is encountered no matching occurs and the current score is passed directly to the neighboring processors. These simple operations implement all the features in the first level of the pattern syntax. The next stage is to return the locations of matching positions to the front-end computer. This requires that the locations and scores of the N matches to the patterns be represented in the first N processors of the array, as opposed to being dispersed throughout the machine. The process is illustrated in Figure 2. The subset of active processors, at which matches start, is enumerated and each processor uses its position within this order as the address of a destination processor, to which it sends both its score and its absolute address in the processor array. For example in Figure 2, processors 1, 9 and 13 all mark the start of a match. Following the enumeration step they carry addresses of 0, 1 and 2 respectively. Scores and absolute coordinates are transferred, using a 'send' operation, to the computed destination processors such that in the figure processor 1 sends its data to processor 0, etc. After this operation the first N processors contain the scores and positions of the N matches to the pattern. These are then read from the machine to the frontend serial computer using an efficient array transfer operation. The titles and sequences of the matching database sequence are then referenced and the search results are displayed for the user. The observed complexity of this stage is approximately O{mlN), where m is the number of matches and N is the number of physical processors. The second level of the pattern syntax provided a way of representing independent pattern blocks that may carry local
379
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
1:
RJones
Processor
0
1
Sequence
D
Disable VPs with scores below the threshold
2".
Enumerate active VPs tor use as a send address
3
Q
0[3i
Rnal Score 1:
2
X
4
5
6
S
T
L
Q
0
0
0
2
X
X
X
7 V
0
X
8 T
9 R
10 Q
11
12
13
14
15
T
K
Q
S
T
0
0
X
X
S
0
0
0
X
X
X
0
X
0
X
3: Compute VP coordinates within the array
1
9
Only the shaded processors are active during these steps
13
The positions of N matches are now In the first N VPs Fig. 2 . Gathering the locations of matches.
Search lor pattern 'QST LIL EDLV
D O S T L Q V T R Q S T K Q S T N L 1: Match block 3 (EDLV)
•
•
2; Match block 2 (LIL)
•
• 1 •
3: Match block 1 (QST)
• «3 '•
IL A C H L
S E D L V F Q S R
1 1
•
1 •
2
4: Unk blocks 2 and 3
5: Link blocks 1 and 2+3
Match 1 Match 2
Q S T
Fig. 3 . Matching pattern blocks linked by variable length gaps.
threshold scores and that may be separated by variable-length spacers. These are first parsed into a series of independent simple patterns, with the relative order and spacing restrictions being recorded. Each basic pattern block is matched to the database independently as described above. Matching positions for each block are flagged and the scores stored in an array on each processor rather than being sent to the front-end. Figure 3 illustrates how these independent matching blocks become linked. First the spacing limits between the end of one block and the start of another are converted into spacing limits between the start positions of adjacent blocks. As with the positions of an individual block, the blocks of a complex pattern are handled
380
in reverse order so that the starts of overall matches are returned at.the end of the procedure. The scores and coordinates of processors that match to the rightmost block (i) are passed in a series of nearest-neighbor communication steps to those processors that fall within the spacing limits for the next lowest pattern block (/ - 1). Any processors within this region that are flagged as matching to block i — 1 are therefore legitmate matches to the pattern segment consisting of blocks i and i — 1. Their score is incremented by that of block /. These processors then continue the procedure by passing their scores and coordinates, as well as the relative coordinates of the associated matches to block j , to those processors that fall
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
4: Send coordinate to send address
Pattern matching on a parallel computer
A Simple Sequence Pattern:
F...D.D..G.X...EL t*cl»
TTHUCS TPTCCI NCUO (CM
CtlBOdull
HCDO HCOO HCUUA HCUUC
CalModuli - Dlctyoic*: IUB ditcoldau* C-lxxh.ll - Olotyoct*. lua discoldaua * P « 2* P ottln - Purpi • M grcblo tpn 2c P ocaln - Pur* • M archIn P«rv«lbuL n beta - Edit 1* froq Parr* 11MM. n beta - b i t l c CKMl P*r**lb-ja o - Tbornback ray C*lclu»-b
CUaoduJi
FVTO
FVCO PVITC RLCII KLM»
sr*
3 •
OMIUTKAA EtHIUTKAA
« 1f * f
uoiurirxA TDtttEIJUA UQIWKU OOSUCLItA
• 1 •
noiuriXA TDTTEE1UXA a WAIT KM
• ctowvrus 7 CTSAADAJOCV 1 ATSADELOL 7 MSOAELAII EHCAKCrifM 11 « XHCOKLDOU If 4 P«tt*rn.
) . c . .o.GJciAcorri |ACG*TI . . . . (ILVI
same protein
Thysldln* Adcoylit*
KIUTt RIHUA KltC A TVHU11A TVKVM TVBYSR
c h U co, • ) - Hi* - lacharlchU coll 1 polyprocvln - 1 — 11 • I roraat v l r p r o t a l n t (tr«n>poBMi Tn7)
bratn
Cilclwfa
PaMM
=j
UttMngSw^i mo* rmrDAOcccoisvuL nwimKicDirrm FSLTDXOCOCTITTIIL ntVTD KDCHCT I >AAI L r ALTD XDCOCTI r m L n v n ROCHC L i I A U L r s LTD KDCOCS i r m L nCVTD KDCMCTICAAIL FKCCCTOCOCK1TS[[L FUIDAOCOGXITPEE L niLDRDUCriEOOIL n i AD ED n c r I u c c L FWVU>GDQ*GY I E V t l L nKYDOOCHCYlDCJICL Ft LYDQOGOCYI DDCt L F. . . D . D . . C . I
CTVTWOCQT CTVHHKLGQT CTVHULCOtl ~ MVHTWLCM CTVmSLCQN MVHTWLCII _ OTVMMLCON ~ MVHTSLGf K _ RAATUIEID KAAntflEID CLTLQWRAJ lUlIATAAD KKTLKCrSDC QAUJCDLCtK DALLXDLCIK
Two domains In the same protein
• 1
21 1
I » t 9
$ 1 SO It
» n
PVHrpACILT cjcuumciir mil KTTttLV VtAKTtTKLV rnfLRTTurv LSIKQGTVIO lAPYlTTWa UMTAOSLL
E V L V V V V r I
OAFGVCVTT GCPCKICT aAKACJCOT CACCVCJtW OAMtVCUA GOGCVGMA GASaSGKST CVPGSOXSA CCSOSatTT
HLHI L RAV f ODLTIWFEF Q C M I VQI Y GYTHLrTCOt. 0 * 0 f 1 HIM Y OltQISTGDM rVDCrOPTIE LT1Q L ION rVOCTDPTIE LTIQ L I OK rVOCTOPTIE LTIQ L 1 0 3 IPDfOOVUD LTKL I OUT DLVTMKXEJI t i u L VTK PCVIYUULH ILHR I LAT
..CL
Fig. 4 . Partial output from a search with a simple sequence pattern.
within the next set of spacing limits. Processors within this region that indicate a match to block / — 2 then increment their score and record the relative coordinates of the previous matching blocks. At each stage only those processors that are legitimate starts for the current aggregate pattern are able to initiate the local searching for matches to the next pattern block. When all the blocks have been linked, the scores and sets of coordinates are returned to the front-end. The relative coordinates of the blocks for an individual match are used to reconstruct the actual spacing when the match is displayed. Execution time for this stage of the algorithm increases linearly with the sum of differences between the upper and lower bounds of the spacings between pattern elements. In the example shown in Figure 3 there is only one processor that can link pattern blocks 2 and 3. However, there are two instances of matches to block 1 that fall within the spacing limits between blocks 1 and 2. Both of these are returned as legitmate matches to the entire pattern. In this way multiple start positions within a sequence can share the same matching blocks with higher coordinates. The current method does not permit a single start position to show ambiguity in blocks with higher coordinates. In other words, each start position returned specifies a unique arrangement of matching blocks. This restriction could be lifted by linking the blocks from both directions, as opposed to the current top-down approach, and then removing duplicate matches. The convenience features that comprise the third level of the syntax described earlier are handled entirely within the parsing routines that reside on the serial front-end computer. The method described forms the core of an interactive search tool, called 'cmsearch'. The program is invoked with the name of a sequence database, which is read from disk and loaded
PMttn:
(C..G.CK.T)
i: utct I- Ontrt l! wncr*
(R
V.Q1
M pwraMM I m i Cl|r«lM
OUCSCItT ie*tcscoT JCLJCICBIT
uj«ucrriDgvLiDC«iAjiiioMurvt
ld«VCWLCi (MVUJMS11 rr«LTMvTfli unnmwnv wuiwivrfli i r t u r a r v
Fig. 5 . (a) A search with a pattern that contains character classes and predefined
names, (b) A search with a pattern that contains a variable
spacer between blocks.
into CM memory. The user can then repeatedly enter patterns at the prompt with the results being dispalyed within a few seconds. The display takes the form of a list of names of matching database entries followed by the matchingregionsof sequence. These regions are displayed in their context, usually with ten characters on either side, which is found to be very useful in the interactive refinement of a pattern. The location of the pattern start within the sequence is displayed as well as the score of the match. With patterns that contain variable spacing between blocks the display is padded to align matches to the individual blocks, using therelativeoffsets between blocks that were recorded during the search. Interpretation of matches to such patterns without this feature is extremely difficult. Results The results described here were obtained with an 8192 processor CM-2 and employed a version of the PIR database that contained 2.1 million amino acids (Protein Identification Resource, Release 25.0, National Biomedical Research Foundation, Washington, DC). In several of the figures the list of matching sequences has been edited for reasons for space.
381
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
7
~] Two domains in the
R.Jones
Pattern: IYGKGGIGKST 1: MICLT* 2: HIKBTP 3t HIAVF
Threshold - 11 (Exact Match)
Search Time versus Length of Pattern (Exact Match)
H l t r o g a n j s a Iron p r o t a l n - c l o s t r l d l u a pastaurljnui m t r o o a n a s a Iron p r o t a l n - K l a b s l a l U poauvonLja l*Uroaana»a Iron p r o t a l n - Axotob«ctar v l n a U n d l l
1.501: 2: 3|
11 11 11
< 7 I
HRQVA :IYGKGCIGXST: TTWLTSGLH TMBQCA :ITCKCCICKSTI TTOILVAALA KAMRQCA :IYCXCCICKIT: TT0NLVHA1A
Pattern: IYGKGGK3KST
1.25-
Threshotd-10, 9, or 8
Additional Matches: KIAir niitni HIXRFT HILVC
«I 5: (l 7t
10 10 10 10
m t r o g a n t s a Iron p r o t a t n - An«b*an« s p . Mltrooan«*a Iron p r o t a l n - Rhlxoblua M l l l o t l Mltrooarusa Iron protaln - fthlzoblua laoiMlnoaaru frxC protaln - Llv«n*ort (H«rch«ntlA polyaorpnt)
11 » • 3
HTSeMROIA HAALAOIA HAAIJtOIA H(IA
Pattern: IYGKGGIGKST
:FICXCCia(JT: iFTCKSCICMTl :FTCXGGICUT: :VTCXCCIGKST:
T3QWTLAAMA TSCWTLAALV T30HTLAALV TSCHllIALA
Threshold-7
Additional Matches : ll «l
10: 11: 12: 13:
I: »: 10: lit 12: 13:
BVAIA1 KIBETU TVDOR1 OOTWT 8VECUA ErBTlA
7 7 7 7
7 7
1• ( 9 29 1 37 (111 ?•
0.50-
0.25-
10
20
30
40
50
Pattern Length
aroM protaln - EMrlcall* nldul«n* Thy*ldlna kln«sa - Fallna harpacvlrus 1 Transforming. proLaln (rji) - Dlctyoatellu
Fig. 7. The effect of the length of a simple pattern on search performance.
uvrA p r o t a l n - E s c h a r l c h U c o l l E l o n o « d o n faccor l - « l p h * A - YaaaL
GPORCKAJIr KSOVHTVSIT KTKYKLV LOARAGDVIS LTLPVCLTTC IGMTTTGHL
(34 a j L C h a i L o t ' l
IICMRCACKST: ACHWVSKAUi 1DCATCICXSL: TAKYLVRAOE lVCCCCVCXSA: LTioLiomr IICStCSCXST: FLRCIHFLCK ITCVSCSdUT: LIHOTLTPIA ITXCSCIDIIRT: IEKFEXEAAE
- output a d l t a d f o r c l i r l t y )
Fig. 6. An example of the effect of pattern threshold scores.
Figure 4 shows an example of the output from the program with a simple pattern as input. The pattern 'F . . . D . D . . G . I . . . EL' is diagnostic of a subset of calcium-binding sites and requires an exact match to all seven specified amino acids. Although too restrictive to match all known calcium-binding proteins, this pattern matched 57 sequences in the database, all of which are known calciumbinding proteins. A number of these are shown in Figure 4 to illustrate the output of the program. Note that in the case of the calmodulin family the pattern has matched to two of the four binding sites present in each sequence. Figure 5(a) illustrates the use of character classes within a pattern, as well as the use of predefined names for both classes and entire pattern blocks. The original pattern consists of a named pattern block, for an ATP-binding site, and five named character classes. During parsing of the pattern these are expanded into an explicit sequence pattern, as shown, and then matched to the database. Matches to the four blocks of the pattern are shown separated by single space columns from the non-matching regions and the flanking sequence. This helps to highlight which positions in a disperse pattern have contributed to a match. In the first matching column of the figure there are instances of all members of the specificed character class
382
0.75-
(aliphatic). An example of specifying spacing limits between pattern blocks is given in Figure 5(b). Matches to blocks of the pattern are shown in aligned columns, but the spacer regions are simply left-justified. The use of threshold scores for matches to a sequence block is illustrated in Figure 6. The original pattern was an 11 amino acid fragment from a single sequence. Three sequences are an exact match to the pattern. Relaxing the stringency of the match, to allow one, two or three mismatches anywhere within the 11 character pattern, results in seven matches. When a threshold of seven is used, allowing up to four mismatches,the number of matches jumps to 34. Among the matching sequences are many examples of known ATP-binding sites from a variety of proteins families. Score thresholds provide a versatile way to find sequences with similar, but not identical, patterns and to help determine which elements of a pattern are the most conserved amongst a diverse set of proteins. The search algorithm is well suited to a massively parallel computer and achieves very high performance for complex patterns. Loading the database from aremotedisk into processor memory takes of the order of 5 s. The elapsed time between entering a pattern at the keyboard and completion of the screen display of the results typically takes between 2 and 10 s, depending on the number of matches that are found. This performance results in a fully interactive search tool, allowing the user to experiment with variations of an original pattern and to see the output immediately. Simple patterns, to which exact matches are required, are searched in time proportional to the pattern length. Figure 7 shows the time taken for the inner loop of the algorithm for a range of pattern lengths. Consensus patterns typically contain
Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
4: Si (t 71
Pattern matching on a parallel computer
patterns from sets of related protein sequences. Proc. Nail. Acad. Sci. (USA), 87, 118-122. Smith.T.F. and Waterman,M.S. (1981) Identification of common molecular subsequences. J. Mol. Biol.. 147, 195-197. Taylor.W.R. (1988) Pattern matching methods in protein sequence comparison and structure prediction. Prot. Eng.. 2, 77-86. Taylor.W.R. (1989) A template based method of pattern matching in protein sequences. Prog. Biophys. Mol. Biol., 54, 159-252. Received on October 17, 1991; accepted on January 12, 1992
Circle No. 9 on Reader Enquiry Card Downloaded from http://bioinformatics.oxfordjournals.org/ at University of Chicago on August 4, 2015
of the order of 10 informative characters and the figure shows that the preloaded database can be searched at a speed of - 6 Mbytes/s for patterns of this length, on an 8192 processor CM-2. Performance scales linearly with increasing numbers of processors such that a 64K processor CM-2 can be predicted to perform the same task at around 48 Mbytes/s for 10 character patterns. In order to search very large databases, the titles of matching database entries are retrieved from file, rather than being stored in the memory. Displaying the search results therefore involves disk access, usually across a network, and is often the limiting stage in the process. Given that the overall time for most searches is only a few seconds, this does not seem to be a serious problem. Discussion High-performance database search tools offer one solution to the problem of searching a body of data that is growing rapidly. The combination of the rich pattern syntax described here and the truly interactive response of the search tool facilitates experimentation with sequence patterns. It allows a user to explore the available sequence data in a manner that has not previously been available. Putative diagnostic patterns can be rapidly tested and the effects of variants examined. Ambiguity in a pattern can be introduced at specific positions through the use of character classes or through the use of threshold scores. The method could be extended to include the use of an amino acid similarity matrix in character matching and to permit the use of amino acid frequency matrices as elements of the syntax. These features would increase the ways in which ambiguity could be expressed in patterns. The performance of the current method is independent of whether or not mismtaches are allowed within a pattern block. This could be used to advantage by sequence comparison methods that employ word neighborhoods to express ambiguity in matching, such as the program BLAST (Altschul etai, 1990). The software described there is available from the author or by anonymous ftp from the Internet node 'think.com'. Connection Machine CM-2 is a trademark of Thinking Machines Corporation. References Altschul.S.F., Gish.W., Miller,W., Myers.E.W. and Lipman.DJ. (1990) Basic local alignment search tool. J. Mol. Biol., 215, 403-410. Bairoch.A. (1990) PROSITE: A Dictionary of Protein Sites and Patterns. University of Geneva. Cohen.F.E., Arbancl.R.M., Kuntz.I.D. and Flctterick.R.D. (1983) Secondary structure assignment for alpha/beta proteins by a combinatorial approach. Biochemistry. 22, 4894-4904. Cohen.F.E., Arbanel.R.M., Kuntz.I.D. and Fletterick.R.D. (1986) Turn prediction in proteins using a pattern-matching approach. Biochemistry, 25, 266-275. Smith,R.F. and Smith.T.F. (1990) Automatic generation of primary sequence
383
