← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Bio/Align/Utilities.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/Align/Utilities.pm
1
 
# $Id: Utilities.pm,v 1.9 2003/06/04 08:36:36 heikki Exp $
 
1
# $Id: Utilities.pm,v 1.20.4.3 2006/10/02 23:10:12 sendu Exp $
2
2
#
3
3
# BioPerl module for Bio::Align::Utilities
4
4
#
17
17
 
18
18
=head1 SYNOPSIS
19
19
 
20
 
use Bio::Align::Utilities qw(aa_to_dna_aln);
21
 
 
22
 
my $dna_aln = aa_to_dna_aln($aaaln,\%dnaseqs);
 
20
  use Bio::Align::Utilities qw(:all);
 
21
  # %dnaseqs is a hash of CDS sequences (spliced)
 
22
 
 
23
 
 
24
  # Even if the protein alignments are local make sure the start/end
 
25
  # stored in the LocatableSeq objects are to the full length protein.
 
26
  # The CoDing Sequence that is passed in should still be the full 
 
27
  # length CDS as the nt alignment will be generated.
 
28
  #
 
29
  my $dna_aln = &aa_to_dna_aln($aa_aln,\%dnaseqs);
 
30
 
 
31
 
 
32
  # generate bootstraps
 
33
  my $replicates = &bootstrap_replicates($aln,$count);
23
34
 
24
35
 
25
36
=head1 DESCRIPTION
42
53
Bioperl modules. Send your comments and suggestions preferably to
43
54
the Bioperl mailing list.  Your participation is much appreciated.
44
55
 
45
 
  bioperl-l@bioperl.org              - General discussion
46
 
  http://bioperl.org/MailList.shtml  - About the mailing lists
 
56
  bioperl-l@bioperl.org                  - General discussion
 
57
  http://bioperl.org/wiki/Mailing_lists  - About the mailing lists
47
58
 
48
59
=head2 Reporting Bugs
49
60
 
50
61
Report bugs to the Bioperl bug tracking system to help us keep track
51
 
of the bugs and their resolution. Bug reports can be submitted via
52
 
email or the web:
 
62
of the bugs and their resolution. Bug reports can be submitted via the
 
63
web:
53
64
 
54
 
  bioperl-bugs@bioperl.org
55
 
  http://bugzilla.bioperl.org/
 
65
  http://bugzilla.open-bio.org/
56
66
 
57
67
=head1 AUTHOR - Jason Stajich
58
68
 
59
69
Email jason@bioperl.org
60
70
 
61
 
=head1 CONTRIBUTORS
62
 
 
63
 
Additional contributors names and emails here
64
 
 
65
71
=head1 APPENDIX
66
72
 
67
73
The rest of the documentation details each of the object methods.
74
80
 
75
81
 
76
82
package Bio::Align::Utilities;
77
 
use vars qw(@ISA @EXPORT @EXPORT_OK);
 
83
use vars qw(@EXPORT @EXPORT_OK $GAP $CODONGAP %EXPORT_TAGS);
78
84
use strict;
79
85
use Carp;
80
86
use Bio::Root::Version;
81
87
require Exporter;
82
88
 
83
 
@ISA = qw(Exporter);
 
89
use base qw(Exporter);
84
90
 
85
91
@EXPORT = qw();
86
 
@EXPORT_OK = qw(aa_to_dna_aln);
87
 
 
88
 
use constant CODONSIZE => 3;
 
92
@EXPORT_OK = qw(aa_to_dna_aln bootstrap_replicates);
 
93
%EXPORT_TAGS = (all =>[@EXPORT, @EXPORT_OK]);
 
94
BEGIN {
 
95
    use constant CODONSIZE => 3;
 
96
    $GAP = '-';
 
97
    $CODONGAP = $GAP x CODONSIZE;
 
98
}
89
99
 
90
100
=head2 aa_to_dna_aln
91
101
 
119
129
        croak('Must provide a valid Bio::Align::AlignI object as the first argument to aa_to_dna_aln, see the documentation for proper usage and the method signature');
120
130
    }
121
131
    my $alnlen = $aln->length;
122
 
    #print "HSP length is $alnlen\n";
123
132
    my $dnaalign = new Bio::SimpleAlign;
 
133
    $aln->map_chars('\.',$GAP);
 
134
 
