Toolbar

Summary

Bio::EnsEMBL::Analysis::Runnable::PSILC;

Package variables

Included modules

Inherit

Bio::EnsEMBL::Analysis::Runnable

Synopsis

  my $blast = Bio::EnsEMBL::Analysis::Runnable::PSILC_BlastP->new 
    (
     '-trans'    => $transcript,
     '-analysis' => $analysis,
    );

  $blast->run;
  my $output = $blast->output;

Description

BLASTP runnable for Bio::EnsEMBL::Analysis::RunnableDB::PSILC
Runs BLASTP for the supplied transcript against a multi species blast database.
Filters the results and returns a hash ref of homologous transcript identifiers

Methods

Methods description

  Arg [1]   : array of Bio::EnsEMBL::Transcript
  Function  : cretes PSILC BLASTP runnable
  Returntype: Bio::EnsEMBL::Analysis::Runnable::PSILC_BlastP
  Exceptions: none

  Arg [1]    : hash ref
  Description: over-ridden output in runnable.pm to allow hash refs rather than array refs
  Returntype : array ref to Bio::EnsEMBL::Transcript objects
  Exceptions : none
  Caller     : general

  Arg [1]   : none
  Function  : runs BLAST on transcript against multi species blast database
  Returntype: Hash containing transcript ids of homologs that pass criteria:
target is non identical to query has at least 80% coverage and 50% sequence identity
  Exceptions: throws if it cannot recognise the transcript identifier

Methods code

sub new {

  my ($class,@args) = @_;
  my $self = $class->SUPER::new(@args);

  $self->{'_trans'} = {};	#gene to test;  

  my( $trans) = $self->_rearrange([qw(
				     TRANS
				    )], @args);
  if ($trans) {
    $self->trans($trans);
  }
  return $self;

}

sub output {

  my ($self, $trans) = @_;
  if ($trans) {
    $self->{'_output'} = $trans;
  }
  return $self->{'_output'};
}



1;

}

sub run {

  my ($self)=@_;
  my $trans = $self->trans;
  $self->run_blast($trans);

}

sub run_blast {

  my ($self,$trans)=@_;
  my %transcript_hash;

  print STDERR "Blast search of ".$trans->dbID.".\n";
  # might want to use minimal blast module
  my $bplitewrapper = Bio::EnsEMBL::Analysis::Tools::BPliteWrapper-> new
    (
     #  -regex => '^\w+\s+(\w+)'
     -query_type => 'dna',
     -database_type => 'dna',
    );
  # make a query object containing the cds

  my $query = Bio::Seq->new(
			    '-display_id' => $trans->dbID,
			    '-seq'        => $trans->translate->seq,
			    );

  my $blast =  Bio::EnsEMBL::Analysis::Runnable::Blast->new 
    ('-query'     => $query,
     '-program'   => 'blastp',
     '-database'  => "$PSILC_BLAST_DB/multi_species.fasta",
     '-threshold' => 1e-6,
     '-parser'    => $bplitewrapper,
     '-options'   => 'V=10 -cpus=1',
     '-analysis'  => $self->analysis,
    );

  $blast->run();
  my $results = $blast->output; 

  # get dnadnaalignfeatures for each HSP
  # include in set if > 50% ID
  # dont include if > 99%id and 80% coverage
  # Coverage based on cds not cdna length

  foreach my $daf (@{$results}) {
    my $coverage = $daf->length/length($query->seq);
    if ($daf->percent_id <= 99 && 
	$coverage > 0.8  && 
	$daf->percent_id > 50 ){    
      # organise the transcripts in to a hash to avoid redundancy and
      # group them by species 
      if ($daf->hseqname =~ /(\w+)_(\w+)/) {
	push @{$transcript_hash{$2}},$1;
      } else {
	$self->throw("Cannot recognise trans_id ".$daf->hseqname."  $@\n");
      }
    }
  }
  $self->output(\%transcript_hash);
}



#######################################
# Containers

}

sub trans {

  my ($self, $trans) = @_;
  if ($trans) {
    unless  ($trans->isa("Bio::EnsEMBL::Transcript")){
      $self->throw("Input isn't a Bio::EnsEMBL::Transcript, it is a $trans\n$@");
    }
    $self->{'_trans'} = $trans;
  }
  return $self->{'_trans'};

}

General documentation

Post questions to the Ensembl development list: ensembl-dev@ebi.ac.uk