Tools For Evo Devo

##############################################################
# =========================================================
# An Evo-Devo Analysis in R: an example with cone snails
# =========================================================
# by Nick Matzke (and whoever else adds to this PhyloWiki page)
# Copyright 2011-infinity
# matzkeATberkeley.edu
# March 2011
#
# Please link/cite if you use this, email me if you have 
#   thoughts/improvements/corrections.
#
##############################################################
#
# Copyright: this lab contains/uses unpublished research data & 
# analyses.
# 
# IB200b students are of course free to re-use the code, others should 
# contact me about re-use of the code or data.
#
# This work was done in collaboration with Monty Slatkin, George
# Oster, and Neil Gong.
# 
###################################################################
#
# Assignment for Tuesday, March 8
#
# In this lab, I will take you through a real analysis I am
# working on that is basically applying ancestral character
# reconstruction methods and model selection to an evo-devo
# question.
#
# George Oster's lab (Oster is a biophysicist) developed a shell-
# growth model for cone snails that runs in Matlab.  (Google cone
# snails to look at some shells.)
#
# Basically, given a series of 24 input parameters, the model grows
# a shell, complete with the pattern, which might be stripes, dots,
# triangles, blank, etc.  The first 19 parameters control various 
# signal and response curves that cells of the cone snail mantle,
# under the model issue and receive, which then control the release
# of pigment.
#
# The last 5 parameters deal with shell shape:
#
# kz, kr, and turn are from Raup's famous shell-shape model, and 
# control the # of turns, the angle, etc. of the shell
#
# N and T are starting conditions for the model (# of cells at the 
# growing edge, and the time point they start at).
#
#
#
# Once they built the model, Oster's lab then took a bunch of living
# species, and figured out the parameters that grow patterns to match
# each of these species.
# 
# Oster then came over to Integrative Biology to "do some evolution
# with this".  Some obvious questions are:
#
# 1. Is there any phylogenetic signal in these parameters, which 
# after all are not direct measurements, but parameters that were
# fit to a computer model by eye.
#
# 1a. Are the parameters statistically independent, or highly 
# correlated with each other?
#
# 2. Is Brownian motion a reasonable model under which to 
# conduct ancestral state reconstruction?
#
# 3. What model of evolution best describes each parameter? 
# Available models are:
# - Brownian motion (BM);
#
# - the Ornstein−Uhlenbeck (OU) stabilizing selection model;
#
# - the lambda model which rescales internal branch lengths by
#   a linear fraction;
#
# - the kappa model (green) which rescales each branch length by a power
#   equal to the kappa parameter, and which becomes a speciational model
#   as kappa approaches 0;
#
# - delta (yellow) which focuses change towards the base or tips;
#
# - early burst (EB, cyan) which has an initial high rate of change that
#   then declines;
#
# - and white noise (white), where observations are produced by a normal
#   distribution with no tree structure, which represents the situation
#   of no phylogenetic signal.
#
# See references under the function fitContinuous.
#
# Assignment: 
# 
# Work through this code, and try to understand what each 
# stage does.
#
# Have a look at the final Likelihood Ratio Test results, and the model
# comparison figure which uses AIC, and email me your verbal summary of
# those results.  E.g., which model is best, how confidence should the 
# reader be in those conclusions, and what are the implications in 
# terms of understanding the evolutionary processes involved,
# and the likely accuracy of using Brownian motion for ancestral
# character reconstruction.
#
# Files to download:
#
# Start by looking at the graphical file Neil_tree.png on the website.
# This has photos of the cone snail species, and the matching shells
# grown in the computer, on a DNA phylogeny.
# 
# The parameters are stored in parameters_njm7.txt, which this script
# takes as input, along with the DNA tree, tr.newick, and the 
# ultrametricized version, tr2.newick.
#
# Just to cut to the end, The final results figure is at:
# tree_DNA_num_ML_inv2.png
#
###################################################################

# Stuff that is handy at the beginning of a script
######################################################
# R, close all devices/windows/ quartz windows
graphics.off()

# Turn off all BS stringsAsFactors silliness
# R tends to read 
options(stringsAsFactors = FALSE)
######################################################

# Working directories:
# One of the first things you want to do, usually, is 
# decide on your working directory

# On my Mac, this is the directory
wd = "/Users/nick/Desktop/_ib200a/ib200b_sp2011/lab10"

# On a PC, you might have to specify paths like this:
#wd = "c:\\Users\\nick\\Desktop\\_ib200a\\ib200b_sp2011\\lab10"

#desc: set working directory
setwd(wd)

#desc: get working directory
getwd()

#desc: list the files in the directory
list.files()

# Sourcing scripts
# Local sourcing:
# source a script (windows machines require \\ instead of /)
# sourcedir = '/_njm/'
#source3 = '_genericR_v1.R'
#source(paste(sourcedir, source3, sep=""))

# You need the APE package to source genericR_v1.R scripts
install.packages("ape")
library(ape)

# Remote sourcing:
source("http://ib.berkeley.edu/courses/ib200b/scripts/_genericR_v1.R")
source("http://ib.berkeley.edu/courses/ib200b/scripts/_conus_analyses_v1.R")
source("http://ib.berkeley.edu/courses/ib200b/scripts/_R_tree_functions_v1.R")

# Setup:
# install these, if needed
library(ape)
library(phangorn)      # helps with "midpoint" function
library(picante)      # for "Kcalc" and "phylosignal", i.e. calculation of the kappa statistic for continuous parameters
library(geiger)     # for e.g. fitContinuous

#########################################
# Plot the original DNA tree
#########################################
trfn = "tr.newick"
tr = read.tree(trfn)
plot(tr)
title("Original DNA tree (digitized from Nam et al. 2009\nusing TreeRogue 0.1)\n(Note: axis scale is in digitization coordinates)")
axisPhylo()
add.scale.bar(x=0, y=23, length=1)