124
135
    foreach my $seq ( $aln->each_seq ) {    
125
 
        my $newseq;         
126
 
        my $dnaseq = $dnaseqs->{$seq->display_id} || croak("cannot find ".
127
 
                                                         $seq->display_id);
128
 
        foreach my $pos ( 1..$alnlen ) {
129
 
            my $loc = $seq->location_from_column($pos);
130
 
            my $dna = ''; 
131
 
            if( !defined $loc || $loc->location_type ne 'EXACT' ) {
132
 
                $dna = '---';
 
136
        my $aa_seqstr = $seq->seq();
 
137
        my $id = $seq->display_id;
 
138
        my $dnaseq = $dnaseqs->{$id} || $aln->throw("cannot find ".
 
139
                                                     $seq->display_id);
 
140
        my $start_offset = ($seq->start - 1) * CODONSIZE;
 
141
 
 
142
        $dnaseq = $dnaseq->seq();
 
143
        my $dnalen = $dnaseqs->{$id}->length;
 
144
        my $nt_seqstr;
 
145
        my $j = 0;
 
146
        for( my $i = 0; $i < $alnlen; $i++ ) {
 
147
            my $char = substr($aa_seqstr,$i + $start_offset,1);     
 
148
            if ( $char eq $GAP || $j >= $dnalen )  { 
 
149
                $nt_seqstr .= $CODONGAP;
133
150
            } else {
134
 
                # To readjust to codon boundaries
135
 
                # end needs to be +1 so we can just multiply by CODONSIZE 
136
 
                # to get this               
137
 
 
138
 
                my ($start,$end) = ((($loc->start - 1)* CODONSIZE) +1,
139
 
                                    ($loc->end)* CODONSIZE);
140
 
                
141
 
                if( $start <=0 || $end > $dnaseq->length() ) {
142
 
                    print STDERR "start is ", $loc->start, " end is ", $loc->end, " while dnaseq length is ", $dnaseq->length(), " and start/end projected are $start,$end \n";
143
 
                    warn("codons don't seem to be matching up for $start,$end");
144
 
                    $dna = '---';                           
145
 
                } else {
146
 
                    $dna = $dnaseq->subseq($start,$end);
147
 
                }
 
151
                $nt_seqstr .= substr($dnaseq,$j,CODONSIZE);
 
152
                $j += CODONSIZE;
148
153
            }
149
 
            $newseq .= $dna;
150
154
        }
151
 
        # funky looking math is to readjust to codon boundaries and deal
152
 
        # with fact that sequence start with 1
153
 
        my $newdna = new Bio::LocatableSeq(-display_id  => $seq->id(),
154
 
                                           -start => (($seq->start - 1) * 
155
 
                                                      CODONSIZE) + 1, 
156
 
                                           -end   => ($seq->end * CODONSIZE),
157
 
                                           -strand => $seq->strand,
158
 
                                           -seq   => $newseq);    
 
155
        $nt_seqstr .= $GAP x (($alnlen * 3) - length($nt_seqstr));
 
156
 
 
157
        my $newdna = new Bio::LocatableSeq(-display_id  => $id,
 
158
                                           -alphabet    => 'dna',
 
159
                                           -start       => $start_offset+1,
 
160
                                           -end         => ($seq->end * 
 
161
                                                            CODONSIZE),
 
162
                                           -strand      => 1,
 
163
                                           -seq         => $nt_seqstr);    
159
164
        $dnaalign->add_seq($newdna);
160
165
    }
161
166
    return $dnaalign;
162
167
}
 
168
 
 
169
=head2 bootstrap_replicates
 
170
 
 
171
 Title   : bootstrap_replicates
 
172
 Usage   : my $alns = &bootstrap_replicates($aln,100);
 
173
 Function: Generate a pseudo-replicate of the data by randomly
 
174
           sampling, with replacement, the columns from an alignment for
 
175
           the non-parametric bootstrap.
 
176
 Returns : Arrayref of L<Bio::SimpleAlign> objects
 
177
 Args    : L<Bio::SimpleAlign> object
 
178
           Number of replicates to generate
 
179
 
 
180
=cut
 
181
 
 
182
sub bootstrap_replicates {
 
183
   my ($aln,$count) = @_;
 
184
   $count ||= 1;
 
185
   my $alen = $aln->length;
 
186
   my (@seqs,@nm);
 
187
   $aln->set_displayname_flat(1);
 
188
   for my $s ( $aln->each_seq ) {
 
189
       push @seqs, $s->seq();
 
190
       push @nm, $s->id;
 
191
   }
 
192
   my (@alns,$i);
 
193
   while( $count-- > 0 ) {
 
194
       my @newseqs;
 
195
       for($i =0; $i < $alen; $i++ ) {
 
196
           my $index = int(rand($alen));
 
197
           my $c = 0;
 
198
           for ( @seqs ) {
 
199
               $newseqs[$c++] .= substr($_,$index,1);
 
200
           }
 
201
       }
 
202
       my $newaln = Bio::SimpleAlign->new();
 
203
       my $i = 0;
 
204
       for my $s ( @newseqs ) {
 
205
 
 
206
           $newaln->add_seq( Bio::LocatableSeq->new
 
207
                             (-start         => 1,
 
208
                              -end           => $alen,
 
209
                              -display_id    => $nm[$i++],
 
210
                              -seq           => $s));
 
211
       }
 
212
       push @alns, $newaln;
 
213
   }
 
214
   return \@alns;
 
215
}
 
216
 
163
217
1;