## DEC (Lagrange) biogeographic analysis using BioGeoBears ## Modified from the example on http://phylo.wikidot.com/biogeobears ## by Cam Webb 2016-03 # install.packages("BioGeoBEARS", dependencies=TRUE, # repos="http://cran.rstudio.com") # 1. Load and prep. library(optimx) library(FD) library(parallel) library(BioGeoBEARS) source("http://phylo.wdfiles.com/local--files/biogeobears/cladoRcpp.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_add_fossils_randomly_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_basics_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_calc_transition_matrices_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_classes_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_detection_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_DNA_cladogenesis_sim_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_extract_Qmat_COOmat_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_generics_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_models_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_on_multiple_trees_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_plots_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_readwrite_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_simulate_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_SSEsim_makePlots_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_SSEsim_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_stochastic_mapping_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_stratified_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_univ_model_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/calc_uppass_probs_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/calc_loglike_sp_v01.R") source("http://phylo.wdfiles.com/local--files/biogeobears/get_stratified_subbranch_top_downpass_likelihoods_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/runBSM_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/stochastic_map_given_inputs.R") source("http://phylo.wdfiles.com/local--files/biogeobears/summarize_BSM_tables_v1.R") source("http://phylo.wdfiles.com/local--files/biogeobears/BioGeoBEARS_traits_v1.R") calc_loglike_sp = compiler::cmpfun(calc_loglike_sp_prebyte) calc_independent_likelihoods_on_each_branch = compiler::cmpfun(calc_independent_likelihoods_on_each_branch_prebyte) # 2. phylogeny # Example Newick file for Hawaiian Psychotria (from Ree & Smith 2008) tr = read.tree("psychotria.new") plot(tr) title("Example Psychotria phylogeny from Ree & Smith (2008)") axisPhylo() # plots timescale # 3. geography # Example geography file for Hawaiian Psychotria (from Ree & Smith 2008) tipranges = getranges_from_LagrangePHYLIP(lgdata_fn="psychotria_geog.phy") tipranges # 4. Intitialize a default model (DEC model) BioGeoBEARS_run_object = define_BioGeoBEARS_run() BioGeoBEARS_run_object$trfn = "psychotria.new" BioGeoBEARS_run_object$geogfn = "psychotria_geog.phy" BioGeoBEARS_run_object$max_range_size = 4 # Other settings BioGeoBEARS_run_object$min_branchlength = 0.000001 BioGeoBEARS_run_object$include_null_range = TRUE BioGeoBEARS_run_object$speedup = TRUE BioGeoBEARS_run_object$use_optimx = TRUE BioGeoBEARS_run_object$num_cores_to_use = 1 BioGeoBEARS_run_object$force_sparse = FALSE BioGeoBEARS_run_object = readfiles_BioGeoBEARS_run(BioGeoBEARS_run_object) BioGeoBEARS_run_object$return_condlikes_table = TRUE BioGeoBEARS_run_object$calc_TTL_loglike_from_condlikes_table = TRUE BioGeoBEARS_run_object$calc_ancprobs = TRUE # 5. View it all and check BioGeoBEARS_run_object BioGeoBEARS_run_object$BioGeoBEARS_model_object BioGeoBEARS_run_object$BioGeoBEARS_model_object@params_table # Run this to check inputs. Read the error messages if you get them! check_BioGeoBEARS_run(BioGeoBEARS_run_object) # 6. Run it! resDEC = bears_optim_run(BioGeoBEARS_run_object) resDEC # 7. Plot it pdf("Psychotria_DEC_M0_unconstrained_v1.pdf", width=6, height=6) analysis_titletxt ="BioGeoBEARS DEC on Psychotria M0_unconstrained" scriptdir = np(system.file("extdata/a_scripts", package="BioGeoBEARS")) plot_BioGeoBEARS_results(resDEC, analysis_titletxt, addl_params=list("j"), plotwhat="text", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) plot_BioGeoBEARS_results(resDEC, analysis_titletxt, addl_params=list("j"), plotwhat="pie", label.offset=0.45, tipcex=0.7, statecex=0.7, splitcex=0.6, titlecex=0.8, plotsplits=TRUE, cornercoords_loc=scriptdir, include_null_range=TRUE, tr=tr, tipranges=tipranges) dev.off()