# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-# Use bowtie to map sequences in the stream against a specified genome.
+# Use bowtie to map sequences in the stream against a specified genome or index.
# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
use warnings;
use strict;
+use Data::Dumper;
use Maasha::Biopieces;
use Maasha::Common;
use Maasha::Fastq;
use Maasha::Fasta;
use Maasha::Calc;
+use constant {
+ SEQ_NAME => 0,
+ SEQ => 1,
+ SCORES => 2,
+};
# >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
]
);
-Maasha::Common::error( qq(both --index_hame and --genome specified) ) if $options->{ "genome" } and $options->{ "index_name" };
+Maasha::Common::error( qq(both --index_name and --genome specified) ) if $options->{ "genome" } and $options->{ "index_name" };
Maasha::Common::error( qq(no --index_name or --genome specified) ) if not $options->{ "genome" } and not $options->{ "index_name" };
$in = Maasha::Biopieces::read_stream( $options->{ "stream_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::Common::error( "Sequence longer than 1024 not allowed") if length( $entry->[ SEQ ] ) > 1024;
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::Common::error( "Sequence longer than 1024 not allowed") if length( $entry->[ SEQ ] ) > 1024;
Maasha::Fasta::put_entry( $entry, $fh_out );
$type = "FASTA";
close $fh_out;
push @args, "-n $options->{ 'mismatches' }";
+push @args, "-v $options->{ 'mismatches' }"; # DANGER: using seed mismatches as alignment mismatches - may work, may not!
push @args, "-f" if $type eq "FASTA";
+push @args, "-p $options->{ 'cpus' }";
+push @args, "--phred64-quals" unless $type eq "FASTA";
if ( defined $options->{ 'max_hits' } ) {
push @args, "-k $options->{ 'max_hits' }";
chomp $line;
@fields = split /\t/, $line;
-
$record = bowtie2biopiece( \@fields );
Maasha::Biopieces::put_record( $record, $out );
# 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 ];
-
+ my ( $record, $s_id, $s_len, $hits );
+
+ $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->{ 'SCORE' } = $entry->[ 6 ] + 1;
+ $record->{ 'ALIGN' } = $entry->[ 7 ] || '.';
+ $record->{ 'S_LEN' } = length $entry->[ 4 ];
$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->{ 'SCORES' } =~ s/(.)/chr( ( ord( $1 ) - 33 ) + 64 )/ge; # convert phred scores to illumina scores
+ $record->{ 'HITS' } = '.';
$record->{ 'REC_TYPE' } = "BOWTIE";
return wantarray ? %{ $record } : $record;