|
0
|
1 =pod
|
|
|
2
|
|
|
3 =head1 NAME
|
|
|
4
|
|
|
5 Bio::EnsEMBL::Hive::RunnableDB::Funcgen::RunBWA
|
|
|
6
|
|
|
7 =head1 DESCRIPTION
|
|
|
8
|
|
|
9 'RunBWA' runs BWA with an input file
|
|
|
10
|
|
|
11 =cut
|
|
|
12
|
|
|
13 package Bio::EnsEMBL::Funcgen::RunnableDB::RunBWA;
|
|
|
14
|
|
|
15 use base ('Bio::EnsEMBL::Funcgen::RunnableDB::Alignment');
|
|
|
16
|
|
|
17 use warnings;
|
|
|
18 use strict;
|
|
|
19 use Bio::EnsEMBL::DBSQL::DBAdaptor;
|
|
|
20 use Bio::EnsEMBL::Funcgen::DBSQL::DBAdaptor;
|
|
|
21 use Bio::EnsEMBL::Utils::Exception qw(throw warning stack_trace_dump);
|
|
|
22
|
|
|
23 sub fetch_input { # fetch parameters...
|
|
|
24 my $self = shift @_;
|
|
|
25
|
|
|
26 $self->SUPER::fetch_input();
|
|
|
27 $self->_input_dir($self->_output_dir());
|
|
|
28
|
|
|
29 my $input_file = $self->param('input_file') || throw "No input file given";
|
|
|
30 $input_file = $self->_input_dir()."/".$input_file;
|
|
|
31 if(! -e $input_file){ throw " Couldn't find $input_file"; }
|
|
|
32 $self->_input_file($input_file);
|
|
|
33
|
|
|
34 my $species = $self->_species();
|
|
|
35 my $gender = $self->_cell_type()->gender();
|
|
|
36 $gender = $gender ? $gender : "male";
|
|
|
37 my $assembly = $self->_assembly();
|
|
|
38 #TODO: Maybe check if the index file is really there? Eventually the bwa output will tell you though
|
|
|
39 my $bwa_index = $self->_work_dir()."/bwa_indexes/".$species."/".$species."_".$gender."_".$assembly."_unmasked.fasta";
|
|
|
40 $self->_bwa_index($bwa_index);
|
|
|
41
|
|
|
42 my $bwa_bin = $self->param('bwa_bin') || $self->_bin_dir().'/bwa';
|
|
|
43 $self->_bwa_bin($bwa_bin);
|
|
|
44
|
|
|
45 return 1;
|
|
|
46 }
|
|
|
47
|
|
|
48 sub run {
|
|
|
49 my $self = shift @_;
|
|
|
50
|
|
|
51 my $input_file = $self->_input_file();
|
|
|
52 my $bwa_index = $self->_bwa_index();
|
|
|
53
|
|
|
54 my $bwa_bin = $self->_bwa_bin();
|
|
|
55
|
|
|
56 #TODO Pass the location of the binary to be sure we'll be running the right version?
|
|
|
57 # my $bwa_cmd = "$bwa_bin aln $bwa_index $input_file";
|
|
|
58 #Allow this to work with paired reads?? Maybe not for the moment...
|
|
|
59 #in that case pass bwa algorithm as parameter...
|
|
|
60 #If using -q make sure we've got the correct Sanger quality scores...
|
|
|
61 # $bwa_cmd .= " | $bwa_bin samse $bwa_index - $input_file";
|
|
|
62 # $bwa_cmd .= " | samtools view -uS - ";
|
|
|
63 # $bwa_cmd .= " | samtools sort - ${input_file}.sorted";
|
|
|
64 # if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd"; }
|
|
|
65
|
|
|
66 #Silent errors seem to be passing... running one command at a time!?
|
|
|
67 my $bwa_cmd = "$bwa_bin aln $bwa_index $input_file > ${input_file}.aln";
|
|
|
68 if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd"; }
|
|
|
69 $bwa_cmd = "$bwa_bin samse $bwa_index ${input_file}.aln $input_file > ${input_file}.aln.sam";
|
|
|
70 if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd"; }
|
|
|
71 if(system("rm ${input_file}.aln") != 0){ warn "Couldn't remove tmp file. Remove it manually."; }
|
|
|
72 $bwa_cmd = $self->_bin_dir()."/samtools view -uS ${input_file}.aln.sam > ${input_file}.aln.bam";
|
|
|
73 if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd"; }
|
|
|
74 if(system("rm ${input_file}.aln.sam") != 0){ warn "Couldn't remove tmp file. Remove it manually."; }
|
|
|
75 $bwa_cmd = $self->_bin_dir()."/samtools sort ${input_file}.aln.bam ${input_file}.sorted";
|
|
|
76 if(system($bwa_cmd) != 0){ throw "Problems running $bwa_cmd"; }
|
|
|
77 if(system("rm ${input_file}.aln.bam") != 0){ warn "Couldn't remove tmp file. Remove it manually."; }
|
|
|
78
|
|
|
79 if(system("rm $input_file") != 0){ warn "Couldn't remove tmp file. Remove it manually."; }
|
|
|
80
|
|
|
81 return 1;
|
|
|
82 }
|
|
|
83
|
|
|
84
|
|
|
85 sub write_output { # Nothing to write
|
|
|
86 my $self = shift @_;
|
|
|
87
|
|
|
88 return 1;
|
|
|
89
|
|
|
90 }
|
|
|
91
|
|
|
92
|
|
|
93 #Private getter / setter to the bwa indexes
|
|
|
94 sub _bwa_index {
|
|
|
95 return $_[0]->_getter_setter('bwa_index',$_[1]);
|
|
|
96 }
|
|
|
97
|
|
|
98 #Private getter / setter to the bwa bin
|
|
|
99 sub _bwa_bin {
|
|
|
100 return $_[0]->_getter_setter('bwa_bin',$_[1]);
|
|
|
101 }
|
|
|
102
|
|
|
103 1;
|
