← Back to branch summary

~ubuntu-branches/ubuntu/trusty/bioperl/trusty

« back to all changes in this revision

Viewing changes to Bio/Assembly/IO/bowtie.pm

  • Committer: Package Import Robot
  • Author(s): Charles Plessy
  • Date: 2013-09-22 13:39:48 UTC
  • mfrom: (3.1.11 sid)
  • Revision ID: package-import@ubuntu.com-20130922133948-c6z62zegjyp7ztou
Tags: 1.6.922-1
http://bugs.debian.org/722910
* New upstream release.
* Replaces and Breaks grinder (<< 0.5.3-3~) because of overlaping contents.
  Closes: #722910
* Stop Replacing and Breaking bioperl ( << 1.6.9 ): not needed anymore. 

Show diffs side-by-side

added added

removed removed

Lines of Context:
Bio/Assembly/IO/bowtie.pm
123
123
# check requirements
124
124
    eval "require Bio::Tools::Run::Bowtie; \$HAVE_BOWTIE = 1";
125
125
    unless ( eval "require IO::Uncompress::Gunzip; \$HAVE_IO_UNCOMPRESS = 1") {
126
 
        Bio::Root::Root->warn("IO::Uncompress::Gunzip is not available; you'll have to do your decompression by hand.");
 
126
        Bio::Root::Root->warn("IO::Uncompress::Gunzip is not available; you'll have to do your decompression by hand.");
127
127
    }
128
128
}
129
129
 
145
145
=cut
146
146
 
147
147
sub new {
148
 
        my $class = shift;
149
 
        my @args = @_;
150
 
        my $self = $class->SUPER::new(@args);
151
 
        $self->_initialize(@args);
152
 
        $self->{'_tempdir'} = $self->tempdir(CLEANUP=>1);
153
 
        my ($file, $index, $no_head, $no_sq) = $self->_rearrange([qw(FILE INDEX NO_HEAD NO_SQ)], @args);
154
 
        $file =~ s/^<//;
155
 
        $self->{'_no_head'} = $no_head;
156
 
        $self->{'_no_sq'} = $no_sq;
157
 
 
158
 
        # get the sequence so Bio::DB::Sam can work with it
159
 
        my $refdb;
160
 
        my $inspector;
161
 
        if (-e $index && -r _ ) {
162
 
                $refdb = ($index =~ m/\.gz[^.]*$/) ? $self->_uncompress($index) : $index;
163
 
                my $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$refdb);
164
 
                $self->throw("'$index' is not a fasta file.")
165
 
                        unless $guesser->guess =~ m/^fasta$/;
166
 
        } elsif ($HAVE_BOWTIE) {
167
 
                $inspector = Bio::Tools::Run::Bowtie->new( -command => 'inspect' );
168
 
                $refdb = $inspector->run($index);
169
 
        } else {
170
 
                $self->throw("Bio::Tools::Run::Bowtie is not available - cannot extract refdb from index.");
171
 
        }
172
 
 
173
 
        my $bam_file = $self->_make_bam($self->_bowtie_to_sam($file, $refdb));
174
 
        my $sam = Bio::Assembly::IO->new( -file => "<$bam_file", -refdb => $refdb , -format => 'sam' );
175
 
 
176
 
        return $sam;
 
148
    my $class = shift;
 
149
    my @args = @_;
 
150
    my $self = $class->SUPER::new(@args);
 
151
    $self->_initialize(@args);
 
152
    $self->{'_tempdir'} = $self->tempdir(CLEANUP=>1);
 
153
    my ($file, $index, $no_head, $no_sq) = $self->_rearrange([qw(FILE INDEX NO_HEAD NO_SQ)], @args);
 
154
    $file =~ s/^<//;
 
155
    $self->{'_no_head'} = $no_head;
 
156
    $self->{'_no_sq'} = $no_sq;
 
157
 
 
158
    # get the sequence so Bio::DB::Sam can work with it
 
