Bio::Tools::Run::BWA - Run wrapper for the BWA short-read assembler *BETA*

 # create an assembly

 # run BWA commands separately

 $bwafac = Bio::Tools::Run::BWA->new( -command => 'aln' );
 # all maq commands
 @all_commands = $bwafac->available_parameters('commands'); 
 @all_commands = $bwafac->available_commands; # alias
 # just for aln
 @aln_params = $bwafac->available_parameters('params');
 @aln_switches = $bwafac->available_parameters('switches');
 @aln_all_options = $bwafac->available_parameters();

  $bwafac = Bio::Tools::Run::BWA->new( -command => 'aln' );
  @filespec = $bwafac->filespec;

 fas
 faq 
 >sai

 $bwafac->run_bwa( -fas => 'my.db.fas', -faq => 'reads.faq',
                   -sai => 'out.sai' );

  bioperl-l@bioperl.org                  - General discussion
http://bioperl.org/wiki/Mailing_lists  - About the mailing lists

  http://redmine.open-bio.org/projects/bioperl/

 Email maj -at- fortinbras -dot- us

 Title   : new
 Usage   : my $obj = new Bio::Tools::Run::BWA();
 Function: Builds a new Bio::Tools::Run::BWA object
 Returns : an instance of Bio::Tools::Run::BWA
 Args    :

 Title   : run
 Usage   : $assembly = $bwafac->run( @args );
 Function: Run the bwa assembly pipeline. 
 Returns : Assembly results (file, IO object or Assembly object)
 Args    : - fastq file containing single-end reads
           - fasta file containing the reference sequence
           - [optional] fastq file containing paired-end reads 

 Title   : run_bwa
 Usage   : $obj->run_bwa( @file_args )
 Function: Run a bwa command as specified during object contruction
 Returns : 
 Args    : a specification of the files to operate on:

 Title   : stdout
 Usage   : $fac->stdout()
 Function: store the output from STDOUT for the run, 
           if no file specified in run_maq()
 Example : 
 Returns : scalar string
 Args    : on set, new value (a scalar or undef, optional)

 Title   : stderr
 Usage   : $fac->stderr()
 Function: store the output from STDERR for the run, 
           if no file is specified in run_maq()
 Example : 
 Returns : scalar string
 Args    : on set, new value (a scalar or undef, optional)

 No-op.

 No-op.

 Convert input fastq and fasta to maq format.

 Title   : _collate_subcmd_args
 Usage   : $args_hash = $self->_collate_subcmd_args
 Function: collate parameters and switches into command-specific
           arg lists for passing to new()
 Returns : hash of named argument lists
 Args    : [optional] composite cmd prefix (scalar string) 
           [default is 'run']

 Title   :   _run
 Usage   :   $factory->_run()
 Function:   Run a bwa assembly pipeline
 Returns :   a text-formatted sam alignment
 Args    :   - single end read file in maq bfq format
             - reference seq file in maq bfa format
             - [optional] paired end read file in maq bfq format

 Title   : available_parameters
 Usage   : @cmds = $fac->available_commands('commands');
 Function: Use to browse available commands, params, or switches
 Returns : array of scalar strings
 Args    : 'commands' : all bwa commands
           'params'   : parameters for this object's command
           'switches' : boolean switches for this object's command
           'filespec' : the filename spec for this object's command
 4Geeks  : Overrides Bio::ParameterBaseI via 
           Bio::Tools::Run::AssemblerBase

# $Id$
package Bio::Tools::Run::BWA;
use strict;
#
# BioPerl module for Bio::Tools::Run::BWA
#
# Please direct questions and support issues to <bioperl-l@bioperl.org>
#
# Cared for by Mark A. Jensen <maj@fortinbras.us>
#
# Copyright Mark A. Jensen
#
# You may distribute this module under the same terms as perl itself

# POD documentation - main docs before the code

# Let the code begin...


