5.3 Restriction.pm: Finding Recognition Sites

The time has now come to write the class that creates an object out of a sequence, enzyme name (s), and the map of the location(s) of the enzyme recognition sites in the sequence.

This module depends a great deal on the module Rebase developed in the previous section, but it is fairly short because it just tries to do a small job. This new Restriction class takes a Rebase object (which has the Rebase database translated to regular expressions), some sequence, and a list of enzymes, and uses the regular expressions to find the recognition sites in the sequence. (Note that it doesn't use inheritance; it simply creates a Rebase object to use.)

In this module, the restriction map (the list of locations where the enzymes have recognition sites in the sequence) is obtained through an auxiliary method map_enzyme that simply lists the locations. Clearly, a more graphical display would be easier and more useful. I'll consider that as the book progresses.

5.3.1 The Restriction.pm Module

Here is the Restriction.pm module:

 package Restriction; # # A class to find locations of restriction enzyme recognition sites in #  DNA sequence data. # use strict; use warnings; use Carp; # Class data and methods {    # A list of all attributes with default values.    # "enzyme" is given as an argument possibly multiple time, set as key to _map hash     my %_attributes = (         _rebase      => { },  # A Rebase.pm hash-based object         # key   = restriction enzyme name         # value = space-separated string of recognition sites => regular expressions         _sequence    => '', # DNA sequence data in raw format (only bases)         _map         => { },# a hash: keys are enzyme names,                              # values are arrays of locations          _enzyme      => '', # space- or comma-separated enzyme names,                              # set as key to _map hash     );     # Global variable to keep count of existing objects     my $_count = 0;     # Return a list of all attributes     sub _all_attributes {         keys %_attributes;     }     # Manage the count of existing objects     sub get_count {         $_count;     }     sub _incr_count {         ++$_count;     }     sub _decr_count {         --$_count;     } } 

This opening block shows that Restriction is also a class with a small and simple set of attributes. One is just the DNA sequence data, _sequence .

The attribute _map is a hash that stores the computed restriction map with a key for each restriction enzyme and a value consisting of an array of the positions in the sequence where recognition sites for that enzyme occur.

The attribute _enzyme is one or more enzyme names separated by spaces or commas.

The remaining attribute _rebase is an object of the class I developed in the last section, the class Rebase . Recall that a Rebase object gives you the ability to get regular expressions for recognition sites of restriction enzymes. When you call Restriction->new you must include as one of its arguments a Rebase object, as is shown in the example given in the POD documentation later in this section.

Since a program can have many restriction maps, I also maintain the count of Restriction objects.

5.3.1.1 Initializing Restriction objects

Next , the new constructor method creates and initializes Restriction objects:

 # The "new" constructor method, called from class, e.g. sub new {     my ($class, %args) = @_;     # Create a new object     my $self = bless {  }, $class;     # Set the attributes for the provided arguments     foreach my $attribute ($self->_all_attributes(  )) {         # E.g. attribute = "_name",  argument = "name"         my($argument) = ($attribute =~ /^_(.*)/);         if (exists $args{$argument}) {             if($argument eq 'enzyme') {                 # permit space or comma separated enzyme names                 $args{$argument} =~ s/,/ /g;             }             $self->{$attribute} = $args{$argument};         }     }     # Check that the correct arguments are given     if( not defined $self->{_rebase} ) {         croak "A Rebase object must be given as an argument";     }elsif( ref($self->{_rebase}) ne 'Rebase' ) {         croak "The argument to rebase is not a Rebase object";     }elsif( not defined $self->{_sequence} ) {         croak "A sequence must be given as an argument";     }     # Calculate the locations for each enzyme, store in _map hash attribute     foreach my $enzyme (split(" ", $self->{_enzyme})) {         $self->map_enzyme($enzyme);     }     $self->_incr_count;     return $self; } # For this simple class I have no AUTOLOAD or DESTROY # No get_rebase method, I don't want to pass around a huge hash # No set mutators: all initialization done by way of "new" constructor # No clone method.  Each sequence and set of enzymes can be easily calculated #  by means of a "new" command. 

5.3.1.2 The methods explained

The new constructor method checks that the proper and required arguments are provided: a sequence, enzyme(s), and a Rebase object.

Then the new constructor performs the required computation by calling the method map_enzyme for each enzyme to determine the (possibly empty) array of locations in which the enzyme can cut the sequence. These enzyme-array key/value pairs are stored in the _map attribute.

As you can see in the next segment of code from Restriction.pm , the method map_enzyme works by getting the regular expressions that have been calculated for the enzyme and then finding the positions where the regular expressions match the sequence.

The get_regular_expressions method accesses the Rebase object that was passed to the Restriction->new constructor method and appears as the $self->{_rebase} attribute. In that object, the attribute _rebase is the hash of the database, which looks for the value of the key $enzyme .

The match_positions method puts a match of a regular expression in a while loop, where it loops once for each location it finds. The offset of the matching sequence is available in the special variable $-[0] after a successful regular-expression match (see the perlvar section of the Perl documentation for more details).

 sub map_enzyme {     my($self, $enzyme) = @_;     my(@positions) = (  );     my(@res) = $self->get_regular_expressions($enzyme);     foreach my $re (@res) {         push @positions, $self->match_positions($re);     }     @{$self->{_map}{$enzyme}} = @positions;     return @positions; } sub get_regular_expressions {     my($self, $enzyme) = @_;     my(%sites) = split(' ', $self->{_rebase}{_rebase}{$enzyme});     # May have duplicate values     return values %sites; } # Find positions of a regular expression in the sequence sub match_positions {     my($self, $regexp) = @_;     my @positions = (  );     # Determine positions of regular expression matches     while ( $self->{_sequence} =~ /$regexp/ig ) {         push @positions, ($-[0] + 1 );     }     return(@positions); } 

Finally, here are some helper methods that report on the attributes of sequence and of the calculated map for a given enzyme (or the entire hash that contains the map of all the enzymes). I didn't use AUTOLOAD here because there are only a handful of attributes that return a few different data types: array, scalar string, and hash.

 sub get_enzyme_map {     my($self, $enzyme) = @_;     @{$self->{_map}{$enzyme}}; } sub get_enzyme_names {     my($self) = @_;     keys %{$self->{_map}}; } sub get_sequence {     my($self) = @_;     $self->{_sequence}; } sub get_map {     my($self) = @_;     %{$self->{_map}}; } 

5.3.1.3 Documentation

Here's the (bare bones) documentation for the class embedded right in the module using the POD documentation language:

 =head1 Restriction Restriction: Given a Rebase object, sequence, and list of restriction enzyme     names, return the locations of the recognition sites in the sequence =head1 Synopsis     use Restriction;     use Rebase;     use strict;     use warnings;     my $rebase = Rebase->new(         dbmfile => 'BIONET',         bionetfile => 'bionet.212'     );     my $restrict = Restriction->new(         rebase => $rebase,         enzyme => 'EcoRI, HindIII',         sequence => 'ACGAATTCCGGAATTCG',     );         print "Locations for EcoRI are ", join(' ',                                      $restrict->get_enzyme_map('EcoRI')), "\n"; =head1 AUTHOR James Tisdall =head1 COPYRIGHT Copyright (c) 2003, James Tisdall =cut 1; 

Saving the Synopsis example in a file and running it (making sure the required bionet.212 file is in place) gives the output:

 EcoRI data in Rebase is GAATTC GAATTC Sequence is ACGAATTCCGGAATTCG Locations for EcoRI are 3 11