Provided by: libbio-perl-perl_1.7.2-3_all bug


       Bio::Restriction::Analysis - cutting sequences with restriction enzymes


         # analyze a DNA sequence for restriction enzymes
         use Bio::Restriction::Analysis;
         use Bio::PrimarySeq;
         use Data::Dumper;

         # get a DNA sequence from somewhere
         my $seq = Bio::PrimarySeq->new
              -primary_id => 'synopsis',
              -molecule => 'dna');

         # now start an analysis.
         # this is using the default set of enzymes
         my $ra = Bio::Restriction::Analysis->new(-seq=>$seq);

         # find unique cutters. This returns a
         # Bio::Restriction::EnzymeCollection object
         my $enzymes = $ra->unique_cutters;
         print "Unique cutters: ", join (', ',
             map {$_->name} $enzymes->unique_cutters), "\n";

         # AluI is one them. Where does it cut?
         # This is will return an array of the sequence strings

         my $enz = 'AluI';
         my @frags = $ra->fragments($enz);
         # how big are the fragments?
         print "AluI fragment lengths: ", join(' & ', map {length $_} @frags), "\n";

         # You can also bypass fragments and call sizes directly:
         # to see all the fragment sizes
         print "All sizes: ", join " ", $ra->sizes($enz), "\n";
         # to see all the fragment sizes sorted by size like on a gel
         print "All sizes, sorted ", join (" ", $ra->sizes($enz, 0, 1)), "\n";

         # how many times does each enzyme cut
         my $cuts = $ra->cuts_by_enzyme('BamHI');
         print "BamHI cuts $cuts times\n";

         # How many enzymes do not cut at all?
         print "There are ", scalar $ra->zero_cutters->each_enzyme,
               " enzymes that do not cut\n";

         # what about enzymes that cut twice?
         my $two_cutters = $ra->cutters(2);
         print join (" ", map {$_->name} $two_cutters->each_enzyme),
             " cut the sequence twice\n";

         # what are all the enzymes that cut, and how often do they cut
         printf "\n%-10s%s\n", 'Enzyme', 'Number of Cuts';
         my $all_cutters = $ra->cutters;
         map {
             printf "%-10s%s\n", $_->name, $ra->cuts_by_enzyme($_->name)
         } $all_cutters->each_enzyme;

         # Finally, we can interact the restriction enzyme object by
         # retrieving it from the collection object see the docs for
         my $enzobj = $enzymes->get_enzyme($enz);


       Bio::Restriction::Analysis describes the results of cutting a DNA sequence with
       restriction enzymes.

       To use this module you can pass a sequence object and optionally a
       Bio::Restriction::EnzymeCollection that contains the enzyme(s) to cut the sequences with.
       There is a default set of enzymes that will be loaded if you do not pass in a

       To cut a sequence, set up a Restriction::Analysis object with a sequence like this:

         use Bio::Restriction::Analysis;
         my $ra = Bio::Restriction::Analysis->new(-seq=>$seqobj);


         my $ra = Bio::Restriction::Analysis->new
             (-seq=>$seqobj, -enzymes=>$enzs);

       Then, to get the fragments for a particular enzyme use this:

         @fragments = $ra->fragments('EcoRI');

       Note that the naming of restriction enzymes is that the last numbers are usually Roman
       numbers (I, II, III, etc). You may want to use something like this:

         # get a reference to an array of unique (single) cutters
         $singles = $re->unique_cutters;
         foreach my $enz ($singles->each_enzyme) {
             @fragments = $re->fragments($enz);
             ... do something here ...

       Note that if your sequence is circular, the first and last fragment will be joined so that
       they are the appropriate length and sequence for further analysis. This fragment will also
       be checked for cuts by the enzyme(s).  However, this will change the start of the

       There are two separate algorithms used depending on whether your enzyme has ambiguity. The
       non-ambiguous algorithm is a lot faster, and if you are using very large sequences you
       should try and use this algorithm. If you have a large sequence (e.g. genome) and want to
       use ambgiuous enzymes you may want to make separate Bio::Restriction::Enzyme objects for
       each of the possible alternatives and make sure that you do not set is_ambiguous!

       This version should correctly deal with overlapping cut sites in both ambiguous and non-
       ambiguous enzymes.

       I have tried to write this module with speed and memory in mind so that it can be
       effectively used for large (e.g. genome sized) sequence. This module only stores the cut
       positions internally, and calculates everything else on an as-needed basis. Therefore when
       you call fragment_maps (for example), there may be another delay while these are


   Mailing Lists
       User feedback is an integral part of the evolution of this and other Bioperl modules. Send
       your comments and suggestions preferably to one of the Bioperl mailing lists. Your
       participation is much appreciated.
                  - General discussion  - About the mailing lists

       Please direct usage questions or support issues to the mailing list:

       rather than to the module maintainer directly. Many experienced and reponsive experts will
       be able look at the problem and quickly address it. Please include a thorough description
       of the problem with code and data examples if at all possible.

   Reporting Bugs
       Report bugs to the Bioperl bug tracking system to help us keep track the bugs and their
       resolution. Bug reports can be submitted via the web:


       Rob Edwards,, Steve Chervitz,


       Heikki Lehvaslaiho, heikki-at-bioperl-dot-org Mark A. Jensen, maj-at-fortinbras-dot-us


       Copyright (c) 2003 Rob Edwards.  Some of this work is Copyright (c) 1997-2002 Steve A.
       Chervitz. All Rights Reserved.

       This module is free software; you can redistribute it and/or modify it under the same
       terms as Perl itself.


       Bio::Restriction::Enzyme, Bio::Restriction::EnzymeCollection


       Methods beginning with a leading underscore are considered private and are intended for
       internal use by this module. They are not considered part of the public interface and are
       described here for documentation purposes only.


        Title     : new
        Function  : Initializes the restriction enzyme object
        Returns   : The Restriction::Analysis object
        Arguments :

                        -enzymes=>Restriction::EnzymeCollection object)
                    -seq requires a Bio::PrimarySeq object
                    -enzymes is optional.
                     If omitted it will use the default set of enzymes

       This is the place to start. Pass in a sequence, and you will be able to get the fragments
       back out.  Several other things are available like the number of zero cutters or single