use IPC::Run;
our $HAVE_IO_UNCOMPRESS;
our $HAVE_SAMTOOLS;
BEGIN {
    eval "require IO::Uncompress::Gunzip; \$HAVE_IO_UNCOMPRESS = 1";
    eval "require Bio::Tools::Run::Samtools; \$HAVE_SAMTOOLS = 1";
}

# Object preamble - inherits from Bio::Root::Root

use Bio::Root::Root;
use Bio::Tools::Run::BWA::Config;
use Bio::Tools::GuessSeqFormat;
use File::Basename qw(fileparse);
use File::Copy;
use Cwd;

use base qw(Bio::Root::Root Bio::Tools::Run::AssemblerBase );

our $program_name = 'bwa'; # name of the executable

# Note:
#  other globals required by Bio::Tools::Run::AssemblerBase are
#  imported from Bio::Tools::Run::BWA::Config

our $qual_param = 'quality_file';
our $use_dash = 1;
our $join = ' ';

# Bio::Assembly::IO::sam now workable!
our $asm_format = 'sam'; 

sub new {
  my ($class,@args) = @_;
  my $self = $class->SUPER::new(@args);
  $self->parameters_changed(1); 
  $self->_register_program_commands( \@program_commands, \%command_prefixes );
  unless (grep /command/, @args) {
      push @args, '-command', 'run';
  }
  $self->_set_program_options(\@args, \@program_params, \@program_switches,
    \%param_translation, $qual_param, $use_dash, $join);
  $self->program_name($program_name) if not defined $self->program_name();
  $self->parameters_changed(1); # set on instantiation, per Bio::ParameterBaseI
  $self->_assembly_format($asm_format);
  return $self;
}

sub run {
  my ($self, $rd1_file, $ref_file, $rd2_file) = @_;

  # Sanity checks
  $self->_check_executable();
  unless ($HAVE_SAMTOOLS) {
      cluck( "Bio::Tools::Run::Samtools is not available. A .sam output alignment will be created, but must be converted to binary SAM (.bam) before it can be passed to Bio::Assembly::IO, as follows: \n\t\$ samtools view -Sb out.sam > out.bam" );
  }
  $rd1_file or $self->throw("Fastq reads file required at arg 1");
  $ref_file or $self->throw("Fasta refseq file required at arg 2");
  for ($rd1_file, $ref_file, $rd2_file) {
      next unless $_;
      if (/\.gz[^.]*$/) {
	  unless ($HAVE_IO_UNCOMPRESS) {
	      croak( "IO::Uncompress::Gunzip not available, can't expand '$_'" );
	  }
	  my ($tfh, $tf) = $self->io->tempfile;
	  my $z = IO::Uncompress::Gunzip->new($_);
	  while (<$z>) { print $tfh $_ }
	  close $tfh;
	  $_ = $tf;
      }
  }

  my $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$rd1_file);

  $guesser->guess eq 'fastq' or $self->throw("Reads file doesn't look like fastq at arg 1");
  $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$ref_file);
  $guesser->guess eq 'fasta' or $self->throw("Refseq file doesn't look like fasta at arg 2");
  if ($rd2_file) {
      $guesser = Bio::Tools::GuessSeqFormat->new(-file=>$rd2_file);
      $guesser->guess eq 'fastq' or $self->throw("Reads file doesn't look like fastq at arg 3");
  }

  #Assemble
  my ($sam_file) = $self->_run($rd1_file, $ref_file, $rd2_file);

  if ($HAVE_SAMTOOLS) {
      my ($nm,$dr,$suf) = fileparse($sam_file, ".sam");
      # goofy kludge for samtools...
      my $pwd = getcwd;
      chdir($dr);
      my $samt = Bio::Tools::Run::Samtools->new( -command => 'view', 
						 -sam_input => 1,
						 -bam_output => 1,
	                                         -refseq => $ref_file);
      my $bam_file = $nm.'.bam';
      $samt->run( -bam => $nm.$suf, -out => $bam_file ) or croak( "Problem converting .sam file");
      $samt = Bio::Tools::Run::Samtools->new( -command => 'sort' );
      $samt->run( -bam => $bam_file, -pfx => $nm.'.srt' ) or croak( "Problem sorting .bam file");
      move( $nm.'.srt.bam', $bam_file );
      $samt = Bio::Tools::Run::Samtools->new( -command => 'index' );
      $samt->run( -bam => $bam_file );
      $bam_file = File::Spec->catfile($dr, $bam_file);
      $sam_file = $bam_file;
      chdir($pwd);
  }
  
  # Export results in desired object type
  my $asm = $self->_export_results($sam_file, -refdb => $ref_file, -keep_asm => 1);
  return $asm;

}

