Create a gist now

Instantly share code, notes, and snippets.

What would you like to do?
Python function for a hacky approximation of the F_Tt factors in the Bryant et al SNAPP method as specialized for use by a homework assignment in the Kelly and Holder "Likelihood methods in Biology" class
def calc_F_Tt(r_Tb, n_Tt, r_Tt, param_vec):
if r_Tt > n_Tt:
return 0.0
Ne = param_vec[ParamOrder.PopSize].value
mu = param_vec[ParamOrder.MutRate].value
tau = param_vec[ParamOrder.DivTime].value
# The probability of coalescence along the Turkish branch goes to 1 as tau gets big or Ne gets small
prob_coalescence = 1 - exp(-tau/(2.0*Ne))
prob_no_coalescence = 1 - prob_coalescence
# Given that there is no coalescence, each branch will have length tau along the Turkish population history
# the following formula would apply to each branch (and are the same as the probabilities for the Spanish population)
expNeg2TauMu = exp(-2*tau*mu)
prob_no_mut_given_no_coal = (1 + expNeg2TauMu)/2.0
prob_mut_given_no_coal = 1 - prob_no_mut_given_no_coal
# Here is some very crude averaging over branch lengths that Bryant et al don't do...
# we'll assume that if there is coalescence it averages out to occurring 1/3 of the way from the present to the
# time of divergence. No real justificaion for this...
hacky_total_branch_length_given_coal = (4/3) * tau # hack mentioned above
expNegHackyBranchLen = exp(-2*hacky_total_branch_length_given_coal*mu)
prob_no_mut_given_coal = (1 + expNegHackyBranchLen)/2.0
prob_mut_given_coal = 1 - prob_no_mut_given_no_coal
# 2/3 of tau is contributed by the common ancestral branch,
# 1/3 of tau is contributed by each terminal (for a total of 4/3)
# so, given that there was coalescence, and only one mutation in the Turkish branch, and
# our new hack. 50% of mutations will be on the ancestral branch of the gene tree,
# and 25% will be on each of the terminal's
prob_mut_common_anc_branch = 0.5*prob_mut_given_coal # hack mentioned above, again
prob_mut_terminal_branch = .25*prob_mut_given_coal # hack mentioned above, again
if (r_Tb == 1 ):
#eqn 10 - eqn 14
if (n_Tt == 2):
#eqn 10 - eqn 12
prob_tree_shape = prob_no_coalescence
if (r_Tt == 1):
# eqn 11 would be the correct form to use
#
# here we'll use the hacky way...
# same # of lineages in state 0 at the beginning and end. This probability will be dominated
# by the scenario of no mutations...
prob_two_mutations = prob_mut_given_no_coal * prob_mut_given_no_coal
prob_zero_mutations = prob_no_mut_given_no_coal * prob_no_mut_given_no_coal
prob_mutation_scenario = prob_two_mutations + prob_zero_mutations
else:
# eqn 10 and 12 are identical
# hacky way -> one lineage mutated, and the other did not...
prob_mutation_scenario = prob_mut_given_no_coal * prob_no_mut_given_no_coal
else:
#eqn 13 and eqn 14 are the same, so only one branch of code here
# the hacky way is to assume there was only one mutation, and it must have been on a terminal...
prob_tree_shape = prob_coalescence
prob_mutation_scenario = prob_no_mut_given_coal*prob_mut_terminal_branch
else:
#eqn 10 - eqn 14
if (n_Tt == 2):
#eqn 5 - eqn 5
prob_tree_shape = prob_no_coalescence
if (r_Tt == r_Tb):
# eqn 5 would be the correct way...
# here is the hack...
prob_mutation_scenario = prob_no_mut_given_no_coal * prob_no_mut_given_no_coal
elif (r_Tt == 1):
# eqn 6 is the correct way...
# hacky way -> one lineage mutated, and the other did not...
prob_mutation_scenario = prob_mut_given_no_coal * prob_no_mut_given_no_coal
else:
#eqn 7 is the correct way.
# if r_Tb is 0 or 2 and r_Tt is 0 or 2, but r_Tt != r_Tb, then both lineages had mutations
prob_mutation_scenario = prob_mut_given_no_coal * prob_mut_given_no_coal
else:
#eqn 8 or eqn 9
prob_tree_shape = prob_coalescence
if (r_Tt == 1):
if (r_Tb == 2):
# eqn 8
# no mutations required when a state-0 lineage splits into two state-0 lineages
prob_mutation_scenario = prob_no_mut_given_coal
else:
# eqn 9
# a mutations before splitting is the easiest way to turn a state-0 lineage splits into two state-1 lineages
prob_mutation_scenario = prob_mut_common_anc_branch
else:
# eqn 9
if (r_Tb == 0):
# eqn 8
# no mutations required when a state-1 lineage splits into two state-1 lineages
prob_mutation_scenario = prob_no_mut_given_coal
else:
# eqn 9
# a mutations before splitting is the easiest way to turn a state-1 lineage splits into two state-0 lineages
prob_mutation_scenario = prob_mut_common_anc_branch
return prob_tree_shape*prob_mutation_scenario
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment