Raw content of Bio::EnsEMBL::Pipeline::GeneDuplication::Finder
<![CDATA[package Bio::EnsEMBL::Pipeline::GeneDuplication::Finder;
use vars qw(@ISA);
use strict;
use Bio::EnsEMBL::Root;
use Bio::Tools::Run::Alignment::Clustalw;
use Bio::EnsEMBL::Pipeline::Runnable::BlastDB;
use Bio::EnsEMBL::Pipeline::Runnable::MinimalBlast;
use Bio::EnsEMBL::Pipeline::SeqFetcher::FetchFromBlastDB;
use Bio::EnsEMBL::Pipeline::GeneDuplication::PAML;
use Bio::EnsEMBL::Pipeline::GeneDuplication::CodonBasedAlignment;
use Bio::EnsEMBL::Pipeline::GeneDuplication::Result;
use Bio::EnsEMBL::Utils::Exception qw(warning throw);
use Bio::EnsEMBL::Utils::Argument qw(rearrange);
my $DEFAULT_DISTANCE_METHOD = 'NeiGojobori';
my $DEFAULT_BLAST_PROGRAM = 'wublastn';
my $DEFAULT_DISTANCE_CUTOFF = 1.000;
@ISA = qw(Bio::EnsEMBL::Root Bio::EnsEMBL::Pipeline::RunnableI);
### Constructor ###
=head2 new
Args : -dbfile - (string) Full path to a fasta formatted file of
nucleotide sequences.
-blastdb - A Bio::EnsEMBL::Pipeline::Runnable::BlastDB
object for which the run method has been invoked.
-query_seq - A Bio::Seq which is comprised of nucleotides.
-blast_program - Manually set the blast program (and
full path) that should be used. Make sure the
program matches the index type chosen. Defaults
to wublastn.
-blast_index_type - The distribution of blast to use. See
Bio::EnsEMBL::Pipeline::Runnable::BlastDB->index_type
documentation. Defaults to 'wu_new'.
-hit_identity - (optional) Hit identity. Defaults to 0.80
-hit_coverage - (optional) Hit coverage. Defults to 0.80
-work_dir - (optional) Dir where working files are
placed. Defaults to /tmp.
-codeml - Full path to codeml executable.
-genetic_code - 0 for Universal, 1 for Mitochondrial. See
the Bio::Seq->translate method for a full
list of options.
-regex_query_species - A regular expression that will parse
some portion of the ids of the query species (and
not the ids of outgroup species). Eg. 'ENSG'.
-regex_outgroup_species - A ref to an array of regular expressions
that will parse the ids of the various outgroup
species (but not the ids of the query species).
E.g. ['ENSMUSG', 'ENSRNOG']
-distance_cutoff - A genetic distance cutoff to apply for
occasions where an outgroup species is not set, or
an outgroup species match was not found.
-distance_method - The method used to calculate the genetic
distance and Ka/Ks ratios. Options are 'NeiGojobori'
and 'ML', ML being the maximum likelihood method
of Yang.
Example : none
Description: Constructs new object
Returntype : Bio::EnsEMBL::Pipeline::GeneDuplication::Finder
Exceptions : Throws if database file not specified or does not exist.
Caller : General
=cut
sub new {
my ($class, @args) = @_;
my $self = bless {}, $class;
my ($query,
$blastdb,
$blast_program,
$blast_index_type,
$work_dir,
$codeml,
$genetic_code,
$regex_query_species,
$regex_outgroup_species,
$identity_cutoff,
$coverage_cutoff,
$distance_cutoff,
$distance_method) = rearrange([qw(QUERY
BLASTDB
BLAST_PROGRAM
BLAST_INDEX_TYPE
WORK_DIR
CODEML_EXECUTABLE
GENETIC_CODE
REGEX_QUERY_SPECIES
REGEX_OUTGROUP_SPECIES
HIT_IDENTITY
HIT_COVERAGE
DISTANCE_CUTOFF
DISTANCE_METHOD
)],@args);
$self->_work_dir($work_dir) if $work_dir;
if ($blastdb && $blastdb->isa("Bio::EnsEMBL::Pipeline::Runnable::BlastDB")){
$self->_blastdb($blastdb);
} else {
throw ("Need a Bio::EnsEMBL::Pipeline::Runnable::BlastDB object.");
}
$self->_query_seq($query) if $query;
$self->_codeml($codeml) if $codeml;
$self->_genetic_code($genetic_code);
$self->_regex_query_species($regex_query_species) if $regex_query_species;
$self->_regex_outgroup_species($regex_outgroup_species) if $regex_outgroup_species;
$self->_identity_cutoff($identity_cutoff) if $identity_cutoff;
$self->_coverage_cutoff($coverage_cutoff) if $coverage_cutoff;
$self->_blast_program($blast_program) if $blast_program;
$self->_identity_cutoff(80) unless $self->_identity_cutoff;
$self->_coverage_cutoff(80) unless $self->_coverage_cutoff;
$self->_distance_method($distance_method) if $distance_method;
$self->_distance_cutoff($distance_cutoff) if $distance_cutoff;
return $self
}
### Public methods ###
=head2 run
Args[1] : [optional] Input sequence (Bio::Seq)
Example : none
Description: Top level method that executes the gene duplicate
finding algorithm.
Returntype : Bio::EnsEMBL::Pipeline::GeneDuplication::Result
Exceptions : Warns if PAML run fails.
Throws if PAML returns multiple result sets (unlikely).
Caller : General
=cut
sub run {
my ($self, $input_seq) = @_;
$self->_query_seq($input_seq)
if $input_seq;
# Derive a list of sequences that look like duplications of
# the query sequence.
my $accepted_hits
= $self->_find_recent_duplications($self->_query_seq);
unless (scalar @$accepted_hits > 1) {
print "No homologous matches were found that satisfied the match criteria.\n";
return 0
}
return $accepted_hits;
}
### Hit Chooser Methods ###
=head2 _find_recent_duplications
Args[1] : Bio::Seq query sequence
Example : none
Description: This is the implementation of the main algorithm of
this module.
Returntype :
Exceptions :
Caller : General
=cut
sub _find_recent_duplications {
my ($self, $query) = @_;
# Perform blast search with this query sequence.
my $bplite_report = $self->_run_blast($query);
# Take blast report and filter for hits with good coverage
# and sufficiently high identity.
my $good_hits = $self->_preliminary_filter($bplite_report);
# Check that good hits have been found.
unless (scalar(keys %$good_hits)){
print "Did not find hits to query sequence.\n";
return []
}
# From the set of promising hits, which will possibly include
# hits to other species, filter using outgroup sequences (if
# any) and the genetic distance between each hit and the query
# sequence.
my $recent_duplications = $self->_phylogenetic_filter($good_hits);
# Return a complex data object that contains the pairwise
# match details for each final accepted hit.
return $recent_duplications;
}
=head2 _run_blast
Args[1] : Bio::Seq query sequence
Example : none
Description:
Returntype :
Exceptions :
Caller : General
=cut
sub _run_blast {
my ($self, $query) = @_;
# We are looking for duplicates of the query gene
# in our blast database.
$self->_query_seq($query) if $query;
my $bplite_report;
eval{
$bplite_report = $self->_blast_obj->run;
};
throw ("Blast did not run successfully. Blast program was [".
$self->_blast_program."]. Index type was [".
$self->_blast_index_type."].")
if $@;
throw ("Blast process did not return a report.")
unless ($bplite_report->isa("Bio::Tools::BPlite"));
return $bplite_report
}
=head2 _preliminary_filter
Args[1] :
Example : none
Description:
Returntype :
Exceptions :
Caller : General
=cut
sub _preliminary_filter {
my ($self, $bplite_report) = @_;
# Process our blast report.
#
# For each blast hit:
# * throw away self matches (if any)
# * filter hits by identity and coverage
# * align WHOLE genes of promising hits, re-filter by
# identity and coverage
# * calculate genetic distance between each subject and the query
my %good_hits;
DISTINCT_HIT:
while (my $sbjct = $bplite_report->nextSbjct){
# Mangle the BPLite::Sbjct object for its own good. Quite
# often the hit ids parsed by BPLite include the whole
# Fasta header description line. This is problematic if
# sequence ids need to be compared or a hash is keyed on
# this sequence id. Here we simply lop the description
# from each Sbjct->name, if there is one.
$sbjct = $self->_fix_sbjct($sbjct);
# It appears that the BPLite::Sbjct object only allows
# HSPs to be accessed once (as this process is closely
# tied to the parsing of the Blast report). Hence, here
# we loop through them all here and store them in an
# array.
my @hsps;
while (my $hsp = $sbjct->nextHSP) {
push (@hsps, $hsp);
}
# Skip hit if it is a match to self.
next DISTINCT_HIT
if ($self->_query_seq->display_id eq $sbjct->name);
# First, filter hits by identity
my $hit_identity = $self->_hit_identity(\@hsps);
next DISTINCT_HIT
unless ($hit_identity >= $self->_identity_cutoff);
# Second, filter hits by coverage
next DISTINCT_HIT
unless ($self->_appraise_hit_coverage($sbjct, \@hsps));
# Third, this has become a serious hit. Make a pairwise
# alignment of the query and hit.
my @seqs = ($self->_fetch_seq($self->_query_seq->display_id),
$self->_fetch_seq($sbjct->name));
throw ("Didnt correctly obtain two sequences for alignment.")
unless scalar @seqs == 2;
my $align = $self->_pairwise_align(\@seqs);
# and filter by coverage and identity for the whole
# aligned sequences.
my ($min_coverage, $identity)
= $self->_calculate_coverage_and_identity($align);
unless ($min_coverage >= $self->_coverage_cutoff &&
$identity >= $self->_identity_cutoff) {
next DISTINCT_HIT;
}
# Forth, calculate the genetic distance while we are here.
my ($ka, $ks, $n, $s)
= $self->_run_pairwise_paml($align);
# Finally, store all this useful information. Perhaps a proper
# object would make this easier to access later?
die "undefined N and S" unless defined $n and defined $s;
$good_hits{$sbjct->name} = {'Ka' => $ka,
'Ks' => $ks,
'N' => $n,
'S' => $s,
'coverage' => $min_coverage,
'identity' => $identity,
'alignment'=> $align};
}
return \%good_hits
}
=head2 _phylogenetic_filter
Args[1] :
Example : none
Description:
Returntype :
Exceptions :
Caller : General
=cut
sub _phylogenetic_filter {
my ($self, $good_hits) = @_;
my %hits_by_species;
HIT:
foreach my $hit_id (keys %$good_hits) {
# Now, partition hits according to their species. The species
# from which the subject is derived is determined by a
# regular expression match to the sequence id.
my $query_regex = $self->_regex_query_species;
if ($hit_id =~ /$query_regex/) {
push(@{$hits_by_species{$self->_regex_query_species}}, $hit_id);
next HIT;
} else {
foreach my $regex (@{$self->{_regex_outgroup_species}}) {
if ($hit_id =~ /$regex/){
push (@{$hits_by_species{$regex}}, $hit_id);
next HIT;
}
}
}
warning("Didnt match hit id to any regex [". $hit_id ."].");
}
# Return now if there are no intra-species hits.
unless (defined @{$hits_by_species{$self->_regex_query_species}} &&
scalar @{$hits_by_species{$self->_regex_query_species}}) {
return []
}
# Sort our hits by their distance to the query sequence.
my %sorted_hits_by_species;
foreach my $species (keys %hits_by_species) {
my @sorted_hits
= sort {$good_hits->{$a}->{Ks} <=> $good_hits->{$b}->{Ks}}
@{$hits_by_species{$species}};
$sorted_hits_by_species{$species} = \@sorted_hits;
}
# Accept all query species hits with a distance less than
# the distance to the most related outgroup species.
$self->outgroup_distance($self->_distance_cutoff);
my $closest_outgroup_distance = $self->_distance_cutoff;
foreach my $regex (@{$self->_regex_outgroup_species}){
next
unless $sorted_hits_by_species{$regex};
my $closest_hit_id = $sorted_hits_by_species{$regex}->[0];
throw("Missing distance value")
unless defined $good_hits->{$closest_hit_id}->{Ks} ||
$good_hits->{$closest_hit_id}->{Ks} == 0;
my $closest_hit_distance = $good_hits->{$closest_hit_id}->{Ks};
if (($closest_hit_distance < $closest_outgroup_distance) &&
($closest_hit_distance > 0)) {
$closest_outgroup_distance = $closest_hit_distance;
}
}
$self->outgroup_distance($closest_outgroup_distance);
my @accepted_hits;
my @query_species_hits = @{$sorted_hits_by_species{$self->_regex_query_species}};
foreach my $hit_id (@query_species_hits) {
if ($good_hits->{$hit_id}->{Ks} <= $closest_outgroup_distance &&
$good_hits->{$hit_id}->{Ks} >= 0) {
push (@accepted_hits, $good_hits->{$hit_id});
} elsif ($good_hits->{$hit_id}->{Ks} >= 0) {
last;
}
}
return \@accepted_hits
}
### Utility Methods ###
=head2 _run_pairwise_paml
Args : An arrayref to a set of aligned Bio::Seq objects.
Example : none
Description: Uses an array of aligned Bio::Seq objects to execute
PAML in pairwise mode.
Returntype : Bio::Tools::Phylo::PAML
Exceptions : none
Caller : $self->run
=cut
sub _run_pairwise_paml {
my ($self, $aligned_seqs) = @_;
# Cope with a bloody PAML quirk. Identical sequences
# cause a fatal error.
if ($aligned_seqs->[0]->seq eq $aligned_seqs->[1]->seq) {
return (0, 0, 0, 0)
}
# The business of running codeml.
my $retries = 0;
my $result;
my $paml; # Object must stay in scope until parsing is complete.
while (! defined $result && $retries < 3) {
$paml = Bio::EnsEMBL::Pipeline::GeneDuplication::PAML->new(
'-work_dir' => $self->_work_dir,
'-executable' => $self->_codeml,
'-aligned_seqs' => $aligned_seqs,
'-runmode' => '-2',
'-seqtype' => '1',
'-model' => '0',
'-nssites' => '0',
'-icode' => ($self->_genetic_code) - 1
);
my $parser = $paml->run_codeml;
# Often a codeml run might look successful, but the Bioperl parser
# will not be able to read the output created. This is usually
# due to a silent codeml failure and a lack of output. Hence,
# here an attempt is made to check the success of the parser.
eval {
$result = $parser->next_result();
};
if ($@ || ! defined $result) {
warning("Problem deriving PAML result object.\n$@")
}
$retries++
}
unless (defined $result) {
throw ("Result has become undefined. This is probably a problem " .
"with Bio::Tools::Phylo::PAML.")
}
my $matrix;
eval {
if ($self->_distance_method eq 'NeiGojobori') {
$matrix = $result->get_NGmatrix()
} else {
$matrix = $result->get_MLmatrix()
}
};
if ($@){
throw ("PAML failed to give a file that could be parsed.\n" .
"This is a PAML/codeml error.\n" .
"The offending aligned sequences are :\n>" .
$aligned_seqs->[0]->display_id . "\n" . $aligned_seqs->[0]->seq . "\n>" .
$aligned_seqs->[1]->display_id . "\n" . $aligned_seqs->[1]->seq . "\n$@");
}
$matrix->[0]->[1]->{N} = 0 unless defined $matrix->[0]->[1]->{N};
$matrix->[0]->[1]->{S} = 0 unless defined $matrix->[0]->[1]->{S};
return ($matrix->[0]->[1]->{dN},
$matrix->[0]->[1]->{dS},
$matrix->[0]->[1]->{N},
$matrix->[0]->[1]->{S})
}
=head2 _pairwise_align
Args[1] :
Example :
Description:
Returntype :
Exceptions :
Caller :
=cut
sub _pairwise_align {
my ($self, $seqs) = @_;
throw ("Pairwise alignment was only expecting two sequences.")
unless ((scalar @$seqs) == 2);
my $cba
= Bio::EnsEMBL::Pipeline::GeneDuplication::CodonBasedAlignment->new(
-genetic_code => 1);
$cba->sequences($seqs);
my $aligned_seqs = $cba->run_alignment;
return $aligned_seqs
}
=head2 _calculate_coverage_and_identity
Args[1] :
Example :
Description:
Returntype : Array (min_coverage, identity)
Exceptions :
Caller :
=cut
sub _calculate_coverage_and_identity {
my ($self, $align) = @_;
my $seq1 = $align->[0]->seq;
my $seq2 = $align->[1]->seq;
my $align_len = length($seq1);
my $seq_length1 = 0;
my $seq_length2 = 0;
my $matches = 0;
my $identical = 0;
for (my $i = 0; $i < $align_len; $i++) {
my $base1 = substr($seq1, $i, 1);
my $base2 = substr($seq2, $i, 1);
if ($base1 ne '-') {
$seq_length1++;
}
if ($base2 ne '-') {
$seq_length2++;
}
if ($base1 ne '-' && $base2 ne '-'){
$matches++;
}
if ($base1 eq $base2 && $base1 ne '-'){
$identical++;
}
}
my $coverage1 = $matches/$seq_length1;
my $coverage2 = $matches/$seq_length2;
my ($min_coverage) = sort {$a <=> $b} ($coverage1, $coverage2);
my $identity = $identical/$matches;
return ($min_coverage, $identity)
}
=head2 _appraise_hit_coverage
Args[1] :
Example :
Description:
Returntype : 0 or 1
Exceptions :
Caller :
=cut
sub _appraise_hit_coverage{
my ($self, $sbjct, $hsps) = @_;
# Select our sequence length as being that of the
# longer sequence.
my $sbjct_length
= $self->_fetch_seq($sbjct->name)->length;
my ($longest_length) = sort {$a <=> $b} ($self->_query_seq->length, $sbjct_length);
# Look at all the hits along the length of the
# query and tally the collective coverage of the hits.
my @query_coverage;
foreach my $hsp (@$hsps) {
for (my $base_position = $hsp->query->start;
$base_position <= $hsp->query->end;
$base_position++){
$query_coverage[$base_position]++;
}
}
my $covered_bases = 0;
foreach my $base (@query_coverage){
$covered_bases++ if $base;
}
# Return true if good coverage exists.
return 1 if ($covered_bases >= $self->_coverage_cutoff * $longest_length);
# Otherwise return false.
return 0;
}
=head2 _fix_sbjct
Args[1] :
Example :
Description: A work-around for a BPLite::Sbjct annoyance. The
sbjct->name object returns the whole fasta description
line for a subject hit. If the input fasta sequence
file includes more than an id on the description line,
this will be passed back every time the name method is
called. This is a real pest is you are trying to
match ids via a regex or use the ids as hash keys.
Returntype :
Exceptions :
Caller :
=cut
sub _fix_sbjct {
my ($self, $sbjct) = @_;
my $sbjct_name = $sbjct->name;
$sbjct_name =~ s/\W*(\S+).*/$1/;
# BAD!
$sbjct->{NAME} = $sbjct_name;
return $sbjct;
}
=head2 _hit_identity
Args[1] :
Example :
Description:
Returntype :
Exceptions :
Caller :
=cut
sub _hit_identity{
my ($self, $hsps) = @_;
my $tally_query_length = 0;
my $tally_matched_bases = 0;
foreach my $hsp (@$hsps) {
$tally_query_length += length($hsp->querySeq);
$tally_matched_bases += $hsp->positive;
}
if ($tally_matched_bases && $tally_query_length){
return $tally_matched_bases/$tally_query_length;
}
return 0
}
=head2 outgroup_distance
Args[1] :
Example :
Description:
Returntype :
Exceptions :
Caller :
=cut
sub outgroup_distance {
my $self = shift;
if (@_) {
$self->{_outgroup_distance} = shift;
}
return $self->{_outgroup_distance};
}
### Sequence fetching ###
=head2 _fetch_seq
Args[1] :
Example :
Description:
Returntype :
Exceptions :
Caller :
=cut
sub _fetch_seq {
my ($self, $id) = @_;
throw ("Cant fetch sequence without an id.")
unless $id;
$self->{_cache} = {}
unless $self->{_cache};
if ($self->{_cache}->{$id}){
return $self->{_cache}->{$id}
}
my $seq = $self->_seq_fetcher->fetch($id);
throw ("Sequence fetch failed for id [$id].")
unless ($seq && $seq->isa("Bio::Seq"));
$self->{_cache}->{$seq->display_id} = $seq;
return $seq
}
=head2 _fetch_seqs
Args : An arrayref of string sequence ids.
Example : none
Description: An alias to $self->_fetch_seq, but handles
multiple sequences.
Returntype : Arrayref of Bio::Seq
Exceptions : none
Caller : $self->run
=cut
sub _fetch_seqs {
my ($self, $seq_ids) = @_;
my @seqs;
foreach my $seq_id (@$seq_ids){
push (@seqs, $self->_fetch_seq($seq_id));
}
return \@seqs;
}
=head2 _force_cache
Args[1] : Bio::Seq
Example : $self->_force_cache($seq);
Description: Allows a sequence to be manually added to the seqfetcher
cache. This is useful for coping with user supplied
sequences (eg. passed as a query sequence) that dont
exist in any database.
Returntype : 1
Exceptions : Warns if sequence already exists in cache. Throws if
a defined sequence isnt supplied.
Caller :
=cut
sub _force_cache {
my ($self, $seq) = @_;
throw ("Trying to add something odd to the sequence cache [$seq].")
unless (defined $seq);
if ($self->{_cache}->{$seq->display_id}){
warning('Sequence [' . $seq->display_id .
'] already exists in cache, but will replace.');
}
$self->{_cache}->{$seq->display_id} = $seq;
return 1
}
=head2 _seq_fetcher
Args : (optional) A seqfetcher of any variety, as long as
it has a 'fetch' method.
Example : none
Description: Holds SeqFetcher object.
Returntype : Bio::EnsEMBL::Pipeline::SeqFetcher::xxx
Exceptions : none
Caller : $self->_candidate_hits, $self->_fetch_seqs
=cut
sub _seq_fetcher {
my $self = shift;
$self->{_seq_fetcher} = shift if @_;
if (! $self->{_seq_fetcher}){
$self->{_seq_fetcher} =
Bio::EnsEMBL::Pipeline::SeqFetcher::FetchFromBlastDB->new(
-db => $self->_blastdb);
}
return $self->{_seq_fetcher};
}
### Blast-related Getter/Setters ###
=head2 _dbfile
Args : (optional) String
Example : none
Description: Holds the filename of the blast database.
Returntype : String.
Exceptions : none
Caller : $self->new, $self->_hit_chooser.
=cut
sub _dbfile {
my $self = shift;
return $self->_blastdb->dbfile
}
=head2 _blast_obj
Args[1] :
Example :
Description:
Returntype :
Exceptions :
Caller :
=cut
sub _blast_obj {
my $self = shift;
if (@_){
$self->{_blast_obj} = shift;
return
}
unless ($self->{_blast_obj}){
# Create a new blast object with our mixed species input database.
$self->{_blast_obj}
= Bio::EnsEMBL::Pipeline::Runnable::MinimalBlast->new(
-program => $self->_blast_program,
-blastdb => $self->_blastdb,
-queryseq => $self->_query_seq,
-options => '-hspmax 200',
-workdir => $self->_work_dir,
-identity_cutoff => $self->_identity_cutoff);
}
return $self->{_blast_obj};
}
=head2 _blastdb
Args[1] :
Example :
Description:
Returntype :
Exceptions :
Caller :
=cut
sub _blastdb {
my $self = shift;
if (@_) {
$self->{_blastdb} = shift;
throw ("Blast database must be a Bio::EnsEMBL::Pipeline::Runnable::BlastDB.")
unless ($self->{_blastdb}->isa("Bio::EnsEMBL::Pipeline::Runnable::BlastDB"));
throw ("Blast database has not been formatted.")
unless $self->{_blastdb}->db_formatted;
throw ("Blast database has been built without the " .
"make_fetchable_index flag set (and this is " .
"a problem because the database can not be " .
"used for sequence fetching).")
unless $self->{_blastdb}->make_fetchable_index
}
throw ("Blast database object not set.")
unless ($self->{_blastdb});
return $self->{_blastdb};
}
=head2 _blast_program
Args[1] :
Example :
Description:
Returntype :
Exceptions :
Caller :
=cut
sub _blast_program {
my $self = shift;
if (@_){
$self->{_blast_program} = shift;
return
}
$self->{_blast_program} = $DEFAULT_BLAST_PROGRAM
unless $self->{_blast_program};
return $self->{_blast_program};
}
=head2 _blast_index_type
Args[1] :
Example :
Description:
Returntype :
Exceptions :
Caller :
=cut
sub _blast_index_type {
my $self = shift;
return $self->_blastdb->index_type;
}
### Getter/Setters ###
=head2 _query_seq
Args : (optional) Bio::Seq
Example : none
Description: Holds the query sequence for which we are searching for duplicates.
Returntype : Bio::Seq
Exceptions : none
Caller : $self->run
=cut
sub _query_seq {
my $self = shift;
if (@_) {
$self->{_query_seq} = shift;
throw ("Query sequence is not a Bio::Seq object [" .
$self->{_query_seq} . "]")
unless $self->{_query_seq}->isa("Bio::Seq");
throw ("Cant add query sequence to sequence cache manually.")
unless $self->_force_cache($self->{_query_seq});
return
}
throw ("Query sequence has not been set.")
unless $self->{_query_seq};
return $self->{_query_seq};
}
=head2 _work_dir
Args : (optional) String.
Example : none
Description: Holds the path to the working directory.
Returntype : String.
Exceptions : none
Caller : $self->new, $self->_hit_chooser, $self->_run_pairwise_paml.
=cut
sub _work_dir {
my $self = shift;
if (@_) {
$self->{_work_dir} = shift;
return
}
throw ("Work directory not set.")
unless $self->{_work_dir};
return $self->{_work_dir};
}
=head2 _identity_cutoff
Args : (optional) an int or a float - a percentage value anyways.
Example : none
Description: Holds identity cutoff percentage value.
Returntype : A float value
Exceptions : none
Caller : $self->new, $self->_hit_chooser.
=cut
sub _identity_cutoff {
my $self = shift;
if (@_) {
$self->{_identity_cutoff} = shift;
$self->{_identity_cutoff} /= 100 if $self->{_identity_cutoff} > 1;
return
}
throw ("Blast match identity cutoff has not been set.")
unless $self->{_identity_cutoff};
return $self->{_identity_cutoff};
}
=head2 _coverage_cutoff
Args : (optional) an int or a float - a percentage value anyways.
Example : none
Description: Holds coverage cutoff percentage value.
Returntype : A float value.
Exceptions : none
Caller : $self->new, $self->_hit_chooser.
=cut
sub _coverage_cutoff {
my $self = shift;
if (@_) {
$self->{_coverage_cutoff} = shift;
$self->{_coverage_cutoff} /= 100 if $self->{_coverage_cutoff} > 1;
return
}
throw ("The coverage cutoff has not been set.")
unless $self->{_coverage_cutoff};
return $self->{_coverage_cutoff};
}
=head2 _distance_cutoff
Args[1] :
Example :
Description:
Returntype :
Exceptions :
Caller :
=cut
sub _distance_cutoff {
my $self = shift;
if (@_) {
$self->{_distance_cutoff} = shift;
return
}
$self->{_distance_cutoff} = $DEFAULT_DISTANCE_CUTOFF
unless $self->{_distance_cutoff};
return $self->{_distance_cutoff};
}
=head2 _regex_query_species
Args : String
Example : $self->_regex_query_species('ENSG')
Description: Holds a regex string that will match the id of
any sequence from the query species.
Returntype : String or 1
Exceptions : Throws when regex is not set.
Caller : $self->new, $self->_hit_chooser
=cut
sub _regex_query_species {
my $self = shift;
if (@_) {
$self->{_regex_query_species} = shift;
return
}
throw ("The regular expression used to match the sequence"
." ids from the query species has not been set.")
unless $self->{_regex_query_species};
return $self->{_regex_query_species};
}
=head2 _regex_outgroup_species
Args : ref to an array of strings
Example : $self->_regex_outgroup_species(['ENSRNO', 'ENSMUS']);
Description: Holds an array of regexs that will allow the sequence
id of all non-query species sequences to be matched.
Returntype : arrayref
Exceptions : Warns if called while unset.
Caller : $self->new, $self->_hit_chooser
=cut
sub _regex_outgroup_species {
my $self = shift;
if (@_) {
$self->{_regex_outgroup_species} = shift;
return
}
warning('No outgroup species regex provided. ' .
'This may or may not be what you intend.')
unless $self->{_regex_outgroup_species};
return $self->{_regex_outgroup_species};
}
=head2 _genetic_code
Args : int
Example : $self->_genetic_code(1);
Description: Holds an integer representing the genetic code. To
choose the correct integer consult the documentation
used by the Bio::Seq->translate method. 1 is universal,
2 is vertebrate mitochondria.
Returntype : int
Exceptions : Warns if called while unset.
Caller : $self->new, $self->run, $self->_run_pairwise_paml
=cut
sub _genetic_code {
my $self = shift;
if (@_) {
$self->{_genetic_code} = shift;
return
}
unless (defined $self->{_genetic_code}) {
throw ('Genetic code unset.')
}
return $self->{_genetic_code};
}
=head2 _distance_method
Args : String
Example :
Description:
Returntype :
Exceptions : Throws if set to an unrecognised string.
Caller :
=cut
sub _distance_method {
my $self = shift;
if (@_) {
$self->{_distance_method} = shift;
unless ($self->{_distance_method} =~ /NeiGojobori/i |
$self->{_distance_method} =~ /ML/i){
throw ("Distance method must be set to either " .
"NeiGojobori or ML, not [".
$self->{_distance_method}."]");
}
}
$self->{_distance_method} = $DEFAULT_DISTANCE_METHOD
unless $self->{_distance_method};
return $self->{_distance_method};
}
=head2 _codeml
Args : [optional] String
Example : $self->_codeml('/path/to/codeml')
Description: Holds the path to the codeml executable
Returntype : String
Exceptions : Throws if a full path is included, but the
file is not executable.
Caller : $self->new, $self->_run_pairwise_paml
=cut
sub _codeml {
my $self = shift;
if (@_) {
$self->{_codeml} = shift;
return
}
$self->{_codeml} = 'codeml'
unless $self->{_codeml};
# If it looks like our executable comes with a full
# path, check that it will work.
throw ("codeml executable not found or not " .
"executable. Trying to execute path " .
"[". $self->{_codeml} ."]")
if ($self->{_codeml} =~ /^\//
& !-x $self->{_codeml});
return $self->{_codeml};
}
return 1
]]>