[Bioperl-l] PAML problem

江 文恺 biology0046 at hotmail.com
Sun Nov 26 13:19:45 UTC 2006


I use scripts at PAML Howto page to caculate M7 or M8 model. but it always 
report errors:
-------------------------------------------------
Error: err: incorrect model for pairwise comparison.check NSsites, alpha, 
aaDist..
-------------------------------------------------
Can't call method "next_result" on an undefined value at pamtest.txt line 
65.
---------------------------------------------------------- 
use Bio::Tools::Run::Phylo::PAML::Codeml;
use Bio::Tools::Run::Alignment::Clustalw;
 
# for projecting alignments from protein to R/DNA space
use Bio::Align::Utilities qw(aa_to_dna_aln);
# for input of the sequence data
use Bio::SeqIO;
use Bio::AlignIO;
 
my $aln_factory = Bio::Tools::Run::Alignment::Clustalw->new;
#my $seqdata = shift || 'cds.fa';
 
my $seqio = new Bio::SeqIO(-file   =>"CG4686.cds",
                           -format => 'fasta');
my %seqs;
my @prots;
# process each sequence
while ( my $seq = $seqio->next_seq ) {
    $seqs{$seq->display_id} = $seq;
    # translate them into protein
    my $protein = $seq->translate();
    my $pseq = $protein->seq();
    if( $pseq =~ /\*/ &&
        $pseq !~ /\*$/ ) {
          warn("provided a CDS sequence with a stop codon, PAML will 
choke!");
          #exit(0);
    }
    # Tcoffee can't handle '*' even if it is trailing
    $pseq =~ s/\*//g;
    $protein->seq($pseq);
    push @prots, $protein;
}
 
if( @prots < 2 ) {
    warn("Need at least 2 CDS sequences to proceed");
    exit(0);
}
 
open(OUT, ">align_output.txt") ||  die("cannot open output align_output for 
writing");
# Align the sequences with clustalw
my $aa_aln = $aln_factory->align(\@prots);
# project the protein alignment back to CDS coordinates
my $dna_aln = aa_to_dna_aln($aa_aln, \%seqs);
 
my @each = $dna_aln->each_seq();

my $kaks_factory = Bio::Tools::Run::Phylo::PAML::Codeml->new
                   ( -params => { 'NSsites' => 7, 'ncatG'=>10,
                                  'seqtype' => 1,
                                } ); #I change the parameters here
 
# set the alignment object
$kaks_factory->alignment($dna_aln);
 
# run the KaKs analysis
my ($rc,$parser) = $kaks_factory->run();
my $result = $parser->next_result;
Seqfile:>D_pseudoobscuraATGATGGACACATTGGACTACATCACCCTCAGCAATCCCGTGAGCAAGGCGGTTATATATTCGGGATCGGCTATCTTCCGAGCCCTTGGGCTGCGTCCGAAACAGTTGGTTCCGAAGGAGACCCAGACAGTGCGTCCGGTAATGACTCAGTACATGTCGCCGGGTGCTGGCTCGCTGCATAATTTGGCCGGCCGGTACTATCATTTCGTGCGACTGGCTGGCCTAGGCGGCGCTTCGGCTATTTTCATGGGAGCCTACTGCAAATACGTCCTGAAGGAGATCAAGAACGAAAAGGAGCAGCTGGACTCGCAGGCCTTTGCCGATGTCGCCAATCGCATACACTTTTTGCATTCGTTTGCACTGATGGCCATGCCACTGGCCCACTATCCCATAATCACGGGCTCTTTGATGACCACTGGCACACTGCTCTTCAGCGGCTGCATGTACTATCGTGCATTGACGGGCGAGAAGCGCTTCCAGCCGTTTGCCACCGTTGGAGGCTTCTGTCTGATAGCCGCTTGGCTGACGCTGATATTT

>D_willistoniATGTCAATAGTTGATACATTGGACTATATAACCCTGGGTAATCCAGTGAGTCGCTTGGTCATATCTTCAACATCGGCATTAATGCGTACAATTGGTCTGCGCCCCAAGCAGGTGCCTATAAAGGATACAGAAATAGCTGGATTGCCACAGCAGGTCCATCATTTCGGAAATCCAAATGGCCCATCCTTATATACAATTGCCGGCTCTCATTATAACTTTATTCGTCTAGCTGGAATTGGCGGAGCATCGGCTATATTTATGGGAGCGTATTGTAAATATTTTCTTAAGGATATCAATGATCCCAAGGAACAGTTGGATTCACAAGCCTTTGCCGATGTGGCCAATCGTATACATTTTCTCCACTCATTTGCATTAATGGCGATGCCCTTGGCTCATTATCCAGTTTTTACTGGCGCCCTTATGACCACCGGCACATTACTTTTCAGTGGGTGTATGTACTATCGCGCTTTAACCGGTGATAAGAGACTTCAACATTATGCCACAATTGGTGGCTTCTGTCTAATGGCTGCTTGGTTATCCTTGGTTTTG

>D_melanogasterATGTCAGTCGCTGACACCATTGAATACGTGACCCTGGGCAATCCTGTCAGCAAGATGGTAGCCTCGTCCGCATCCGCCCTGCTCCGCACGCTTGGTCTGCGTCCCAAGAAGGTGCCGGTGCAGGAGACGAGTATGGCGGTGCTCCCTGCCGCCCACAGCTACGCGCATTCACACGGATCACTCTATCGACTGGCCGGCTGCCATTACCACTTCATTCGGCTGGCCGGGATCGTCGGCGCGTCGGCCATCTTTATGGGCGCCTACTGCAAGTACGTCCTGAAGGACGTCAGCGATCCCAAGGAGCAGGTGGACTCGCAGGCCTTTGCTGATGTGGCCAATCGCATCCACTTTCTGCACTCCTTTGCCATGATGGCCATGCCTCTGGCCCACTATCCCGTATTCACTGGCACTTTGATGATTACGGGCATGATGCTTTTCAGCGGCTGCATGTACTACCGCGCTTTGACTGGCGAGAAGCGTCTGCAACCGTACGCGACCGTCGGAGGATTCTGCCTGATGGCCGCGTGGCTGTCGCTGGTCCTG


_________________________________________________________________
与联机的朋友进行交流,请使用 MSN Messenger:  http://messenger.msn.com/cn  




More information about the Bioperl-l mailing list