Applying a discrete-time, non-stationary nucleotide model

We fit a discrete-time Markov nucleotide model. This corresponds to a Barry and Hartigan 1987 model.

from cogent3 import get_app

loader = get_app("load_aligned", format_name="fasta", moltype="dna")
aln = loader("data/primate_brca1.fasta")
model = get_app("model", "BH", tree="data/primate_brca1.tree")
result = model(aln)
result

BH
key	lnL	nfp	DLC	unique_Q
'BH'	-6941.6028	132	True	True

Note

DLC stands for diagonal largest in column and the value is a check on the identifiability of the model. unique_Q is another identifiability check, but it not applicable to a discrete-time model and so remains as None.

Looking at the likelihood function, we see these maximum likelihood estimated values

result.lf

BH

log-likelihood = -6941.6028

number of free parameters = 132

Edge motif motif2 params
edge	motif	motif2	psubs
Galago	T	T	0.88
Galago	T	C	0.06
Galago	T	A	0.04
Galago	T	G	0.02
Galago	C	T	0.11
...	...	...	...
edge.3	A	G	0.01
edge.3	G	T	0.00
edge.3	G	C	0.00
edge.3	G	A	0.00
edge.3	G	G	0.99

Motif params
A	C	G	T
0.38	0.17	0.21	0.24

Get a tree with branch lengths as paralinear

This is the only possible length metric for a discrete-time process.

tree = result.tree
fig = tree.get_figure()
fig.scale_bar = "top right"
fig.show(width=500, height=500)

Getting parameter estimates

For a discrete-time model, aside from the root motif probabilities, everything is edge specific. But note that the tabular_result has different keys from the continuous-time case, as demonstrated below.

tabulator = get_app("tabulate_stats")
stats = tabulator(result)
stats

2x tabular_result('edge motif motif2 params': Table, 'motif params': Table)

stats["edge motif motif2 params"]

edge motif motif2 params
edge	motif	motif2	psubs
Galago	T	T	0.88
Galago	T	C	0.06
Galago	T	A	0.04
Galago	T	G	0.02
Galago	C	T	0.11
...	...	...	...
edge.3	A	G	0.01
edge.3	G	T	0.00
edge.3	G	C	0.00
edge.3	G	A	0.00
edge.3	G	G	0.99

Top 5 and bottom 5 rows from 176 rows x 4 columns