159
    my $refdb;
 
160
    my $inspector;
 
161
    if (-e $index && -r _ ) {
 
162
        $refdb = ($index =~ m/\.gz[^.]*$/) ? $self->_uncompress($index) : $index;
 
163
        my $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$refdb);
 
164
        $self->throw("'$index' is not a fasta file.")
 
165
            unless $guesser->guess =~ m/^fasta$/;
 
166
    } elsif ($HAVE_BOWTIE) {
 
167
        $inspector = Bio::Tools::Run::Bowtie->new( -command => 'inspect' );
 
168
        $refdb = $inspector->run($index);
 
169
    } else {
 
170
        $self->throw("Bio::Tools::Run::Bowtie is not available - cannot extract refdb from index.");
 
171
    }
 
172
 
 
173
    my $bam_file = $self->_make_bam($self->_bowtie_to_sam($file, $refdb));
 
174
    my $sam = Bio::Assembly::IO->new( -file => "<$bam_file", -refdb => $refdb , -format => 'sam' );
 
175
 
 
176
    return $sam;
177
177
}
178
178
 
179
179
sub _bowtie_to_sam {
180
 
        my ($self, $file, $refdb) = @_;
181
 
 
182
 
        $self->throw("'$file' does not exist or is not readable.")
183
 
                unless ( -e $file && -r _ );
184
 
 
185
 
        if ($file =~ m/\.gz[^.]*$/) {
186
 
                $file = $self->_uncompress($file);
187
 
                $self->close;
188
 
                open (my $fh,$file);
189
 
                $self->file($file);
190
 
                $self->_fh($fh);
191
 
        }
192
 
 
193
 
        my $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$file);
194
 
        $self->throw("'$file' is not a bowtie formatted file.") unless $guesser->guess =~ m/^bowtie$/;
195
 
 
196
 
        my %SQ;
197
 
        my $mapq = 255;
198
 
        my $in_pair;
199
 
        my @mate_line;
200
 
        my $mlen;
201
 
 
202
 
        # create temp file for working
203
 
        my ($sam_tmp_h, $sam_tmp_f) = $self->tempfile( -dir => $self->{'_tempdir'}, -suffix => '.sam' );
