Skip to content

Commit

Permalink
add support for parsing PAML/codeml M2a_rel model results
Browse files Browse the repository at this point in the history
  • Loading branch information
@brandoninvergo @peterjc
brandoninvergo authored and peterjc committed Apr 6, 2016
1 parent ead7c44 commit 4031b63
Show file tree
Hide file tree
Showing 8 changed files with 635 additions and 8 deletions.
5 changes: 3 additions & 2 deletions Bio/Phylo/PAML/_parse_codeml.py
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
# Copyright (C) 2011 by Brandon Invergo (b.invergo@gmail.com)
# Copyright (C) 2011, 2016 by Brandon Invergo (b.invergo@gmail.com)
# This code is part of the Biopython distribution and governed by its
# license. Please see the LICENSE file that should have been included
# as part of this package.
Expand Down Expand Up @@ -105,7 +105,8 @@ def parse_nssites(lines, results, multi_models, multi_genes):
"PositiveSelection": 2,
"discrete": 3,
"beta": 7,
"beta&w>1": 8}[siteclass_model]
"beta&w>1": 8,
"M2a_rel": 22}[siteclass_model]
if multi_genes:
genes = results["genes"]
current_gene = None
Expand Down
25 changes: 25 additions & 0 deletions Tests/PAML/Control_files/codeml/m2a_rel.ctl
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,25 @@
seqfile = ../alignment.phylip
outfile = temp.out
treefile = ../species.tree
verbose = 0
CodonFreq = 2
cleandata = 1
fix_blength = 1
NSsites = 22
fix_omega = 0
clock = 0
runmode = 0
fix_kappa = 0
fix_alpha = 1
alpha = 0
method = 0
Malpha = 0
RateAncestor = 0
icode = 0
seqtype = 1
omega = 1
getSE = 0
noisy = 0
Mgene = 0
kappa = 4.54006
model = 0
145 changes: 145 additions & 0 deletions Tests/PAML/Results/codeml/m2a_rel/m2a_rel-4_6.out
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,145 @@
CODONML (in paml version 4.6, August 2012) ../Alignments/alignment.phylip
Model: One dN/dS ratio for branches,
Codon frequency model: F3x4
Site-class models: M2a_rel (3 categories)
ns = 5 ls = 74

Codon usage in sequences
--------------------------------------------------------------------------------------------------
Phe TTT 0 0 0 0 0 | Ser TCT 1 1 1 1 1 | Tyr TAT 0 0 1 0 0 | Cys TGT 3 2 3 3 3
TTC 1 1 1 1 1 | TCC 1 1 1 1 1 | TAC 1 1 0 1 1 | TGC 0 1 0 0 0
Leu TTA 0 0 0 0 0 | TCA 1 1 1 1 1 | *** TAA 0 0 0 0 0 | *** TGA 0 0 0 0 0
TTG 1 1 1 1 1 | TCG 0 0 0 0 0 | TAG 0 0 0 0 0 | Trp TGG 0 0 0 0 0
--------------------------------------------------------------------------------------------------
Leu CTT 0 0 0 0 0 | Pro CCT 0 0 0 0 0 | His CAT 0 0 0 0 0 | Arg CGT 0 0 0 0 0
CTC 2 2 2 2 2 | CCC 1 1 1 1 1 | CAC 0 0 0 0 0 | CGC 0 0 0 0 0
CTA 1 1 1 1 1 | CCA 3 3 3 3 3 | Gln CAA 0 0 0 0 0 | CGA 1 1 1 1 1
CTG 3 3 3 3 3 | CCG 0 0 0 0 0 | CAG 1 1 1 1 1 | CGG 0 0 0 0 0
--------------------------------------------------------------------------------------------------
Ile ATT 2 2 2 2 2 | Thr ACT 0 0 0 0 0 | Asn AAT 2 2 2 2 2 | Ser AGT 0 0 0 0 0
ATC 2 2 2 2 2 | ACC 0 0 0 0 0 | AAC 0 0 0 0 0 | AGC 0 0 0 0 0
ATA 0 0 0 0 0 | ACA 2 2 2 2 2 | Lys AAA 5 5 5 5 5 | Arg AGA 2 2 2 2 2
Met ATG 3 3 3 3 3 | ACG 0 0 0 0 0 | AAG 5 5 5 5 5 | AGG 0 0 0 0 0
--------------------------------------------------------------------------------------------------
Val GTT 3 3 3 3 1 | Ala GCT 0 0 0 0 0 | Asp GAT 2 2 2 2 2 | Gly GGT 0 0 0 0 0
GTC 0 0 0 0 1 | GCC 0 0 0 0 0 | GAC 4 4 4 4 4 | GGC 3 3 3 3 3
GTA 3 3 3 3 4 | GCA 0 0 0 0 0 | Glu GAA 8 8 8 8 8 | GGA 1 1 1 1 1
GTG 2 2 2 2 2 | GCG 0 0 0 0 0 | GAG 4 4 4 4 4 | GGG 0 0 0 0 0
--------------------------------------------------------------------------------------------------

