bp_bin/read_454

   1 #!/usr/bin/env perl
   2
   3 # Copyright (C) 2007-2010 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24 # Read 454 entries from fasta and quality file.
  25
  26 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  27
  28
  29 use warnings;
  30 use strict;
  31 use Data::Dumper;
  32 use Maasha::Biopieces;
  33 use Maasha::Filesys;
  34 use Maasha::Fasta;
  35 use Maasha::Fastq;
  36
  37
  38 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  39
  40
  41 my ( $options, $in, $out, $record, $data_in, $qual_in, $num, $fasta, $qual, @seqs, @scores, $i );
  42
  43 $options = Maasha::Biopieces::parse_options(
  44     [
  45         { long => 'data_in', short => 'i', type => 'file!',  mandatory => 'yes', default => undef, allowed => undef, disallowed => undef },
  46         { long => 'qual_in', short => 'q', type => 'file!',  mandatory => 'yes', default => undef, allowed => undef, disallowed => undef },
  47         { long => 'num',     short => 'n', type => 'uint',   mandatory => 'no',  default => undef, allowed => undef, disallowed => '0' },
  48     ]
  49 );
  50
  51 $in  = Maasha::Biopieces::read_stream( $options->{ "stream_in" } );
  52 $out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } );
  53
  54 while ( $record = Maasha::Biopieces::get_record( $in ) ) {
  55     Maasha::Biopieces::put_record( $record, $out );
  56 }
  57
  58 if ( $options->{ 'data_in' } )
  59 {
  60     $data_in = Maasha::Filesys::file_read_open( $options->{ 'data_in' } );
  61     $qual_in = Maasha::Filesys::file_read_open( $options->{ 'qual_in' } );
  62
  63     $num = 1;
  64
  65     while ( $fasta = Maasha::Fasta::get_entry( $data_in ) )
  66     {
  67         $qual = get_qual( $qual_in );
  68
  69         check_names( $fasta, $qual );
  70         check_lengths( $fasta, $qual );
  71
  72         $record = {
  73             SEQ_NAME    => $fasta->[ 0 ],
  74             SEQ         => $fasta->[ 1 ],
  75             SEQ_LEN     => length $fasta->[ 1 ],
  76             SCORES      => $qual->[ 1 ],
  77         };
  78
  79         Maasha::Biopieces::put_record( $record, $out );
  80
  81         last if $options->{ "num" } and $num == $options->{ "num" };
  82
  83         $num++;
  84     }
  85
  86     close $data_in;
  87     close $qual_in;
  88 }
  89
  90 Maasha::Biopieces::close_stream( $in );
  91 Maasha::Biopieces::close_stream( $out );
  92
  93
  94 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  95
  96
  97 sub get_qual
  98 {
  99     # Martin A. Hansen, May 2010.
 100
 101     # Get the next 454 quality entry from an open file.
 102
 103     my ( $fh,   # file handle
 104        ) = @_;
 105
 106     # Returns list.
 107
 108     my ( $block, $name, $qual, $scores, $entry );
 109
 110     local $/ = "\n>";
 111
 112     while ( $block = <$fh> )
 113     {
 114         chomp $block;
 115
 116         last if $block !~ /^\s+$/;
 117     }
 118
 119     return if not defined $block;
 120
 121     $block =~ /^>?(.+)\n/m;
 122     $name = $1;
 123     $qual = $';
 124
 125     local $/ = "\n";
 126
 127     chomp $qual;
 128
 129     $name =~ tr/\r//d;
 130     $qual =~ tr/ \n\r/;;;/;
 131     $qual =~ s/;;/;/g;
 132
 133     $scores = Maasha::Fastq::dec_str2phred_str( $qual );
 134
 135     $entry = [ $name, $scores ];
 136
 137     return wantarray ? @{ $entry } : $entry;
 138 }
 139
 140
 141 sub check_names
 142 {
 143     # Martin A. Hansen, May 2010.
 144
 145     # Check if the names of the fasta and qual entries match
 146     # and raise an error if they don't.
 147
 148     my ( $fasta,    # fasta entry
 149          $qual,     # qual entry
 150        ) = @_;
 151
 152     # Returns nothing.
 153
 154     my ( $f_name, $q_name );
 155
 156     if ( $fasta->[ 0 ] =~ /^([^ ]+)/ ) {
 157         $f_name = $1;
 158     }
 159
 160     if ( $qual->[ 0 ] =~ /^([^ ]+)/ ) {
 161         $q_name = $1;
 162     }
 163
 164     Maasha::Common::error( qq(names don't match "$fasta->[ 0 ]" ne "$qual->[ 0 ]") ) if $f_name ne $q_name;
 165 }
 166
 167
 168 sub check_lengths
 169 {
 170     # Martin A. Hansen, April 2011.
 171
 172     # Check if the lengths of the fasta and qual strings are the same
 173     # and raise an error if not.
 174
 175     my ( $fasta,    # fasta entry
 176          $qual,     # qual entry
 177        ) = @_;
 178
 179     # Returns nothing.
 180
 181     my ( $f_len, $q_len );
 182
 183     $f_len = length $fasta->[ 1 ];
 184     $q_len = length $qual->[ 1 ];
 185
 186     Maasha::Common::error( qq(lengths don't match "$f_len" != "$q_len") ) if $f_len != $q_len;
 187 }
 188
 189
 190 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 191
 192
 193 BEGIN
 194 {
 195     Maasha::Biopieces::status_set();
 196 }
 197
 198
 199 END
 200 {
 201     Maasha::Biopieces::status_log();
 202 }
 203
 204
 205 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 206
 207
 208 __END__