← Back to branch summary

~ubuntu-branches/ubuntu/raring/bioperl/raring

« back to all changes in this revision

Viewing changes to Bio/Restriction/Analysis.pm

  • Committer: Bazaar Package Importer
  • Author(s): Charles Plessy
  • Date: 2008-03-18 14:44:57 UTC
  • mfrom: (4 hardy)
  • mto: This revision was merged to the branch mainline in revision 6.
  • Revision ID: james.westby@ubuntu.com-20080318144457-1jjoztrvqwf0gruk
http://bugs.debian.org/448890
* debian/control:
  - Removed MIA Matt Hope (dopey) from the Uploaders field.
    Thank you for your work, Matt. I hope you are doing well.
  - Downgraded some recommended package to the 'Suggests' priority,
    according to the following discussion on Upstream's mail list.
    http://bioperl.org/pipermail/bioperl-l/2008-March/027379.html
    (Closes: #448890)
* debian/copyright converted to machine-readable format.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Restriction/Analysis.pm
1
 
#------------------------------------------------------------------
 
1
# $Id: Analysis.pm,v 1.19.4.1 2006/10/02 23:10:23 sendu Exp $
2
2
#
3
3
# BioPerl module Bio::Restriction::Analysis
4
4
#
5
5
# Cared for by Rob Edwards <redwards@utmem.edu>
6
6
#
7
7
# You may distribute this module under the same terms as perl itself
8
 
#------------------------------------------------------------------
9
8
 
10
9
## POD Documentation:
11
10
 
22
21
  use Data::Dumper;
23
22
 
24
23
  # get a DNA sequence from somewhere
25
 
  my $seq=new Bio::PrimarySeq
26
 
      (-seq =>'AGCTTAATTCATTAGCTCTGACTGCAACGGGCAATATGTCTC'.
27
 
       'TGTGTGGATTAAAAAAAGAGTGAGCTTCTGATAGCAGC',
 
24
  my $seq = new Bio::PrimarySeq
 
25
      (-seq =>'AGCTTAATTCATTAGCTCTGACTGCAACGGGCAATATGTCTC',
28
26
       -primary_id => 'synopsis',
29
27
       -molecule => 'dna');
30
28
 
31
29
  # now start an analysis.
32
30
  # this is using the default set of enzymes
33
 
  my $ra=Bio::Restriction::Analysis->new(-seq=>$seq);
 
31
  my $ra = Bio::Restriction::Analysis->new(-seq=>$seq);
34
32
 
35
33
  # find unique cutters. This returns a
36
34
  # Bio::Restriction::EnzymeCollection object
37
 
  my $enzymes=$ra->unique_cutters;
 
35
  my $enzymes = $ra->unique_cutters;
38
36
  print "Unique cutters: ", join (', ', 
39
37
      map {$_->name} $enzymes->unique_cutters), "\n";
40
38
 
42
40
  # This is will return an array of the sequence strings
43
41
 
44
42
  my $enz = 'AluI';
45
 
  my @frags=$ra->fragments($enz);
 
43
  my @frags = $ra->fragments($enz);
46
44
  # how big are the fragments?
47
45
  print "AluI fragment lengths: ", join(' & ', map {length $_} @frags), "\n";
48
46
 
53
51
  print "All sizes, sorted ", join (" ", $ra->sizes($enz, 0, 1)), "\n";
54
52
 
55
53
  # how many times does each enzyme cut
56
 
  my $cuts=$ra->cuts_by_enzyme('BamHI');
 
54
  my $cuts = $ra->cuts_by_enzyme('BamHI');
57
55
  print "BamHI cuts $cuts times\n";
58
56
 
59
57
  # How many enzymes do not cut at all?
61
59
        " enzymes that do not cut\n";
62
60
 
63
61
  # what about enzymes that cut twice?
64
 
  my $two_cutters=$ra->cutters(2);
 
62
  my $two_cutters = $ra->cutters(2);
65
63
  print join (" ", map {$_->name} $two_cutters->each_enzyme),
66
64
      " cut the sequence twice\n";
67
65
 
68
66
  # what are all the enzymes that cut, and how often do they cut
69
67
  printf "\n%-10s%s\n", 'Enzyme', 'Number of Cuts';
70
 
  my $all_cutters=$ra->cutters;
 
68
  my $all_cutters = $ra->cutters;
71
69
  map {
72
70
      printf "%-10s%s\n", $_->name, $ra->cuts_by_enzyme($_->name)
73
71
  } $all_cutters->each_enzyme;
75
73
  # Finally, we can interact the restriction enzyme object by
76
74
  # retrieving it from the collection object see the docs for
77
75
  # Bio::Restriction::Enzyme.pm
78
 
  my $enzobj=$enzymes->get_enzyme($enz);
 
76
  my $enzobj = $enzymes->get_enzyme($enz);
79
77
 
80
78
 
81
79
=head1 DESCRIPTION
92
90
like this:
93
91
 
94
92
  use Bio::Restriction::Analysis;
95
 
  my $ra=Bio::Restriction::Analysis->new(-seq=>$seqobj);
 
93
  my $ra = Bio::Restriction::Analysis->new(-seq=>$seqobj);
96
94
 
97
95
or
98
96
 
99
 
  my $ra=Bio::Restriction::Analysis->new
 
97
  my $ra = Bio::Restriction::Analysis->new
100
98
      (-seq=>$seqobj, -enzymes=>$enzs);
101
99
 
102
100
Then, to get the fragments for a particular enzyme use this:
103
101
 
104
 
  @fragments=$ra->fragments('EcoRI');
 
102
  @fragments = $ra->fragments('EcoRI');
105
103
 
106
104
Note that the naming of restriction enzymes is that the last numbers
107
105
are usually Roman numbers (I, II, III, etc). You may want to use
110
108
  # get a reference to an array of unique (single) cutters
111
109
  $singles = $re->unique_cutters;
112
110
  foreach my $enz ($singles->each_enzyme) {
113
 
      @fragments=$re->fragments($enz);
 
111
      @fragments = $re->fragments($enz);
114
112
      ... do something here ...
115
113
  }
116
114
 
126
124
this algorithm. If you have a large sequence (e.g. genome) and 
127
125
want to use ambgiuous enzymes you may want to make seperate
128
126
Bio::Restriction::Enzyme objects for each of the possible
129
 
alternatives and make sure that you don't set is_ambiguous!
 
127
alternatives and make sure that you do not set is_ambiguous!
130
128
 
131
129
This version should correctly deal with overlapping cut sites
132
130
in both ambiguous and non-ambiguous enzymes.
146
144
Bioperl modules. Send your comments and suggestions preferably to one
147
145
of the Bioperl mailing lists. Your participation is much appreciated.
148
146
 
149
 
    bioperl-l@bioperl.org             - General discussion
150
 
    http://bioperl.org/MailList.shtml - About the mailing lists
 
147
  bioperl-l@bioperl.org                  - General discussion
 
148
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
151
149
 
152
150
=head2 Reporting Bugs
153
151
 
154
152
Report bugs to the Bioperl bug tracking system to help us keep track
155
 
the bugs and their resolution. Bug reports can be submitted via email
156
 
or the web:
 
153
the bugs and their resolution. Bug reports can be submitted via the
 
154
web:
157
155
 
158
 
     bioperl-bugs@bio.perl.org
159
 
     http://bugzilla.bioperl.org/
 
156
  http://bugzilla.open-bio.org/
160
157
 
161
158
=head1 AUTHOR
162
159
 
165
162
 
166
163
=head1 CONTRIBUTORS
167
164
 
168
 
Heikki Lehvaslaiho, heikki@ebi.ac.uk
 
165
Heikki Lehvaslaiho, heikki-at-bioperl-dot-org
169
166
 
170
167
=head1 COPYRIGHT
171
168
 
191
188
 
192
189
package Bio::Restriction::Analysis;
193
190
use Bio::Restriction::EnzymeCollection;
194
 
use Bio::Root::Root;
195
191
use strict;
196
192
use Data::Dumper;
197
193
 
198
 
use vars qw (@ISA);
199
 
@ISA = qw(Bio::Root::Root);
 
194
use vars qw ();
 
195
use base qw(Bio::Root::Root);
200
196
 
201
197
=head1 new
202
198
 
481
477
  use Bio::Restriction::EnzymeCollection;
482
478
  use Bio::PrimarySeq;
483
479
 
484
 
  my $seq=new Bio::PrimarySeq
 
480
  my $seq = new Bio::PrimarySeq
485
481
      (-seq =>'AGCTTAATTCATTAGCTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATCCAAAAAAGAGTGAGCTTCTGAT',
486
482
       -primary_id => 'synopsis',
487
483
       -molecule => 'dna');
488
484
 
489
 
  my $ra=Bio::Restriction::Analysis->new(-seq=>$seq);
 
485
  my $ra = Bio::Restriction::Analysis->new(-seq=>$seq);
490
486
 
491
487
  my @gel;
492
488
  my @bam_maps = $ra->fragment_maps('BamHI');
518
514
    } elsif ($enz->isa("Bio::Restriction::EnzymeI")) {
519
515
        @cut_positions=@{$self->{'_cut_positions'}->{$enz->name}};
520
516
    } elsif ($enz->isa("Bio::Restriction::EnzymeCollection")) {
521
 
        $self->cuts('multiple', $enz);
 
517
        $self->cut('multiple', $enz);
522
518
        @cut_positions=@{$self->{'_cut_positions'}->{'multiple_digest'}};
523
519
    }
524
520
 
549
545
        $start=$stop+1;
550
546
    }