Codon position x base (3x4) table for each sequence.

#1: Homo_sapie
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.17568 C:0.20270 A:0.36486 G:0.25676
Average T:0.20270 C:0.16216 A:0.36937 G:0.26577

#2: Pan_troglo
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.16216 C:0.21622 A:0.36486 G:0.25676
Average T:0.19820 C:0.16667 A:0.36937 G:0.26577

#3: Gorilla_go
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.18919 C:0.18919 A:0.36486 G:0.25676
Average T:0.20721 C:0.15766 A:0.36937 G:0.26577

#4: Pongo_pygm
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.17568 C:0.20270 A:0.36486 G:0.25676
Average T:0.20270 C:0.16216 A:0.36937 G:0.26577

#5: Macaca_mul
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.14865 C:0.21622 A:0.37838 G:0.25676
Average T:0.19369 C:0.16667 A:0.37387 G:0.26577

Sums of codon usage counts
------------------------------------------------------------------------------
Phe F TTT 0 | Ser S TCT 5 | Tyr Y TAT 1 | Cys C TGT 14
TTC 5 | TCC 5 | TAC 4 | TGC 1
Leu L TTA 0 | TCA 5 | *** * TAA 0 | *** * TGA 0
TTG 5 | TCG 0 | TAG 0 | Trp W TGG 0
------------------------------------------------------------------------------
Leu L CTT 0 | Pro P CCT 0 | His H CAT 0 | Arg R CGT 0
CTC 10 | CCC 5 | CAC 0 | CGC 0
CTA 5 | CCA 15 | Gln Q CAA 0 | CGA 5
CTG 15 | CCG 0 | CAG 5 | CGG 0
------------------------------------------------------------------------------
Ile I ATT 10 | Thr T ACT 0 | Asn N AAT 10 | Ser S AGT 0
ATC 10 | ACC 0 | AAC 0 | AGC 0
ATA 0 | ACA 10 | Lys K AAA 25 | Arg R AGA 10
Met M ATG 15 | ACG 0 | AAG 25 | AGG 0
------------------------------------------------------------------------------
Val V GTT 13 | Ala A GCT 0 | Asp D GAT 10 | Gly G GGT 0
GTC 1 | GCC 0 | GAC 20 | GGC 15
GTA 16 | GCA 0 | Glu E GAA 40 | GGA 5
GTG 10 | GCG 0 | GAG 20 | GGG 0
------------------------------------------------------------------------------


Codon position x base (3x4) table, overall

position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.17027 C:0.20541 A:0.36757 G:0.25676
Average T:0.20090 C:0.16306 A:0.37027 G:0.26577


Nei & Gojobori 1986. dN/dS (dN, dS)
(Note: This matrix is not used in later ML. analysis.
Use runmode = -2 for ML pairwise comparison.)

Homo_sapie
Pan_troglo -1.0000 (0.0000 0.0207)
Gorilla_go -1.0000 (0.0000 0.0207)-1.0000 (0.0000 0.0421)
Pongo_pygm -1.0000 (0.0000 0.0000)-1.0000 (0.0000 0.0207)-1.0000 (0.0000 0.0207)
Macaca_mul -1.0000 (0.0000 0.0421)-1.0000 (0.0000 0.0640)-1.0000 (0.0000 0.0640)-1.0000 (0.0000 0.0421)


TREE # 1: (((1, 2), 3), 4, 5); MP score: 4
check convergence..
lnL(ntime: 7 np: 12): -308.031579 +0.000000
6..7 7..8 8..1 8..2 7..3 6..4 6..5
0.000004 0.000004 0.000004 0.013438 0.013431 0.000004 0.027839 1.518248 0.999999 0.000000 0.000001 0.000001

Note: Branch length is defined as number of nucleotide substitutions per codon (not per neucleotide site).

tree length = 0.05472

(((1: 0.000004, 2: 0.013438): 0.000004, 3: 0.013431): 0.000004, 4: 0.000004, 5: 0.027839);

(((Homo_sapie: 0.000004, Pan_troglo: 0.013438): 0.000004, Gorilla_go: 0.013431): 0.000004, Pongo_pygm: 0.000004, Macaca_mul: 0.027839);

Detailed output identifying parameters

kappa (ts/tv) = 1.51825

Parameters in M22 (M2a_rel):


dN/dS (w) for site classes (K=3)

p: 1.00000 0.00000 0.00000
w: 0.00000 1.00000 0.00000

dN & dS for each branch

