Raw content of Bio::EnsEMBL::Analysis::Runnable::Funcgen::Splitter
<![CDATA[# Ensembl module for Bio::EnsEMBL::Analysis::Runnable::Funcgen::Splitter
#
# Copyright (c) 2007 Ensembl
#
=head1 NAME
Bio::EnsEMBL::Analysis::Runnable::Funcgen::Splitter
=head1 SYNOPSIS
my $runnable = Bio::EnsEMBL::Analysis::Runnable::Funcgen::Splitter->new
(
-analysis => $analysis,
-query => 'slice',
-program => 'program.pl',
);
$runnable->run;
my @features = @{$runnable->output};
=head1 DESCRIPTION
Splitter expects to run the command line version of the Splitter program
(http://zlab.bu.edu/splitter) and predicts features which can be stored in
the annotated_feature table in the eFG database
=head1 LICENCE
This code is distributed under an Apache style licence. Please see
/info/about/code_licence.html for details.
=head1 AUTHOR
Stefan Graf, Ensembl Functional Genomics -
=head1 CONTACT
Post questions to the Ensembl development list: ensembl-dev@ebi.ac.uk
=cut
package Bio::EnsEMBL::Analysis::Runnable::Funcgen::Splitter;
use strict;
use warnings;
use Data::Dumper;
use Bio::EnsEMBL::Analysis::Runnable;
use Bio::EnsEMBL::Analysis::Runnable::Funcgen;
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use Bio::EnsEMBL::Utils::Argument qw( rearrange );
use vars qw(@ISA);
@ISA = qw(Bio::EnsEMBL::Analysis::Runnable::Funcgen);
=head2 write_infile
Arg [1] : Bio::EnsEMBL::Analysis::Runnable::Splitter
Arg [2] : filename
Description :
Returntype :
Exceptions :
Example :
=cut
sub write_infile {
my ($self, $filename) = @_;
if (! $filename) {
$filename = $self->infile();
}
# everything needs to happen in a particular subdirectory of the workdir,
# filenames are set by the Splitter program
(my $workdir = $filename) =~ s/\.dat$//;
eval { system("mkdir -p $workdir") };
throw ("Couldn't create directory $workdir: $@") if ($@);
$self->workdir($workdir);
$filename =~ s,.+/,,;
# determine both number of result sets (replicates) and features
my $noa = scalar(keys %{$self->result_features});
warn('no. of result sets: ', $noa);
my $nof = scalar(@{(values %{$self->result_features})[0]});
warn('no. of result fts: ', $nof);
# dump features
open(F, "> $workdir/$filename")
or throw("Can't open file $workdir/$filename.");
#print Dumper $self->result_features;
for (my $i=0; $i<$nof; $i++) {
my $coord=0;
foreach my $rset (values %{$self->result_features}) {
printf F "chr%s\t%d\t%d\t", ${$rset}[$i][0], ${$rset}[$i][1], ${$rset}[$i][2] unless ($coord);
$coord=1;
print F "\t".${$rset}[$i][3];
}
print F "\n";
}
close F;
$self->infile($workdir.'/'.$filename);
return "$workdir/$filename";
}
=head2 run_analysis
Arg [1] : Bio::EnsEMBL::Analysis::Runnable::Splitter
Arg [2] : string, program name
Usage :
Description :
Returns :
Exceptions :
=cut
### step 1
# ./splitter.cgi "id=SplitterJob.$(gdate -I)" 'submit=Step by step' 'signalfile=samplesplittersignal' 'sorted=yes' 'norm=qnorm' 'combine=median' 'splitterstep=10' 'splitterratio=2' 'cutoff=2.5sd' 'maxgap=100' 'minrun=5'
### step 2
# ./splitter.cgi "id=SplitterJob.$(gdate -I)" 'submit=Step 2' 'sorted=yes' 'norm=none' 'combine=median'
### step 3
# ./splitter.cgi "id=SplitterJob.$(gdate -I)" 'submit=Step 3' 'last=pipeline.median' 'splitterfrom=0.0149' 'splitterto=0.4884' 'splitterstep=10' 'splitterratio=2' 'cutoff=2.5sd' 'custom=0.4884' '2.5sd=0.4884' '1pct=0.395' '5pct=0.25' '2sd=0.3937' 'maxgap=100' 'minrun=5'
### One shot (all in one)
# ./splitter.cgi "id=SplitterJob.TEST" 'submit=One shot' 'signalfile=samplesplittersignal' 'sorted=yes' 'norm=none' 'combine=median' 'splitterstep=10' 'splitterratio=2' 'cutoff=2.5sd' 'maxgap=100' 'minrun=5'
#available Splitter parameters (for details see http://zlab.bu.edu/splitter)
# -Input already sorted by genomic positions of probes
# 'sorted=yes|no'
# -Perform normalization between replicates
# 'norm=none|qnorm|rnorm'
# -Combine replicates
# 'combine=mean|median'
# -Signal cutoff for the combined intensities>=
# 'cutoff=splitter|5pct|1pct|2sd|2.5sd'
# 'splitterfrom='
# 'splitterto='
# 'splitterstep=10'
# 'splitterratio=2'
# -Determine how the probes are clustered
# Gap (Maxgap) <= n base pairs
# 'maxgap=100'
# Clustering (Minrun) >= n probes
# 'minrun=4'
sub run_analysis {
my ($self, $program) = @_;
if(!$program){
$program = $self->program;
}
throw($program." is not executable Splitter::run_analysis ")
unless($program && -x $program);
(my $id = $self->workdir) =~ s,^.+/,,;
my $command = $self->program." \'id=$id\' \'submit=One shot\' \'signalfile=".$self->infile."\' " .
$self->analysis->parameters;
warn("Running analysis " . $command . "\n");
eval { system( $command ) };
throw("FAILED to run $command: ", $@) if ($@);
opendir(DIR, $self->workdir())
or throw("Can't open workdir ".$self->workdir());
my @resultfiles = grep { /\.scored\.txt$/ } readdir DIR;
closedir(DIR);
throw("More than one resultfile in ".$self->workdir())
if (scalar(@resultfiles) > 1);
# set resultsfile
$self->resultsfile($self->workdir.'/'.$resultfiles[0]);
}
=head2 infile
Arg [1] : Bio::EnsEMBL::Analysis::Runnable::Splitter
Arg [2] : filename (string)
Description : will hold a given filename or if one is requested but none
defined it will use the create_filename method to create a filename
Returntype : string, filename
Exceptions : none
Example :
=cut
sub infile{
my ($self, $filename) = @_;
if($filename){
$self->{'infile'} = $filename;
}
if(!$self->{'infile'}){
$self->{'infile'} = $self->create_filename($self->analysis->logic_name, 'dat');
}
return $self->{'infile'};
}
=head2 parse_results
Arg [1] : Bio::EnsEMBL::Analysis::Runnable::Chipotle
Arg [2] : filename (string)
Decription : open and parse resultsfile
Returntype : none
Exceptions : throws
Example :
=cut
sub parse_results{
my ($self, $resultsfile) = @_;
if (! defined $resultsfile) {
$resultsfile = $self->resultsfile;
}
throw ('No resultsfile defined in instance!')
if(! defined $resultsfile);
throw("parse_results: results file ".$resultsfile." does not exist.")
if (! -e $resultsfile);
throw("parse_results: can't open file ".$resultsfile.": ". $!)
unless (open(F, $resultsfile));
my @output = ();
while (<F>) {
next unless (s/^chr//);
chomp;
my @ft = split;
my $coord = shift(@ft);
my @coord = split(/[:-]/, $coord);
push(@output, [ @coord, $ft[0] ]);
}
throw("No features to annotate on slice ".$self->query->seq_region_name."!") if (scalar(@output) == 0);
$self->output(\@output);
#print Dumper $self->output();
throw("parse_results: can't close file ".$resultsfile.".")
unless (close(F));
}
1;
]]>