Methods to set parameters

        Title    : seq
        Usage    : $ranalysis->seq($newval);
        Function : get/set method for the  sequence to be cut
        Example  : $re->seq($seq);
        Returns  : value of seq
        Args     : A Bio::PrimarySeqI dna object (optional)

        Title    : enzymes
        Usage    : $re->enzymes($newval)
        Function : gets/Set the restriction enzyme enzymes
        Example  : $re->enzymes('EcoRI')
        Returns  : reference to the collection
        Args     : an array of Bio::Restriction::EnzymeCollection and/or
                   Bio::Restriction::Enzyme objects

       The default object for this method is Bio::Restriction::EnzymeCollection.  However, you
       can also pass it a list of Bio::Restriction::Enzyme objects - even mixed with Collection
       objects.  They will all be stored into one collection.

Perform the analysis

        Title    : cut
        Usage    : $re->cut()
        Function : Cut the sequence with the enzymes
        Example  : $re->cut(); $re->cut('single'); or $re->cut('multiple', $enzymecollection);
        Returns  : $self
        Args     : 'single' (optional), 'multiple' with enzyme collection.

       An explicit cut method is needed to pass arguments to it.

       There are two varieties of cut. Single is the default, and need not be explicitly called.
       This cuts the sequence with each enzyme separately.

       Multiple cuts a sequence with more than one enzyme. You must pass it a
       Bio::Restriction::EnzymeCollection object of the set of enzymes that you want to use in
       the double digest. The results will be stored as an enzyme named "multiple_digest", so you
       can use all the retrieval methods to get the data.

       If you want to use the default setting there is no need to call cut directly. Every method
       in the class that needs output checks the object's internal status and recalculates the
       cuts if needed.

       Note: cut doesn't now re-initialize everything before figuring out cuts. This is so that
       you can do multiple digests, or add more data or whatever. You'll have to use new to reset

       See also the comments in above about ambiguous and non-ambiguous sequences.

        Title     : multiple_digest
        Function  : perform a multiple digest on a sequence
        Returns   : $self so you can go and get any of the other methods
        Arguments : An enzyme collection

        Multiple digests can use 1 or more enzymes, and the data is stored
        in as if it were an enzyme called multiple_digest. You can then
        retrieve information about multiple digests from any of the other

        You can use this method in place of $re->cut('multiple', $enz_coll);