# Get the parameters file
params_orig = read_table_good("parameters_njm7.txt")

# identify the types of params
uniform1_params_names = c("s1g", "s2eg", "s2ig", "s3g", "sf2ia", "i1s1g", "i1s2eg", "i1s2ig", "i1s3g", "i1sf2ia", "i2s1g", "i2s2eg", "i2s2ig", "i2s3g", "i2sf2ia")
continuous_params_names = c("s1n", "s1t", "s2en", "s2et", "s2in", "s2it", "s3n", "s3t", "sf2ea", "sf2ed", "sf2id", "tf2eb1", "tf2ec1", "tf2eb2", "tf2ec2", "tf2ib1", "tf2ic1", "tf2ib2", "tf2ic2")
shell_params_names = c("kz", "kr", "turn")
starting_params_names = c("N", "T")
complex_params_names = c("s3n_std", "s3t_spd")
clad_chars_names = c("stripes", "triangles", "dots", "color", "conical", "food")
complex1_names = names(params_orig)[(ncol(params_orig)-25*2) : (ncol(params_orig)-26)]
complex2_names = names(params_orig)[(ncol(params_orig)-25) : ncol(params_orig)]
nonsucky_names = c(continuous_params_names, shell_params_names, starting_params_names, complex_params_names, clad_chars_names)
all_cont_names = c(continuous_params_names, shell_params_names, starting_params_names)

# remove the uniform params
cols_tokeep = names(params_orig) %in% uniform1_params_names == FALSE

all_params = params_orig[, cols_tokeep]
row.names(all_params) = all_params$species
(all_params)

###################################
# Drop the tips we're not using
###################################
all_params_orig = all_params
#n1 = xy3$nodenames[xy3$nodenames != ""]
n1 = tr$tip.label
n2 = all_params$species
n1_to_keep = n1 %in% n2
print(cbind(n1, n1_to_keep))

# drop the tips we're not using
tr2 = drop.tip(tr, n1[n1_to_keep == FALSE])

####################################################################
# reorder the parameters to match the order of tips in the tree
# (this is CRUCIAL!!)
####################################################################
# tree tip order, from bottom to top:
(labels1 = tr2$tip.label)

# order of rows in char_dtf
(labels2 = row.names(all_params))

# reorder char_dtf to match tip.labels
new_order = match(labels1, labels2)
all_params = all_params_orig[new_order, ]

params = all_params

# Cut params down to just continuous parameters for correlation analysis
params = params[, all_cont_names]

####################################

#pairs(params[, -1])
cps = cor(params)
cps = as.data.frame(cps)
names_cps = names(cps)
#dim(params[, -1])

# Rank pairwise correlations
pairwise_cors_list = c()
for (i in 2:nrow(cps))
    {
    if ((i+1) < ncol(cps))
        {
        for (j in (i+1):ncol(cps))
            {
            tmprow = c(names_cps[i], names_cps[j], cps[i,j], abs(cps[i, j]))
            pairwise_cors_list = rbind(pairwise_cors_list, tmprow)
            }
        }
    }
row.names(pairwise_cors_list) = NULL
pairwise_cors_list = as.data.frame(pairwise_cors_list)
names(pairwise_cors_list) = c("x", "y", "cor", "abs_cor")
#rank(pairwise_cors_list$abs_cor)
order_cors = order(pairwise_cors_list$abs_cor, decreasing=TRUE)

pairwise_cors_list = pairwise_cors_list[order_cors, ]

print("Top correlations:")
(pairwise_cors_list[1:20, ])

######################################################
# Phylogenetically Independent Contrasts
######################################################
# we are going to use an ultrametric tree:
chtr2 = read.tree("chtr2.newick")
#chronogram(tr2, scale = 1, expo = 2, minEdgeLength = 1e-06)
write.tree(chtr2, file="chtr2.newick")
plot(chtr2)
title("DNA tree, ultrametricized with APE's chronogram() function: NPRS\n(nonparametric rate smoothing), basically a least squares fit after Sanderson 1997")
axisPhylo()
add.scale.bar(x=0, y=0.9, length=0.1)

pic_params = c()
for (i in 1:ncol(params))
    {
    x = pic(params[,i], chtr2)
    pic_params = cbind(pic_params, x)
    }

pic_params = as.data.frame(pic_params)
names(pic_params) = names(params)

#pairs(pic_params)
cor(pic_params)

# Plot all params vs. all params
tmpparams = params[, all_cont_names]
tmpparams = dfnums_to_numeric(tmpparams)
tmpparams[tmpparams == "?"] = NA
pairs.panels(tmpparams, smooth=FALSE)
title("Regular correlations")
pairs.panels(pic_params, smooth=FALSE)
title("PIC correlations")

# Assemble the tables of correlations for params and PICs of params
res1_txt_table = c()
respic_txt_table = c()
for (i in 1:ncol(params))
    {
    for (j in i:ncol(params))
        {
        if (i == j)
            {
            next()
            }
        cat("x=", names(params)[i], "        vs.        y=", names(params)[j], "\n", sep="")

        res1 = cor.test(params[, i], params[, j])
        respic = cor.test(pic_params[, i], pic_params[, j])
        #print()
        #print()

        lm_res1 = lm(params[, j] ~ params[, i])
        lm_respic = lm(pic_params[, j] ~ pic_params[, i])

        res1_txt = c(names(params)[i], names(params)[j], res1$statistic, res1$parameter, res1$p.value, res1$estimate, res1$null.value, res1$alternative, res1$method, res1$data.name, res1$conf.int[1], res1$conf.int[2], lm_res1$coefficients[2], lm_res1$coefficients[1])

        respic_txt = c(names(params)[i], names(params)[j], respic$statistic, respic$parameter, respic$p.value, respic$estimate, respic$null.value, respic$alternative, respic$method, respic$data.name, respic$conf.int[1], respic$conf.int[2], lm_respic$coefficients[2], lm_respic$coefficients[1])

        res1_txt_table = rbind(res1_txt_table, res1_txt)
        respic_txt_table = rbind(respic_txt_table, respic_txt)
        }
    }

