view release on metacpan or  search on metacpan

lib/Bio/Phylo/Models/Substitution/Dna.pm  view on Meta::CPAN

=item set_pi

Setter for base frequencies.

 Type    : method
 Title   : get_pi
 Usage   : $model->set_pi((0.2, 0.2, 0.3, 0.3));
 Function: Set base frequencies.
 Returns : A Bio::Phylo::Models::Substitution::Dna object.
 Args    : array of four base frequencies (A, C, G, T)
 Comments: Base frequencies must sum to one

=cut

sub set_pi {
    my ( $self, $pi ) = @_;
    ref $pi eq 'ARRAY' or throw 'BadArgs' => "Not an array ref!";
    my $total = 0;
    $total += $_ for @{$pi};
    my $epsilon = 0.000001;
    abs(1 - $total) < $epsilon or throw 'BadArgs' => 'Frequencies must sum to one';
    $self->{'_pi'} = $pi;
    return $self;
}

=item set_median

Setter for median for gamma-modeled rate categories.

 Type    : method
 Title   : set_median
 Usage   : $model->set_median(1);
 Function: Setter for median.
 Returns : A Bio::Phylo::Models::Substitution::Dna object.
 Args    : scalar

=cut

sub set_median {
    my $self = shift;
    $self->{'_median'} = !!shift;
    return $self;
}

=item modeltest

Performing a modeltest using the package 'phangorn' in
R (Schliep, Bioinformatics (2011) 27 (4): 592-593) from an
DNA alignment. If no tree is given as argument, a neighbor-joining
tree is generated from the alignment to perform model testing.
Selects the model with the minimum AIC.

 Type    : method
 Title   : modeltest
 Usage   : $model->modeltest(-matrix=>$matrix);
 Function: Determine DNA substitution model from alignment.
 Returns : An object which is subclass of Bio::Phylo::Models::Substitution::Dna.
 Args    : -matrix: A Bio::Phylo::Matrices::Matrix object
           Optional:
           -tree: A Bio::Phylo::Forest::Tree object
           -timeout: Timeout in seconds to prevent getting stuck in an R process.
 Comments: Prerequisites: Statistics::R, R, and the R package phangorn.

=cut

sub modeltest {
	my ($self, %args) = @_;

	my $matrix = $args{'-matrix'};
	my $tree = $args{'-tree'};
	my $timeout = $args{'-timeout'};

	my $model;

	if ( looks_like_class 'Statistics::R' ) {

		eval {
			# phangorn needs files as input
			my ($fasta_fh, $fasta) = tempfile();
			print $fasta_fh unparse('-phylo'=>$matrix, '-format'=>'fasta');
			close $fasta_fh;

			# instanciate R and lcheck if phangorn is installed
			my $R = Statistics::R->new;
			$R->timeout($timeout) if $timeout;
			$R->run(q[options(device=NULL)]);
			$R->run(q[package <- require("phangorn")]);

			if ( ! $R->get(q[package]) eq "TRUE") {
				$logger->warn("R library phangorn must be installed to run modeltest");
				return $model;
			}

			# read data
			$R->run(qq[data <- read.FASTA("$fasta")]);

			# remove temp file
			cleanup();

			if ( $tree ) {
				# make copy of tree since it will be pruned
				my $current_tree = parse('-format'=>'newick', '-string'=>$tree->to_newick)->first;
				# prune out taxa from tree that are not present in the data
				my @taxon_names = map {$_->get_name} @{ $matrix->get_entities };
				$logger->debug('pruning input tree');
				$current_tree->keep_tips(\@taxon_names);
				$logger->debug('pruned input tree: ' . $current_tree->to_newick);

				if ( ! $current_tree or scalar( @{ $current_tree->get_terminals } ) < 3 ) {
					$logger->warn('pruned tree has too few tip labels, determining substitution model using NJ tree');
					$R->run(q[test <- modelTest(phyDat(data))]);
				}
				else {
					my $newick = $current_tree->to_newick;

					$R->run(qq[tree <- read.tree(text="$newick")]);
					# call modelTest
					$logger->debug("calling modelTest from R package phangorn");
					$R->run(q[test <- modelTest(phyDat(data), tree=tree)]);
				}
			}