Take to be the total instantaneous mutation rate away from a base, so each specified alternative has rate . Under the Jukes-Cantor model, the uniform component is stationary and the three contrast components decay as , givingAn implementation is:
def jukes(a, b, t, mu):
decay = np.exp(-4*mu*t/3)
return .25 + .75*decay if a == b else .25*(1-decay)
print(f"{jukes(2,2,1,.5):.3f} {jukes(2,3,1,.5):.3f}")The two probabilities are and . The probabilities in any row sum to one, at the transition matrix is the identity, and at long times every entry tends to .
The Jukes-Cantor model is also parameterized by a rate for each specified replacement, in which case the exponential is and . The total-rate convention is used consistently below; it yields an optimum within the stated search interval, whereas the corresponding per-replacement rate is about .
Initialize all eight leaves, not seven: every observed base contributes a one-hot vector in the order A, C, G, T. The reference to leaf rows in the prose is inconsistent with its eight-leaf tree.
def initialise_prob(leafList):
alphabet = {base: j for j, base in enumerate("ACGT")}
result = np.zeros((2*len(leafList)-1, 4), dtype=float)
for i, base in enumerate(leafList):
result[i, alphabet[base]] = 1.
return result
leaves = ["G", "A", "T", "G", "G", "G", "C", "G"]
prob = initialise_prob(leaves)
print(prob[:len(leaves)])The initialized leaf rows are:
0 0 1 0
1 0 0 0
0 0 0 1
0 0 1 0
0 0 1 0
0 0 1 0
0 1 0 0
0 0 1 0The remaining seven internal-node rows initially contain zeros. Their eventual entries are conditional likelihoods of descendant observations given the ancestral base, rather than normalized posterior probabilities of that ancestral base.
For node with children , the phylogenetic pruning algorithm combines the two conditionally independent descendant subtrees:The likelihood function is under a uniform ancestral prior. In the supplied numbering, children have smaller indices than their parents, so increasing node order is a valid bottom-up traversal:
def felsenstein(prob, mu, L, left, right, left_time, right_time):
# All four entries are overwritten at every internal node.
for node in range(L, 2*L-1):
Tl = np.array([[jukes(a,b,left_time[node],mu)
for b in range(4)] for a in range(4)])
Tr = np.array([[jukes(a,b,right_time[node],mu)
for b in range(4)] for a in range(4)])
prob[node] = (Tl @ prob[left[node]]) * (Tr @ prob[right[node]])
return float(np.log(prob[-1].mean()))
print(felsenstein(prob, .5, len(leaves), left, right, left_time, right_time))At , the log-likelihood is . Internal rows must remain unnormalized in this implementation, because their magnitudes contribute to the tree likelihood function; normalization without retaining the removed factors would lose that information. The final root normalization required for a posterior probability is a separate step in part (e).
Evaluate the log-likelihood function on the specified interval, retaining the rate with the largest value. A grid spacing of finds ; refining locally avoids relying on a coarse-grid location when reporting three significant figures:
def evaluate(mu):
p = initialise_prob(leaves)
return felsenstein(p, mu, len(leaves), left, right, left_time, right_time)
coarse = np.linspace(.1, 1., 901)
centre = coarse[np.argmax([evaluate(mu) for mu in coarse])]
fine = np.linspace(max(.1, centre-.001), min(1., centre+.001), 2001)
mu_max = float(fine[np.argmax([evaluate(mu) for mu in fine])])
print(f"mu = {mu_max:.3g}; log-likelihood = {evaluate(mu_max):.8f}")The maximum-likelihood estimate is to three significant figures, with refined value about and log-likelihood . Under the alternative per-replacement convention this same fitted transition matrix has ; these are equivalent parameterizations, not competing fitted trees.
Recompute the root conditional likelihoods at the fitted rate, then multiply by the uniform root prior and normalize. The factors cancel between numerator and denominator:
prob = initialise_prob(leaves)
felsenstein(prob, mu_max, len(leaves), left, right, left_time, right_time)
root_posterior = prob[-1] / prob[-1].sum()
print(dict(zip("ACGT", root_posterior)))At the refined optimum, the posterior distribution isThus G is the most probable ancestral base, but it is not inferred with certainty. This conditions on the fitted rate and does not additionally average over uncertainty in that rate.
Articles by others on the same topic
There are currently no matching articles.