res1_txt_table = as.data.frame(res1_txt_table, row.names=1:nrow(res1_txt_table))
names(res1_txt_table) = (c("x", "y", "statistic_t", "parameter_df", "pval", "cor_estimate", "null", "alternative", "method", "data.name", "conf.int1", "conf.int2", "slope", "intercept"))

respic_txt_table = as.data.frame(respic_txt_table, row.names=1:nrow(respic_txt_table))
names(respic_txt_table) = (c("x", "y", "statistic_t", "parameter_df", "pval", "cor_estimate", "null", "alternative", "method", "data.name", "conf.int1", "conf.int2", "slope", "intercept"))

res1_txt_table = dfnums_to_numeric(res1_txt_table)
respic_txt_table = dfnums_to_numeric(respic_txt_table)

res1_txt_table$pval_bonf = res1_txt_table$pval * nrow(res1_txt_table)
respic_txt_table$pval_bonf = respic_txt_table$pval * nrow(respic_txt_table)

res1_txt_table = res1_txt_table[order(respic_txt_table$pval_bonf),]
respic_txt_table = respic_txt_table[order(respic_txt_table$pval_bonf),]

(res1_txt_table[1:20,])
(respic_txt_table[1:20,])

# sort by ascending p-value
res1_txt_table_sort = res1_txt_table[order(res1_txt_table$pval), ]
respic_txt_table_sort = respic_txt_table[order(respic_txt_table$pval), ]

(res1_txt_table_sort[1:20,])
(respic_txt_table_sort[1:20,])

x = res1_txt_table$cor_estimate
y = respic_txt_table$cor_estimate
xlab = "r (regular)"
ylab="r (PIC)"
xlim = c(-1, 1)
ylim = c(-1, 1)
titlestart_txt = "correlation (r) estimates for pairwise correlations,\nnormal vs. phylogenetically independent contrasts;\n"
basic_xy_plot(x, y, titlestart_txt, xlab, ylab, xlim, ylim)

x = res1_txt_table$pval
y = respic_txt_table$pval
xlab="p-val (regular)"
ylab="p-val (PIC)"
xlim = c(0, 1)
ylim = c(0, 1)
titlestart_txt = "p-values for pairwise correlations,\nnormal vs. phylogenetically independent contrasts\n"
basic_xy_plot(x, y, titlestart_txt, xlab, ylab, xlim, ylim)

# Write the tables to text
write.table(res1_txt_table, "res1_txt_table.txt")
write.table(respic_txt_table, "respic_txt_table.txt")
############################################################

#####################################################
# Assess the parsimony stats of the discrete characters
#####################################################

colnames_to_use = c(complex_params_names, clad_chars_names)
char_dtf_orig = all_params[, colnames_to_use]
char_dtf = char_dtf_orig

# Make the final figure comparing the trees

# quite fast
#(parstats = parsim_stats_fast(tr2, dfnums_to_numeric(char_dtf)))
(parstats = parsim_stats_fast(tr2, dfnums_to_numeric(char_dtf)))
write_table_good(parstats, "parstats_v8.txt")

# slower
#(parstats = parsim_stats(tr2, char_dtf))

# Now, resample (takes a few seconds)
num_resamps = 100
parstats_null = array(data=NA, dim=c(ncol(char_dtf)+1, 6, num_resamps+1))
for (i in 1:num_resamps)
    {
    tmpsamp = resample_dtf_matrix(char_dtf)
    parstats_null[, , i] = as.matrix(parsim_stats_fast(tr2, tmpsamp))
    }
parstats_null[, , num_resamps+1] = as.matrix(parstats)

summary_parstats_means = parstats
summary_parstats_sds = parstats
summary_parstats_pvals = parstats

# do histograms of the six different distances
# set outer margins & number of subplots
par(oma=c(1,1,5,1), mfrow=c(2,3))
for (i in 1:nrow(parstats))
    {
    for (j in 1:ncol(parstats))
        {
        x = parstats_null[i, j, ]
        summary_parstats_means[i, j] = mean(x)
        summary_parstats_sds[i, j] = sd(x)

        observed_val = parstats_null[i, j, num_resamps+1]
        pval = round(empirical_pval_fast(x, observed_val), 4)

        # redo pval depending on whether value should be low or high
        if (j > 3)
            {
            pval = round(1-pval, 4)
            }

        summary_parstats_pvals[i, j] = pval

        tmp_title = paste("non-parametric p-value = ", pval, sep="")
        h = hist(x, main=tmp_title, breaks=50, cex.main=1, xlab=names(parstats)[j])

        # add the arrow
        arrowtop = 0.3*max(h$counts, na.rm=TRUE)
        arrowbot = 0.03*max(h$counts, na.rm=TRUE)
        arrows(observed_val, arrowtop, observed_val, arrowbot, col="blue", lwd=4)    
        }
    titletxt = paste("Parsimony stats for ", row.names(parstats)[i], sep="")
    mtext(titletxt, outer=TRUE)
    }

write_table_good(summary_parstats_means, "parstats_means_v8.txt")
write_table_good(summary_parstats_sds, "parstats_sds_v8.txt")
write_table_good(summary_parstats_pvals, "parstats_pvals_v8.txt")

#######################################################
# Ancestral character reconstruction --  continuous, plots
#######################################################

par(mfrow=c(1,2))