branch t N S dN/dS dN dS N*dN S*dS

6..7 0.000 167.7 54.3 0.0000 0.0000 0.0000 0.0 0.0
7..8 0.000 167.7 54.3 0.0000 0.0000 0.0000 0.0 0.0
8..1 0.000 167.7 54.3 0.0000 0.0000 0.0000 0.0 0.0
8..2 0.013 167.7 54.3 0.0000 0.0000 0.0183 0.0 1.0
7..3 0.013 167.7 54.3 0.0000 0.0000 0.0183 0.0 1.0
6..4 0.000 167.7 54.3 0.0000 0.0000 0.0000 0.0 0.0
6..5 0.028 167.7 54.3 0.0000 0.0000 0.0379 0.0 2.1


Naive Empirical Bayes (NEB) analysis
Time used: 0:05
145 changes: 145 additions & 0 deletions Tests/PAML/Results/codeml/m2a_rel/m2a_rel-4_7.out
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,145 @@
CODONML (in paml version 4.7, January 2013) ../Alignments/alignment.phylip
Model: One dN/dS ratio for branches,
Codon frequency model: F3x4
Site-class models: M2a_rel (3 categories)
ns = 5 ls = 74

Codon usage in sequences
--------------------------------------------------------------------------------------------------
Phe TTT 0 0 0 0 0 | Ser TCT 1 1 1 1 1 | Tyr TAT 0 0 1 0 0 | Cys TGT 3 2 3 3 3
TTC 1 1 1 1 1 | TCC 1 1 1 1 1 | TAC 1 1 0 1 1 | TGC 0 1 0 0 0
Leu TTA 0 0 0 0 0 | TCA 1 1 1 1 1 | *** TAA 0 0 0 0 0 | *** TGA 0 0 0 0 0
TTG 1 1 1 1 1 | TCG 0 0 0 0 0 | TAG 0 0 0 0 0 | Trp TGG 0 0 0 0 0
--------------------------------------------------------------------------------------------------
Leu CTT 0 0 0 0 0 | Pro CCT 0 0 0 0 0 | His CAT 0 0 0 0 0 | Arg CGT 0 0 0 0 0
CTC 2 2 2 2 2 | CCC 1 1 1 1 1 | CAC 0 0 0 0 0 | CGC 0 0 0 0 0
CTA 1 1 1 1 1 | CCA 3 3 3 3 3 | Gln CAA 0 0 0 0 0 | CGA 1 1 1 1 1
CTG 3 3 3 3 3 | CCG 0 0 0 0 0 | CAG 1 1 1 1 1 | CGG 0 0 0 0 0
--------------------------------------------------------------------------------------------------
Ile ATT 2 2 2 2 2 | Thr ACT 0 0 0 0 0 | Asn AAT 2 2 2 2 2 | Ser AGT 0 0 0 0 0
ATC 2 2 2 2 2 | ACC 0 0 0 0 0 | AAC 0 0 0 0 0 | AGC 0 0 0 0 0
ATA 0 0 0 0 0 | ACA 2 2 2 2 2 | Lys AAA 5 5 5 5 5 | Arg AGA 2 2 2 2 2
Met ATG 3 3 3 3 3 | ACG 0 0 0 0 0 | AAG 5 5 5 5 5 | AGG 0 0 0 0 0
--------------------------------------------------------------------------------------------------
Val GTT 3 3 3 3 1 | Ala GCT 0 0 0 0 0 | Asp GAT 2 2 2 2 2 | Gly GGT 0 0 0 0 0
GTC 0 0 0 0 1 | GCC 0 0 0 0 0 | GAC 4 4 4 4 4 | GGC 3 3 3 3 3
GTA 3 3 3 3 4 | GCA 0 0 0 0 0 | Glu GAA 8 8 8 8 8 | GGA 1 1 1 1 1
GTG 2 2 2 2 2 | GCG 0 0 0 0 0 | GAG 4 4 4 4 4 | GGG 0 0 0 0 0
--------------------------------------------------------------------------------------------------

Codon position x base (3x4) table for each sequence.

#1: Homo_sapie
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.17568 C:0.20270 A:0.36486 G:0.25676
Average T:0.20270 C:0.16216 A:0.36937 G:0.26577

#2: Pan_troglo
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.16216 C:0.21622 A:0.36486 G:0.25676
Average T:0.19820 C:0.16667 A:0.36937 G:0.26577

#3: Gorilla_go
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.18919 C:0.18919 A:0.36486 G:0.25676
Average T:0.20721 C:0.15766 A:0.36937 G:0.26577

#4: Pongo_pygm
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.17568 C:0.20270 A:0.36486 G:0.25676
Average T:0.20270 C:0.16216 A:0.36937 G:0.26577

