R 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

gistfile1.r

calc.F_Tt <- function (r_Tb, n_Tt, r_Tt, param.vec) {
    if (r_Tt > n_Tt) {
        return (0.0);
    }
    # unnumbered eqn on the bottom of page 5
    Ne <- param.vec[2];
    mu <- param.vec[3];
    tau <- param.vec[4];

    # 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 = .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
            } else if (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)
}