# character mapping
# uniform1_params_names = c("s1g", "s2eg", "s2ig", "s3g", "sf2ia")
# continuous_params_names = c("s1n", "s1t", "s2en", "s2et", "s2in", "s2it", "s3n", "s3t", "sf2ea", "sf2ed", "sf2id", "tf2eb1", "tf2ec1", "tf2eb2", "tf2ec2", "tf2ib1", "tf2ic1", "tf2ib2", "tf2ic2")
# shell_params_names = c("kz", "kr", "turn")
# starting_params_names = c("N", "T")
# complex_params_names = c("s3n_std", "s3t_spd")
# clad_chars_names

# continuous parameter characters
colnames_to_use = c(continuous_params_names, shell_params_names, starting_params_names)
char_dtf_orig = all_params[, colnames_to_use]
char_dtf = char_dtf_orig

#################
# reorder char_dtf
################
# tree tip order, from bottom to top:
#(labels1 = chtr2$tip.label)
# order of rows in char_dtf
#(labels2 = row.names(char_dtf))
# reorder char_dtf to match tip.labels
#new_order = match(labels1, labels2)
#char_dtf = char_dtf[new_order, ]

par(mfrow=c(6,4))
for (colnum in 1:ncol(char_dtf))
    {
    hist(as.numeric(char_dtf[,colnum]), 17, main=title(names(char_dtf)[colnum]) )
    }

par(mfrow=c(1,1))

###########################################################################
# Compare different models for the evolution of this character
###########################################################################
runthis4 = TRUE
if (runthis4 == TRUE)
    {
    # models_list = c("BM", "OU", "lambda", "kappa", "delta", "EB", "white", "trend")
    # cut "trend" model, doesn't matter for ultrametric trees (same as BM model, evidently)
    models_list = c("BM", "OU", "lambda", "kappa", "delta", "EB", "white")
    model_fits = NULL

    # GET THE VARIOUS PARAMETERS ALSO

    # beta (typically the rate)
    # alpha (central tendency on BM walk that is proportional to alpha)

    # Lambda model
    # multiplies all internal branches of the tree by lambda, leaving tip
    # branches as their original length
    # 
    # kappa raises all branch lengths to the power kappa. As kappa approaches
    # zero, the model becomes speciational.
    #
    # (apparently, kappa is saved as "lamda" in the kappa model)

    # delta (delta>1, evo tipwards, < 1, evo at base)

    # EB (a = proportion of initial rate represented by the end rate (?))

    # White Noise (draw from mean/variance)
    # mean = mean
    # nv = variance (?)

    # THESE BOUNDS PRODUCE PATHOLOGICAL RESULTS (on my first dataset)
    bounds_list = NULL
    # bounds_list$BM = list(beta=c(1,100))
    # bounds_list$OU = list(alpha=c(0.0001, 10), beta=c(1, 100))
    # bounds_list$lambda = list(beta=c(1,100))
    # bounds_list$kappa = list(beta=c(1,100))
    # bounds_list$EB = list(beta=c(1,100))
    # bounds_list$white = NULL

    # THESE BOUNDS ARE THE DEFAULTS
    bounds_list = NULL

    # db = default bounds
    #db_beta = c(1e-08, 20)
    #db_lambda = c(1e-07, 1)
    #db_kappa = c(1e-06, 1)
    #db_delta = c(1e-05, 2.999999)
    #db_alpha = c(1e-07, 50)
    #db_a = c(-3, 0)
    #db_nv = c(1e-10, 100)
    #db_mu = c(-100, 100)

    # YOU SHOULD ALWAYS CHECK IF YOUR ESTIMATES HIT THE BOUNDS, AND
    # THEN RE-ADJUST
    bounds_list = NULL

    # db = default bounds
    db_beta = c(1e-08, 1e08)
    db_alpha = c(1e-07, 50)
    db_lambda = c(1e-07, 1)
    db_kappa = c(1e-06, 1)
    db_delta = c(1e-05, 2.9999)
    db_a = c(-3, 0)
    db_nv = c(1e-10, 1000000)
    #db_mu = c(-100, 100)

    bounds_list$BM = list(beta=db_beta)
    bounds_list$OU = list(beta=db_beta, alpha=db_alpha)
    bounds_list$lambda = list(beta=db_beta, lambda=db_lambda)
    bounds_list$kappa = list(beta=db_beta, lambda=db_kappa)
    bounds_list$EB = list(beta=db_beta, a=db_a)
    bounds_list$white = list(nv=db_nv)
    # bounds_list$model = list(mu=db_mu)

    for (i in 1:length(models_list))
        {
        # Chose the right bounds
        cmdstr = paste("tmp_bounds = bounds_list$", models_list[i], sep="")
        eval(parse(text = cmdstr))

        model_fit = fitContinuous(chtr2, char_dtf, model=models_list[i], bounds = tmp_bounds)
        cmdstr = paste("model_fits$", models_list[i], " = model_fit", sep="")
        eval(parse(text = cmdstr))
        }
    # Save R object to a file
    save(model_fits, file="model_fits.Rdata")
    } else {
    load(file="model_fits.Rdata")
    }

model_fits_orig = model_fits
model_fits = model_fits_orig

# Get log likelihoods and AICs
results_list = c("lnl", "aic", "aicc", "k", "beta", "alpha", "lambda", "delta", "a", "mean", "nv")

# fill in blank results
for (item in results_list)
    {
    cmdstr = paste(item, "_matrix = matrix(NA, nrow=length(models_list), ncol=length(all_cont_names))", sep="")
    eval(parse(text = cmdstr))    
    }

# put in the data from the model results
for (i in 1:length(models_list))
    {
    for (j in 1:length(all_cont_names))
        {
        for (item in results_list)
            {
            cmdstr = paste("tmp_item = model_fits$", models_list[i], "$", all_cont_names[j], "$", item, sep="")
            eval(parse(text = cmdstr))
            if (is.null(tmp_item))
                {
                cmdstr = paste(item, "_matrix[i,j] = NA", sep="")
                eval(parse(text = cmdstr))    
                } else {
                cmdstr = paste(item, "_matrix[i,j] = tmp_item", sep="")
                eval(parse(text = cmdstr))
                }
            }
        }
    }

