[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