sub run_bwa {
    my ($self, @args) = @_;
    # _translate_params will provide an array of command/parameters/switches
    # -- these are set at object construction
    # to set up the run, need to add the files to the call
    # -- provide these as arguments to this function

    my $cmd = $self->command if $self->can('command');
    $self->throw("No maq command specified for the object") unless $cmd;
    # setup files necessary for this command
    my $filespec = $command_files{$cmd};
    $self->throw("No command-line file specification is defined for command '$cmd'; check Bio::Tools::Run::Maq::Config") unless $filespec;

    # parse args based on filespec
    # require named args
    $self->throw("Named args are required") unless !(@args % 2);
    s/^-// for @args;
    my %args = @args; 
    # validate
    my @req = map { 
	my $s = $_; 
	$s =~ s/^[012]?[<>]//;
	$s =~ s/[^a-zA-Z0-9_]//g; 
	$s
    } grep !/[#]/, @$filespec;
    !defined($args{$_}) && $self->throw("Required filearg '$_' not specified") for @req;
    # set up redirects
    my ($in, $out, $err);
    for (@$filespec) {
	m/^1?>(.*)/ && do {
	    defined($args{$1}) && ( open($out,">", $args{$1}) or $self->throw("Open for write error : $!"));
	    next;
	};
	m/^2>#?(.*)/ && do {
	    defined($args{$1}) && (open($err, ">", $args{$1}) or $self->throw("Open for write error : $!"));
	    next;
	};
	m/^<#?(.*)/ && do {
	    defined($args{$1}) && (open($in, "<", $args{$1}) or $self->throw("Open for read error : $!"));
	    next;
	}
    }
    my $dum;
    $in || ($in = \$dum);
    $out || ($out = \$self->{'stdout'});
    $err || ($err = \$self->{'stderr'});
    
    # Get program executable
    my $exe = $self->executable;
    # Get command-line options
    my $options = $self->_translate_params();
    # Get file specs sans redirects in correct order
    my @specs = map { 
	my $s = $_; 
	$s =~ s/[^a-zA-Z0-9_]//g; 
	$s
    } grep !/[<>]/, @$filespec;
    my @files = @args{@specs};
    # expand arrayrefs
    my $l = $#files;
    for (0..$l) {
	splice(@files, $_, 1, @{$files[$_]}) if (ref($files[$_]) eq 'ARRAY');
    }
    @files = map { defined $_ ? $_ : () } @files; # squish undefs
    my @ipc_args = ( $exe, @$options, @files );

    eval {
	IPC::Run::run(\@ipc_args, $in, $out, $err) or
	    die ("There was a problem running $exe : $!");
    };
    if ($@) {
	$self->throw("$exe call crashed: $@");
    }

    # return arguments as specified on call
    return @args;
}

sub stdout {
    my $self = shift;
    
    return $self->{'stdout'} = shift if @_;
    return $self->{'stdout'};
}

sub stderr {
    my $self = shift;
    
    return $self->{'stderr'} = shift if @_;
    return $self->{'stderr'};
}

sub _check_sequence_input {
    return 1;
}

sub _check_optional_quality_input {
    return 1;
}

sub _prepare_input_sequences {
    shift->throw_not_implemented;
}