# Apply the row & column names
for (item in results_list)
    {
    cmdstr = paste(item, "_matrix = adf(", item, "_matrix)", sep="")
    eval(parse(text = cmdstr))

    cmdstr = paste("names(", item, "_matrix) = names(char_dtf)", sep="")
    eval(parse(text = cmdstr))

    cmdstr = paste("row.names(", item, "_matrix) = models_list", sep="")
    eval(parse(text = cmdstr))

    # write to txt file
    cmdstr = paste("write_table_good(", item, "_matrix, '", item, "_matrix.txt')", sep="")
    eval(parse(text = cmdstr))
    }

lnl_matrix
aic_matrix
aicc_matrix
k_matrix
beta_matrix
alpha_matrix
lambda_matrix
delta_matrix
a_matrix
mean_matrix
nv_matrix

# Compare the models using AIC
# Another way: use AIC
# For each col, subtract the min

diff_aic_matrix = aic_matrix

for (colnum in 1:ncol(aic_matrix))
    {
    diff_aic_matrix[, colnum] = aic_matrix[, colnum] - min(aic_matrix[, colnum], na.rm=TRUE)
    }

cat("Table of difference-from-best-AIC values; 0 = best model\n")
print(diff_aic_matrix, digits=2)

# Citation for this calculation:
# http://www.brianomeara.info/tutorials/aic
# 
# or:
#
# 
aic_wi = exp(1) ^ (-0.5 * diff_aic_matrix)
colsum_aic_wi = apply(aic_wi, 2, sum, na.rm=TRUE)

# calculate the relative weights
aic_weights = aic_wi
for (colnum in 1:ncol(aic_wi))
    {
    aic_weights[,colnum] = aic_wi[,colnum] / colsum_aic_wi[colnum]
    }

# 
cat("These are the AIC weights, calculated according to Burnham & Anderson (2002).\n")
print(aic_weights, digits=2)

# flip aic_weights for graphing
aic_weights[nrow(aic_weights):1, ] = aic_weights
row.names(aic_weights) = row.names(aic_weights)[nrow(aic_weights):1]

#####################################
# Plot AIC weights, with labels
#####################################
# What are the best models under AIC
best_model_list = rownames(aic_weights)

best_models = rep(NA, length(colnames(aic_weights)))
for (colnum in 1:length(colnames(aic_weights)))
    {
    TF_row = aic_weights[,colnum] == max(aic_weights[,colnum])
    best_models[colnum] = best_model_list[TF_row]
    }
print(best_models)

# Get results
#pdf(file = "conus_cors_all_vs_all_v3.pdf", width=15, height=15)
doPDFs2 = TRUE
if (doPDFs2 == TRUE)
    {
    pdffn = "Figure_SX_model_compare_v8.pdf"
    pdf(file=pdffn, width=8, height=10, paper="USr")
    }

par(oma=c(1,1,1,1), mar=c(15.6, 4.1, 4.1, 7.1), mfrow = c(1,1))
par(xpd=TRUE)
tmp_colors = c("white", "cyan", "yellow", "green2", "orange", "red", "brown")
tmp_xnames = rep("", length(colnames(aic_weights)))
legend_text = rev(c("Brownian", "O-U", "lambda", "kappa", "delta", "early burst", "white noise"))
tmp_legend = list(x=38, y=0.9, pt.cex=4)

# don't plot y-axis here
par(yaxt="n")
bar_x_positions = barplot(aic_weights, names.arg=tmp_xnames, col=tmp_colors, ylab="Proportion of AIC weight", xlab="shell growth parameter", legend.text=legend_text, args.legend=tmp_legend)

# plot y-axis here
par(yaxt="s")
ticks = c(0, 0.25, 0.5, 0.75, 1)
axis(side=2, at=ticks, line=0, labels=ticks, lty = "solid", lwd = 1)

## Create plot and get bar midpoints in 'mp'
#bar_x_positions = barplot(1:length(colnames(aic_weights)))

## Set up x axis with tick marks alone
#axis(1, at = bar_x_positions, labels = FALSE)

# Create some text labels
tmp_x_labels = colnames(aic_weights)

# Plot x axis labels at mp
text(x=bar_x_positions-0.1, y=-0.015, adj=0, srt=315, labels = tmp_x_labels, xpd=TRUE)

# Plot the best models across the top
text(x=bar_x_positions-0.1, y=1.015, adj=0, srt=45, labels = best_models, xpd=TRUE)

toptxt = "(Best model under AIC)"
text(x=15, y=1.14, labels=toptxt)

caption = "Figure SX: Akaike weights of 7 models for the evolution of shell parameters.  The models \ncompared are: Brownian motion (BM, brown); the Ornstein-Uhlenbeck (OU, red) \nstabilizing selection model; the lambda model (orange) which rescales internal branch \nlengths by a linear fraction; the kappa model (green) which rescales each branch length \nby a power equal to the kappa parameter, and which becomes a speciational model as \nkappa approaches 0; delta (yellow) which focuses change towards the base or tips; \nearly burst (EB, cyan) which has an initial high rate of change that then declines; and white \nnoise (white), where observations are produced by a normal distribution with no tree  \nstructure, which represents the situation of no phylogenetic signal. Brownian motion (BM) \nhas the highest AIC weight for 63% of the shell parameters (15/24).  White noise \n(no phylogenetic signal) is superior for 25% (6/24) shell parameters."

text(x=-5.4, y=-0.56, labels=caption, adj=0, xpd=TRUE)

##########################################################

# Turn off the PDF and open
if (doPDFs2 == TRUE)
    {
    dev.off()
    cmdstr = paste("open ", pdffn, sep="")
    system(cmdstr)
    }

##########################################################
# Do LRT on the different models
##########################################################
# Compare all models to BM
pchisq_table = matrix(NA, nrow=5, ncol=ncol(lnl_matrix))