Query the results of the analysis

         Title    : positions
         Function : Retrieve the positions that an enzyme cuts at
         Returns  : An array of the positions that an enzyme cuts at
                  : or an empty array if the enzyme doesn't cut
         Arguments: An enzyme name to retrieve the positions for
         Comments : The cut occurs after the base specified.

         Title    : fragments
         Function : Retrieve the fragments that we cut
         Returns  : An array of the fragments retrieved.
         Arguments: An enzyme name to retrieve the fragments for

       For example this code will retrieve the fragments for all enzymes that cut your sequence

         my $all_cutters = $analysis->cutters;
         foreach my $enz ($$all_cutters->each_enzyme}) {

         Title     : fragment_maps
         Function  : Retrieves fragment sequences with start and end
                     points. Useful for feature construction.

         Returns   : An array containing a hash reference for each fragment,
                     containing the start point, end point and DNA
                     sequence. The hash keys are 'start', 'end' and
                     'seq'. Returns an empty array if not defined.

         Arguments : An enzyme name, enzyme object,
                     or enzyme collection to retrieve the fragments for.

       If passes an enzyme collection it will return the result of a multiple digest. This : will
       also cause the special enzyme 'multiple_digest' to be created so you can get : other
       information about this multiple digest. (TMTOWTDI).

       There is a minor problem with this and $self->fragments that I haven't got a good answer
       for (at the moment). If the sequence is not cut, do we return undef, or the whole

       For linear fragments it would be good to return the whole sequence. For circular fragments
       I am not sure.

       At the moment it returns the whole sequence with start of 1 and end of length of the
       sequence.  For example:

         use Bio::Restriction::Analysis;
         use Bio::Restriction::EnzymeCollection;
         use Bio::PrimarySeq;

         my $seq = Bio::PrimarySeq->new
              -primary_id => 'synopsis',
              -molecule => 'dna');

         my $ra = Bio::Restriction::Analysis->new(-seq=>$seq);

         my @gel;
         my @bam_maps = $ra->fragment_maps('BamHI');
         foreach my $i (@bam_maps) {
            my $start = $i->{start};
            my $end = $i->{end};
            my $sequence = $i->{seq};
            push @gel, "$start--$sequence--$end";
            @gel = sort {length $b <=> length $a} @gel;
         print join("\n", @gel) . "\n";

         Title    : sizes
         Function : Retrieves an array with the sizes of the fragments
         Returns  : Array that has the sizes of the fragments ordered from
                    largest to smallest like they would appear in a gel.
         Arguments: An enzyme name to retrieve the sizes for is required and
                    kilobases to the nearest 0.1 kb, else it will be in
                    bp. If the optional third entry is set the results will
                    be sorted.

       This is designed to make it easy to see what fragments you should get on a gel!

       You should be able to do these:

         # to see all the fragment sizes,
         print join "\n", $re->sizes($enz), "\n";
         # to see all the fragment sizes sorted
         print join "\n", $re->sizes($enz, 0, 1), "\n";
         # to see all the fragment sizes in kb sorted
         print join "\n", $re->sizes($enz, 1, 1), "\n";

How many times does enzymes X cut?

        Title     : cuts_by_enzyme
        Function  : Return the number of cuts for an enzyme
        Returns   : An integer with the number of times each enzyme cuts.
                    Returns 0 if doesn't cut or undef if not defined
        Arguments : An enzyme name string

Which enzymes cut the sequence N times?

        Title     : cutters
        Function  : Find enzymes that cut a given number of times
        Returns   : a Bio::Restriction::EnzymeCollection
        Arguments : 1. exact time or lower limit,
                       non-negative integer, optional
                    2. upper limit, non-negative integer,
                       larger or equalthan first, optional

       If no arguments are given, the method returns all enzymes that do cut the sequence. The
       argument zero, '0', is same as method zero_cutters().  The argument one, '1', corresponds
       to unique_cutters.  If either of the limits is larger than number of cuts any enzyme cuts
       the sequence, the that limit is automagically lowered. The method max_cuts() gives the
       largest number of cuts.

       See Also : unique_cutters, zero_cutters, max_cuts

        Title     : unique_cutters
        Function  : A special case if cutters() where enzymes only cut once
        Returns   : a Bio::Restriction::EnzymeCollection
        Arguments : -

       See also:  cutters, zero_cutters

        Title     : zero_cutters
        Function  : A special case if cutters() where enzymes don't cut the sequence
        Returns   : a Bio::Restriction::EnzymeCollection
        Arguments : -

       See also:  cutters, unique_cutters

        Title     : max_cuts
        Function  : Find the most number of cuts
        Returns   : The number of times the enzyme that cuts most cuts.
        Arguments : None

       This is not a very practical method, but if you are curious...

Internal methods

        Title     : _cuts
        Function  : Figures out which enzymes we know about and cuts the sequence.
        Returns   : Nothing.
        Arguments : None.
        Comments  : An internal method. This will figure out where the sequence
                    should be cut, and provide the appropriate results.

        Title     : _enzyme_sites
        Function  : An internal method to figure out the two sides of an enzyme
        Returns   : The sequence before the cut and the sequence after the cut
        Arguments : A Bio::Restriction::Enzyme object,
                    $comp : boolean, calculate based on $enz->complementary_cut()
                            if true, $enz->cut() if false
        Status    : NOW DEPRECATED - maj

         Title    : _non_pal_enz
         Function : Analyses non_palindromic enzymes for cuts in both ways
                    (in fact, delivers only minus strand cut positions in the
                     plus strand coordinates/maj)
         Returns  : A reference to an array of cut positions
         Arguments: The sequence to check and the enzyme object
         NOW DEPRECATED/maj

        Title     : _ambig_cuts
        Function  : An internal method to localize the cuts in the sequence
        Returns   : A reference to an array of cut positions
        Arguments : The separated enzyme site, the target sequence, and the enzyme object
        Comments  : This is a slow implementation but works for ambiguous sequences.
                    Whenever possible, _nonambig_cuts should be used as it is a lot faster.

        Title     : _nonambig_cuts
        Function  : Figures out which enzymes we know about and cuts the sequence.
        Returns   : Nothing.
        Arguments : The separated enzyme site, the target sequence, and the enzyme object

       An internal method. This will figure out where the sequence should be cut, and provide the
       appropriate results.  This is a much faster implementation because it doesn't use a
       regexp, but it can not deal with ambiguous sequences

        Title   : _make_cuts
        Usage   : $an->_make_cuts( $target_sequence, $enzyme, $complement_q )
        Function: Returns an array of cut sites on target seq, using enzyme
                  on the plus strand ($complement_q = 0) or minus strand
                  ($complement_q = 1); follows Enzyme objects in
        Returns : array of scalar integers
        Args    : sequence string, B:R:Enzyme object, boolean

        Title     : _multiple_cuts
        Function  : Figures out multiple digests
        Returns   : An array of the cut sites for multiply digested DNA
        Arguments : A Bio::Restriction::EnzymeCollection object
        Comments  : Double digests is one subset of this, but you can use
                    as many enzymes as you want.

        Title     : _circular
        Function  : Identifies cuts at the join of the end of the target with
                    the beginning of the target
        Returns   : array of scalar integers ( cut sites near join, if any )
        Arguments : scalar string (target sequence), Bio::Restriction::Enzyme obj

        Title     : _expanded_string
        Function  : Expand nucleotide ambiguity codes to their representative letters
        Returns   : The full length string
        Arguments : The string to be expanded.

       Stolen from the original