File Coverage

blib/lib/FASTX/Reader.pm
Criterion Covered Total %
statement 198 234 84.6
branch 74 110 66.3
condition 26 45 57.7
subroutine 14 14 100.0
pod 5 5 100.0
total 317 408 77.4


line stmt bran cond sub pod time code
1             package FASTX::Reader;
2 25     25   1670674 use 5.012;
  25         321  
3 25     25   139 use warnings;
  25         45  
  25         846  
4 25     25   141 use Carp qw(confess);
  25         49  
  25         1344  
5 25     25   10265 use Data::Dumper;
  25         110285  
  25         2344  
6             $FASTX::Reader::VERSION = '0.90';
7             require Exporter;
8             our @ISA = qw(Exporter);
9             #ABSTRACT: A lightweight module to parse FASTA and FASTQ files, supporting compressed files and paired-ends.
10              
11 25     25   176 use constant GZIP_SIGNATURE => pack('C3', 0x1f, 0x8b, 0x08);
  25         48  
  25         57581  
12              
13              
14             sub new {
15             # Instantiate object
16 33     33 1 5114 my ($class, $args) = @_;
17              
18 33 100       156 if (defined $args->{loadseqs}) {
19 2 100 66     19 if ($args->{loadseqs} eq 'name' or $args->{loadseqs} eq 'names' ) {
    50 33        
      33        
20 1         3 $args->{loadseqs} = 'name';
21             } elsif ($args->{loadseqs} eq 'seq' or $args->{loadseqs} eq 'seqs' or $args->{loadseqs} == 1) {
22 1         3 $args->{loadseqs} = 'seq';
23             } else {
24 0         0 confess("attribute should be 'name' or 'seq' to specify the key of the hash.");
25             }
26             }
27             my $self = {
28             filename => $args->{filename},
29             loadseqs => $args->{loadseqs},
30 33         186 };
31              
32              
33 33         120 my $object = bless $self, $class;
34              
35             # Initialize auxiliary array for getRead
36 33         231 $object->{aux} = [undef];
37 33         137 $object->{compressed} = 0;
38              
39             # Check if a filename was provided and not {{STDIN}}
40             # uncoverable branch false
41              
42 33 100 66     299 if (defined $self->{filename} and $self->{filename} ne '{{STDIN}}') {
43 31 100       2803 open my $fh, '<', $self->{filename} or confess "Unable to read file ", $self->{filename}, "\n";
44 27         534 read( $fh, my $magic_byte, 4 );
45 27         364 close $fh;
46              
47             # See: __BioX::Seq::Stream__ for GZIP (and other) compressed file reader
48 27 100       619 if (substr($magic_byte,0,3) eq GZIP_SIGNATURE) {
    50          
49 16         45 $object->{compressed} = 1;
50 16         120 our $GZIP_BIN = _which('pigz', 'gzip');
51 16         5855 close $fh;
52 16 50       58 if (! defined $GZIP_BIN) {
53 0         0 require IO::Uncompress::Gunzip;
54 0         0 $fh = IO::Uncompress::Gunzip->new($self->{filename}, MultiStream => 1);
55             } else {
56 16 50       33114 open $fh, '-|', "$GZIP_BIN -dc $self->{filename}" or confess "Error opening gzip file ", $self->{filename}, ": $!\n";
57             }
58             } elsif (-B $self->{filename}) {
59             # BINARY FILE NOT SUPPORTED?
60 0         0 close $fh;
61 0         0 $self->{fh} = undef;
62 0         0 $self->{status} = 1;
63 0         0 $self->{message} = 'Binary file not supported';
64             } else {
65 11         39 close $fh;
66 10 50   10   64 open $fh, '<:encoding(utf8)', $self->{filename} or confess "Unable to read file ", $self->{filename}, ": ", $!, "\n";
  10         20  
  10         65  
  11         312  
67             }
68 27         112967 $object->{fh} = $fh;
69             } else {
70 2         6 $self->{fh} = \*STDIN;
71 2 50       7 if ($self->{loadseqs}) {
72 0         0 confess("Load sequences not supported for STDIN");
73             }
74             }
75              
76 29 100       346 if ($self->{loadseqs}) {
77 2         5 _load_seqs($self);
78              
79             }
80              
81 29         930 return $object;
82             }
83              
84              
85              
86             sub getRead {
87 37     37 1 4120 my $self = shift;
88             #my ($fh, $aux) = @_;
89             #@::::::: :::
90 37         118 my $fh = $self->{fh};
91              
92 37 50 33     247 return undef if (defined $self->{status} and $self->{status} == 0);
93              
94 37         134 my $aux = $self->{aux};
95 37         89 my $return;
96 37 50       214 @$aux = [undef, 0] if (!(@$aux)); # remove deprecated 'defined'
97 37 100       294 return if ($aux->[1]);
98             # uncoverable branch true
99 34 100       143 if (!defined($aux->[0])) {
100 28         2043 while (<$fh>) {
101 28         170 chomp;
102 28 100 100     351 if (substr($_, 0, 1) eq '>' || substr($_, 0, 1) eq '@') {
103 27         118 $aux->[0] = $_;
104 27         71 last;
105             }
106             }
107 28 100       136 if (!defined($aux->[0])) {
108 1         4 $aux->[1] = 1;
109 1         3 return;
110             }
111             }
112              
113             #my $comm = /^.\S+\s+(.*)/? $1 : ''; # retain "comment"
114             # Comments can have more spaces: xx
115 33 50       670 my ($name, $comm) = /^.(\S+)(?:\s+)(.+)/ ? ($1, $2) :
    100          
116             /^.(\S+)/ ? ($1, '') : ('', '');
117              
118 33         106 my $seq = '';
119 33         59 my $c;
120 33         63 $aux->[0] = undef;
121 33         188 while (<$fh>) {
122 83         158 chomp;
123 83         195 $c = substr($_, 0, 1);
124 83 100 66     493 last if ($c eq '>' || $c eq '@' || $c eq '+');
      100        
125 53         221 $seq .= $_;
126             }
127 33         87 $aux->[0] = $_;
128 33 100       126 $aux->[1] = 1 if (!defined($aux->[0]));
129 33         155 $return->{name} = $name;
130 33         157 $return->{comment} = $comm;
131 33         148 $return->{seq} = $seq;
132 33         96 $self->{counter}++;
133 33 100       129 return $return if ($c ne '+');
134 22         51 my $qual = '';
135 22         118 while (<$fh>) {
136 22         49 chomp;
137 22         68 $qual .= $_;
138 22 50       84 if (length($qual) >= length($seq)) {
139 22         42 $aux->[0] = undef;
140 22         64 $return->{name} = $name;
141 22         51 $return->{seq} = $seq;
142 22         52 $return->{comment} = $comm;
143 22         81 $return->{qual} = $qual;
144             #$self->{counter}+=100;
145 22         120 return $return;
146             }
147             }
148 0         0 $aux->[1] = 1;
149 0         0 $return->{name} = $name;
150 0         0 $return->{seq} = $seq;
151 0         0 $return->{comment} = $comm;
152 0         0 $self->{counter}++;
153 0         0 return $return;
154              
155             }
156              
157              
158             sub getFastqRead {
159 7     7 1 2887 my $self = shift;
160 7         22 my $seq_object = undef;
161              
162 7 50 66     56 return undef if (defined $self->{status} and $self->{status} == 0);
163              
164 7         24 $self->{status} = 1;
165 7         495 my $header = readline($self->{fh});
166 7         72 my $seq = readline($self->{fh});
167 7         27 my $check = readline($self->{fh});
168 7         24 my $qual = readline($self->{fh});
169              
170              
171             # Check 4 lines were found (FASTQ)
172              
173 7 100       30 unless (defined $qual) {
174 1 50       5 if (defined $header) {
175 0         0 $self->{message} = "Unknown format: FASTQ truncated at " . $header . "?";
176 0         0 $self->{status} = 0;
177             }
178 1         4 return undef;
179             }
180              
181             # Fast format control: header and separator
182 6 100 66     74 if ( (substr($header, 0, 1) eq '@') and (substr($check, 0, 1) eq '+') ) {
183 5         17 chomp($header);
184 5         13 chomp($seq);
185 5         9 chomp($qual);
186             # Also control sequence integrity
187 5 100 66     76 if ($seq=~/^[ACGTNacgtn]+$/ and length($seq) == length($qual) ) {
188 4         32 my ($name, $comments) = split /\s+/, substr($header, 1);
189 4         18 $seq_object->{name} = $name;
190 4         15 $seq_object->{comments} = $comments;
191 4         15 $seq_object->{seq} = $seq;
192 4         15 $seq_object->{qual} = $qual;
193 4         16 $self->{counter}++;
194              
195             } else {
196             # Return error (corrupted FASTQ)
197 1         3 $self->{message} = "Unknown format: expecting FASTQ (corrupted?)";
198 1         2 $self->{status} = 0;
199              
200             }
201             } else {
202             # Return error (not a FASTQ)
203 1         3 $self->{message} = "Unknown format: expecting FASTQ but @ header not found";
204              
205 1 50       4 if (substr($header, 0,1 ) eq '>' ) {
206             # HINT: is FASTA?
207 1         4 $self->{message} .= " (might be FASTA instead)";
208             }
209 1         3 $self->{status} = 0;
210             }
211              
212 6         26 return $seq_object;
213             }
214              
215              
216             sub getIlluminaRead {
217 16     16 1 9191 my $self = shift;
218 16         32 my $seq_object = undef;
219              
220 16 50 66     90 return undef if (defined $self->{status} and $self->{status} == 0);
221              
222 16         35 $self->{status} = 1;
223 16         587 my $header = readline($self->{fh});
224 16         101 my $seq = readline($self->{fh});
225 16         41 my $check = readline($self->{fh});
226 16         54 my $qual = readline($self->{fh});
227              
228              
229             # Check 4 lines were found (FASTQ)
230              
231 16 100       47 unless (defined $qual) {
232 2 50       8 if (defined $header) {
233 0         0 $self->{message} = "Unknown format: FASTQ truncated at " . $header . "?";
234 0         0 $self->{status} = 0;
235             }
236 2         9 return undef;
237             }
238              
239             # Fast format control: header and separator
240 14 50 33     103 if ( (substr($header, 0, 1) eq '@') and (substr($check, 0, 1) eq '+') ) {
241 14         28 chomp($header);
242 14         21 chomp($seq);
243 14         24 chomp($qual);
244             # Also control sequence integrity
245 14 50 33     152 if ($seq=~/^[ACGTNacgtn]+$/ and length($seq) == length($qual) ) {
246 14         83 my ($name, $comments) = split /\s+/, substr($header, 1);
247             #@::::::: :::
248 14         65 my ($instrument, $run, $flowcell, $lane, $tile, $x, $y, $umi) = split /:/, $name;
249 14         31 my ($read, $filtered, $ctrl, $index1, $index2);
250              
251 14 50       29 if (not defined $y) {
252 0         0 $self->{message} = "Unknown format: not Illumina naming: ::::::";
253 0         0 $self->{status} = 0;
254             }
255              
256 14         38 $seq_object->{name} = $name;
257 14         31 $seq_object->{comments} = $comments;
258 14         29 $seq_object->{seq} = $seq;
259 14         28 $seq_object->{qual} = $qual;
260              
261 14         33 $seq_object->{instrument} = $instrument;
262 14         34 $seq_object->{run} = $run;
263 14         24 $seq_object->{flowcell} = $flowcell;
264 14         36 $seq_object->{lane} = $lane;
265 14         29 $seq_object->{tile} = $tile;
266 14         28 $seq_object->{x} = $x;
267 14         23 $seq_object->{y} = $y;
268 14         25 $seq_object->{umi} = $umi;
269 14         20 $seq_object->{instrument} = $instrument;
270 14         20 $seq_object->{run} = $run;
271 14         18 $seq_object->{flowcell} = $flowcell;
272              
273 14 100       29 if (defined $comments) {
274 7         41 ($read, $filtered, $ctrl, $index1, $index2) = split /[:+]/, $comments;
275 7 50       16 if ( defined $ctrl ) {
276 7         19 $seq_object->{read} = $read;
277 7 50       17 $filtered eq 'N' ? $seq_object->{filtered} = 0 : $seq_object->{filtered} = 1;
278 7         10 $seq_object->{control} = $ctrl;
279 7 50       67 if ($read eq '1') {
280 7         23 $seq_object->{index} = $index1;
281 7         20 $seq_object->{paired_index} = $index2;
282             } else {
283 0         0 $seq_object->{paired_index} = $index1;
284 0         0 $seq_object->{index} = $index2;
285             }
286             }
287             }
288 14         28 $self->{counter}++;
289              
290             } else {
291             # Return error (corrupted FASTQ)
292 0         0 $self->{message} = "Unknown format: expecting FASTQ (corrupted?)";
293 0         0 $self->{status} = 0;
294              
295             }
296             } else {
297             # Return error (not a FASTQ)
298 0         0 $self->{message} = "Unknown format: expecting FASTQ but @ header not found";
299              
300 0 0       0 if (substr($header, 0,1 ) eq '>' ) {
301             # HINT: is FASTA?
302 0         0 $self->{message} .= " (might be FASTA instead)";
303             }
304 0         0 $self->{status} = 0;
305             }
306              
307 14         43 return $seq_object;
308             }
309              
310             sub getFileFormat {
311 4     4 1 1266 my $self = shift;
312 4         9 my ($filename) = shift;
313 4 50       11 return 0 if (not defined $filename);
314            
315 4 50       149 open my $fh, '<', $filename or confess "Unable to read file ", $filename, "\n";
316 4         69 read( $fh, my $magic_byte, 4 );
317 4         59 close $fh;
318              
319 4 50       19 if (substr($magic_byte,0,3) eq GZIP_SIGNATURE) {
320            
321 0 0       0 if (! defined $self->{GZIP_BIN}) {
322 0         0 require IO::Uncompress::Gunzip;
323 0         0 $fh = IO::Uncompress::Gunzip->new($filename, MultiStream => 1);
324             } else {
325 0 0       0 open $fh, '-|', "$self->{GZIP_BIN} -dc $filename" or confess "Error opening gzip file ", $filename, ": $!\n";
326             }
327             } else {
328 4   33 1   113 open $fh, '<:encoding(utf8)', "$filename" || confess "Unable to read $filename\n$!\n";
  1         6  
  1         2  
  1         5  
329             }
330 4         11257 my $first = readline($fh);
331 4 100       55 if (substr($first, 0,1) eq '>') {
    100          
332             #should be FASTA
333 1         24 return 'fasta';
334             } elsif (substr($first, 0, 1) eq '@') {
335             #should be fastq
336 2         7 readline($fh);
337 2         5 my $sep = readline($fh);
338 2 50       9 if ( substr($sep, 0, 1) eq '+' ) {
339             #second check for fastq
340 2         36 return 'fastq';
341             }
342             } else {
343             #unknown format
344 1         33 return undef;
345             }
346             }
347             sub _load_seqs {
348 2     2   5 my ($self) = @_;
349 2 50       7 return 0 unless (defined $self->{loadseqs});
350              
351 2         3 my $seqs = undef;
352 2         7 while (my $s = $self->getRead() ) {
353 6         15 my ($name, $seq) = ($s->{name}, $s->{seq});
354 6 100       11 if ($self->{loadseqs} eq 'name') {
355 3         15 $seqs->{$name} = $seq;
356             } else {
357 3         15 $seqs->{$seq} = $name;
358             }
359              
360             }
361 2         5 $self->{seqs} = $seqs;
362             }
363              
364              
365             sub _which {
366 16 50   16   136 return undef if ($^O eq 'MSWin32');
367 16         66 my $has_which = eval { require File::Which; File::Which->import(); 1};
  16         4659  
  16         11053  
  16         52  
368 16 50       102 if ($has_which) {
369 16         70 foreach my $cmd (@_) {
370 32 100       4123 return which($cmd) if (which($cmd));
371             }
372             } else {
373 0         0 foreach my $cmd (@_) {
374 0         0 `which $cmd 2> /dev/null`;
375 0 0       0 return $cmd if (not $?);
376             }
377             }
378 0         0 return undef;
379             }
380             1;
381              
382             __END__