for (colnum in 1:ncol(lnl_matrix))
    {
    tmp_lnls = lnl_matrix[, colnum]

    # Calculate the log-likelihood differences from 
    # the Brownian motion model
    #
    # You can't compare white noise, since it's not nested
    tmp_lnls_diffs = tmp_lnls[2:(length(tmp_lnls)-1)] - tmp_lnls[1]

    # take twice the difference
    tmp_lnls_twice_diffs = 2 * tmp_lnls_diffs
    dim(tmp_lnls_twice_diffs) = c(1, length(tmp_lnls_twice_diffs))

    # Do chi-squared test with 1 degree of freedom
    tmp_pchisq = apply(tmp_lnls_twice_diffs, 1, pchisq, 1, lower.tail=FALSE)

    pchisq_table[,colnum] = tmp_pchisq

    #print(tmp_pchisq)
    }

pchisq_table = adf(pchisq_table)
names(pchisq_table) = names(lnl_matrix)
row.names(pchisq_table) = row.names(lnl_matrix)[2:(nrow(lnl_matrix)-1)]

print(pchisq_table, digits=2)

write_table_good(pchisq_table, "pchisq_table.txt")

# 11 / 120 times there is a model significantly better than BM
pchisq_table < 0.05
sum(pchisq_table < 0.05)

# 5/ 24 params reject BM significantly:
# sf2id, tf2ib1, tf2ib2, kr, turn
BM_rejected = apply(pchisq_table < 0.05, 2, sum)

# All false with Bonferroni correction, but that's probably not fair
pchisq_table < (0.05 / (nrow(pchisq_table) * ncol(pchisq_table)) )

##########################################################

#x = dtf_to_phyDat(char_dtf)
txt_suffix = ""

# ACE methods:
# a character specifying the method used for estimation. Three choices are possible: "ML", "pic", or "GLS".

# ACE with ML method
# BM = brownian motion
num_internal_nodes = chtr2$Nnode

PICcont_means = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
PICcont_lower025 = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
PICcont_upper975 = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )

for (i in 1:ncol(char_dtf))
    {
    PICcont = ace(char_dtf[,i], chtr2, type="continuous", method="pic", CI=TRUE, model="BM")
    PICcont_means[, i] = PICcont$ace
    PICcont_lower025[, i] = PICcont$CI95[, 1]
    PICcont_upper975[, i] = PICcont$CI95[, 2]

    }
PICcont_means = adf(PICcont_means)
names(PICcont_means) = names(char_dtf)
tree_dtf = prt(chtr2)
row.names(PICcont_means) = tree_dtf$node[ (length(chtr2$tip.label)+1): length(tree_dtf$node)]
PICcont_lower025 = adf(PICcont_lower025)
names(PICcont_lower025) = names(char_dtf)
row.names(PICcont_lower025) = row.names(PICcont_means)
PICcont_upper975 = adf(PICcont_upper975)
names(PICcont_upper975) = names(char_dtf)
row.names(PICcont_upper975) = row.names(PICcont_means)

write_table_good_w_rownames(PICcont_means, "PICcont_means_v8.txt")
write_table_good_w_rownames(PICcont_lower025, "PICcont_lower025_v8.txt")
write_table_good_w_rownames(PICcont_upper975, "PICcont_upper975_v8.txt")

# plot tip & node values to check them

# plot tip & node values to check them
txt1 = "PIC"
txt2 = txt_suffix
for (i in 1:ncol(char_dtf))
    {
    recon_txt = paste("ACE technique: ", txt1, txt2, sep="")

    tmptipvals = char_dtf[,i]
    tmpnodevals = PICcont_means[,i]
    maxval = max(c(tmptipvals, tmpnodevals))
    texttxt = paste(recon_txt, "\nParameter: ", names(char_dtf)[i], "\nMaximum value = ", maxval, sep="") 

    plot_cont_vals_on_tree(chtr2, tmptipvals, tmpnodevals, texttxt)
    }

MLcont_means = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
MLcont_lower025 = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
MLcont_upper975 = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
Blombergs_Kstats = c()
PIC.variance.obs = c()
PIC.variance.rnd.mean = c()
PIC.variance.P = c()
PIC.variance.Z = c()

for (i in 1:ncol(char_dtf))
    {
    MLcont = ace(char_dtf[,i], chtr2, type="continuous", method="ML", CI=TRUE, model="BM")
    phylosignals = phylosignal(char_dtf[,i], chtr2)
    Blombergs_Kstats = c(Blombergs_Kstats, phylosignals$K)
    PIC.variance.obs = c(PIC.variance.obs, phylosignals$PIC.variance.obs)
    PIC.variance.rnd.mean = c(PIC.variance.rnd.mean, phylosignals$PIC.variance.rnd.mean)
    PIC.variance.P = c(PIC.variance.P, phylosignals$PIC.variance.P)
    PIC.variance.Z = c(PIC.variance.Z, phylosignals$PIC.variance.Z)

    MLcont_means[, i] = MLcont$ace
    MLcont_lower025[, i] = MLcont$CI95[, 1]
    MLcont_upper975[, i] = MLcont$CI95[, 2]
    }
MLcont_means = adf(MLcont_means)
names(MLcont_means) = names(char_dtf)
tree_dtf = prt(chtr2)
row.names(MLcont_means) = tree_dtf$node[ (length(chtr2$tip.label)+1): length(tree_dtf$node)]
MLcont_lower025 = adf(MLcont_lower025)
names(MLcont_lower025) = names(char_dtf)
row.names(MLcont_lower025) = row.names(MLcont_means)
MLcont_upper975 = adf(MLcont_upper975)
names(MLcont_upper975) = names(char_dtf)
row.names(MLcont_upper975) = row.names(MLcont_means)

