#!/usr/bin/env perl # Copyright (C) 2007-2009 Martin A. Hansen. # This program is free software; you can redistribute it and/or # modify it under the terms of the GNU General Public License # as published by the Free Software Foundation; either version 2 # of the License, or (at your option) any later version. # This program is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # GNU General Public License for more details. # You should have received a copy of the GNU General Public License # along with this program; if not, write to the Free Software # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. # http://www.gnu.org/copyleft/gpl.html # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< # Use bowtie to map sequences in the stream against a specified genome. # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< use warnings; use strict; use Maasha::Biopieces; use Maasha::Common; use Maasha::Fastq; use Maasha::Fasta; use Maasha::Calc; # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< my ( $options, $in, $out, $index, $tmp_dir, $tmp_in, $tmp_out, $fh_in, $fh_out, $record, $entry, $line, @fields, $type, $args ); $options = Maasha::Biopieces::parse_options( [ { long => 'genome', short => 'g', type => 'genome', mandatory => 'yes', default => undef, allowed => undef, disallowed => undef }, ] ); $in = Maasha::Biopieces::read_stream( $options->{ "stream_in" } ); $out = Maasha::Biopieces::write_stream( $options->{ "stream_out" } ); $index = "$ENV{ 'BP_DATA' }/genomes/$options->{ 'genome' }/bowtie/$options->{ 'genome' }"; $tmp_dir = Maasha::Biopieces::get_tmpdir(); $tmp_in = "$tmp_dir/bowtie.seq"; $tmp_out = "$tmp_dir/bowtie.out"; $fh_out = Maasha::Filesys::file_write_open( $tmp_in ); while ( $record = Maasha::Biopieces::get_record( $in ) ) { if ( $entry = Maasha::Fastq::biopiece2fastq( $record ) ) { Maasha::Common::error( "Mixed FASTA and FASTQ entries in stream" ) if defined $type and $type ne "FASTQ"; Maasha::Fastq::put_entry( $entry, $fh_out ); $type = "FASTQ"; } elsif ( $entry = Maasha::Fasta::biopiece2fasta( $record ) ) { Maasha::Common::error( "Mixed FASTA and FASTQ entries in stream" ) if defined $type and $type ne "FASTA"; Maasha::Fasta::put_entry( $entry, $fh_out ); $type = "FASTA"; } Maasha::Biopieces::put_record( $record, $out ); } close $fh_out; $args = ""; $args = "-f" if $type eq "FASTA"; if ( $options->{ 'verbose' } ) { Maasha::Common::run( "bowtie", "$args $index $tmp_in $tmp_out" ); } else { Maasha::Common::run( "bowtie", "$args $index $tmp_in $tmp_out > /dev/null 2>&1" ); } unlink $tmp_in; $fh_in = Maasha::Filesys::file_read_open( $tmp_out ); while ( $line = <$fh_in> ) { chomp $line; @fields = split /\t/, $line; $record = bowtie2biopiece( \@fields ); Maasha::Biopieces::put_record( $record, $out ); } close $fh_out; unlink $tmp_out; Maasha::Biopieces::close_stream( $in ); Maasha::Biopieces::close_stream( $out ); # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SUBROUTINES <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< sub bowtie2biopiece { # Martin A. Hansen, July 2009 # Convert a bowtie entry to a Biopiece record. my ( $entry, # bowtie entry ) = @_; # Returns a hash. my ( $record, @scores ); $record->{ 'Q_ID' } = $entry->[ 0 ]; $record->{ 'STRAND' } = $entry->[ 1 ]; $record->{ 'S_ID' } = $entry->[ 2 ]; $record->{ 'S_BEG' } = $entry->[ 3 ]; $record->{ 'SEQ' } = $entry->[ 4 ]; $record->{ 'SCORES' } = $entry->[ 5 ]; $record->{ 'MISMATCH' } = $entry->[ 6 ]; $record->{ 'SEQ_LEN' } = length $entry->[ 4 ]; $record->{ 'S_END' } = $record->{ 'S_BEG' } + $record->{ 'SEQ_LEN' } - 1; $record->{ 'SCORES' } =~ s/(.)/ord( $1 ) - 33 . ";"/ge; # http://maq.sourceforge.net/fastq.shtml $record->{ 'SCORE_MEAN' } = sprintf( "%.2f", Maasha::Calc::mean( [ split /;/, $record->{ 'SCORES' } ] ) ); $record->{ 'REC_TYPE' } = "BOWTIE"; return wantarray ? %{ $record } : $record; } # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< BEGIN { Maasha::Biopieces::status_set(); } END { Maasha::Biopieces::status_log(); } # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< __END__