# something like this
my $query = new Bio::Seq(-file => $queryfile,
-format => 'Fasta');
my $hmm = Bio::EnsEMBL::Pipeline::Runnable::Protein::Hmmpfam->new
('-query' => $query,
'-program' => 'hmmpfam' or '/usr/local/pubseq/bin/hmmpfam',
'-database' => '/data/Pfam_ls');
$hmm->run;
my @results = @{$hmm->output};
sub parse_results
{ my ($self) = @_;
my (@hits, $id);
my $f = $self->resultsfile;
if (-e $f and not -z $f) {
my $fh;
open ($fh, "<$f") or throw ("Error opening $f");
while (<$fh>) {
chomp;
if (/^Query sequence:\s+(\S+)/) {
$id = $1;
next;
}
if (my ($hid,
$start,
$end,
$hstart,
$hend,
$score,
$evalue) = /^(\S+)\s+\d+\/\d+\s+(\d+)\s+(\d+)\s+\S+\s+(\d+)\s+(\d+)\s+\S+\s+(\S+)\s+(\S+)/) {
$evalue = sprintf ("%.3e", $evalue);
my $fp = $self->create_protein_feature($start,
$end,
$score,
$id,
$hstart,
$hend,
$hid,
$self->analysis,
$evalue,
0);
push @hits, $fp;
}
}
}
$self->output(\@hits);
}
1; } |
sub run_analysis
{ my ($self) = @_;
my $options = "";
if (defined($self->options)) {
$options .= $self->options;
}
if ($options !~ /\-\-acc/) {
$options .= ' --acc';
}
if ($options !~ /\-\-cpu/) {
$options .= ' --cpu 1';
}
my $cmd = $self->program
. ' ' . $options
. ' ' . $self->database
. ' ' . $self->queryfile
.' > '. $self->resultsfile;
print STDERR "Running:$cmd\n";
throw ("Error running ".$self->program." on ".
$self->queryfile." against ".$self->database)unless
((system ($cmd)) == 0); } |