# Bonferroni correction
PIC.variance.P_bonferroni = PIC.variance.P * length(Blombergs_Kstats)
phylosignal_table = cbind(Blombergs_Kstats, PIC.variance.obs, PIC.variance.rnd.mean, PIC.variance.Z, PIC.variance.P, PIC.variance.P_bonferroni)
row.names(phylosignal_table) = names(char_dtf)
(phylosignal_table)

write_table_good_w_rownames(MLcont_means, "MLcont_means_v8.txt")
write_table_good_w_rownames(MLcont_lower025, "MLcont_lower025_v8.txt")
write_table_good_w_rownames(MLcont_upper975, "MLcont_upper975_v8.txt")
write_table_good_w_rownames(phylosignal_table, "phylosignal_table_v8.txt")

# plot tip & node values to check them
txt1 = "ML"
txt2 = txt_suffix
for (i in 1:ncol(char_dtf))
    {
    recon_txt = paste("ACE technique: ", txt1, txt2, sep="")

    tmptipvals = char_dtf[,i]
    tmpnodevals = PICcont_means[,i]
    maxval = max(c(tmptipvals, tmpnodevals))
    texttxt = paste(recon_txt, "\nParameter: ", names(char_dtf)[i], "\nMaximum value = ", maxval, sep="") 

    plot_cont_vals_on_tree(chtr2, tmptipvals, tmpnodevals, texttxt)
    }

GLScont_means = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
GLScont_lower025 = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
GLScont_upper975 = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
for (i in 1:ncol(char_dtf))
    {
    GLScont = ace(char_dtf[,i], chtr2, type="continuous", method="GLS", CI=TRUE, corStruct = corBrownian(1, chtr2))
    GLScont_means[, i] = GLScont$ace
    GLScont_lower025[, i] = GLScont$CI95[, 1]
    GLScont_upper975[, i] = GLScont$CI95[, 2]
    }
GLScont_means = adf(GLScont_means)
names(GLScont_means) = names(char_dtf)
tree_dtf = prt(chtr2)
row.names(GLScont_means) = tree_dtf$node[ (length(chtr2$tip.label)+1): length(tree_dtf$node)]
GLScont_lower025 = adf(GLScont_lower025)
names(GLScont_lower025) = names(char_dtf)
row.names(GLScont_lower025) = row.names(GLScont_means)
GLScont_upper975 = adf(GLScont_upper975)
names(GLScont_upper975) = names(char_dtf)
row.names(GLScont_upper975) = row.names(GLScont_means)

write_table_good_w_rownames(GLScont_means, "GLScont_means_v8.txt")
write_table_good_w_rownames(GLScont_lower025, "GLScont_lower025_v8.txt")
write_table_good_w_rownames(GLScont_upper975, "GLScont_upper975_v8.txt")

# write the tree out to a table for Excel
tmptree_table = prt(chtr2)
tmptree_table$daughter_nds = as.character(tmptree_table$daughter_nds)
#tmptree_table = tmptree_table[, !(names(tmptree_table) == "daughter_nds")]
write_table_good(tmptree_table, "chtr2_v8.txt")

# plot tip & node values to check them
txt1 = "GLS"
txt2 = txt_suffix
for (i in 1:ncol(char_dtf))
    {
    recon_txt = paste("ACE technique: ", txt1, txt2, sep="")

    tmptipvals = char_dtf[,i]
    tmpnodevals = PICcont_means[,i]
    maxval = max(c(tmptipvals, tmpnodevals))
    texttxt = paste(recon_txt, "\nParameter: ", names(char_dtf)[i], "\nMaximum value = ", maxval, sep="") 

    plot_cont_vals_on_tree(chtr2, tmptipvals, tmpnodevals, texttxt)
    }

# Compare the estimates against each other with pairs plots
for (i in 1:ncol(char_dtf))
    {
    tmpcols = cbind(PICcont_means[,i], MLcont_means[,i], GLScont_means[,i])
    tmpcols = adf(tmpcols)
    names(tmpcols) = c("ACE_pic", "ACE_ML", "ACE_GLS")

    pairs.panels(tmpcols)

    titletxt = paste("param #", i, ": ", names(char_dtf)[i], sep="")
    title(titletxt, line=-1.35, outer=TRUE)
    }

#######################################################
# Reconstruct the ancestors -- discrete, plots
# discrete parameter characters
#######################################################
colnames_to_use = c(complex_params_names, clad_chars_names)
char_dtf = all_params[, colnames_to_use]
num_uniq = apply(char_dtf, 2, function(x) length(unique(x)))

# parsimony reconstruction
#
# https://stat.ethz.ch/pipermail/r-sig-phylo/2010-January/000565.html
# Hi Santiago,
# 
# Look at the ace function in the Ape package - to get a max parsimony  
# reconstruction set branch lengths =1 and run a maximum likelihood  
# estimation.  You can find further information here: http://www.r-phylo.org/wiki/HowTo/Ancestral_State_Reconstruction
# 

#x = dtf_to_phyDat(char_dtf)

# make a parsimony tree: 
pars_tr = chtr2
# set all branch lengths to 1
pars_tr$edge.length[!is.na(pars_tr$edge.length)] = 1
#ace(char_dtf[,1], pars_tr, type="discrete", method="ML", model="ER")

