[Bioperl-l] Error calling alignment method
Lorenzo Carretero Paulet
locarpau at upvnet.upv.es
Wed May 25 18:41:22 UTC 2011
El 25/05/11 19:58, Lorenzo Carretero Paulet escribió:
>
>
> Hi all,
> I'm trying to create the following subroutine but I get the error
> message:
> "Can't call method "alignment" on an undefined value at
> /Users/marioafares/Documents/workspace/PlantEvolGen/test.pl line 80,
> <GEN1> line 2."
> which indicates $dna_aln is undefined. However, I manage to print it
> using Dumper, so I guess it is actually defined.
> Can anyone see where the problem is?
> I attach the two files I'm using in the test.
> Thanks in advance,
> Lorenzo
>
> #!/usr/local/bin/perl -w
> use 5.010;
> use strict;
> use Data::Dumper;
> # definition of the environmental variable CLUSTALDIR
> BEGIN {$ENV{CLUSTALDIR} =
> '/Applications/Bioinformatics/clustalw-2.1-macosx/'}
> use Bio::Tools::Run::Alignment::Clustalw;
> use Bio::Align::Utilities qw(aa_to_dna_aln);
>
> BEGIN {$ENV{PAMLDIR} = '/Applications/Bioinformatics/paml44/bin/'}
> use Bio::Tools::Run::Phylo::PAML::Codeml;
>
> my $sequencesfilenameAA =
>
> "/Users/marioafares/Documents/SequenceDatabase/plaza_public_02_Apr27/plaza_public_02/BLAST_Parsed_results/PerSpecies/test_vs_test.par.aa.1.fas";
>
> my $sequencesfilenameNT =
>
> "/Users/marioafares/Documents/SequenceDatabase/plaza_public_02_Apr27/plaza_public_02/BLAST_Parsed_results/PerSpecies/test_vs_test.par.nt.1.fas";
>
> my $format = 'fasta';
>
> GettingBioperlAlignmentAAtoDNAplusPAMLcalculation
> ($sequencesfilenameAA, $sequencesfilenameNT, $format);
>
>
> sub GettingBioperlAlignmentAAtoDNAplusPAMLcalculation
> {
> my ( $sequencesfilenameAA, $sequencesfilenameNT, $format, $ktuple,
> $matrix ) = @_;
> my %hashNTseqs = ();
> my $likelihood;
> my $Ks;
> my $Ka;
> my $omega;
> my $factory = Bio::Tools::Run::Alignment::Clustalw->new (); #use
> default parameters in here
> my $pep_aln = $factory->align($sequencesfilenameAA);
>
>
> my $inseq = Bio::SeqIO->new(-file => "<$sequencesfilenameNT",
> -format => $format );
> while (my $seq = $inseq->next_seq)
> {
> my $seq_id = $seq->display_id();
> $hashNTseqs{$seq_id} = $seq;
> }
> my $dna_aln = aa_to_dna_aln($pep_aln, \%hashNTseqs);
> say Dumper \%hashNTseqs;
> say Dumper $dna_aln;
> ###########################################################
> #my $codeml_runs = 5;
> my $ks_factory = new Bio::Tools::Run::Phylo::PAML::Codeml
> ( -params => {
> #'verbose' => 0,
> #'noisy' => 0,
> 'runmode' => -2,
> 'seqtype' => 1,
> #'model' => 0,
> #'NSsites' => 0,
> #'icode' => 0,
> #'fix_alpha' => 0,
> #'fix_kappa' => 0,
> #'RateAncestor' => 0,
> 'CodonFreq' => 2,
> 'cleandata' => 1,
> 'ndata' => 1
> },
>
> );
>
> my $kaks_factory->alignment($dna_aln);
>
> say "\nCalculating Ks-values of current cluster ...";
>
> # The estimation of the Ks-values is repeated $codeml_runs
> times ...
> # for(my $k=1;$k <= $codeml_runs;$k++)
> # {
> # print "\nCodeml-Run $k:\n\n";
>
> # Ka and Ks-vlaues are calculated using codeml
> my ($rc,$parser) = $kaks_factory->run();
> #$kaks_factory->cleanup();
> # If the calculation was succsessful ...
> # if($rc != 0)
> # {
> my $result = $parser->next_result;#not sure what it
> does here
> #my $NGmatrix = $result->get_MLmatrix();
> my $MLmatrix = $result->get_MLmatrix();
> $likelihood = $MLmatrix->[0]->[1]->{'lnL'};
> $Ks = $MLmatrix->[0]->[1]->{'dS'};
> $Ka = $MLmatrix->[0]->[1]->{'dA'};
> $omega = $MLmatrix->[0]->[1]->{'omega'};
> print " likelihood = $likelihood, Ka = $Ka, Ks =
> $Ks, Ka/Ks = $omega\n";
> # }
> # else # If an error occured during the Ks-calculation ...
> # {
> # return (-1);
> # }
> # }
> return ( $likelihood, $Ks, $Ka, $omega );
> }
>
>
>
> _______________________________________________
> Bioperl-l mailing list
> Bioperl-l at lists.open-bio.org
> http://lists.open-bio.org/mailman/listinfo/bioperl-l
Similarly, when I try to run the subroutine from a manually created
codon alignment as:
my $alignment =
"/Users/marioafares/Documents/SequenceDatabase/plaza_public_02_Apr27/plaza_public_02/BLAST_Parsed_results/PerSpecies/Alyrata_vs_Alyrata.par.cds.aln.fas";
my $format = 'fasta';
#GettingBioperlAlignmentAAtoDNAplusPAMLcalculation
($sequencesfilenameAA, $sequencesfilenameNT, $format);
GetPAMLcalculation($alignment, $format);
sub GetPAMLcalculation
{
my ( $alignment, $format ) = @_;
my $alignio = Bio::AlignIO->new(-format => $format,
-file => $alignment);
my $dna_aln = $alignio->next_aln;
my $kaks_factory = new Bio::Tools::Run::Phylo::PAML::Codeml
( -params => {
#'verbose' => 0,
#'noisy' => 0,
'runmode' => -2,
'seqtype' => 1,
#'model' => 0,
#'NSsites' => 0,
#'icode' => 0,
#'fix_alpha' => 0,
#'fix_kappa' => 0,
#'RateAncestor' => 0,
'CodonFreq' => 2,
'cleandata' => 1,
'ndata' => 1
},
);
$kaks_factory->alignment($dna_aln);
say "\nCalculating Ks-values of current cluster ...";
# The estimation of the Ks-values is repeated $codeml_runs times ...
# for(my $k=1;$k <= $codeml_runs;$k++)
# {
# print "\nCodeml-Run $k:\n\n";
# Ka and Ks-vlaues are calculated using codeml
my ($rc,$parser) = $kaks_factory->run();
#$kaks_factory->cleanup();
# If the calculation was succsessful ...
# if($rc != 0)
# {
my $result = $parser->next_result;#not sure what it does here
#my $NGmatrix = $result->get_MLmatrix();
my $MLmatrix = $result->get_MLmatrix();
my @otus = $result->get_seqs();
my $likelihood = $MLmatrix->[0]->[1]->{'lnL'};
my $Ks = $MLmatrix->[0]->[1]->{'dS'};
my $Ka = $MLmatrix->[0]->[1]->{'dA'};
my $omega = $MLmatrix->[0]->[1]->{'omega'};
print " likelihood = $likelihood, Ka = $Ka, Ks = $Ks,
Ka/Ks = $omega\n";
# }
# else # If an error occured during the Ks-calculation ...
# {
# return (-1);
# }
# }
return ( $likelihood, $Ks, $Ka, $omega );
}
I also get the following error message, which seems to indicate that the
alignment is not recognized as seqtype 1 by codeml (i.e. codon alignment).
--------------------- WARNING ---------------------
MSG: There was an error - see error_string for the program output
---------------------------------------------------
------------- EXCEPTION Bio::Root::NotImplemented -------------
MSG: Unknown format of PAML output did not see seqtype
STACK Bio::Tools::Phylo::PAML::_parse_summary
/Library/Perl//5.10.0/Bio/Tools/Phylo/PAML.pm:461
STACK Bio::Tools::Phylo::PAML::next_result
/Library/Perl//5.10.0/Bio/Tools/Phylo/PAML.pm:270
STACK main::GetPAMLcalculation
/Users/marioafares/Documents/workspace/PlantEvolGen/test.pl:78
STACK toplevel
/Users/marioafares/Documents/workspace/PlantEvolGen/test.pl:36
---------------------------------------------------------------
--
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
Lorenzo Carretero Paulet
Institute for Plant Molecular and Cell Biology - IBMCP (CSIC-UPV)
Integrative Systems Biology Group
C/ Ingeniero Fausto Elio s/n.
46022 Valencia, Spain
Phone: +34 963879934
Fax: +34 963877859
e-mail: locarpau at upvnet.upv.es
*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*-*
More information about the Bioperl-l
mailing list