#5: Macaca_mul
position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.14865 C:0.21622 A:0.37838 G:0.25676
Average T:0.19369 C:0.16667 A:0.37387 G:0.26577

Sums of codon usage counts
------------------------------------------------------------------------------
Phe F TTT 0 | Ser S TCT 5 | Tyr Y TAT 1 | Cys C TGT 14
TTC 5 | TCC 5 | TAC 4 | TGC 1
Leu L TTA 0 | TCA 5 | *** * TAA 0 | *** * TGA 0
TTG 5 | TCG 0 | TAG 0 | Trp W TGG 0
------------------------------------------------------------------------------
Leu L CTT 0 | Pro P CCT 0 | His H CAT 0 | Arg R CGT 0
CTC 10 | CCC 5 | CAC 0 | CGC 0
CTA 5 | CCA 15 | Gln Q CAA 0 | CGA 5
CTG 15 | CCG 0 | CAG 5 | CGG 0
------------------------------------------------------------------------------
Ile I ATT 10 | Thr T ACT 0 | Asn N AAT 10 | Ser S AGT 0
ATC 10 | ACC 0 | AAC 0 | AGC 0
ATA 0 | ACA 10 | Lys K AAA 25 | Arg R AGA 10
Met M ATG 15 | ACG 0 | AAG 25 | AGG 0
------------------------------------------------------------------------------
Val V GTT 13 | Ala A GCT 0 | Asp D GAT 10 | Gly G GGT 0
GTC 1 | GCC 0 | GAC 20 | GGC 15
GTA 16 | GCA 0 | Glu E GAA 40 | GGA 5
GTG 10 | GCG 0 | GAG 20 | GGG 0
------------------------------------------------------------------------------


Codon position x base (3x4) table, overall

position 1: T:0.12162 C:0.16216 A:0.31081 G:0.40541
position 2: T:0.31081 C:0.12162 A:0.43243 G:0.13514
position 3: T:0.17027 C:0.20541 A:0.36757 G:0.25676
Average T:0.20090 C:0.16306 A:0.37027 G:0.26577


Nei & Gojobori 1986. dN/dS (dN, dS)
(Note: This matrix is not used in later ML. analysis.
Use runmode = -2 for ML pairwise comparison.)

Homo_sapie
Pan_troglo -1.0000 (0.0000 0.0207)
Gorilla_go -1.0000 (0.0000 0.0207)-1.0000 (0.0000 0.0421)
Pongo_pygm -1.0000 (0.0000 0.0000)-1.0000 (0.0000 0.0207)-1.0000 (0.0000 0.0207)
Macaca_mul -1.0000 (0.0000 0.0421)-1.0000 (0.0000 0.0640)-1.0000 (0.0000 0.0640)-1.0000 (0.0000 0.0421)


TREE # 1: (((1, 2), 3), 4, 5); MP score: 4
check convergence..
lnL(ntime: 7 np: 12): -308.031579 +0.000000
6..7 7..8 8..1 8..2 7..3 6..4 6..5
0.000004 0.000004 0.000004 0.013438 0.013431 0.000004 0.027839 1.518244 1.000000 0.000000 0.000001 0.000001

Note: Branch length is defined as number of nucleotide substitutions per codon (not per neucleotide site).

tree length = 0.05472

(((1: 0.00000, 2: 0.01344): 0.00000, 3: 0.01343): 0.00000, 4: 0.00000, 5: 0.02784);

(((Homo_sapie: 0.00000, Pan_troglo: 0.01344): 0.00000, Gorilla_go: 0.01343): 0.00000, Pongo_pygm: 0.00000, Macaca_mul: 0.02784);

Detailed output identifying parameters

kappa (ts/tv) = 1.51824

Parameters in M22 (M2a_rel):


dN/dS (w) for site classes (K=3)

p: 1.00000 0.00000 0.00000
w: 0.00000 1.00000 0.00000

dN & dS for each branch

branch t N S dN/dS dN dS N*dN S*dS

6..7 0.000 167.7 54.3 0.0000 0.0000 0.0000 0.0 0.0
7..8 0.000 167.7 54.3 0.0000 0.0000 0.0000 0.0 0.0
8..1 0.000 167.7 54.3 0.0000 0.0000 0.0000 0.0 0.0
8..2 0.013 167.7 54.3 0.0000 0.0000 0.0183 0.0 1.0
7..3 0.013 167.7 54.3 0.0000 0.0000 0.0183 0.0 1.0
6..4 0.000 167.7 54.3 0.0000 0.0000 0.0000 0.0 0.0
6..5 0.028 167.7 54.3 0.0000 0.0000 0.0379 0.0 2.1


Naive Empirical Bayes (NEB) analysis
Time used: 0:05
Loading

0 comments on commit 4031b63

Please sign in to comment.