par(mfrow=c(1,2))
ACE_pars_estimates = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
for (i in 1:ncol(char_dtf))
    {
    (ACE_pars = ace(char_dtf[,i], pars_tr, type="discrete", method="ML", model="ER", ip=0.3))

    # convert to parsimony
    maxpar_node_vals = ACE_pars$lik.anc
    maxpar_node_maxes = apply(maxpar_node_vals, 1, max)

    maxpar_node_vals = adf(maxpar_node_vals)
    maxpar_node_MP = array(NA, dim=c(num_internal_nodes, 1))
    for (j in 1:nrow(maxpar_node_vals))
        {
        tmpTFs = maxpar_node_vals[j,] == maxpar_node_maxes[j]
        names(tmpTFs)
        colnum = which(tmpTFs)
        if (length(colnum) > 1)
            {
            tmplist = c()
            for (k in 1:length(colnum))
                {
                tmplist[k] = names(maxpar_node_vals)[colnum[k]]
                }
            tmpstr = paste(tmplist, collapse="&")
            maxpar_node_MP[j] = tmpstr
            } else {
            maxpar_node_MP[j] = names(maxpar_node_vals)[colnum]
            }
        }

    ACE_pars_estimates[, i] = maxpar_node_MP

    # plot results
    tree_ht = get_max_height_tree(chtr2)
    label_offset = 0.05*tree_ht
    plot(chtr2, show.node.label=FALSE, cex=0.9, no.margin = TRUE, x.lim=1.6*tree_ht, label.offset=label_offset)

    if (num_uniq[i] == 2)
        {
        colors = c("yellow", "blue")
        tmpcols = char_dtf[,i]
        tmpcols[tmpcols == 0] = colors[1]
        tmpcols[tmpcols == 1] = colors[2]
        }
    if (num_uniq[i] == 3)
        {
        colors = c("yellow", "blue", "red")
        tmpcols = char_dtf[,i]
        tmpcols[tmpcols == 0] = colors[1]
        tmpcols[tmpcols == 1] = colors[2]
        tmpcols[tmpcols == 2] = colors[3]
        }

    tiplabels(pch = 19, col = tmpcols, cex = 2)  #adj = 0, 
    legend(0, 3, legend=as.character(0:max(char_dtf[,i])), fill=colors)

    nodelabels(thermo=ACE_pars$lik.anc, cex=1, piecol=colors)

    titletxt = paste("ML, parsimony-like reconstruction of\n", names(char_dtf)[i], " character", sep="")
    title(titletxt, cex=0.5)

    }

ACE_pars_estimates = adf(ACE_pars_estimates)
names(ACE_pars_estimates) = names(char_dtf)
pars_tr_dtf = prt(pars_tr)
row.names(ACE_pars_estimates) = tree_dtf$node[ (length(pars_tr$tip.label)+1): length(pars_tr_dtf$node)]
#row.names(ACE_pars_estimates) = row.names(char_dtf)
(ACE_pars_estimates)

write_table_good(ACE_pars_estimates, "ACE_pars_estimates_v8.txt")

# Do ML reconstruction, with varying branch lengths (i.e. not branches = 1 like "parsimony")
par(mfrow=c(1,2))
ACE_ML_estimates = array(NA, dim=c(num_internal_nodes, ncol(char_dtf)) )
for (i in 1:ncol(char_dtf))
    {
    (ACE_ML = ace(char_dtf[,i], chtr2, type="discrete", method="ML", model="ER", ip=0.3))

    # convert to the single best ML estimate (plurality vote)
    maxML_node_vals = ACE_ML$lik.anc
    maxML_node_maxes = apply(maxML_node_vals, 1, max)

    maxML_node_vals = adf(maxML_node_vals)
    maxML_node_best = array(NA, dim=c(num_internal_nodes, 1))
    for (j in 1:nrow(maxML_node_vals))
        {
        tmpTFs = round(maxML_node_vals[j,], 3) == round(maxML_node_maxes[j], 3)
        names(tmpTFs)
        colnum = which(tmpTFs)
        if (length(colnum) > 1)
            {
            tmplist = c()
            for (k in 1:length(colnum))
                {
                tmplist[k] = names(maxML_node_vals)[colnum[k]]
                }
            tmpstr = paste(tmplist, collapse="&")
            maxML_node_best[j] = tmpstr
            } else {
            maxML_node_best[j] = names(maxML_node_vals)[colnum]
            }
        }

    #cbind(maxML_node_vals, maxML_node_maxes, maxML_node_best)

    ACE_ML_estimates[, i] = maxML_node_best

    # plot results    
    tree_ht = get_max_height_tree(chtr2)
    label_offset = 0.05*tree_ht
    plot(chtr2, show.node.label=FALSE, cex=0.9, no.margin = TRUE, x.lim=1.6*tree_ht, label.offset=label_offset)

    if (num_uniq[i] == 2)
        {
        colors = c("yellow", "blue")
        tmpcols = char_dtf[,i]
        tmpcols[tmpcols == 0] = colors[1]
        tmpcols[tmpcols == 1] = colors[2]
        }
    if (num_uniq[i] == 3)
        {
        colors = c("yellow", "blue", "red")
        tmpcols = char_dtf[,i]
        tmpcols[tmpcols == 0] = colors[1]
        tmpcols[tmpcols == 1] = colors[2]
        tmpcols[tmpcols == 2] = colors[3]
        }
    tiplabels(pch = 19, col = tmpcols, cex = 2)  #adj = 0, 
    legend(0, 3, legend=as.character(0:max(char_dtf[,i])), fill=colors)

    nodelabels(thermo=ACE_ML$lik.anc, cex=1, piecol=colors)

    titletxt = paste("ML, equal rates reconstruction of\n", names(char_dtf)[i], " character", sep="")
    title(titletxt, cex=0.5)
    }

ACE_ML_estimates = adf(ACE_ML_estimates)
names(ACE_ML_estimates) = names(char_dtf)
pars_tr_dtf = prt(pars_tr)
row.names(ACE_ML_estimates) = tree_dtf$node[ (length(pars_tr$tip.label)+1): length(pars_tr_dtf$node)]
#row.names(ACE_ML_estimates) = row.names(char_dtf)
(ACE_ML_estimates)

write_table_good(ACE_ML_estimates, "ACE_ML_estimates_v8.txt")

# Write out the trees we generates, also...
write.tree(tr, "tr.newick")
write.tree(tr2, "tr2.newick")
write.tree(chtr2, "chtr2.newick")