551
547
    $stop=$self->{'_seq'}->length;
552
 
    $seq{$start}=$self->{'_seq'}->subseq($start, $stop);
553
 
    $stop{$start}=$stop;
 
548
    if ($start > $stop) {
 
549
        # borderline case. The enzyme cleaved at the end of the sequence
 
550
        # what do I do now?
 
551
    }
 
552
    else {
 
553
         $seq{$start}=$self->{'_seq'}->subseq($start, $stop);
 
554
         $stop{$start}=$stop;
 
555
    }
554
556
 
555
557
    if ($self->{'_seq'}->is_circular) {
556
558
        # join the first and last fragments
701
703
    $end = $self->{'maximum_cuts'} if $end > $self->{'maximum_cuts'};
702
704
    my $set = new Bio::Restriction::EnzymeCollection(-empty => 1);
703
705
 
 
706
    #return an empty set if nothing cuts
 
707
    return $set unless $self->{'maximum_cuts'};
 
708
 
704
709
    for (my $i=$start; $i<=$end; $i++) {
705
710
        $set->enzymes( @{$self->{_number_of_cuts_by_cuts}->{$i}} )
706
711
            if defined $self->{_number_of_cuts_by_cuts}->{$i};
862
867
    my $site=$enz->cut;
863
868
    # split it into the two fragments for the sequence before and after.
864
869
    $site=0 unless defined $site;
865
 
    # shouldn't happen, but what if the cut site is outside of the sequence
866
 
    if ($site < 0 || $site > length($enz->string)) {
867
 
       $self->throw("This is (probably) not your fault.\nGot a cut site of $site and a sequence of ".$enz->string);
868
 
    }
 
870
    # The following should not be an exception, both Type I and Type III
 
871
         # enzymes cut outside of their recognition sequences
 
872
    #if ($site < 0 || $site > length($enz->string)) {
 
873
    #   $self->throw("This is (probably) not your fault.\nGot a cut site of $site and a     # sequence of ".$enz->string);
 
874
    # }
869
875
 
870
876
    # the default values just stop an error from an undefined
871
877
    # string. But they don't affect the split.