[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