| line |
stmt |
bran |
cond |
sub |
pod |
time |
code |
|
1
|
|
|
|
|
|
|
# |
|
2
|
|
|
|
|
|
|
# BioPerl module for Bio::Tools::Run::Bowtie |
|
3
|
|
|
|
|
|
|
# |
|
4
|
|
|
|
|
|
|
# Please direct questions and support issues to |
|
5
|
|
|
|
|
|
|
# |
|
6
|
|
|
|
|
|
|
# Cared for by Dan Kortschak |
|
7
|
|
|
|
|
|
|
# |
|
8
|
|
|
|
|
|
|
# Copyright Dan Kortschak and Mark A. Jensen |
|
9
|
|
|
|
|
|
|
# |
|
10
|
|
|
|
|
|
|
# You may distribute this module under the same terms as perl itself |
|
11
|
|
|
|
|
|
|
|
|
12
|
|
|
|
|
|
|
# POD documentation - main docs before the code |
|
13
|
|
|
|
|
|
|
|
|
14
|
|
|
|
|
|
|
=head1 NAME |
|
15
|
|
|
|
|
|
|
|
|
16
|
|
|
|
|
|
|
Bio::Tools::Run::Bowtie - Run wrapper for the Bowtie short-read assembler *BETA* |
|
17
|
|
|
|
|
|
|
|
|
18
|
|
|
|
|
|
|
=head1 SYNOPSIS |
|
19
|
|
|
|
|
|
|
|
|
20
|
|
|
|
|
|
|
# create an index |
|
21
|
|
|
|
|
|
|
$bowtie_build = Bio::Tools::Run::Bowtie->new(); |
|
22
|
|
|
|
|
|
|
$index = $bowtie_fac->run( 'reference.fasta', 'index_base' ); |
|
23
|
|
|
|
|
|
|
|
|
24
|
|
|
|
|
|
|
# or with named args... |
|
25
|
|
|
|
|
|
|
|
|
26
|
|
|
|
|
|
|
$index = $bowtie_fac->run( -ref => 'reference.fasta', -ind => 'index_base' ); |
|
27
|
|
|
|
|
|
|
|
|
28
|
|
|
|
|
|
|
# get the base name of the last index from an index builder |
|
29
|
|
|
|
|
|
|
$index = $bowtie_fac->result; |
|
30
|
|
|
|
|
|
|
|
|
31
|
|
|
|
|
|
|
# create an assembly |
|
32
|
|
|
|
|
|
|
$bowtie_fac = Bio::Tools::Run::Bowtie->new(); |
|
33
|
|
|
|
|
|
|
$bowtie_fac->want('Bio::Assembly::Scaffold'); |
|
34
|
|
|
|
|
|
|
$bowtie_assy = $bowtie_fac->run( 'reads.fastq', 'index_base' ); |
|
35
|
|
|
|
|
|
|
|
|
36
|
|
|
|
|
|
|
# if IO::Uncompress::Gunzip is available and with named args... |
|
37
|
|
|
|
|
|
|
$bowtie_assy = $bowtie_fac->run( -seq => 'reads.fastq.gz', -ind => 'index_base' ); |
|
38
|
|
|
|
|
|
|
|
|
39
|
|
|
|
|
|
|
# paired-end |
|
40
|
|
|
|
|
|
|
$bowtie_fac = Bio::Tools::Run::Bowtie->new(-command => 'paired', |
|
41
|
|
|
|
|
|
|
-want => 'Bio::Assembly::Scaffold'); |
|
42
|
|
|
|
|
|
|
$bowtie_assy = $bowtie_fac->run( 'reads.fastq', 'index_base', 'paired-reads.fastq' ); |
|
43
|
|
|
|
|
|
|
|
|
44
|
|
|
|
|
|
|
# be more strict |
|
45
|
|
|
|
|
|
|
$bowtie_fac->set_parameters( -max_qual_mismatch => 50 ); |
|
46
|
|
|
|
|
|
|
|
|
47
|
|
|
|
|
|
|
# create a Bio::Assembly::Scaffold object |
|
48
|
|
|
|
|
|
|
$bowtie_assy = $bowtie_fac->run( 'reads.fastq', 'index_base', 'paired-reads.fastq' ); |
|
49
|
|
|
|
|
|
|
|
|
50
|
|
|
|
|
|
|
# print consensus sequences from assembly object |
|
51
|
|
|
|
|
|
|
for $contig ($bowtie_assy->all_contigs) { |
|
52
|
|
|
|
|
|
|
print $contig->get_consensus_sequence->seq,"\n"; |
|
53
|
|
|
|
|
|
|
} |
|
54
|
|
|
|
|
|
|
|
|
55
|
|
|
|
|
|
|
# get the file object of the last assembly |
|
56
|
|
|
|
|
|
|
$io = $bowtie_fac->result( -want => 'Bio::Root::IO' ); |
|
57
|
|
|
|
|
|
|
|
|
58
|
|
|
|
|
|
|
# get a merged SeqFeature::Collection of all hits |
|
59
|
|
|
|
|
|
|
# - currently only available with SAM format |
|
60
|
|
|
|
|
|
|
$io = $bowtie_fac->result( -want => 'Bio::SeqFeature::Collection' ); |
|
61
|
|
|
|
|
|
|
|
|
62
|
|
|
|
|
|
|
#... or the file name directly |
|
63
|
|
|
|
|
|
|
$filename = $bowtie_fac->result( -want => 'raw' ); |
|
64
|
|
|
|
|
|
|
|
|
65
|
|
|
|
|
|
|
=head1 DESCRIPTION |
|
66
|
|
|
|
|
|
|
|
|
67
|
|
|
|
|
|
|
This module provides a wrapper interface for Ben Langmead and Col |
|
68
|
|
|
|
|
|
|
Trapnell's ultrafast memory-efficient short read aligner C |
|
69
|
|
|
|
|
|
|
(see L for manuals and downloads). |
|
70
|
|
|
|
|
|
|
|
|
71
|
|
|
|
|
|
|
|
|
72
|
|
|
|
|
|
|
=head1 OPTIONS |
|
73
|
|
|
|
|
|
|
|
|
74
|
|
|
|
|
|
|
C is complex, with many command-line options. This module attempts to |
|
75
|
|
|
|
|
|
|
provide and options comprehensively. You can browse the choices like so: |
|
76
|
|
|
|
|
|
|
|
|
77
|
|
|
|
|
|
|
$bowtiefac = Bio::Tools::Run::Bowtie->new( -command => 'single' ); |
|
78
|
|
|
|
|
|
|
# all bowtie commands |
|
79
|
|
|
|
|
|
|
@all_commands = $bowtiefac->available_parameters('commands'); |
|
80
|
|
|
|
|
|
|
@all_commands = $bowtiefac->available_commands; # alias |
|
81
|
|
|
|
|
|
|
# just for single |
|
82
|
|
|
|
|
|
|
@assemble_params = $bowtiefac->available_parameters('params'); |
|
83
|
|
|
|
|
|
|
@assemble_switches = $bowtiefac->available_parameters('switches'); |
|
84
|
|
|
|
|
|
|
@assemble_all_options = $bowtiefac->available_parameters(); |
|
85
|
|
|
|
|
|
|
|
|
86
|
|
|
|
|
|
|
Reasonably mnemonic names have been assigned to the single-letter |
|
87
|
|
|
|
|
|
|
command line options. These are the names returned by |
|
88
|
|
|
|
|
|
|
C, and can be used in the factory constructor |
|
89
|
|
|
|
|
|
|
like typical BioPerl named parameters. |
|
90
|
|
|
|
|
|
|
|
|
91
|
|
|
|
|
|
|
As a number of options are mutually exclusive, and the interpretation of |
|
92
|
|
|
|
|
|
|
intent is based on last-pass option reaching bowtie with potentially unpredicted |
|
93
|
|
|
|
|
|
|
results. This module will prevent inconsistent switches and parameters |
|
94
|
|
|
|
|
|
|
from being passed. |
|
95
|
|
|
|
|
|
|
|
|
96
|
|
|
|
|
|
|
See L for details of bowtie |
|
97
|
|
|
|
|
|
|
options. |
|
98
|
|
|
|
|
|
|
|
|
99
|
|
|
|
|
|
|
=head1 FILES |
|
100
|
|
|
|
|
|
|
|
|
101
|
|
|
|
|
|
|
When a command requires filenames, these are provided to the C method, not |
|
102
|
|
|
|
|
|
|
the constructor (C). To see the set of files required by a command, use |
|
103
|
|
|
|
|
|
|
C or the alias C: |
|
104
|
|
|
|
|
|
|
|
|
105
|
|
|
|
|
|
|
$bowtiefac = Bio::Tools::Run::Bowtie->new( -command => 'paired' ); |
|
106
|
|
|
|
|
|
|
@filespec = $bowtiefac->filespec; |
|
107
|
|
|
|
|
|
|
|
|
108
|
|
|
|
|
|
|
This example returns the following array: |
|
109
|
|
|
|
|
|
|
|
|
110
|
|
|
|
|
|
|
ind |
|
111
|
|
|
|
|
|
|
seq |
|
112
|
|
|
|
|
|
|
seq2 |
|
113
|
|
|
|
|
|
|
#out |
|
114
|
|
|
|
|
|
|
|
|
115
|
|
|
|
|
|
|
This indicates that ind (C index file base name), seq (fasta/fastq),and seq2 |
|
116
|
|
|
|
|
|
|
(fasta/fastq) files MUST be specified, and that the out file MAY be specified. Use |
|
117
|
|
|
|
|
|
|
these in the C call like so: |
|
118
|
|
|
|
|
|
|
|
|
119
|
|
|
|
|
|
|
$bowtiefac->run( -ind => 'index_base', -seq => 'seq-a.fq', |
|
120
|
|
|
|
|
|
|
-seq2 => 'seq-b.fq', -out => 'align.out' ); |
|
121
|
|
|
|
|
|
|
|
|
122
|
|
|
|
|
|
|
Note that named parameters in this form allow you to specify the location of the outfile; |
|
123
|
|
|
|
|
|
|
without named parameters, the outfile is located in a tempdir and does not persist beyond |
|
124
|
|
|
|
|
|
|
the life of the object - with the exception of index creation. |
|
125
|
|
|
|
|
|
|
|
|
126
|
|
|
|
|
|
|
The object will store the programs STDOUT and STDERR output for you in the C |
|
127
|
|
|
|
|
|
|
and C attributes: |
|
128
|
|
|
|
|
|
|
|
|
129
|
|
|
|
|
|
|
handle_map_warning($bowtiefac) if ($bowtiefac->stderr =~ /warning/); |
|
130
|
|
|
|
|
|
|
|
|
131
|
|
|
|
|
|
|
=head1 FEEDBACK |
|
132
|
|
|
|
|
|
|
|
|
133
|
|
|
|
|
|
|
=head2 Mailing Lists |
|
134
|
|
|
|
|
|
|
|
|
135
|
|
|
|
|
|
|
User feedback is an integral part of the evolution of this and other |
|
136
|
|
|
|
|
|
|
Bioperl modules. Send your comments and suggestions preferably to |
|
137
|
|
|
|
|
|
|
the Bioperl mailing list. Your participation is much appreciated. |
|
138
|
|
|
|
|
|
|
|
|
139
|
|
|
|
|
|
|
bioperl-l@bioperl.org - General discussion |
|
140
|
|
|
|
|
|
|
http://bioperl.org/wiki/Mailing_lists - About the mailing lists |
|
141
|
|
|
|
|
|
|
|
|
142
|
|
|
|
|
|
|
=head2 Support |
|
143
|
|
|
|
|
|
|
|
|
144
|
|
|
|
|
|
|
Please direct usage questions or support issues to the mailing list: |
|
145
|
|
|
|
|
|
|
|
|
146
|
|
|
|
|
|
|
L |
|
147
|
|
|
|
|
|
|
|
|
148
|
|
|
|
|
|
|
Rather than to the module maintainer directly. Many experienced and |
|
149
|
|
|
|
|
|
|
reponsive experts will be able look at the problem and quickly |
|
150
|
|
|
|
|
|
|
address it. Please include a thorough description of the problem |
|
151
|
|
|
|
|
|
|
with code and data examples if at all possible. |
|
152
|
|
|
|
|
|
|
|
|
153
|
|
|
|
|
|
|
=head2 Reporting Bugs |
|
154
|
|
|
|
|
|
|
|
|
155
|
|
|
|
|
|
|
Report bugs to the Bioperl bug tracking system to help us keep track |
|
156
|
|
|
|
|
|
|
of the bugs and their resolution. Bug reports can be submitted via |
|
157
|
|
|
|
|
|
|
the web: |
|
158
|
|
|
|
|
|
|
|
|
159
|
|
|
|
|
|
|
http://redmine.open-bio.org/projects/bioperl/ |
|
160
|
|
|
|
|
|
|
|
|
161
|
|
|
|
|
|
|
=head1 AUTHOR - Dan Kortschak |
|
162
|
|
|
|
|
|
|
|
|
163
|
|
|
|
|
|
|
Email dan.kortschak adelaide.edu.au |
|
164
|
|
|
|
|
|
|
|
|
165
|
|
|
|
|
|
|
=head1 CONTRIBUTORS |
|
166
|
|
|
|
|
|
|
|
|
167
|
|
|
|
|
|
|
Mark A. Jensen (maj -at- fortinbras -dot- us) |
|
168
|
|
|
|
|
|
|
|
|
169
|
|
|
|
|
|
|
=head1 APPENDIX |
|
170
|
|
|
|
|
|
|
|
|
171
|
|
|
|
|
|
|
The rest of the documentation details each of the object methods. |
|
172
|
|
|
|
|
|
|
Internal methods are usually preceded with a _ |
|
173
|
|
|
|
|
|
|
|
|
174
|
|
|
|
|
|
|
=cut |
|
175
|
|
|
|
|
|
|
|
|
176
|
|
|
|
|
|
|
# Let the code begin... |
|
177
|
|
|
|
|
|
|
|
|
178
|
|
|
|
|
|
|
|
|
179
|
|
|
|
|
|
|
package Bio::Tools::Run::Bowtie; |
|
180
|
1
|
|
|
1
|
|
124184
|
use strict; |
|
|
1
|
|
|
|
|
2
|
|
|
|
1
|
|
|
|
|
31
|
|
|
181
|
|
|
|
|
|
|
our $HAVE_IO_UNCOMPRESS; |
|
182
|
|
|
|
|
|
|
|
|
183
|
|
|
|
|
|
|
BEGIN { |
|
184
|
1
|
|
|
1
|
|
49
|
eval 'require IO::Uncompress::Gunzip; $HAVE_IO_UNCOMPRESS = 1'; |
|
185
|
|
|
|
|
|
|
} |
|
186
|
|
|
|
|
|
|
|
|
187
|
1
|
|
|
1
|
|
6
|
use IPC::Run; |
|
|
1
|
|
|
|
|
1
|
|
|
|
1
|
|
|
|
|
36
|
|
|
188
|
|
|
|
|
|
|
|
|
189
|
|
|
|
|
|
|
# Object preamble - inherits from Bio::Root::Root |
|
190
|
|
|
|
|
|
|
|
|
191
|
1
|
|
|
1
|
|
374
|
use lib '../../..'; |
|
|
1
|
|
|
|
|
447
|
|
|
|
1
|
|
|
|
|
4
|
|
|
192
|
1
|
|
|
1
|
|
494
|
use Bio::Tools::Run::Bowtie::Config; |
|
|
1
|
|
|
|
|
3
|
|
|
|
1
|
|
|
|
|
125
|
|
|
193
|
1
|
|
|
1
|
|
367
|
use Bio::Tools::Run::WrapperBase; |
|
|
1
|
|
|
|
|
2
|
|
|
|
1
|
|
|
|
|
8
|
|
|
194
|
1
|
|
|
1
|
|
515
|
use Bio::Tools::Run::WrapperBase::CommandExts; |
|
|
1
|
|
|
|
|
2
|
|
|
|
1
|
|
|
|
|
7
|
|
|
195
|
1
|
|
|
1
|
|
504
|
use Bio::Tools::GuessSeqFormat; |
|
|
1
|
|
|
|
|
2552
|
|
|
|
1
|
|
|
|
|
9
|
|
|
196
|
1
|
|
|
1
|
|
381
|
use Bio::Tools::Run::Samtools; |
|
|
1
|
|
|
|
|
1
|
|
|
|
1
|
|
|
|
|
6
|
|
|
197
|
1
|
|
|
1
|
|
509
|
use Bio::Seq; |
|
|
1
|
|
|
|
|
25386
|
|
|
|
1
|
|
|
|
|
8
|
|
|
198
|
1
|
|
|
1
|
|
25
|
use File::Basename; |
|
|
1
|
|
|
|
|
1
|
|
|
|
1
|
|
|
|
|
55
|
|
|
199
|
|
|
|
|
|
|
|
|
200
|
1
|
|
|
1
|
|
3
|
use base qw( Bio::Tools::Run::WrapperBase Bio::Tools::Run::AssemblerBase ); |
|
|
1
|
|
|
|
|
1
|
|
|
|
1
|
|
|
|
|
455
|
|
|
201
|
|
|
|
|
|
|
|
|
202
|
|
|
|
|
|
|
## bowtie |
|
203
|
|
|
|
|
|
|
our $program_name = '*bowtie'; |
|
204
|
|
|
|
|
|
|
our $default_cmd = 'single'; |
|
205
|
|
|
|
|
|
|
|
|
206
|
|
|
|
|
|
|
our $asm_format; # this is determined dynamically |
|
207
|
|
|
|
|
|
|
|
|
208
|
|
|
|
|
|
|
# Note: |
|
209
|
|
|
|
|
|
|
# other globals required by Bio::Tools::Run::AssemblerBase are |
|
210
|
|
|
|
|
|
|
# imported from Bio::Tools::Run::Bowtie::Config |
|
211
|
|
|
|
|
|
|
|
|
212
|
|
|
|
|
|
|
our $qual_param = undef; |
|
213
|
|
|
|
|
|
|
our $use_dash = 'mixed'; |
|
214
|
|
|
|
|
|
|
our $join = ' '; |
|
215
|
|
|
|
|
|
|
|
|
216
|
|
|
|
|
|
|
=head2 new() |
|
217
|
|
|
|
|
|
|
|
|
218
|
|
|
|
|
|
|
Title : new |
|
219
|
|
|
|
|
|
|
Usage : my $obj = new Bio::Tools::Run::Bowtie(); |
|
220
|
|
|
|
|
|
|
Function: Builds a new Bio::Tools::Run::Bowtie object |
|
221
|
|
|
|
|
|
|
Returns : an instance of Bio::Tools::Run::Bowtie |
|
222
|
|
|
|
|
|
|
Args : |
|
223
|
|
|
|
|
|
|
|
|
224
|
|
|
|
|
|
|
=cut |
|
225
|
|
|
|
|
|
|
|
|
226
|
|
|
|
|
|
|
sub new { |
|
227
|
2
|
|
|
2
|
1
|
83
|
my ($class,@args) = @_; |
|
228
|
2
|
50
|
|
|
|
18
|
unless (grep /command/, @args) { |
|
229
|
0
|
|
|
|
|
0
|
push @args, '-command', $default_cmd; |
|
230
|
|
|
|
|
|
|
} |
|
231
|
|
|
|
|
|
|
#default to SAM output if no other format specified and we are running an alignment |
|
232
|
2
|
|
|
|
|
12
|
my %args=@args; |
|
233
|
2
|
50
|
|
|
|
20
|
if ($args{'-command'} =~ m/(?:single|paired|crossbow)/) { |
|
234
|
2
|
50
|
|
|
|
19
|
unless (grep /(?:sam_format|concise|quiet|refout|refidx)/, @args) { |
|
235
|
2
|
|
|
|
|
6
|
push @args, ('-sam_format', 1); |
|
236
|
|
|
|
|
|
|
} |
|
237
|
|
|
|
|
|
|
} |
|
238
|
2
|
|
|
|
|
42
|
my $self = $class->SUPER::new(@args); |
|
239
|
2
|
|
|
|
|
9
|
foreach (keys %command_executables) { |
|
240
|
10
|
|
|
|
|
33
|
$self->executables($_, $self->_find_executable($command_executables{$_})); |
|
241
|
|
|
|
|
|
|
} |
|
242
|
2
|
|
|
|
|
10
|
my ($want) = $self->_rearrange([qw(WANT)],@args); |
|
243
|
2
|
|
|
|
|
69
|
$self->want($want); |
|
244
|
2
|
|
|
|
|
6
|
$asm_format = $self->_assembly_format; |
|
245
|
2
|
|
|
|
|
5
|
$self->parameters_changed(1); # set on instantiation, per Bio::ParameterBaseI |
|
246
|
2
|
|
|
|
|
14
|
return $self; |
|
247
|
|
|
|
|
|
|
} |
|
248
|
|
|
|
|
|
|
|
|
249
|
|
|
|
|
|
|
=head2 run() |
|
250
|
|
|
|
|
|
|
|
|
251
|
|
|
|
|
|
|
Title : run |
|
252
|
|
|
|
|
|
|
Usage : $assembly = $bowtie_assembler->run($read1_fastq_file, |
|
253
|
|
|
|
|
|
|
$index_location, |
|
254
|
|
|
|
|
|
|
$read2_fastq_file); |
|
255
|
|
|
|
|
|
|
$assembly = $bowtie_assembler->run(%params); |
|
256
|
|
|
|
|
|
|
Function: Run the bowtie assembly pipeline. |
|
257
|
|
|
|
|
|
|
Returns : Assembly results (file, IO object or Assembly object) |
|
258
|
|
|
|
|
|
|
Args : - fastq file containing single-end reads |
|
259
|
|
|
|
|
|
|
- name of the base of the bowtie index |
|
260
|
|
|
|
|
|
|
- [optional] fastq file containing paired-end reads |
|
261
|
|
|
|
|
|
|
Named params are also available with args: |
|
262
|
|
|
|
|
|
|
-seq, -seq2, -ind (bowtie index), -ref (fasta reference) and -out |
|
263
|
|
|
|
|
|
|
Note : gzipped inputs are allowed if IO::Uncompress::Gunzip |
|
264
|
|
|
|
|
|
|
is available |
|
265
|
|
|
|
|
|
|
The behaviour for locating indexes follows the definition in |
|
266
|
|
|
|
|
|
|
the bowtie manual - you may use the environment variable |
|
267
|
|
|
|
|
|
|
BOWTIE_INDEXES to specify the index path or use an 'indexes' |
|
268
|
|
|
|
|
|
|
directory under the directory where the bowtie executable |
|
269
|
|
|
|
|
|
|
is located |
|
270
|
|
|
|
|
|
|
|
|
271
|
|
|
|
|
|
|
=cut |
|
272
|
|
|
|
|
|
|
|
|
273
|
|
|
|
|
|
|
sub run { |
|
274
|
0
|
|
|
0
|
1
|
0
|
my $self = shift; |
|
275
|
|
|
|
|
|
|
|
|
276
|
0
|
|
|
|
|
0
|
my ($arg1, $arg2, $arg3); # these are useless names because the different |
|
277
|
|
|
|
|
|
|
# programs take very different arguments |
|
278
|
0
|
|
|
|
|
0
|
my ($index, $seq, $seq2, $ref, $out); # these are the meaningful names that are used |
|
279
|
|
|
|
|
|
|
# with named args |
|
280
|
|
|
|
|
|
|
|
|
281
|
0
|
0
|
|
|
|
0
|
if (!(@_ % 2)) { |
|
282
|
0
|
|
|
|
|
0
|
my %args = @_; |
|
283
|
0
|
0
|
|
|
|
0
|
if ((grep /^-\w+/, keys %args) == keys %args) { |
|
|
|
0
|
|
|
|
|
|
|
284
|
0
|
|
|
|
|
0
|
($index, $seq, $seq2, $ref, $out) = |
|
285
|
|
|
|
|
|
|
$self->_rearrange([qw( IND SEQ SEQ2 REF OUT )], @_); |
|
286
|
|
|
|
|
|
|
} elsif (grep /^-\w+/, keys %args) { |
|
287
|
0
|
|
|
|
|
0
|
$self->throw("Badly formed named args: ".join(' ',@_)); |
|
288
|
|
|
|
|
|
|
} else { |
|
289
|
0
|
|
|
|
|
0
|
($arg1, $arg2) = @_; |
|
290
|
|
|
|
|
|
|
} |
|
291
|
|
|
|
|
|
|
} else { |
|
292
|
0
|
0
|
|
|
|
0
|
if (grep /^-\w+/, @_) { |
|
293
|
0
|
|
|
|
|
0
|
$self->throw("Badly formed named args: ".join(' ',@_)); |
|
294
|
|
|
|
|
|
|
} else { |
|
295
|
0
|
|
|
|
|
0
|
($arg1, $arg2, $arg3) = @_; |
|
296
|
|
|
|
|
|
|
} |
|
297
|
|
|
|
|
|
|
} |
|
298
|
|
|
|
|
|
|
|
|
299
|
|
|
|
|
|
|
# Sanity checks |
|
300
|
0
|
|
|
|
|
0
|
$self->_check_executable(); |
|
301
|
0
|
0
|
|
|
|
0
|
my $cmd = $self->command if $self->can('command'); |
|
302
|
|
|
|
|
|
|
|
|
303
|
0
|
|
|
|
|
0
|
for ($cmd) { |
|
304
|
0
|
0
|
|
|
|
0
|
m/(?:single|paired|crossbow)/ && do { |
|
305
|
0
|
|
0
|
|
|
0
|
$seq ||= $arg1; |
|
306
|
0
|
|
0
|
|
|
0
|
$index ||= $arg2; |
|
307
|
0
|
|
0
|
|
|
0
|
$seq2 ||= $arg3; |
|
308
|
0
|
0
|
|
|
|
0
|
$seq or $self->throw("Fasta/fastq/raw read(s) file/Bio::Seq required at arg 1/-seq"); |
|
309
|
0
|
0
|
|
|
|
0
|
$index or $self->throw("Bowtie index base required at arg 2/-index"); |
|
310
|
|
|
|
|
|
|
|
|
311
|
|
|
|
|
|
|
# expand gzipped files as nec. |
|
312
|
0
|
|
|
|
|
0
|
for ($seq, $seq2) { |
|
313
|
0
|
0
|
|
|
|
0
|
next unless $_; |
|
314
|
0
|
0
|
|
|
|
0
|
if (/\.gz[^.]*$/) { |
|
315
|
0
|
0
|
|
|
|
0
|
unless ($HAVE_IO_UNCOMPRESS) { |
|
316
|
0
|
|
|
|
|
0
|
croak( "IO::Uncompress::Gunzip not available, can't expand '$_'" ); |
|
317
|
|
|
|
|
|
|
} |
|
318
|
0
|
|
|
|
|
0
|
my ($tfh, $tf) = $self->io->tempfile; |
|
319
|
0
|
|
|
|
|
0
|
my $z = IO::Uncompress::Gunzip->new($_); |
|
320
|
0
|
|
|
|
|
0
|
while (<$z>) { print $tfh $_ } |
|
|
0
|
|
|
|
|
0
|
|
|
321
|
0
|
|
|
|
|
0
|
close $tfh; |
|
322
|
0
|
|
|
|
|
0
|
$_ = $tf; |
|
323
|
|
|
|
|
|
|
} |
|
324
|
|
|
|
|
|
|
} |
|
325
|
|
|
|
|
|
|
|
|
326
|
|
|
|
|
|
|
# confirm index files exist |
|
327
|
|
|
|
|
|
|
$self->_validate_file_input( -ind => $index ) or |
|
328
|
|
|
|
|
|
|
($self->_validate_file_input( -ind => $self->io->catfile(dirname($self->executable),'indexes',$index)) and |
|
329
|
|
|
|
|
|
|
$index = $self->io->catfile(dirname($self->executable),'indexes',$index)) or |
|
330
|
|
|
|
|
|
|
($self->_validate_file_input( -ind => $self->io->catfile($ENV{BOWTIE_INDEXES},$index)) and |
|
331
|
0
|
0
|
0
|
|
|
0
|
$index = $self->io->catfile($ENV{BOWTIE_INDEXES},$index)) or |
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
332
|
|
|
|
|
|
|
$self->throw("Incorrect filetype (expecting bowtie index) or absent file arg 2/-index"); |
|
333
|
|
|
|
|
|
|
|
|
334
|
|
|
|
|
|
|
# bowtie prepare the multiple input types |
|
335
|
0
|
|
|
|
|
0
|
$seq = $self->_prepare_input_sequences($seq); |
|
336
|
0
|
0
|
|
|
|
0
|
if ($cmd =~ m/^p/) { |
|
337
|
0
|
0
|
|
|
|
0
|
$seq2 && ($seq2 = $self->_prepare_input_sequences($seq2)); |
|
338
|
|
|
|
|
|
|
} else { |
|
339
|
0
|
0
|
|
|
|
0
|
$seq2 && $self->throw("Second sequence input not wanted for command: $cmd"); |
|
340
|
|
|
|
|
|
|
} |
|
341
|
|
|
|
|
|
|
|
|
342
|
|
|
|
|
|
|
# Assemble |
|
343
|
0
|
|
|
|
|
0
|
my $format = $self->_assembly_format; |
|
344
|
0
|
|
|
|
|
0
|
my $suffix = '.'.$format; |
|
345
|
|
|
|
|
|
|
|
|
346
|
0
|
0
|
|
|
|
0
|
if ($out) { |
|
347
|
0
|
|
|
|
|
0
|
$out .= $suffix; |
|
348
|
|
|
|
|
|
|
} else { |
|
349
|
0
|
|
|
|
|
0
|
my ($bowtieh, $bowtief) = $self->io->tempfile( -dir => $self->tempdir(), -suffix => $suffix ); |
|
350
|
0
|
|
|
|
|
0
|
$bowtieh->close; |
|
351
|
0
|
|
|
|
|
0
|
$out = $bowtief; |
|
352
|
|
|
|
|
|
|
} |
|
353
|
|
|
|
|
|
|
|
|
354
|
0
|
|
|
|
|
0
|
my %params = ( -ind => $index, -seq => $seq, -seq2 => $seq2, -out => $out ); |
|
355
|
|
|
|
|
|
|
map { |
|
356
|
0
|
0
|
|
|
|
0
|
delete $params{$_} unless defined $params{$_} |
|
|
0
|
|
|
|
|
0
|
|
|
357
|
|
|
|
|
|
|
} keys %params; |
|
358
|
0
|
|
|
|
|
0
|
$self->_run(%params); |
|
359
|
|
|
|
|
|
|
|
|
360
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'index'} = $index; |
|
361
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'file_name'} = $out; |
|
362
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'format'} = $format; |
|
363
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'file'} = Bio::Root::IO->new( -file => $out ); |
|
364
|
|
|
|
|
|
|
|
|
365
|
0
|
|
|
|
|
0
|
return $self->result; |
|
366
|
|
|
|
|
|
|
}; |
|
367
|
|
|
|
|
|
|
|
|
368
|
0
|
0
|
|
|
|
0
|
m/build/ && do { |
|
369
|
0
|
|
0
|
|
|
0
|
$ref ||= $arg1; |
|
370
|
0
|
|
0
|
|
|
0
|
$index ||= $arg2; |
|
371
|
0
|
0
|
|
|
|
0
|
$ref or $self->throw("Fasta read(s) file/Bio::Seq required at arg 1/-ref"); |
|
372
|
0
|
|
0
|
|
|
0
|
$index ||= $self->io->tempdir(CLEANUP => 1).'/index'; # we want a new one each time |
|
373
|
0
|
0
|
|
|
|
0
|
$arg3 && $self->throw("Second sequence input not wanted for command: $cmd"); |
|
374
|
|
|
|
|
|
|
|
|
375
|
0
|
|
|
|
|
0
|
my $format = $self->_assembly_format; |
|
376
|
|
|
|
|
|
|
|
|
377
|
|
|
|
|
|
|
# expand gzipped file as nec. |
|
378
|
0
|
0
|
|
|
|
0
|
if ($ref =~ (m/\.gz[^.]*$/)) { |
|
379
|
0
|
0
|
|
|
|
0
|
unless ($HAVE_IO_UNCOMPRESS) { |
|
380
|
0
|
|
|
|
|
0
|
croak( "IO::Uncompress::Gunzip not available, can't expand '$_'" ); |
|
381
|
|
|
|
|
|
|
} |
|
382
|
0
|
|
|
|
|
0
|
my ($tfh, $tf) = $self->io->tempfile; |
|
383
|
0
|
|
|
|
|
0
|
my $z = IO::Uncompress::Gunzip->new($_); |
|
384
|
0
|
|
|
|
|
0
|
while (<$z>) { print $tfh $_ } |
|
|
0
|
|
|
|
|
0
|
|
|
385
|
0
|
|
|
|
|
0
|
close $tfh; |
|
386
|
0
|
|
|
|
|
0
|
$ref = $tf; |
|
387
|
|
|
|
|
|
|
} |
|
388
|
|
|
|
|
|
|
|
|
389
|
|
|
|
|
|
|
# bowtie prepare the two input types for the first argument |
|
390
|
0
|
|
|
|
|
0
|
$ref = $self->_prepare_input_sequences($ref); |
|
391
|
|
|
|
|
|
|
|
|
392
|
|
|
|
|
|
|
# Build index |
|
393
|
0
|
|
|
|
|
0
|
$self->_run( -ref => $ref, -out => $index ); |
|
394
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'format'} = $format; |
|
395
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'file_name'} = $index; |
|
396
|
|
|
|
|
|
|
|
|
397
|
0
|
|
|
|
|
0
|
return $index; |
|
398
|
|
|
|
|
|
|
}; |
|
399
|
|
|
|
|
|
|
|
|
400
|
0
|
0
|
|
|
|
0
|
m/inspect/ && do { |
|
401
|
0
|
|
0
|
|
|
0
|
$index ||= $arg1; |
|
402
|
0
|
|
0
|
|
|
0
|
$out ||= $arg2; |
|
403
|
0
|
0
|
|
|
|
0
|
$index or $self->throw("Bowtie index required at arg 1"); |
|
404
|
|
|
|
|
|
|
|
|
405
|
|
|
|
|
|
|
$self->_validate_file_input( -ind => $index ) or |
|
406
|
|
|
|
|
|
|
($self->_validate_file_input( -ind => $self->io->catfile(dirname($self->executable),'indexes',$index)) and |
|
407
|
|
|
|
|
|
|
$index = $self->io->catfile(dirname($self->executable),'indexes',$index)) or |
|
408
|
|
|
|
|
|
|
($self->_validate_file_input( -ind => $self->io->catfile($ENV{BOWTIE_INDEXES},$index)) and |
|
409
|
0
|
0
|
0
|
|
|
0
|
$index = $self->io->catfile($ENV{BOWTIE_INDEXES},$index)) or |
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
410
|
|
|
|
|
|
|
$self->throw("'$index' doesn't look like a bowtie index or index component is missing at arg 1/-ind"); |
|
411
|
0
|
0
|
|
|
|
0
|
$arg3 && $self->throw("Second sequence input not wanted for command: $cmd"); |
|
412
|
|
|
|
|
|
|
|
|
413
|
|
|
|
|
|
|
# Inspect index |
|
414
|
0
|
|
|
|
|
0
|
my $format = $self->_assembly_format; |
|
415
|
0
|
|
|
|
|
0
|
my $suffix = '.'.$format; |
|
416
|
|
|
|
|
|
|
|
|
417
|
0
|
0
|
|
|
|
0
|
if ($out) { |
|
418
|
0
|
|
|
|
|
0
|
$out .= $suffix; |
|
419
|
|
|
|
|
|
|
} else { |
|
420
|
0
|
|
|
|
|
0
|
my ($desch, $descf) = $self->io->tempfile( -dir => $self->tempdir(), -suffix => $suffix ); |
|
421
|
0
|
|
|
|
|
0
|
$desch->close; |
|
422
|
0
|
|
|
|
|
0
|
$out = $descf; |
|
423
|
|
|
|
|
|
|
} |
|
424
|
|
|
|
|
|
|
|
|
425
|
0
|
|
|
|
|
0
|
$self->_run( -ind => $index, -out => $out ); |
|
426
|
|
|
|
|
|
|
|
|
427
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'file_name'} = $out; |
|
428
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'format'} = $format; |
|
429
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'file'} = Bio::Root::IO->new( -file => $out ); |
|
430
|
|
|
|
|
|
|
|
|
431
|
0
|
|
|
|
|
0
|
return $self->result; |
|
432
|
|
|
|
|
|
|
} |
|
433
|
|
|
|
|
|
|
} |
|
434
|
|
|
|
|
|
|
} |
|
435
|
|
|
|
|
|
|
|
|
436
|
|
|
|
|
|
|
=head2 want() |
|
437
|
|
|
|
|
|
|
|
|
438
|
|
|
|
|
|
|
Title : want |
|
439
|
|
|
|
|
|
|
Usage : $bowtiefac->want( $class ) |
|
440
|
|
|
|
|
|
|
Function: make factory return $class, or raw (scalar) results in file |
|
441
|
|
|
|
|
|
|
Returns : return wanted type |
|
442
|
|
|
|
|
|
|
Args : [optional] string indicating class or raw of wanted result |
|
443
|
|
|
|
|
|
|
|
|
444
|
|
|
|
|
|
|
=cut |
|
445
|
|
|
|
|
|
|
|
|
446
|
|
|
|
|
|
|
sub want { |
|
447
|
2
|
|
|
2
|
1
|
2
|
my $self = shift; |
|
448
|
2
|
50
|
|
|
|
9
|
return $self->{'_want'} = shift if @_; |
|
449
|
0
|
|
|
|
|
0
|
return $self->{'_want'}; |
|
450
|
|
|
|
|
|
|
} |
|
451
|
|
|
|
|
|
|
|
|
452
|
|
|
|
|
|
|
=head2 result() |
|
453
|
|
|
|
|
|
|
|
|
454
|
|
|
|
|
|
|
Title : result |
|
455
|
|
|
|
|
|
|
Usage : $bowtiefac->result( [-want => $type|$format] ) |
|
456
|
|
|
|
|
|
|
Function: return result in wanted format |
|
457
|
|
|
|
|
|
|
Returns : results |
|
458
|
|
|
|
|
|
|
Args : [optional] hashref of wanted type |
|
459
|
|
|
|
|
|
|
|
|
460
|
|
|
|
|
|
|
=cut |
|
461
|
|
|
|
|
|
|
|
|
462
|
|
|
|
|
|
|
sub result { |
|
463
|
0
|
|
|
0
|
1
|
0
|
my ($self, @args) = @_; |
|
464
|
|
|
|
|
|
|
|
|
465
|
0
|
0
|
|
|
|
0
|
my $want = $self->want ? $self->want : $self->want($self->_rearrange([qw(WANT)],@args)); |
|
466
|
0
|
0
|
|
|
|
0
|
my $cmd = $self->command if $self->can('command'); |
|
467
|
0
|
|
|
|
|
0
|
my $format = $self->{'_result'}->{'format'}; |
|
468
|
|
|
|
|
|
|
|
|
469
|
0
|
0
|
0
|
|
|
0
|
return $self->{'_result'}->{'format'} if (defined $want && $want eq 'format'); |
|
470
|
0
|
0
|
0
|
|
|
0
|
return $self->{'_result'}->{'file_name'} if (!$want || $want eq 'raw' || $cmd eq 'build'); |
|
|
|
|
0
|
|
|
|
|
|
471
|
0
|
0
|
|
|
|
0
|
return $self->{'_result'}->{'file'} if ($want =~ m/^Bio::Root::IO/); |
|
472
|
|
|
|
|
|
|
|
|
473
|
0
|
|
|
|
|
0
|
for ($cmd) { |
|
474
|
0
|
0
|
|
|
|
0
|
m/(?:single|paired|crossbow)/ && do { |
|
475
|
0
|
|
|
|
|
0
|
my $scaffold; |
|
476
|
0
|
|
|
|
|
0
|
for ($format) { |
|
477
|
0
|
0
|
0
|
|
|
0
|
m/^bowtie/i && $want =~ m/^Bio::Assembly::Scaffold/ && do { |
|
478
|
0
|
0
|
0
|
|
|
0
|
unless (defined $self->{'_result'}->{'object'} && |
|
479
|
|
|
|
|
|
|
ref($self->{'_result'}->{'object'}) =~ m/^Bio::Assembly::Scaffold/) { |
|
480
|
|
|
|
|
|
|
$self->{'_result'}->{'object'} = |
|
481
|
|
|
|
|
|
|
$self->_export_results( $self->{'_result'}->{'file_name'}, |
|
482
|
0
|
|
|
|
|
0
|
-index => $self->{'_result'}->{'index'}, |
|
483
|
|
|
|
|
|
|
-keep_asm => 1 ); |
|
484
|
|
|
|
|
|
|
} |
|
485
|
0
|
|
|
|
|
0
|
last; |
|
486
|
|
|
|
|
|
|
}; |
|
487
|
0
|
0
|
0
|
|
|
0
|
m/^bowtie/i && $want =~ m/^Bio::SeqFeature::Collection/ && do { |
|
488
|
0
|
|
|
|
|
0
|
$self->warn("Don't know how to create a $want object for $cmd with bowtie format - try SAM format."); |
|
489
|
0
|
|
|
|
|
0
|
last; |
|
490
|
|
|
|
|
|
|
}; |
|
491
|
0
|
0
|
0
|
|
|
0
|
m/^sam/i && $want =~ m/^Bio::Assembly::Scaffold/ && do { |
|
492
|
0
|
0
|
0
|
|
|
0
|
unless (defined $self->{'_result'}->{'object'} && |
|
493
|
|
|
|
|
|
|
ref($self->{'_result'}->{'object'}) =~ m/^Bio::Assembly::Scaffold/) { |
|
494
|
0
|
|
|
|
|
0
|
my $bamf = $self->_make_bam($self->{'_result'}->{'file_name'}, 1); |
|
495
|
0
|
|
|
|
|
0
|
my $inspector = Bio::Tools::Run::Bowtie->new( -command => 'inspect' ); |
|
496
|
0
|
|
|
|
|
0
|
my $refdb = $inspector->run($self->{'_result'}->{'index'}); |
|
497
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'object'} = |
|
498
|
|
|
|
|
|
|
$self->_export_results($bamf, -refdb => $refdb, -keep_asm => 1 ); |
|
499
|
|
|
|
|
|
|
} |
|
500
|
0
|
|
|
|
|
0
|
last; |
|
501
|
|
|
|
|
|
|
}; |
|
502
|
0
|
0
|
0
|
|
|
0
|
m/^sam/i && $want =~ m/^Bio::SeqFeature::Collection/ && do { |
|
503
|
0
|
0
|
0
|
|
|
0
|
unless (defined $self->{'_result'}->{'object'} && |
|
504
|
|
|
|
|
|
|
ref($self->{'_result'}->{'object'}) =~ m/^Bio::Assembly::Scaffold/) { |
|
505
|
0
|
|
|
|
|
0
|
my $bamf = $self->_make_bam($self->{'_result'}->{'file_name'}, 0); |
|
506
|
0
|
|
|
|
|
0
|
my $convert = Bio::Tools::Run::BEDTools->new( -command => 'bam_to_bed' ); |
|
507
|
0
|
|
|
|
|
0
|
my $bedf = $convert->run( -bed => $bamf ); |
|
508
|
0
|
|
|
|
|
0
|
my $merge = Bio::Tools::Run::BEDTools->new( -command => 'merge' ); |
|
509
|
0
|
|
|
|
|
0
|
$merge->run($self->{'_result'}->{'index'}); |
|
510
|
0
|
|
|
|
|
0
|
$self->{'_result'}->{'object'} = $merge->result( -want => $want ); |
|
511
|
|
|
|
|
|
|
} |
|
512
|
0
|
|
|
|
|
0
|
last; |
|
513
|
|
|
|
|
|
|
}; |
|
514
|
0
|
|
|
|
|
0
|
do { |
|
515
|
0
|
|
|
|
|
0
|
$self->warn("Don't know how to create a $want object for $cmd."); |
|
516
|
0
|
|
|
|
|
0
|
return; |
|
517
|
|
|
|
|
|
|
} |
|
518
|
|
|
|
|
|
|
}; |
|
519
|
0
|
|
|
|
|
0
|
last; |
|
520
|
|
|
|
|
|
|
}; |
|
521
|
0
|
0
|
|
|
|
0
|
m/inspect/ && do { |
|
522
|
0
|
|
|
|
|
0
|
for ($want) { |
|
523
|
0
|
0
|
0
|
|
|
0
|
m/^Bio::SeqIO/ && $format eq 'fasta' && do { |
|
524
|
0
|
0
|
0
|
|
|
0
|
unless (defined $self->{'_result'}->{'object'} && |
|
525
|
|
|
|
|
|
|
ref($self->{'_result'}->{'object'}) =~ m/^Bio::SeqIO/) { |
|
526
|
|
|
|
|
|
|
$self->{'_result'}->{'object'} = |
|
527
|
0
|
|
|
|
|
0
|
Bio::SeqIO->new(-file => $self->{'_result'}->{'file'}, |
|
528
|
|
|
|
|
|
|
-format => 'fasta'); |
|
529
|
|
|
|
|
|
|
} |
|
530
|
0
|
|
|
|
|
0
|
last; |
|
531
|
|
|
|
|
|
|
}; |
|
532
|
0
|
0
|
0
|
|
|
0
|
m/^Bio::SeqIO/ && $format ne 'fasta' && do { |
|
533
|
0
|
|
|
|
|
0
|
$self->warn("Don't know how to create a $want object for names only - try -want => 'Bio::Root::IO'."); |
|
534
|
0
|
|
|
|
|
0
|
return; |
|
535
|
|
|
|
|
|
|
}; |
|
536
|
0
|
|
|
|
|
0
|
do { |
|
537
|
0
|
|
|
|
|
0
|
$self->warn("Don't know how to create a $want object for $cmd."); |
|
538
|
0
|
|
|
|
|
0
|
return; |
|
539
|
|
|
|
|
|
|
} |
|
540
|
|
|
|
|
|
|
} |
|
541
|
|
|
|
|
|
|
} |
|
542
|
|
|
|
|
|
|
} |
|
543
|
|
|
|
|
|
|
|
|
544
|
0
|
|
|
|
|
0
|
return $self->{'_result'}->{'object'}; |
|
545
|
|
|
|
|
|
|
} |
|
546
|
|
|
|
|
|
|
|
|
547
|
|
|
|
|
|
|
=head2 _determine_format() |
|
548
|
|
|
|
|
|
|
|
|
549
|
|
|
|
|
|
|
Title : _determine_format |
|
550
|
|
|
|
|
|
|
Usage : $bowtiefac->_determine_format |
|
551
|
|
|
|
|
|
|
Function: determine the format of output for current options |
|
552
|
|
|
|
|
|
|
Returns : format of bowtie output |
|
553
|
|
|
|
|
|
|
Args : |
|
554
|
|
|
|
|
|
|
|
|
555
|
|
|
|
|
|
|
=cut |
|
556
|
|
|
|
|
|
|
|
|
557
|
|
|
|
|
|
|
sub _determine_format { |
|
558
|
2
|
|
|
2
|
|
4
|
my ($self) = shift; |
|
559
|
|
|
|
|
|
|
|
|
560
|
2
|
50
|
|
|
|
66
|
my $cmd = $self->command if $self->can('command'); |
|
561
|
2
|
|
|
|
|
13
|
for ($cmd) { |
|
562
|
2
|
50
|
|
|
|
6
|
m/build/ && do { |
|
563
|
0
|
|
|
|
|
0
|
return 'ebwt'; |
|
564
|
|
|
|
|
|
|
}; |
|
565
|
2
|
50
|
|
|
|
5
|
m/inspect/ && do { |
|
566
|
0
|
0
|
|
|
|
0
|
$self->{'_summary'} && return 'text'; |
|
567
|
0
|
0
|
|
|
|
0
|
return $self->{'_names_only'} ? 'text' : 'fasta'; |
|
568
|
|
|
|
|
|
|
}; |
|
569
|
2
|
50
|
|
|
|
17
|
m/(?:single|paired|crossbow)/ && do { |
|
570
|
2
|
|
|
|
|
4
|
my $format = 'bowtie'; # this is our default position |
|
571
|
2
|
|
|
|
|
7
|
for (keys %format_lookup) { |
|
572
|
10
|
100
|
|
|
|
22
|
$format = $format_lookup{$_} if $self->{'_'.$_}; |
|
573
|
|
|
|
|
|
|
} |
|
574
|
2
|
|
|
|
|
6
|
return $format; |
|
575
|
|
|
|
|
|
|
} |
|
576
|
|
|
|
|
|
|
} |
|
577
|
|
|
|
|
|
|
} |
|
578
|
|
|
|
|
|
|
|
|
579
|
|
|
|
|
|
|
=head2 _make_bam() |
|
580
|
|
|
|
|
|
|
|
|
581
|
|
|
|
|
|
|
Title : _make_bam |
|
582
|
|
|
|
|
|
|
Usage : $bowtiefac->_make_bam( $file, $sort ) |
|
583
|
|
|
|
|
|
|
Function: make a sorted BAM format file from SAM file |
|
584
|
|
|
|
|
|
|
Returns : sorted BAM file name |
|
585
|
|
|
|
|
|
|
Args : SAM file name and boolean flag to select sorted BAM format |
|
586
|
|
|
|
|
|
|
|
|
587
|
|
|
|
|
|
|
=cut |
|
588
|
|
|
|
|
|
|
|
|
589
|
|
|
|
|
|
|
sub _make_bam { |
|
590
|
0
|
|
|
0
|
|
0
|
my ($self, $file, $sort) = @_; |
|
591
|
|
|
|
|
|
|
|
|
592
|
0
|
0
|
0
|
|
|
0
|
$self->throw("'$file' does not exist or is not readable") |
|
593
|
|
|
|
|
|
|
unless ( -e $file && -r _ ); |
|
594
|
|
|
|
|
|
|
|
|
595
|
|
|
|
|
|
|
# make a sorted bam file from a sam file input |
|
596
|
0
|
|
|
|
|
0
|
my ($bamh, $bamf) = $self->io->tempfile( -dir => $self->tempdir(), -suffix => '.bam' ); |
|
597
|
0
|
|
|
|
|
0
|
$bamh->close; |
|
598
|
|
|
|
|
|
|
|
|
599
|
0
|
|
|
|
|
0
|
my $samt = Bio::Tools::Run::Samtools->new( -command => 'view', |
|
600
|
|
|
|
|
|
|
-sam_input => 1, |
|
601
|
|
|
|
|
|
|
-bam_output => 1 ); |
|
602
|
|
|
|
|
|
|
|
|
603
|
0
|
|
|
|
|
0
|
$samt->run( -bam => $file, -out => $bamf ); |
|
604
|
|
|
|
|
|
|
|
|
605
|
0
|
0
|
|
|
|
0
|
if ($sort) { |
|
606
|
0
|
|
|
|
|
0
|
my ($srth, $srtf) = $self->io->tempfile( -dir => $self->io->tempdir(CLEANUP=>1), -suffix => '.srt' ); |
|
607
|
|
|
|
|
|
|
# shared tempdir, so make new - otherwise it is scrubbed during Bio::DB::Sam |
|
608
|
0
|
|
|
|
|
0
|
$srth->close; |
|
609
|
|
|
|
|
|
|
|
|
610
|
0
|
|
|
|
|
0
|
$samt = Bio::Tools::Run::Samtools->new( -command => 'sort' ); |
|
611
|
0
|
|
|
|
|
0
|
$samt->run( -bam => $bamf, -pfx => $srtf); |
|
612
|
|
|
|
|
|
|
|
|
613
|
0
|
|
|
|
|
0
|
return $srtf.'.bam'; |
|
614
|
|
|
|
|
|
|
} else { |
|
615
|
0
|
|
|
|
|
0
|
return $bamf; |
|
616
|
|
|
|
|
|
|
} |
|
617
|
|
|
|
|
|
|
} |
|
618
|
|
|
|
|
|
|
|
|
619
|
|
|
|
|
|
|
=head2 _validate_file_input() |
|
620
|
|
|
|
|
|
|
|
|
621
|
|
|
|
|
|
|
Title : _validate_file_input |
|
622
|
|
|
|
|
|
|
Usage : $bowtiefac->_validate_file_input( -type => $file ) |
|
623
|
|
|
|
|
|
|
Function: validate file type for file spec |
|
624
|
|
|
|
|
|
|
Returns : file type if valid type for file spec |
|
625
|
|
|
|
|
|
|
Args : hash of filespec => file_name |
|
626
|
|
|
|
|
|
|
|
|
627
|
|
|
|
|
|
|
=cut |
|
628
|
|
|
|
|
|
|
|
|
629
|
|
|
|
|
|
|
sub _validate_file_input { |
|
630
|
0
|
|
|
0
|
|
0
|
my ($self, @args) = @_; |
|
631
|
0
|
|
|
|
|
0
|
my (%args); |
|
632
|
0
|
0
|
|
|
|
0
|
if (grep (/^-/, @args)) { # named parms |
|
633
|
0
|
0
|
|
|
|
0
|
$self->throw("Wrong number of args - requires one named arg") if (@args > 2); |
|
634
|
0
|
|
|
|
|
0
|
s/^-// for @args; |
|
635
|
0
|
|
|
|
|
0
|
%args = @args; |
|
636
|
|
|
|
|
|
|
} else { |
|
637
|
0
|
|
|
|
|
0
|
$self->throw("Must provide named filespec"); |
|
638
|
|
|
|
|
|
|
} |
|
639
|
|
|
|
|
|
|
|
|
640
|
0
|
|
|
|
|
0
|
for (keys %args) { |
|
641
|
0
|
0
|
|
|
|
0
|
m/^seq|seq2|ref$/ && do { |
|
642
|
0
|
0
|
0
|
|
|
0
|
return unless ( -e $args{$_} && -r _ ); |
|
643
|
0
|
|
|
|
|
0
|
my $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$args{$_}); |
|
644
|
0
|
0
|
|
|
|
0
|
return $guesser->guess if grep {$guesser->guess =~ m/$_/} @{$accepted_types{$_}}; |
|
|
0
|
|
|
|
|
0
|
|
|
|
0
|
|
|
|
|
0
|
|
|
645
|
|
|
|
|
|
|
}; |
|
646
|
0
|
0
|
|
|
|
0
|
m/^ind$/ && do { |
|
647
|
|
|
|
|
|
|
return 'ebwt' if |
|
648
|
|
|
|
|
|
|
(-e $args{$_}.'.1.ebwt' && -e $args{$_}.'.2.ebwt' && -e $args{$_}.'.3.ebwt' && |
|
649
|
0
|
0
|
0
|
|
|
0
|
-e $args{$_}.'.4.ebwt' && -e $args{$_}.'.rev.1.ebwt' && -e $args{$_}.'.rev.2.ebwt'); |
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
|
|
|
0
|
|
|
|
|
|
650
|
|
|
|
|
|
|
} |
|
651
|
|
|
|
|
|
|
} |
|
652
|
0
|
|
|
|
|
0
|
return; |
|
653
|
|
|
|
|
|
|
} |
|
654
|
|
|
|
|
|
|
|
|
655
|
|
|
|
|
|
|
=head1 Bio::Tools::Run::AssemblerBase overrides |
|
656
|
|
|
|
|
|
|
|
|
657
|
|
|
|
|
|
|
=head2 _assembly_format() |
|
658
|
|
|
|
|
|
|
|
|
659
|
|
|
|
|
|
|
Title : _assembly_format |
|
660
|
|
|
|
|
|
|
Usage : $bowtiefac->_determine_format |
|
661
|
|
|
|
|
|
|
Function: set the format of output for current options |
|
662
|
|
|
|
|
|
|
Returns : format of bowtie output |
|
663
|
|
|
|
|
|
|
Args : |
|
664
|
|
|
|
|
|
|
|
|
665
|
|
|
|
|
|
|
=cut |
|
666
|
|
|
|
|
|
|
|
|
667
|
|
|
|
|
|
|
sub _assembly_format { |
|
668
|
2
|
|
|
2
|
|
3
|
my $self = shift; |
|
669
|
|
|
|
|
|
|
|
|
670
|
2
|
|
|
|
|
6
|
my $format = $self->_determine_format; |
|
671
|
2
|
|
|
|
|
13
|
return $self->SUPER::_assembly_format($format); |
|
672
|
|
|
|
|
|
|
} |
|
673
|
|
|
|
|
|
|
|
|
674
|
|
|
|
|
|
|
|
|
675
|
|
|
|
|
|
|
=head2 _check_sequence_input() |
|
676
|
|
|
|
|
|
|
|
|
677
|
|
|
|
|
|
|
No-op. |
|
678
|
|
|
|
|
|
|
|
|
679
|
|
|
|
|
|
|
=cut |
|
680
|
|
|
|
|
|
|
|
|
681
|
|
|
|
|
|
|
sub _check_sequence_input { |
|
682
|
0
|
|
|
0
|
|
0
|
return 1; |
|
683
|
|
|
|
|
|
|
} |
|
684
|
|
|
|
|
|
|
|
|
685
|
|
|
|
|
|
|
=head2 _check_optional_quality_input() |
|
686
|
|
|
|
|
|
|
|
|
687
|
|
|
|
|
|
|
No-op. |
|
688
|
|
|
|
|
|
|
|
|
689
|
|
|
|
|
|
|
=cut |
|
690
|
|
|
|
|
|
|
|
|
691
|
|
|
|
|
|
|
sub _check_optional_quality_input { |
|
692
|
0
|
|
|
0
|
|
0
|
return 1; |
|
693
|
|
|
|
|
|
|
} |
|
694
|
|
|
|
|
|
|
|
|
695
|
|
|
|
|
|
|
=head2 _prepare_input_sequences() |
|
696
|
|
|
|
|
|
|
|
|
697
|
|
|
|
|
|
|
Prepare and check input sequences for bowtie. |
|
698
|
|
|
|
|
|
|
|
|
699
|
|
|
|
|
|
|
=cut |
|
700
|
|
|
|
|
|
|
|
|
701
|
|
|
|
|
|
|
sub _prepare_input_sequences { |
|
702
|
|
|
|
|
|
|
|
|
703
|
0
|
|
|
0
|
|
0
|
my ($self, @args) = @_; |
|
704
|
0
|
|
|
|
|
0
|
my (%args, $read); |
|
705
|
0
|
0
|
|
|
|
0
|
if (grep (/^-/, @args)) { # named parms |
|
706
|
0
|
0
|
|
|
|
0
|
$self->throw("Input args not an even number") unless !(@args % 2); |
|
707
|
0
|
|
|
|
|
0
|
%args = @args; |
|
708
|
0
|
|
|
|
|
0
|
($read) = @args{qw( -sequence )}; |
|
709
|
|
|
|
|
|
|
} else { |
|
710
|
0
|
|
|
|
|
0
|
($read) = @args; |
|
711
|
|
|
|
|
|
|
} |
|
712
|
|
|
|
|
|
|
|
|
713
|
|
|
|
|
|
|
# Could use the AssemblerBase routine for this, except that would not permit |
|
714
|
|
|
|
|
|
|
# an array of strings |
|
715
|
0
|
0
|
|
|
|
0
|
if ($self->inline) { # expect inline data |
|
716
|
0
|
0
|
0
|
|
|
0
|
if (UNIVERSAL::isa($read,'can') && $read->isa("Bio::PrimarySeqI")) { # we have a Bio::*Seq* |
|
717
|
0
|
|
|
|
|
0
|
$read=$read->seq(); |
|
718
|
|
|
|
|
|
|
} else { # we have something else |
|
719
|
0
|
0
|
|
|
|
0
|
if (ref($read) =~ /ARRAY/i) { |
|
720
|
0
|
|
|
|
|
0
|
my @ts; |
|
721
|
0
|
|
|
|
|
0
|
foreach my $seq (@$read) { |
|
722
|
0
|
0
|
|
|
|
0
|
if ($seq->isa("Bio::PrimarySeqI")) { |
|
723
|
0
|
|
|
|
|
0
|
$seq=$seq->seq(); |
|
724
|
|
|
|
|
|
|
} else { |
|
725
|
0
|
0
|
|
|
|
0
|
next if $read=~m/[[^:alpha:]]/; |
|
726
|
|
|
|
|
|
|
} |
|
727
|
0
|
|
|
|
|
0
|
push @ts,$seq; |
|
728
|
|
|
|
|
|
|
} |
|
729
|
0
|
0
|
|
|
|
0
|
$self->throw("bowtie requires at least one sequence read") unless (@ts); |
|
730
|
0
|
0
|
|
|
|
0
|
if (@ts>1) { |
|
731
|
0
|
|
|
|
|
0
|
$read="'".join(',',@ts)."'"; |
|
732
|
|
|
|
|
|
|
} else { |
|
733
|
0
|
|
|
|
|
0
|
($read)=@ts; |
|
734
|
|
|
|
|
|
|
} |
|
735
|
|
|
|
|
|
|
} else { #must be a string... fail if non-alpha |
|
736
|
0
|
0
|
|
|
|
0
|
$self->throw("bowtie requires at least one valid sequence read") if $read=~m/[[^:alpha:]]/; |
|
737
|
|
|
|
|
|
|
} |
|
738
|
|
|
|
|
|
|
} |
|
739
|
|
|
|
|
|
|
} else { # expect file(s) - so test whether it's/they're appropriate |
|
740
|
|
|
|
|
|
|
# and make a comma-separated list of filenames |
|
741
|
0
|
0
|
|
|
|
0
|
my @ts = (ref($read) =~ /ARRAY/i) ? @$read : ($read); |
|
742
|
0
|
|
|
|
|
0
|
for my $file (@ts) { |
|
743
|
0
|
0
|
|
|
|
0
|
if ( -e $file ) { |
|
744
|
0
|
0
|
|
|
|
0
|
my $cmd = $self->command if $self->can('command'); |
|
745
|
0
|
|
|
|
|
0
|
my $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$file); |
|
746
|
0
|
|
|
|
|
0
|
for ($guesser->guess) { |
|
747
|
0
|
0
|
|
|
|
0
|
m/^fasta$/ && do { |
|
748
|
0
|
0
|
|
|
|
0
|
$cmd =~ m/^b/ && last; |
|
749
|
0
|
0
|
0
|
|
|
0
|
($cmd =~ m/^c/ or $self->fastq or $self->raw) and $self->throw("Fasta reads file '$file' inappropriate"); |
|
|
|
|
0
|
|
|
|
|
|
750
|
0
|
|
|
|
|
0
|
$self->fasta(1); |
|
751
|
0
|
|
|
|
|
0
|
last; |
|
752
|
|
|
|
|
|
|
}; |
|
753
|
0
|
0
|
|
|
|
0
|
m/^fastq$/ && do { |
|
754
|
0
|
0
|
0
|
|
|
0
|
($cmd =~ m/^[cb]/ or $self->fasta or $self->raw) and $self->throw("Fastq reads file '$file' inappropriate"); |
|
|
|
|
0
|
|
|
|
|
|
755
|
0
|
|
|
|
|
0
|
$self->fastq(1); |
|
756
|
0
|
|
|
|
|
0
|
last; |
|
757
|
|
|
|
|
|
|
}; |
|
758
|
0
|
0
|
|
|
|
0
|
m/^tab$/ && do { |
|
759
|
0
|
0
|
|
|
|
0
|
$cmd =~ m/^c/ or $self->throw("Crossbow reads file '$file' inappropriate"); # this is unrecoverable since the object has default program defined |
|
760
|
0
|
|
|
|
|
0
|
last; |
|
761
|
|
|
|
|
|
|
}; |
|
762
|
0
|
0
|
|
|
|
0
|
m/^raw$/ && do { |
|
763
|
0
|
0
|
0
|
|
|
0
|
($cmd =~ m/^[cb]/ or $self->fasta or $self->fastq) and $self->throw("Raw reads file '$file' inappropriate"); |
|
|
|
|
0
|
|
|
|
|
|
764
|
0
|
|
|
|
|
0
|
$self->raw(1); |
|
765
|
0
|
|
|
|
|
0
|
last; |
|
766
|
|
|
|
|
|
|
}; |
|
767
|
0
|
|
|
|
|
0
|
do { |
|
768
|
0
|
|
|
|
|
0
|
$self->throw("File '$file' not a recognised bowtie input filetype"); |
|
769
|
|
|
|
|
|
|
} |
|
770
|
|
|
|
|
|
|
} |
|
771
|
|
|
|
|
|
|
} else { |
|
772
|
0
|
|
|
|
|
0
|
$self->throw("Sequence read file '$file' does not exist"); |
|
773
|
|
|
|
|
|
|
} |
|
774
|
|
|
|
|
|
|
} |
|
775
|
0
|
0
|
|
|
|
0
|
if (@ts>1) { |
|
776
|
0
|
|
|
|
|
0
|
$read="'".join(',',@ts)."'"; |
|
777
|
|
|
|
|
|
|
} else { |
|
778
|
0
|
|
|
|
|
0
|
($read)=@ts; |
|
779
|
|
|
|
|
|
|
} |
|
780
|
|
|
|
|
|
|
} |
|
781
|
|
|
|
|
|
|
|
|
782
|
0
|
|
|
|
|
0
|
return $read; |
|
783
|
|
|
|
|
|
|
} |
|
784
|
|
|
|
|
|
|
|
|
785
|
|
|
|
|
|
|
=head2 set_parameters() |
|
786
|
|
|
|
|
|
|
|
|
787
|
|
|
|
|
|
|
Title : set_parameters |
|
788
|
|
|
|
|
|
|
Usage : $bowtiefac->set_parameters(%params); |
|
789
|
|
|
|
|
|
|
Function: sets the parameters listed in the hash or array, |
|
790
|
|
|
|
|
|
|
maintaining sane options. |
|
791
|
|
|
|
|
|
|
Returns : true on success |
|
792
|
|
|
|
|
|
|
Args : [optional] hash or array of parameter/values. |
|
793
|
|
|
|
|
|
|
Note : This will unset conflicts and set required options, |
|
794
|
|
|
|
|
|
|
but will not prevent non-sane requests in the arguments |
|
795
|
|
|
|
|
|
|
|
|
796
|
|
|
|
|
|
|
=cut |
|
797
|
|
|
|
|
|
|
|
|
798
|
|
|
|
|
|
|
sub set_parameters { |
|
799
|
8
|
|
|
8
|
1
|
864
|
my ($self, @args) = @_; |
|
800
|
|
|
|
|
|
|
|
|
801
|
|
|
|
|
|
|
# Mutually exclusive switches/params prevented from being set to |
|
802
|
|
|
|
|
|
|
# avoid confusion resulting from setting incompatible switches. |
|
803
|
|
|
|
|
|
|
|
|
804
|
8
|
50
|
|
|
|
22
|
$self->throw("Input args not an even number") if (@args % 2); |
|
805
|
8
|
|
|
|
|
50
|
my %args = @args; |
|
806
|
|
|
|
|
|
|
|
|
807
|
|
|
|
|
|
|
|
|
808
|
8
|
|
|
|
|
26
|
foreach (keys %args) { |
|
809
|
88
|
|
|
|
|
50
|
my @added; |
|
810
|
|
|
|
|
|
|
my @removed; |
|
811
|
88
|
|
|
|
|
108
|
s/^-//; |
|
812
|
88
|
|
|
|
|
56
|
foreach my $conflict (@{$incompat_params{$_}}) { |
|
|
88
|
|
|
|
|
118
|
|
|
813
|
66
|
50
|
|
|
|
78
|
return if grep /$conflict/, @added; |
|
814
|
66
|
|
|
|
|
65
|
delete $args{'-'.$conflict}; |
|
815
|
66
|
100
|
|
|
|
92
|
$args{'-'.$conflict} = undef if $self->{'_'.$conflict}; |
|
816
|
66
|
|
|
|
|
62
|
push @removed, $conflict; |
|
817
|
|
|
|
|
|
|
} |
|
818
|
88
|
|
|
|
|
47
|
foreach my $requirement (@{$corequisite_switches{$_}}) { |
|
|
88
|
|
|
|
|
132
|
|
|
819
|
12
|
50
|
|
|
|
29
|
return if grep /$requirement/, @removed; |
|
820
|
12
|
50
|
|
|
|
17
|
$args{'-'.$requirement}=1 if $args{$_}; |
|
821
|
12
|
|
|
|
|
17
|
push @added, $requirement; |
|
822
|
|
|
|
|
|
|
} |
|
823
|
|
|
|
|
|
|
} |
|
824
|
|
|
|
|
|
|
|
|
825
|
8
|
|
|
|
|
50
|
return $self->SUPER::set_parameters(%args); |
|
826
|
|
|
|
|
|
|
} |
|
827
|
|
|
|
|
|
|
|
|
828
|
|
|
|
|
|
|
=head2 version() |
|
829
|
|
|
|
|
|
|
|
|
830
|
|
|
|
|
|
|
Title : version |
|
831
|
|
|
|
|
|
|
Usage : $version = $bowtiefac->version() |
|
832
|
|
|
|
|
|
|
Function: Returns the program version (if available) |
|
833
|
|
|
|
|
|
|
Returns : string representing location and version of the program |
|
834
|
|
|
|
|
|
|
|
|
835
|
|
|
|
|
|
|
=cut |
|
836
|
|
|
|
|
|
|
|
|
837
|
|
|
|
|
|
|
sub version{ |
|
838
|
0
|
|
|
0
|
1
|
|
my ($self) = @_; |
|
839
|
|
|
|
|
|
|
|
|
840
|
0
|
0
|
|
|
|
|
my $cmd = $self->command if $self->can('command'); |
|
841
|
|
|
|
|
|
|
|
|
842
|
0
|
0
|
|
|
|
|
defined $cmd or $self->throw("No command defined - cannot determine program executable"); |
|
843
|
|
|
|
|
|
|
|
|
844
|
0
|
|
|
|
|
|
my ($in, $out, $err); |
|
845
|
0
|
|
|
|
|
|
my $dum; |
|
846
|
0
|
|
|
|
|
|
$in = \$dum; |
|
847
|
0
|
|
|
|
|
|
$out = \$self->{'stdout'}; |
|
848
|
0
|
|
|
|
|
|
$err = \$self->{'stderr'}; |
|
849
|
|
|
|
|
|
|
|
|
850
|
|
|
|
|
|
|
# Get program executable |
|
851
|
0
|
|
|
|
|
|
my $exe = $self->executable; |
|
852
|
|
|
|
|
|
|
# Get version switch from switches, translate and dash it |
|
853
|
0
|
|
|
|
|
|
my $version_switch = $param_translation{"$command_prefixes{$cmd}|version"}; |
|
854
|
0
|
|
|
|
|
|
$version_switch = $self->_dash_switch( $version_switch ); |
|
855
|
|
|
|
|
|
|
|
|
856
|
0
|
|
|
|
|
|
my @ipc_args = ( $exe, $version_switch ); |
|
857
|
|
|
|
|
|
|
|
|
858
|
0
|
|
|
|
|
|
eval { |
|
859
|
0
|
0
|
|
|
|
|
IPC::Run::run(\@ipc_args, $in, $out, $err) or |
|
860
|
|
|
|
|
|
|
die ("There was a problem running $exe : $!"); |
|
861
|
|
|
|
|
|
|
}; |
|
862
|
0
|
0
|
|
|
|
|
if ($@) { |
|
863
|
0
|
|
|
|
|
|
$self->throw("$exe call crashed: $@"); |
|
864
|
|
|
|
|
|
|
} |
|
865
|
|
|
|
|
|
|
|
|
866
|
0
|
|
|
|
|
|
my @details = split("\n",$self->stdout); |
|
867
|
0
|
|
|
|
|
|
(my $version) = grep /$exe version [[:graph:]]*$/, @details; |
|
868
|
0
|
|
|
|
|
|
$version =~ s/version //; |
|
869
|
0
|
|
|
|
|
|
(my $addressing) = grep /-bit$/, @details; |
|
870
|
|
|
|
|
|
|
|
|
871
|
0
|
|
|
|
|
|
return $version.' '.$addressing; |
|
872
|
|
|
|
|
|
|
} |
|
873
|
|
|
|
|
|
|
|
|
874
|
0
|
|
|
0
|
0
|
|
sub available_commands { shift->available_parameters('commands') }; |
|
875
|
|
|
|
|
|
|
|
|
876
|
0
|
|
|
0
|
0
|
|
sub filespec { shift->available_parameters('filespec') }; |
|
877
|
|
|
|
|
|
|
|
|
878
|
|
|
|
|
|
|
1; |