diff variant_effect_predictor/Bio/EnsEMBL/Funcgen/RunnableDB/RunBWA.pm @ 0:1f6dce3d34e0

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/variant_effect_predictor/Bio/EnsEMBL/Funcgen/RunnableDB/RunBWA.pm	Thu Apr 11 02:01:53 2013 -0400
@@ -0,0 +1,103 @@
+=pod 
+
+=head1 NAME
+
+Bio::EnsEMBL::Hive::RunnableDB::Funcgen::RunBWA
+
+=head1 DESCRIPTION
+
+'RunBWA' runs BWA with an input file
+
+=cut
+
+package Bio::EnsEMBL::Funcgen::RunnableDB::RunBWA;
+
+use base ('Bio::EnsEMBL::Funcgen::RunnableDB::Alignment');
+
+use warnings;
+use strict;
+use Bio::EnsEMBL::DBSQL::DBAdaptor;
+use Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor; 
+use Bio::EnsEMBL::Utils::Exception qw(throw warning stack_trace_dump);
+
+sub fetch_input {   # fetch parameters...
+  my $self = shift @_;
+  
+  $self->SUPER::fetch_input();
+  $self->_input_dir($self->_output_dir());
+
+  my $input_file = $self->param('input_file') || throw "No input file given";
+  $input_file = $self->_input_dir()."/".$input_file;
+  if(! -e $input_file){ throw " Couldn't find $input_file"; }
+  $self->_input_file($input_file);
+
+  my $species = $self->_species();
+  my $gender = $self->_cell_type()->gender();
+  $gender = $gender ? $gender : "male";
+  my $assembly = $self->_assembly();
+  #TODO: Maybe check if the index file is really there? Eventually the bwa output will tell you though
+  my $bwa_index = $self->_work_dir()."/bwa_indexes/".$species."/".$species."_".$gender."_".$assembly."_unmasked.fasta";
+  $self->_bwa_index($bwa_index);
+
+  my $bwa_bin = $self->param('bwa_bin') || $self->_bin_dir().'/bwa';
+  $self->_bwa_bin($bwa_bin);
+
+  return 1;
+}
+
+sub run {   
+  my $self = shift @_;
+
+  my $input_file = $self->_input_file();
+  my $bwa_index = $self->_bwa_index();
+
+  my $bwa_bin = $self->_bwa_bin();
+
+  #TODO Pass the location of the binary to be sure we'll be running the right version?
+#  my $bwa_cmd = "$bwa_bin aln $bwa_index $input_file";
+  #Allow this to work with paired reads?? Maybe not for the moment...
+  #in that case pass bwa algorithm as parameter...
+  #If using -q make sure we've got the correct Sanger quality scores...
+#  $bwa_cmd .= " | $bwa_bin samse $bwa_index - $input_file";
+#  $bwa_cmd .= " | samtools view -uS - ";
+#  $bwa_cmd .= " | samtools sort - ${input_file}.sorted";
+#  if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd";  }
+
+  #Silent errors seem to be passing... running one command at a time!?
+  my $bwa_cmd = "$bwa_bin aln $bwa_index $input_file > ${input_file}.aln";
+  if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd";  }
+  $bwa_cmd = "$bwa_bin samse $bwa_index ${input_file}.aln $input_file > ${input_file}.aln.sam";
+  if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd";  }
+  if(system("rm ${input_file}.aln") != 0){ warn "Couldn't remove tmp file. Remove it manually."; }
+  $bwa_cmd = $self->_bin_dir()."/samtools view -uS ${input_file}.aln.sam > ${input_file}.aln.bam";
+  if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd";  }
+  if(system("rm ${input_file}.aln.sam") != 0){ warn "Couldn't remove tmp file. Remove it manually."; }
+  $bwa_cmd = $self->_bin_dir()."/samtools sort ${input_file}.aln.bam ${input_file}.sorted";
+  if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd";  }
+  if(system("rm ${input_file}.aln.bam") != 0){ warn "Couldn't remove tmp file. Remove it manually."; }
+
+  if(system("rm $input_file") != 0){ warn "Couldn't remove tmp file. Remove it manually."; }
+
+  return 1;
+}
+
+
+sub write_output {  # Nothing to write
+  my $self = shift @_;
+
+  return 1;
+
+}
+
+
+#Private getter / setter to the bwa indexes
+sub _bwa_index {
+  return $_[0]->_getter_setter('bwa_index',$_[1]);
+}
+
+#Private getter / setter to the bwa bin
+sub _bwa_bin {
+  return $_[0]->_getter_setter('bwa_bin',$_[1]);
+}
+
+1;