In a bioinformatics programming course I am currently teaching (using Python, and approached from a mostly CS rather than a biology background) I assigned a project which involves finding sites where restriction enzymes bind to a DNA sequence. The project is heavily influenced by Chapter 9 of the now classic text "Beginning Perl for Bioinformatics" by James Tisdall. I am using the REBASE database, in particular REBASE Format #9. In this format, each restriction enzyme has a recognition sequence which looks like e.g. C^YCGRG
where ^
gives information about the exact cut-point and upper-case letters such as Y
and R
other than ACGT
are ambiguity codes. The point of the project is to turn these into regular expressions where cleavage information is discarded and the ambiguity codes are expanded into character classes, and then use these regular expressions to search for matches.
All fairly straightforward, but one very observant student noticed that some enzymes had multiple sequences, for example M.BssHIII
has a total of 5 sequences: ACGCGT,CCGCGG,RGCGCY,RCCGGY,GCGCGC
. It is easy to see that these correspond to a total of 11 sequences (after the Rs and Ys are expanded), far fewer than the 144 sequences which would match VSSSSB
(which is the least ambiguous single sequence capable of matching the above).
My question is, what is going on here? Am I safe in simply replacing ,
by |
(for "or") in the regular expression, or is there something more subtle going on? Only 8 of the nearly 5000 enzymes have such multiple sequences, and neither the Tisdall book nor even the description of the format (both of which seem to assume single sequences) shed any light on the situation.