204
 
        
205
 
        while (my $line = $self->_readline) {
206
 
                chomp($line);
207
 
                my ($qname,$strand,$rname,$pos,$seq,$qual,$m,$details)=split("\cI",$line);
208
 
                $SQ{$rname} = 1;
209
 
                
210
 
                my $paired_f =  ($qname =~ m#/[12]#) ? 0x03 : 0;
211
 
                my $strand_f = ($strand eq '-') ? 0x10 : 0;
212
 
                my $op_strand_f = ($strand eq '+' && $paired_f) ? 0x20 : 0;
213
 
                my $first_f =  ($qname =~ m#/1#) ? 0x40 : 0;
214
 
                my $second_f = ($qname =~ m#/2#) ? 0x80 : 0;
215
 
                my $flag = $paired_f | $strand_f | $op_strand_f | $first_f | $second_f;
216
 
 
217
 
                $pos++;
218
 
                my $len = length $seq;
219
 
                die unless $len == length $qual;
220
 
                my $cigar = $len.'M';
221
 
                my @detail = split(',',$details);
222
 
                my $dist = 'NM:i:'.scalar @detail;
223
 
                
224
 
                my @mismatch;
225
 
                my $last_pos = 0;
226
 
                for (@detail) {
227
 
                        m/(\d+):(\w)>\w/;
228
 
                        my $err = ($1-$last_pos);
229
 
                        $last_pos = $1+1;
230
 
                        push @mismatch,($err,$2);
231
 
                }
232
 
                push @mismatch, $len-$last_pos;
233
 
                @mismatch = reverse @mismatch if $strand eq '-';
234
 
                my $mismatch = join('',('MD:Z:',@mismatch));
235
 
 
236
 
                if ($paired_f) {
237
 
                        my $mrnm = '=';
238
 
                        if ($in_pair) {
239
 
                                my $mpos = $mate_line[3];
240
 
                                $mate_line[7] = $pos;
241
 
                                my $isize = $mpos-$pos-$len;
242
 
                                $mate_line[8] = -$isize;
243
 
                                print $sam_tmp_h join("\t",@mate_line),"\n";
244
 
                                print $sam_tmp_h join("\t",$qname, $flag, $rname, $pos, $mapq, $cigar, $mrnm, $mpos, $isize, $seq, $qual, $mismatch, $dist),"\n";
245
 
                                $in_pair = 0;
246
 
                        } else {
247
 
                                $mlen = $len;
248
 
                                @mate_line = ($qname, $flag, $rname, $pos, $mapq, $cigar, $mrnm, undef, undef, $seq, $qual, $mismatch, $dist);
249
 
                                $in_pair = 1;
250
 
                        }
251
 
                } else {
252
 
                        my $mrnm = '*';
253
 
                        my $mpos = 0;
254
 
                        my $isize = 0;
255
 
                        print $sam_tmp_h join("\t",$qname, $flag, $rname, $pos, $mapq, $cigar, $mrnm, $mpos, $isize, $seq, $qual, $mismatch, $dist),"\n";
256
 
                }
257
 
        }
258
 
 
259
 
        $sam_tmp_h->close;
260
 
        
261
 
        return $sam_tmp_f if $self->{'_no_head'};
262
 
 
263
 
        my ($samh, $samf) = $self->tempfile( -dir => $self->{'_tempdir'}, -suffix => '.sam' );
264
 
 
265
 
        # print header
266
 
        print $samh $HD;
267
 
        
268
 
        # print sequence dictionary
269
 
        unless ($self->{'_no_sq'}) {
270
 
                my $db  = Bio::SeqIO->new( -file => $refdb, -format => 'fasta' );
271
 
                while ( my $seq = $db->next_seq() ) {
272
 
                        $SQ{$seq->id} = $seq->length if $SQ{$seq->id};
273
 
                }
274
 
        
275
 
                map { print $samh join("\t", ('@SQ', "SN:$_", "LN:$SQ{$_}")), "\n" } keys %SQ;
276
 
        }
277
 
        
278
 
        # print program
279
 
        print $samh $PG;
280
 
        
281
 
        # print alignments
282
 
        open($sam_tmp_h, $sam_tmp_f) or
283
 
                $self->throw("Can not open '$sam_tmp_f' for reading: $!");
284
 
 
285
 
        print $samh $_ while (<$sam_tmp_h>);
286
 
        
287
 
        close($sam_tmp_h);
288
 
        $samh->close;
289
 
        
290
 
        return $samf;
 
180
    my ($self, $file, $refdb) = @_;
 
181
 
 
182
    $self->throw("'$file' does not exist or is not readable.")
 
183
        unless ( -e $file && -r _ );
 
184
 
 
185
    if ($file =~ m/\.gz[^.]*$/) {
 
186
        $file = $self->_uncompress($file);
 
187
        $self->close;
 
188
        open (my $fh,$file);
 
189
        $self->file($file);
 
190
        $self->_fh($fh);
 
191
    }
 
192
 
 
193
    my $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$file);
 
194
    $self->throw("'$file' is not a bowtie formatted file.") unless $guesser->guess =~ m/^bowtie$/;
 
195
 
 
196
    my %SQ;
 
197
    my $mapq = 255;
 
198
    my $in_pair;
 
199
    my @mate_line;
 
200
    my $mlen;
 
201
 
 
202
    # create temp file for working
 
203
    my ($sam_tmp_h, $sam_tmp_f) = $self->tempfile( -dir => $self->{'_tempdir'}, -suffix => '.sam' );
 
204
    
 
205
    while (my $line = $self->_readline) {
 
206
        chomp($line);
 
207
        my ($qname,$strand,$rname,$pos,$seq,$qual,$m,$details)=split("\cI",$line);
 
208
        $SQ{$rname} = 1;
 
209
        
 
210
        my $paired_f =  ($qname =~ m#/[12]#) ? 0x03 : 0;
 
211
        my $strand_f = ($strand eq '-') ? 0x10 : 0;
 
212
        my $op_strand_f = ($strand eq '+' && $paired_f) ? 0x20 : 0;
 
213
        my $first_f =  ($qname =~ m#/1#) ? 0x40 : 0;
 
214
        my $second_f = ($qname =~ m#/2#) ? 0x80 : 0;
 
215
        my $flag = $paired_f | $strand_f | $op_strand_f | $first_f | $second_f;
 
216
 
 
217
        $pos++;
 
218
        my $len = length $seq;
 
219
        die unless $len == length $qual;
 
220
        my $cigar = $len.'M';
 
221
        my @detail = split(',',$details);
 
222
        my $dist = 'NM:i:'.scalar @detail;
 
223
 
 
224
        my @mismatch;
 
225
        my $last_pos = 0;
 
226
        for (@detail) {
 
227
            m/(\d+):(\w)>\w/;
 
228
            my $err = ($1-$last_pos);
 
229
            $last_pos = $1+1;
 
230
            push @mismatch,($err,$2);
 
231
        }
 
232
        push @mismatch, $len-$last_pos;
 
233
        @mismatch = reverse @mismatch if $strand eq '-';
 
234
        my $mismatch = join('',('MD:Z:',@mismatch));
 
235
 
 
236
        if ($paired_f) {
 
237
            my $mrnm = '=';
 
238
            if ($in_pair) {
 
239
                my $mpos = $mate_line[3];
 
240
                $mate_line[7] = $pos;
 
241
                my $isize = $mpos-$pos-$len;
 
242
                $mate_line[8] = -$isize;
 
243
                print $sam_tmp_h join("\t",@mate_line),"\n";
 
244
                print $sam_tmp_h join("\t",$qname, $flag, $rname, $pos, $mapq, $cigar, $mrnm, $mpos, $isize, $seq, $qual, $mismatch, $dist),"\n";
 
245
                $in_pair = 0;
 
246
            } else {
 
247
                $mlen = $len;
 
248
                @mate_line = ($qname, $flag, $rname, $pos, $mapq, $cigar, $mrnm, undef, undef, $seq, $qual, $mismatch, $dist);
 
249
                $in_pair = 1;
 
250
            }
 
251
        } else {
 
252
            my $mrnm = '*';
 
253
            my $mpos = 0;
 
254
            my $isize = 0;
 
255
            print $sam_tmp_h join("\t",$qname, $flag, $rname, $pos, $mapq, $cigar, $mrnm, $mpos, $isize, $seq, $qual, $mismatch, $dist),"\n";
 
256
        }
 
257
    }
 
258
 
 
259
    $sam_tmp_h->close;
 
260
 
 
261
    return $sam_tmp_f if $self->{'_no_head'};
 
262
 
 
263
    my ($samh, $samf) = $self->tempfile( -dir => $self->{'_tempdir'}, -suffix => '.sam' );
 
264
 
 
265
    # print header
 
266
    print $samh $HD;
 
267
 
 
268
    # print sequence dictionary
 
269
    unless ($self->{'_no_sq'}) {
 
270
        my $db  = Bio::SeqIO->new( -file => $refdb, -format => 'fasta' );
 
271
        while ( my $seq = $db->next_seq() ) {
 
272
            $SQ{$seq->id} = $seq->length if $SQ{$seq->id};
 
273
        }
 
274
    
 
275
        map { print $samh join("\t", ('@SQ', "SN:$_", "LN:$SQ{$_}")), "\n" } keys %SQ;
 
276
    }
 
277
    
 
278
    # print program
 
279
    print $samh $PG;
 
280
    
 
281
    # print alignments
 
282
    open($sam_tmp_h, $sam_tmp_f) or
 
283
        $self->throw("Can not open '$sam_tmp_f' for reading: $!");
 
284
 
 
285
    print $samh $_ while (<$sam_tmp_h>);
 
286
    
 
287
    close($sam_tmp_h);
 
288
    $samh->close;
 
289
    
 
290
    return $samf;
291
291
}
292
292
 
293
293
sub _uncompress {
294
 
        my ($self, $file) = @_;
295
 
 
296
 
        unless ($HAVE_IO_UNCOMPRESS) {
297
 
                croak( "IO::Uncompress::Gunzip not available, can't expand '$file'" );
298
 
        }
299
 
        my ($tfh, $tf) = $self->tempfile( -dir => $self->{'_tempdir'} );
300
 
        my $z = IO::Uncompress::Gunzip->new($file);
301
 
        while (my $block = $z->getline) { print $tfh $block }
302
 
        close $tfh;
303
 
        $file = $tf;
304
 
 
305
 
        return $file
 
294
    my ($self, $file) = @_;
 
295
 
 
296
    unless ($HAVE_IO_UNCOMPRESS) {
 
297
        croak( "IO::Uncompress::Gunzip not available, can't expand '$file'" );
 
298
    }
 
299
    my ($tfh, $tf) = $self->tempfile( -dir => $self->{'_tempdir'} );
 
300
    my $z = IO::Uncompress::Gunzip->new($file);
 
301
    while (my $block = $z->getline) { print $tfh $block }
 
302
    close $tfh;
 
303
    $file = $tf;
 
304
 
 
305
    return $file
306
306
}
307
307
 
308
308
sub _make_bam {
309
 
        my ($self, $file) = @_;
310
 
        
311
 
        $self->throw("'$file' does not exist or is not readable")
312
 
                unless ( -e $file && -r _ );
313
 
 
314
 
        # make a sorted bam file from a sam file input
315
 
        my ($bamh, $bamf) = $self->tempfile( -dir => $self->{'_tempdir'}, -suffix => '.bam' );
316
 
        my ($srth, $srtf) = $self->tempfile( -dir => $self->{'_tempdir'}, -suffix => '.srt' );
317
 
        $_->close for ($bamh, $srth);
318
 
        
319
 
        my $samt = Bio::Tools::Run::Samtools->new( -command => 'view',
320
 
                                                   -sam_input => 1,
321
 
                                                   -bam_output => 1 );
322
 
 
323
 
        $samt->run( -bam => $file, -out => $bamf );
324
 
 
325
 
        $samt = Bio::Tools::Run::Samtools->new( -command => 'sort' );
326
 
 
327
 
        $samt->run( -bam => $bamf, -pfx => $srtf);
328
 
 
329
 
        return $srtf.'.bam'
 
309
    my ($self, $file) = @_;
 
310
    
 
311
    $self->throw("'$file' does not exist or is not readable")
 
312
        unless ( -e $file && -r _ );
 
313
 
 
314
    # make a sorted bam file from a sam file input
 
315
    my ($bamh, $bamf) = $self->tempfile( -dir => $self->{'_tempdir'}, -suffix => '.bam' );
 
316
    my ($srth, $srtf) = $self->tempfile( -dir => $self->{'_tempdir'}, -suffix => '.srt' );
 
317
    $_->close for ($bamh, $srth);
 
318
    
 
319
    my $samt = Bio::Tools::Run::Samtools->new( -command => 'view',
 
320
                                               -sam_input => 1,
 
321
                                               -bam_output => 1 );
 
322
 
 
323
    $samt->run( -bam => $file, -out => $bamf );
 
324
 
 
325
    $samt = Bio::Tools::Run::Samtools->new( -command => 'sort' );
 
326
 
 
327
    $samt->run( -bam => $bamf, -pfx => $srtf);
 
328
 
 
329
    return $srtf.'.bam'
330
330
}
331
331
 
332
332
1;