sub _collate_subcmd_args {
    my $self = shift;
    my $cmd = shift;
    my %ret;
    # default command is 'run'
    $cmd ||= 'run';
    my @subcmds = @{$composite_commands{$cmd}};
    my %subcmds;
    my $cur_options = $self->{'_options'};

    # collate
    foreach my $subcmd (@subcmds) {
	# find the composite cmd form of the argument in 
	# the current params and switches
	# e.g., map_max_mismatches
	my @params = grep /^${subcmd}_/, @{$$cur_options{'_params'}};
	my @switches = grep /^${subcmd}_/, @{$$cur_options{'_switches'}};
	$ret{$subcmd} = [];
	# create an argument list suitable for passing to new() of
	# the subcommand factory...
	foreach my $opt (@params, @switches) {
	    my $subopt = $opt; 
	    $subopt =~ s/^${subcmd}_//; 
	    push(@{$ret{$subcmd}}, '-'.$subopt => $self->$opt) if defined $self->$opt;
	}
    }
    return \%ret;
}

sub _run {
  my ($self, $rd1_file, $ref_file, $rd2_file) = @_;
  my ($cmd, $filespec, @ipc_args);
  # Get program executable
  my $exe = $self->executable;
  my $paired = $rd1_file && $rd2_file;
  
  my $tdir = $self->tempdir();
  my ($saih, $saif) = $self->io->tempfile(-template => 'saiXXXX', -dir=>$tdir);
  my ($sai2h, $sai2f) = $self->io->tempfile(-template => 'saiXXXX', -dir=>$tdir)
      if $paired;
  my ($samh, $samf) = $self->io->tempfile(-template => 'saiXXXX', -dir=>$tdir);
  $_->close for ($saih, $samh);
  $sai2h->close if $paired;

  my $subcmd_args = $self->_collate_subcmd_args();
  # index the fasta file (bwa's, not samtools', index...)
  my $bwa = Bio::Tools::Run::BWA->new(
      -command => 'index',
      @{$subcmd_args->{idx}}
      );
  $bwa->run_bwa( -fas => $ref_file );
  # map reads to reference seqs
  $bwa = Bio::Tools::Run::BWA->new(
      -command => 'aln',
      @{$subcmd_args->{aln}}
      );

  $bwa->run_bwa( -fas => $ref_file, -faq => $rd1_file, -sai => $saif );
  # do paired run if nec
  $bwa->run_bwa( -fas => $ref_file, -faq => $rd2_file, -sai => $sai2f ) 
      if $paired;
  # assemble reads
  $bwa = Bio::Tools::Run::BWA->new(
      -command => ($paired ? 'sampe' : 'samse'),
      @{$subcmd_args->{ ($paired ? 'smp' : 'sms' ) }}
      );
  if ($paired) {
      $bwa->run_bwa( -fas => $ref_file, -sai1 => $saif, -faq1 => $rd1_file,
		     -sai2 => $sai2f, -faq2 => $rd2_file, -sam => $samf );
  }
  else {
      $bwa->run_bwa( -fas => $ref_file, -sai => $saif, -faq => $rd1_file,
		     -sam => $samf );
  }
  # note this returns a text-sam file-- needs conversion for B:A:IO::sam...
  # conversion done in run(), if Bio::Tools::Run::Samtools available.
  return $samf;

}

sub available_parameters {
    my $self = shift;
    my $subset = shift;
    for ($subset) { # get commands
	!defined && do { # delegate
	    return $self->SUPER::available_parameters($subset);
	};	    
	m/^c/i && do {
	    return grep !/^run$/, @program_commands;
	};
	m/^f/i && do { # get file spec
	    return @{$command_files{$self->command}};
	};
	do { #else delegate...
	    return $self->SUPER::available_parameters($subset);
	};
    }
}

sub available_commands { shift->available_parameters('commands') };

sub filespec { shift->available_parameters('filespec') };
	
1;