Phyllotis Phylogeny
MrBayes and RAxML
Options
Set Directory and Environment
setwd("D:/Phyllotis/PhyProj/PhyloProj/MrBayes/MrBayes")
Sys.setenv(PATH = paste(Sys.getenv("PATH"), "D:\\Phyllotis\\PhyProj\\PhyloProj\\MrBayes\\MrBayes;", sep = ":"))
#Sys.getenv("PATH")Load Libraries
suppressMessages(library(ape))
suppressMessages(library(ips))
suppressMessages(library(phangorn))
suppressMessages(library(geiger))
suppressMessages(library(dendextend))
suppressMessages(library(ggplot2))
suppressMessages(library(ggtree))
suppressMessages(library(phylobase))
suppressMessages(library(dplyr))Load Alignment
cytb = read.dna("Phyllotis.cytb.050817.fas", format = "fasta")
write.nexus.data(cytb , file = "cytb.nex",
format = "dna",
datablock = TRUE,
interleaved = FALSE,
gap = "-",
missing = "N")Explore Models of Evolution
Nucleotide substitution models
cytb1 = as.phyDat(cytb)
mt = modelTest(cytb1)## negative edges length changed to 0!## [1] "JC+I"
## [1] "JC+G"
## [1] "JC+G+I"
## [1] "F81+I"
## [1] "F81+G"
## [1] "F81+G+I"
## [1] "K80+I"
## [1] "K80+G"
## [1] "K80+G+I"
## [1] "HKY+I"
## [1] "HKY+G"
## [1] "HKY+G+I"
## [1] "SYM+I"
## [1] "SYM+G"
## [1] "SYM+G+I"
## [1] "GTR+I"
## [1] "GTR+G"
## [1] "GTR+G+I"env = attr(mt, "env")
ls(envir=env)## [1] "data" "F81" "F81+G" "F81+G+I"
## [5] "F81+I" "GTR" "GTR+G" "GTR+G+I"
## [9] "GTR+I" "HKY" "HKY+G" "HKY+G+I"
## [13] "HKY+I" "JC" "JC+G" "JC+G+I"
## [17] "JC+I" "K80" "K80+G" "K80+G+I"
## [21] "K80+I" "SYM" "SYM+G" "SYM+G+I"
## [25] "SYM+I" "tree_F81" "tree_F81+G" "tree_F81+G+I"
## [29] "tree_F81+I" "tree_GTR" "tree_GTR+G" "tree_GTR+G+I"
## [33] "tree_GTR+I" "tree_HKY" "tree_HKY+G" "tree_HKY+G+I"
## [37] "tree_HKY+I" "tree_JC" "tree_JC+G" "tree_JC+G+I"
## [41] "tree_JC+I" "tree_K80" "tree_K80+G" "tree_K80+G+I"
## [45] "tree_K80+I" "tree_SYM" "tree_SYM+G" "tree_SYM+G+I"
## [49] "tree_SYM+I"arrange(mt, BIC) #Models summary, sorted by BIC## Model df logLik AIC AICw AICc AICcw
## 1 GTR+G+I 591 -15453.21 32088.43 1.000000e+00 33356.08 1.000000e+00
## 2 HKY+G+I 587 -15492.01 32158.01 7.752190e-16 33399.58 3.576480e-10
## 3 GTR+G 590 -15514.30 32208.61 8.015366e-27 33469.69 2.137810e-25
## 4 HKY+G 586 -15602.69 32377.37 1.806541e-63 33612.49 2.090632e-56
## 5 SYM+G+I 588 -15825.64 32827.28 3.641691e-161 34075.32 6.596723e-157
## 6 SYM+G 587 -15892.62 32959.25 8.023903e-190 34200.81 3.701835e-184
## 7 K80+G+I 584 -15959.78 33087.57 1.095162e-217 34309.89 7.621189e-208
## 8 HKY+I 586 -16002.65 33177.30 3.578947e-237 34412.43 4.141761e-230
## 9 GTR+I 590 -15993.44 33166.88 6.565363e-235 34427.96 1.751074e-233
## 10 K80+G 583 -16020.27 33206.54 1.603505e-243 34422.51 2.675479e-232
## 11 SYM+I 587 -16249.55 33673.10 0.000000e+00 34914.66 0.000000e+00
## 12 K80+I 583 -16454.55 34075.11 0.000000e+00 35291.08 0.000000e+00
## 13 F81+G+I 586 -17135.75 35443.50 0.000000e+00 36678.62 0.000000e+00
## 14 F81+G 585 -17190.12 35550.25 0.000000e+00 36778.96 0.000000e+00
## 15 JC+G+I 583 -17495.58 36157.15 0.000000e+00 37373.12 0.000000e+00
## 16 JC+G 582 -17542.70 36249.41 0.000000e+00 37459.06 0.000000e+00
## 17 F81+I 585 -17654.71 36479.43 0.000000e+00 37708.13 0.000000e+00
## 18 JC+I 582 -17929.68 37023.36 0.000000e+00 38233.01 0.000000e+00
## 19 GTR 589 -18067.12 37312.23 0.000000e+00 38566.78 0.000000e+00
## 20 SYM 586 -18133.37 37438.74 0.000000e+00 38673.86 0.000000e+00
## 21 HKY 585 -18178.27 37526.53 0.000000e+00 38755.24 0.000000e+00
## 22 K80 582 -18389.43 37942.85 0.000000e+00 39152.50 0.000000e+00
## 23 F81 584 -19532.17 40232.35 0.000000e+00 41454.67 0.000000e+00
## 24 JC 581 -19793.86 40749.73 0.000000e+00 41953.08 0.000000e+00
## BIC
## 1 35068.42
## 2 35117.84
## 3 35183.55
## 4 35332.15
## 5 35792.14
## 6 35919.07
## 7 36032.26
## 8 36132.08
## 9 36141.83
## 10 36146.19
## 11 36632.92
## 12 37014.76
## 13 38398.28
## 14 38499.99
## 15 39096.80
## 16 39184.02
## 17 39429.16
## 18 39957.97
## 19 40282.14
## 20 40393.52
## 21 40476.27
## 22 40877.46
## 23 43177.04
## 24 43679.30Mr. Bayes
(Bayesian Tree)
Call MrBayes
mrbayes(cytb,
file = "cytb.nex",
lset = mrbayes.lset(nst = 6, rates = "invgamma"), ## GTR + Gamma + I
#ngammacat = 4, ## Default rate categories Gamma distribution
#prset, ## Default priors
mcmc = mrbayes.mcmc(ngen = 3500000),
#samplefreq = 100,
#printfreq = 100,
#nchains = 4,
#temp = 0.2,
#savebrlens = "yes",
burnin = 1000000,
contype = "allcompat", ## Strict consensus tree
run = TRUE)Plot Bayesian Tree
mrbayes.tree = read.nexus("cytb.nex.con.tre")
BaseTree = root(mrbayes.tree, outgroup = c(1, 2, 3))
p1 = ggtree(BaseTree,
ladderize=TRUE) +
geom_tiplab(cex = 1.75) +
geom_treescale(x=0.3, y=0.5, width=0.005,
offset=-1, fontsize = 3, color = "red") +
ggtitle("05-12-17 Phyllotis.cytb.050817.fas GTR+G+I MrBayes g3.5m b1m")
p1
Number of Trees
run1 = read.nexus("cytb.nex.run1.t")
run1## 7001 phylogenetic treesrun2 = read.nexus("cytb.nex.run2.t")
run2## 7001 phylogenetic treesBayesian Support Values
tree.mrbayes = read.beast("MrBayesSupp2.nex") #Previously extracted using FigTree
ggtree(tree.mrbayes, ladderize=TRUE) +
geom_text2(aes(subset=!isTip,
label=round(prob, 3),
hjust=-.3)) +
geom_tiplab(cex = 1.75) +
geom_treescale(x=0.3, y=0.5, width=0.005,
offset=-1, fontsize = 3, color = "red") +
ggtitle("05-12-17 Posteriors Phyllotis.cytb.050817.fas GTR+G+I MrBayes g3.5m b1m") 
Change Directory and Environment
setwd("D:/RAxML/standardRAxML/RAxML-8.0.3")
Sys.setenv(PATH = paste(Sys.getenv("PATH"), "D:\\RAxML\\standardRAxML\\RAxML-8.0.3", sep = ":"))
Sys.getenv("PATH")Load Alignment
cytb = read.dna("D:/RAxML/standardRAxML/RAxML-8.0.3/Phyllotis.cytb.050817.fas", format = "fasta")
write.nexus.data(cytb , file = "cytb.nex",
format = "dna",
datablock = TRUE,
interleaved = FALSE,
gap = "-",
missing = "N")Call RAxML
exec = "D:/RAxML/standardRAxML/RAxML-8.0.3/raxmlHPC.exe"
cytb.tr = raxml(cytb,
m = "GTRGAMMA",
f = "a",
N = 2,
outgroup = c("Graomys_domorum_AF159291",
"Graomys_griseoflavus_AF159290",
"Graomys_centralis_AF281207"),
p = 1976,
x = 1976,
k = TRUE,
threads = 2,
file = "cytb.nex",
exec = exec)Plot RAxML Best Tree
RAxML.tree = read.tree("D:/RAxML/standardRAxML/RAxML-8.0.3/RAxML_bestTree.cytb.nex")
p1 = ggtree(RAxML.tree,
ladderize=TRUE) +
geom_tiplab(cex = 1.75) +
geom_treescale(x=0.3, y=0.5, width=0.005,
offset=-1, fontsize = 3, color = "red") +
ggtitle("05-15-17 Phyllotis.cytb.050817.fas GTRGAMMA RAxML p1976")
p1 
RAxML Support Values
RAxML.Stree = read.tree("D:/RAxML/standardRAxML/RAxML-8.0.3/RAxML_bipartitions.cytb.nex")
p1 = ggtree(RAxML.Stree,
ladderize=TRUE) +
geom_tiplab(cex = 1.75) +
geom_text2(aes(subset=!isTip,
label=label,
hjust=-.3)) +
geom_treescale(x=0.3, y=0.5, width=0.005,
offset=-1, fontsize = 3, color = "red") +
ggtitle("05-15-17 Clade Support Phyllotis.cytb.050817.fas GTRGAMMA RAxML p1976")
p1 
Compare Trees
RAxML vs MrBayes
chronotr2 = chronopl(RAxML.tree,
lambda = 0,
age.min = 1) #just to convert format
dend = as.phylo(chronotr2)
dend = as.hclust(reorder(multi2di(dend)))
dend = as.dendrogram(dend)
chronotr = chronopl(mrbayes.tree,
lambda = 0,
age.min = 1) #just to convert format
dend2 = as.phylo(chronotr)
dend2 = as.hclust(reorder(multi2di(dend2)))
dend2 = as.dendrogram(dend2)
#dend_diff(dend, dend2)
dends = dendlist(dend, dend2)
x = dends %>% dendextend::untangle(method = "step2side")
entanglement(x)## [1] 0.0005502955#x %>% plot(main = paste("entanglement =", round(entanglement(x), 2)))
x %>% set("rank_branches") %>%
tanglegram(main_left = "RAxML",
cex_main_left = 3,
main_right = "MrBayes",
cex_main_right = 3,
columns_width = c(5,0.5,5),
lwd = 2,
lab.cex = 0.5,
dLeaf_left = -0.1,
dLeaf_right = 0.1,
margin_outer = 4,
margin_inner = 10,
color_lines = c(1:14),
common_subtrees_color_branches = FALSE)
Phyllotis xanthopygus Complex
tree = RAxML.tree
tree = groupClade(tree, node=c(505, 369))
p1 = ggtree(tree, ladderize=TRUE, aes(color=group)) +
geom_cladelabel(node=505,
label= "Altiplano Group",
align=T, color='darkgreen', angle = 270,
offset = -0.01, hjust = 0.5, geom = "text", fill = "gray") +
geom_cladelabel(node=369,
label= "North-South Group",
align=T, color='blue', angle = 270,
offset = -0.01, hjust = 0.5, geom = "text", fill = "gray") +
scale_color_manual(values=c("black", "darkgreen", "blue")) +
geom_point2(aes(subset=368), color="red", size=3) +
#geom_hilight(node=506, fill="yellow", alpha=0.4) +
ggtitle("Phyllotis xanthopygus Complex (red node)")
p1 
Altiplano Group
Ap.tree = as(RAxML.tree, "phylo4")
Ap.tree = subset(Ap.tree, node.subtree=505)
Ap.tree = as(Ap.tree, "phylo")
p2 = ggtree(Ap.tree,
ladderize=T) +
geom_tiplab(cex = 3)
p2 
North-South Group
Ap.tree = as(RAxML.tree, "phylo4")
Ap.tree = subset(Ap.tree, node.subtree=369)
Ap.tree = as(Ap.tree, "phylo")
p2 = ggtree(Ap.tree,
ladderize=T) +
geom_tiplab(cex = 3)
p2 
sessionInfo()## R version 3.3.3 (2017-03-06)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] dplyr_0.5.0 phylobase_0.8.4 ggtree_1.6.11
## [4] ggplot2_2.2.1.9000 dendextend_1.5.2 geiger_2.0.6
## [7] phangorn_2.2.0 ips_0.0-10 ape_4.1
## [10] knitr_1.15.1
##
## loaded via a namespace (and not attached):
## [1] viridis_0.4.0 httr_1.2.1 tidyr_0.6.2
## [4] jsonlite_1.4 viridisLite_0.2.0 foreach_1.4.3
## [7] bold_0.4.0 assertthat_0.2.0 stats4_3.3.3
## [10] taxize_0.8.4 robustbase_0.92-7 progress_1.1.2
## [13] backports_1.0.5 lattice_0.20-34 quadprog_1.5-5
## [16] uuid_0.1-2 digest_0.6.12 colorspace_1.3-2
## [19] htmltools_0.3.6 Matrix_1.2-8 plyr_1.8.4
## [22] XML_3.98-1.6 rncl_0.8.2 mvtnorm_1.0-6
## [25] scales_0.4.1 whisker_0.3-2 tibble_1.3.0
## [28] nnet_7.3-12 lazyeval_0.2.0 magrittr_1.5
## [31] mclust_5.2.3 evaluate_0.10 nlme_3.1-131
## [34] MASS_7.3-45 xml2_1.1.1 class_7.3-14
## [37] prettyunits_1.0.2 tools_3.3.3 data.table_1.10.4
## [40] trimcluster_0.1-2 stringr_1.2.0 kernlab_0.9-25
## [43] munsell_0.4.3 cluster_2.0.5 fpc_2.1-10
## [46] ade4_1.7-6 RNeXML_2.0.7 grid_3.3.3
## [49] iterators_1.0.8 subplex_1.2-2 igraph_1.0.1
## [52] labeling_0.3 rmarkdown_1.5 gtable_0.2.0
## [55] codetools_0.2-15 deSolve_1.14 flexmix_2.3-14
## [58] DBI_0.6-1 reshape_0.8.6 reshape2_1.4.2
## [61] R6_2.2.0 gridExtra_2.2.1 prabclus_2.2-6
## [64] fastmatch_1.1-0 rprojroot_1.2 modeltools_0.2-21
## [67] stringi_1.1.5 parallel_3.3.3 Rcpp_0.12.10
## [70] DEoptimR_1.0-8 diptest_0.75-7 coda_0.19-1