Data analysis workflow:
Correlated evolution of conspicuous coloration and burrowing in crayfish
Dylan Padilla, School of Life Sciences, Arizona State University, Tempe, AZ 85287, USA.
Library
library(ape)
library(caper)
library(EloChoice)
library(geiger)
library(irr)
library(phytools)
library(plotrix)
library(png)
library(reshape)
library(reshape2)
library(xtable)R.version> _
> platform aarch64-apple-darwin20
> arch aarch64
> os darwin20
> system aarch64, darwin20
> status
> major 4
> minor 3.3
> year 2024
> month 02
> day 29
> svn rev 86002
> language R
> version.string R version 4.3.3 (2024-02-29)
> nickname Angel Food CakesessionInfo()> R version 4.3.3 (2024-02-29)
> Platform: aarch64-apple-darwin20 (64-bit)
> Running under: macOS Sonoma 14.5
>
> Matrix products: default
> BLAS: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRblas.0.dylib
> LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
>
> locale:
> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> time zone: America/Phoenix
> tzcode source: internal
>
> attached base packages:
> [1] stats graphics grDevices utils datasets methods base
>
> other attached packages:
> [1] rmarkdown_2.26 xtable_1.8-4 reshape2_1.4.4 reshape_0.8.9 png_0.1-8
> [6] plotrix_3.8-4 irr_0.84.1 lpSolve_5.6.20 geiger_2.0.11 phytools_2.1-1
> [11] maps_3.4.2 EloChoice_0.29.4 caper_1.0.3 mvtnorm_1.2-4 MASS_7.3-60.0.1
> [16] ape_5.8
>
> loaded via a namespace (and not attached):
> [1] sass_0.4.9 generics_0.1.3 stringi_1.8.3 lattice_0.22-5
> [5] digest_0.6.35 magrittr_2.0.3 evaluate_0.23 grid_4.3.3
> [9] iterators_1.0.14 fastmap_1.1.1 jsonlite_1.8.8 foreach_1.5.2
> [13] doParallel_1.0.17 plyr_1.8.9 Matrix_1.6-5 deSolve_1.40
> [17] optimParallel_1.0-2 combinat_0.0-8 subplex_1.8 jquerylib_0.1.4
> [21] codetools_0.2-19 numDeriv_2016.8-1.1 mnormt_2.1.1 Rdpack_2.6
> [25] cli_3.6.2 rlang_1.1.3 psychotools_0.7-3 expm_0.999-9
> [29] scatterplot3d_0.3-44 rbibutils_2.2.16 cachem_1.0.8 yaml_2.3.8
> [33] tools_4.3.3 parallel_4.3.3 coda_0.19-4.1 fastmatch_1.1-4
> [37] R6_2.5.1 lifecycle_1.0.4 stringr_1.5.1 pkgconfig_2.0.3
> [41] bslib_0.7.0 glue_1.7.0 phangorn_2.11.1 Rcpp_1.0.12
> [45] clusterGeneration_1.3.8 highr_0.10 xfun_0.43 knitr_1.45
> [49] htmltools_0.5.8.1 nlme_3.1-164 igraph_2.0.3 compiler_4.3.3
> [53] quadprog_1.5-8## Loading data
craydata <- read.csv("col_data.csv", row.names = 1)
#head(craydata)
craydata_updated <- read.csv("col_data_updated.csv", row.names = 1)
head(craydata_updated)> family stern_phylogeny traditional_hobbs_burrow_category
> Cambarus_harti Cambaridea 1 primary_burrower
> Cherax_quadricarinatus Parastacidea 1 tertiary_burrower
> Euastacus_fleckeri Parastacidea 1 tertiary_burrower
> Parastacus_defossus Parastacidea 1 primary_burrower
> Cambarus_gentryi Cambaridea 1 primary_burrower
> Procambarus_pygmaeus Cambaridea 1 secondary_burrower
> combined_burrow_category conspic_color conspicous_body_color_present
> Cambarus_harti terrestrial_burrower yes yes
> Cherax_quadricarinatus aquatic_burrower yes yes
> Euastacus_fleckeri aquatic_burrower yes yes
> Parastacus_defossus terrestrial_burrower no no
> Cambarus_gentryi terrestrial_burrower yes yes
> Procambarus_pygmaeus terrestrial_burrower no no
> body_color_zg body_color_Observer1 body_color_Observer2 body_color_Observer3
> Cambarus_harti blue blue blue blue
> Cherax_quadricarinatus blue blue blue blue
> Euastacus_fleckeri blue blue blue blue
> Parastacus_defossus cryptic cryptic blue blue
> Cambarus_gentryi blue blue blue blue
> Procambarus_pygmaeus cryptic blue blue blue
> body_color_Observer4 body_color_Observer5 body_color_Observer6
> Cambarus_harti blue blue blue
> Cherax_quadricarinatus blue blue blue
> Euastacus_fleckeri blue blue blue
> Parastacus_defossus blue blue blue
> Cambarus_gentryi cryptic blue blue
> Procambarus_pygmaeus cryptic blue blue
> body_color_Observer7 body_color_Observer8 body_color_Observer9
> Cambarus_harti blue blue blue
> Cherax_quadricarinatus blue blue blue
> Euastacus_fleckeri blue blue blue
> Parastacus_defossus blue blue blue
> Cambarus_gentryi blue blue blue
> Procambarus_pygmaeus blue blue blue
> body_color_Observer10 conspicious_claw_color_present claw_color body_color_red
> Cambarus_harti blue yes blue 0
> Cherax_quadricarinatus blue yes blue 0
> Euastacus_fleckeri blue yes blue 0
> Parastacus_defossus blue no no 0
> Cambarus_gentryi blue yes blue 0
> Procambarus_pygmaeus blue no no 0
> body_color_blue body_color_orange body_color_purple body_color_unpigmented
> Cambarus_harti 1 0 0 0
> Cherax_quadricarinatus 1 0 0 0
> Euastacus_fleckeri 1 0 0 0
> Parastacus_defossus 0 0 0 0
> Cambarus_gentryi 1 0 0 0
> Procambarus_pygmaeus 0 0 0 0
> claw_color_red claw_color_blue claw_color_orange claw_color_purple
> Cambarus_harti 0 1 0 0
> Cherax_quadricarinatus 0 1 0 0
> Euastacus_fleckeri 0 1 0 0
> Parastacus_defossus 0 0 0 0
> Cambarus_gentryi 0 1 0 0
> Procambarus_pygmaeus 0 0 0 0length(unique(craydata_updated$body_color_zg))> [1] 17length(unique(craydata_updated$body_color_Observer10))> [1] 14craytree <- read.tree("stern.2017.ml.new.tree")
craytree>
> Phylogenetic tree with 466 tips and 465 internal nodes.
>
> Tip labels:
> Samastacus_spinifrons, Virilastacus_araucanius, Virilastacus_retamali, Virilastacus_rucapihuelensis, Parastacus_nicoleti, Parastacus_pugnax, ...
>
> Rooted; includes branch lengths.## Creating a table for reliability analysis
names(craydata_updated)> [1] "family" "stern_phylogeny"
> [3] "traditional_hobbs_burrow_category" "combined_burrow_category"
> [5] "conspic_color" "conspicous_body_color_present"
> [7] "body_color_zg" "body_color_Observer1"
> [9] "body_color_Observer2" "body_color_Observer3"
> [11] "body_color_Observer4" "body_color_Observer5"
> [13] "body_color_Observer6" "body_color_Observer7"
> [15] "body_color_Observer8" "body_color_Observer9"
> [17] "body_color_Observer10" "conspicious_claw_color_present"
> [19] "claw_color" "body_color_red"
> [21] "body_color_blue" "body_color_orange"
> [23] "body_color_purple" "body_color_unpigmented"
> [25] "claw_color_red" "claw_color_blue"
> [27] "claw_color_orange" "claw_color_purple"head(craydata_updated)> family stern_phylogeny traditional_hobbs_burrow_category
> Cambarus_harti Cambaridea 1 primary_burrower
> Cherax_quadricarinatus Parastacidea 1 tertiary_burrower
> Euastacus_fleckeri Parastacidea 1 tertiary_burrower
> Parastacus_defossus Parastacidea 1 primary_burrower
> Cambarus_gentryi Cambaridea 1 primary_burrower
> Procambarus_pygmaeus Cambaridea 1 secondary_burrower
> combined_burrow_category conspic_color conspicous_body_color_present
> Cambarus_harti terrestrial_burrower yes yes
> Cherax_quadricarinatus aquatic_burrower yes yes
> Euastacus_fleckeri aquatic_burrower yes yes
> Parastacus_defossus terrestrial_burrower no no
> Cambarus_gentryi terrestrial_burrower yes yes
> Procambarus_pygmaeus terrestrial_burrower no no
> body_color_zg body_color_Observer1 body_color_Observer2 body_color_Observer3
> Cambarus_harti blue blue blue blue
> Cherax_quadricarinatus blue blue blue blue
> Euastacus_fleckeri blue blue blue blue
> Parastacus_defossus cryptic cryptic blue blue
> Cambarus_gentryi blue blue blue blue
> Procambarus_pygmaeus cryptic blue blue blue
> body_color_Observer4 body_color_Observer5 body_color_Observer6
> Cambarus_harti blue blue blue
> Cherax_quadricarinatus blue blue blue
> Euastacus_fleckeri blue blue blue
> Parastacus_defossus blue blue blue
> Cambarus_gentryi cryptic blue blue
> Procambarus_pygmaeus cryptic blue blue
> body_color_Observer7 body_color_Observer8 body_color_Observer9
> Cambarus_harti blue blue blue
> Cherax_quadricarinatus blue blue blue
> Euastacus_fleckeri blue blue blue
> Parastacus_defossus blue blue blue
> Cambarus_gentryi blue blue blue
> Procambarus_pygmaeus blue blue blue
> body_color_Observer10 conspicious_claw_color_present claw_color body_color_red
> Cambarus_harti blue yes blue 0
> Cherax_quadricarinatus blue yes blue 0
> Euastacus_fleckeri blue yes blue 0
> Parastacus_defossus blue no no 0
> Cambarus_gentryi blue yes blue 0
> Procambarus_pygmaeus blue no no 0
> body_color_blue body_color_orange body_color_purple body_color_unpigmented
> Cambarus_harti 1 0 0 0
> Cherax_quadricarinatus 1 0 0 0
> Euastacus_fleckeri 1 0 0 0
> Parastacus_defossus 0 0 0 0
> Cambarus_gentryi 1 0 0 0
> Procambarus_pygmaeus 0 0 0 0
> claw_color_red claw_color_blue claw_color_orange claw_color_purple
> Cambarus_harti 0 1 0 0
> Cherax_quadricarinatus 0 1 0 0
> Euastacus_fleckeri 0 1 0 0
> Parastacus_defossus 0 0 0 0
> Cambarus_gentryi 0 1 0 0
> Procambarus_pygmaeus 0 0 0 0rel.dat <- cbind(rownames(craydata_updated), craydata_updated$combined_burrow_category, craydata_updated[8:17])
colnames(rel.dat)[c(1, 2)] <- c("species", "combined_burrow_category")
head(rel.dat)> species combined_burrow_category body_color_Observer1
> Cambarus_harti Cambarus_harti terrestrial_burrower blue
> Cherax_quadricarinatus Cherax_quadricarinatus aquatic_burrower blue
> Euastacus_fleckeri Euastacus_fleckeri aquatic_burrower blue
> Parastacus_defossus Parastacus_defossus terrestrial_burrower cryptic
> Cambarus_gentryi Cambarus_gentryi terrestrial_burrower blue
> Procambarus_pygmaeus Procambarus_pygmaeus terrestrial_burrower blue
> body_color_Observer2 body_color_Observer3 body_color_Observer4
> Cambarus_harti blue blue blue
> Cherax_quadricarinatus blue blue blue
> Euastacus_fleckeri blue blue blue
> Parastacus_defossus blue blue blue
> Cambarus_gentryi blue blue cryptic
> Procambarus_pygmaeus blue blue cryptic
> body_color_Observer5 body_color_Observer6 body_color_Observer7
> Cambarus_harti blue blue blue
> Cherax_quadricarinatus blue blue blue
> Euastacus_fleckeri blue blue blue
> Parastacus_defossus blue blue blue
> Cambarus_gentryi blue blue blue
> Procambarus_pygmaeus blue blue blue
> body_color_Observer8 body_color_Observer9 body_color_Observer10
> Cambarus_harti blue blue blue
> Cherax_quadricarinatus blue blue blue
> Euastacus_fleckeri blue blue blue
> Parastacus_defossus blue blue blue
> Cambarus_gentryi blue blue blue
> Procambarus_pygmaeus blue blue bluemat <- matrix(NA, nrow = nrow(rel.dat), ncol = 10, byrow = FALSE)
for(j in 1:10){
for(i in 1:nrow(rel.dat)){
dat <- rel.dat[3:ncol(rel.dat)]
obs <- dat[i, j]
if(obs == "cryptic" | obs == "unpigmented"){
mat[i, j] <- "cryptic"
}
else{
mat[i, j] <- "conspicuous"
}
}
}
colnames(mat) <- names(dat)
head(mat)> body_color_Observer1 body_color_Observer2 body_color_Observer3 body_color_Observer4
> [1,] "conspicuous" "conspicuous" "conspicuous" "conspicuous"
> [2,] "conspicuous" "conspicuous" "conspicuous" "conspicuous"
> [3,] "conspicuous" "conspicuous" "conspicuous" "conspicuous"
> [4,] "cryptic" "conspicuous" "conspicuous" "conspicuous"
> [5,] "conspicuous" "conspicuous" "conspicuous" "cryptic"
> [6,] "conspicuous" "conspicuous" "conspicuous" "cryptic"
> body_color_Observer5 body_color_Observer6 body_color_Observer7 body_color_Observer8
> [1,] "conspicuous" "conspicuous" "conspicuous" "conspicuous"
> [2,] "conspicuous" "conspicuous" "conspicuous" "conspicuous"
> [3,] "conspicuous" "conspicuous" "conspicuous" "conspicuous"
> [4,] "conspicuous" "conspicuous" "conspicuous" "conspicuous"
> [5,] "conspicuous" "conspicuous" "conspicuous" "conspicuous"
> [6,] "conspicuous" "conspicuous" "conspicuous" "conspicuous"
> body_color_Observer9 body_color_Observer10
> [1,] "conspicuous" "conspicuous"
> [2,] "conspicuous" "conspicuous"
> [3,] "conspicuous" "conspicuous"
> [4,] "conspicuous" "conspicuous"
> [5,] "conspicuous" "conspicuous"
> [6,] "conspicuous" "conspicuous"## Computing Fleiss Kappa
kappam.fleiss(mat, detail = TRUE)> Fleiss' Kappa for m Raters
>
> Subjects = 401
> Raters = 10
> Kappa = 0.821
>
> z = 110
> p-value = 0
>
> Kappa z p.value
> conspicuous 0.821 110.293 0.000
> cryptic 0.821 110.293 0.000rel.tab <- cbind(rel.dat[1:2], mat)
names(rel.tab)> [1] "species" "combined_burrow_category" "body_color_Observer1"
> [4] "body_color_Observer2" "body_color_Observer3" "body_color_Observer4"
> [7] "body_color_Observer5" "body_color_Observer6" "body_color_Observer7"
> [10] "body_color_Observer8" "body_color_Observer9" "body_color_Observer10"rownames(rel.tab) <- rel.tab$species
dim(rel.tab)> [1] 401 12mlt <- melt(rel.tab, id.vars = c("species", "combined_burrow_category"), measure.vars = names(rel.tab[3:12]), direction = "long", value.name = "color", variable.name = "raterID")
mlt <- mlt[1:nrow(rel.tab), ]
rownames(mlt) <- mlt$species
head(mlt)> species combined_burrow_category raterID color
> Cambarus_harti Cambarus_harti terrestrial_burrower body_color_Observer1 conspicuous
> Cherax_quadricarinatus Cherax_quadricarinatus aquatic_burrower body_color_Observer1 conspicuous
> Euastacus_fleckeri Euastacus_fleckeri aquatic_burrower body_color_Observer1 conspicuous
> Parastacus_defossus Parastacus_defossus terrestrial_burrower body_color_Observer1 cryptic
> Cambarus_gentryi Cambarus_gentryi terrestrial_burrower body_color_Observer1 conspicuous
> Procambarus_pygmaeus Procambarus_pygmaeus terrestrial_burrower body_color_Observer1 conspicuoustail(mlt)> species combined_burrow_category raterID color
> Procambarus_leitheuseri Procambarus_leitheuseri aquatic_burrower body_color_Observer1 cryptic
> Procambarus_milleri Procambarus_milleri aquatic_burrower body_color_Observer1 cryptic
> Procambarus_orcinus Procambarus_orcinus aquatic_burrower body_color_Observer1 cryptic
> Procambarus_pallidus Procambarus_pallidus aquatic_burrower body_color_Observer1 cryptic
> Procambarus_xochitlanae Procambarus_xochitlanae aquatic_burrower body_color_Observer1 cryptic
> Troglocambarus_maclanei Troglocambarus_maclanei aquatic_burrower body_color_Observer1 crypticmlt$combined_burrow_category[mlt$combined_burrow_category == "terrestrial_burrower"] <- "semi-terrestrial"
mlt$combined_burrow_category[mlt$combined_burrow_category == "aquatic_burrower"] <- "aquatic"
mlt$combined_burrow_category <- as.factor(mlt$combined_burrow_category)
levels(mlt$combined_burrow_category)> [1] "aquatic" "semi-terrestrial"## Pruning the data and tree
chk <- name.check(craytree, craydata)
summary(chk)> 62 taxa are present in the tree but not the data:
> Astacus_leptodactylus,
> Austropotamobius_italicus_carinthiacus,
> Austropotamobius_italicus_carsicus,
> Austropotamobius_italicus_italicus,
> Austropotamobius_torrentium,
> Cambarellus_blacki,
> ....
> 2 taxa are present in the data but not the tree:
> Austropotamobius_torrentium_torrentium,
> Lacunicambarus_miltus
>
> To see complete list of mis-matched taxa, print object.pruned.tree <- drop.tip(craytree, chk$tree_not_data)
pruned.data <- craydata[!(rownames(craydata) %in% chk$data_not_tree),, drop = FALSE]
name.check(pruned.tree, pruned.data)> [1] "OK"bodycol <- setNames(as.factor(pruned.data$body_color), rownames(pruned.data))
head(bodycol)> Cambarus_gentryi Cambarus_harti Cambarus_striatus
> blue blue blue
> Cherax_destructor_albidus Cherax_destructor_destructor Cherax_quadricarinatus
> blue blue blue
> Levels: blue blue_orange_red blue_purple blue_red no orange orange_red purple red unpigmentedburrowing <- setNames(as.factor(pruned.data$combined_burrow_category), rownames(pruned.data))
head(burrowing)> Cambarus_gentryi Cambarus_harti Cambarus_striatus
> terrestrial_burrower terrestrial_burrower terrestrial_burrower
> Cherax_destructor_albidus Cherax_destructor_destructor Cherax_quadricarinatus
> terrestrial_burrower terrestrial_burrower aquatic_burrower
> Levels: aquatic_burrower terrestrial_burrower## Pruning the data and tree (updated)
chk <- name.check(craytree, mlt)
summary(chk)> 67 taxa are present in the tree but not the data:
> Astacus_leptodactylus,
> Austropotamobius_italicus_carinthiacus,
> Austropotamobius_italicus_carsicus,
> Austropotamobius_italicus_italicus,
> Austropotamobius_torrentium,
> Cambarellus_blacki,
> ....
> 2 taxa are present in the data but not the tree:
> Austropotamobius_torrentium_torrentium,
> Lacunicambarus_miltus
>
> To see complete list of mis-matched taxa, print object.pruned.tree <- drop.tip(craytree, chk$tree_not_data)
pruned.data <- mlt[!(rownames(mlt) %in% chk$data_not_tree),, drop = FALSE]
name.check(pruned.tree, pruned.data)> [1] "OK"bodycol <- setNames(as.factor(pruned.data$color), rownames(pruned.data))
head(bodycol)> Cambarus_harti Cherax_quadricarinatus Euastacus_fleckeri Parastacus_defossus
> conspicuous conspicuous conspicuous cryptic
> Cambarus_gentryi Procambarus_pygmaeus
> conspicuous conspicuous
> Levels: conspicuous crypticburrowing <- setNames(as.factor(pruned.data$combined_burrow_category), rownames(pruned.data))
head(burrowing)> Cambarus_harti Cherax_quadricarinatus Euastacus_fleckeri Parastacus_defossus
> semi-terrestrial aquatic aquatic semi-terrestrial
> Cambarus_gentryi Procambarus_pygmaeus
> semi-terrestrial semi-terrestrial
> Levels: aquatic semi-terrestrial## Fitting models of correlated evolution
interdependent.model <- fitPagel(pruned.tree, bodycol, burrowing)
print(interdependent.model)>
> Pagel's binary character correlation test:
>
> Assumes "ARD" substitution model for both characters
>
> Independent model rate matrix:
> conspicuous|aquatic conspicuous|semi-terrestrial cryptic|aquatic
> conspicuous|aquatic -0.042822 0.004817 0.038005
> conspicuous|semi-terrestrial 0.003155 -0.041160 0.000000
> cryptic|aquatic 0.009393 0.000000 -0.014209
> cryptic|semi-terrestrial 0.000000 0.009393 0.003155
> cryptic|semi-terrestrial
> conspicuous|aquatic 0.000000
> conspicuous|semi-terrestrial 0.038005
> cryptic|aquatic 0.004817
> cryptic|semi-terrestrial -0.012547
>
> Dependent (x & y) model rate matrix:
> conspicuous|aquatic conspicuous|semi-terrestrial cryptic|aquatic
> conspicuous|aquatic -0.079303 0.020697 0.058606
> conspicuous|semi-terrestrial 0.000000 -0.013729 0.000000
> cryptic|aquatic 0.008131 0.000000 -0.009887
> cryptic|semi-terrestrial 0.000000 0.004680 0.005129
> cryptic|semi-terrestrial
> conspicuous|aquatic 0.000000
> conspicuous|semi-terrestrial 0.013729
> cryptic|aquatic 0.001755
> cryptic|semi-terrestrial -0.009809
>
> Model fit:
> log-likelihood AIC
> independent -339.7840 687.5680
> dependent -323.8532 663.7064
>
> Hypothesis test result:
> likelihood-ratio: 31.8617
> p-value: 2.04178e-06
>
> Model fitting method used was fitMkdependent.bodycol <- fitPagel(pruned.tree, bodycol, burrowing, dep.var = "y")
print(dependent.bodycol)>
> Pagel's binary character correlation test:
>
> Assumes "ARD" substitution model for both characters
>
> Independent model rate matrix:
> conspicuous|aquatic conspicuous|semi-terrestrial cryptic|aquatic
> conspicuous|aquatic -0.042822 0.004817 0.038005
> conspicuous|semi-terrestrial 0.003155 -0.041160 0.000000
> cryptic|aquatic 0.009393 0.000000 -0.014209
> cryptic|semi-terrestrial 0.000000 0.009393 0.003155
> cryptic|semi-terrestrial
> conspicuous|aquatic 0.000000
> conspicuous|semi-terrestrial 0.038005
> cryptic|aquatic 0.004817
> cryptic|semi-terrestrial -0.012547
>
> Dependent (y only) model rate matrix:
> conspicuous|aquatic conspicuous|semi-terrestrial cryptic|aquatic
> conspicuous|aquatic -0.054462 0.033937 0.020525
> conspicuous|semi-terrestrial 0.000000 -0.020525 0.000000
> cryptic|aquatic 0.006429 0.000000 -0.008112
> cryptic|semi-terrestrial 0.000000 0.006429 0.005186
> cryptic|semi-terrestrial
> conspicuous|aquatic 0.000000
> conspicuous|semi-terrestrial 0.020525
> cryptic|aquatic 0.001682
> cryptic|semi-terrestrial -0.011616
>
> Model fit:
> log-likelihood AIC
> independent -339.7840 687.5680
> dependent -327.2957 666.5914
>
> Hypothesis test result:
> likelihood-ratio: 24.9766
> p-value: 3.77049e-06
>
> Model fitting method used was fitMkdependent.burrowing <- fitPagel(pruned.tree, bodycol, burrowing, dep.var = "x")
print(dependent.burrowing)>
> Pagel's binary character correlation test:
>
> Assumes "ARD" substitution model for both characters
>
> Independent model rate matrix:
> conspicuous|aquatic conspicuous|semi-terrestrial cryptic|aquatic
> conspicuous|aquatic -0.042822 0.004817 0.038005
> conspicuous|semi-terrestrial 0.003155 -0.041160 0.000000
> cryptic|aquatic 0.009393 0.000000 -0.014209
> cryptic|semi-terrestrial 0.000000 0.009393 0.003155
> cryptic|semi-terrestrial
> conspicuous|aquatic 0.000000
> conspicuous|semi-terrestrial 0.038005
> cryptic|aquatic 0.004817
> cryptic|semi-terrestrial -0.012547
>
> Dependent (x only) model rate matrix:
> conspicuous|aquatic conspicuous|semi-terrestrial cryptic|aquatic
> conspicuous|aquatic -0.073189 0.004861 0.068328
> conspicuous|semi-terrestrial 0.003147 -0.207910 0.000000
> cryptic|aquatic 0.007214 0.000000 -0.012074
> cryptic|semi-terrestrial 0.000000 0.083837 0.003147
> cryptic|semi-terrestrial
> conspicuous|aquatic 0.000000
> conspicuous|semi-terrestrial 0.204763
> cryptic|aquatic 0.004861
> cryptic|semi-terrestrial -0.086984
>
> Model fit:
> log-likelihood AIC
> independent -339.7840 687.5680
> dependent -331.0367 674.0734
>
> Hypothesis test result:
> likelihood-ratio: 17.4946
> p-value: 0.000158891
>
> Model fitting method used was fitMkanova(dependent.burrowing, dependent.bodycol, interdependent.model)> log(L) d.f. AIC weight
> independent -339.7840 4 687.5680 5.301580e-06
> dependent.burrowing -331.0367 6 674.0734 4.515618e-03
> dependent.bodycol -327.2957 6 666.5914 1.902910e-01
> interdependent.model -323.8532 8 663.7064 8.051881e-01##xtable(anova(dependent.burrowing, dependent.bodycol, interdependent.model), digits = 3)
## Plotting model
plot(interdependent.model, fit.signif = 4, cex.main = 1,
cex.sub = 0.8, cex.traits = 0.7, color = TRUE,
cex.rates = 0.7, lwd = 1, lwd.by.rate = TRUE)
## Trees facing to each other
layout(matrix(c(1, 1, 2, 2,
1, 1, 3, 3,
1, 1, 3, 3,
1, 1, 2, 0,
1, 1, 2, 2), nrow = 5, ncol = 4, byrow = TRUE))
c <- pruned.data$species[pruned.data$color == "cryptic"]
cons <- pruned.data$species[pruned.data$color == "conspicuous"]
## paint the edges
tt <- paintBranches(pruned.tree, edge = sapply(cons, match, pruned.tree$tip.label),
state = "conspicuous", anc.state = "cryptic")
cols1 <- setNames(c("orange", "black"), c("conspicuous", "cryptic"))
## Body color tree
plot(tt, cols1, ftype = "off", lwd = 2, offset = 0.4)
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
h <- max(nodeHeights(tt))
points(pp$xx[1:Ntip(tt)] + 0.02*h, pp$yy[1:Ntip(tt)], pch = 16, col = cols1[bodycol[tt$tip.label]], cex = 0.6)
legend("topleft", legend = c("conspicuous", "cryptic"), pch = 21, pt.bg = cols1,
pt.cex = 1.5, bty = "n", col = "transparent", cex = 1.5)
cryptic <- readPNG("imgs/cryptic.png")
conspicuous <- readPNG("imgs/conspicuous.png")
rasterImage(cryptic, 65, 300, 170, 380)
rasterImage(conspicuous, 60, 213, 168, 295)
## Burrowing tree
a <- pruned.data$species[pruned.data$combined_burrow_category == "aquatic"]
t <- pruned.data$species[pruned.data$combined_burrow_category == "semi-terrestrial"]
## paint the edges
bt <- paintBranches(pruned.tree, edge = sapply(a, match, pruned.tree$tip.label),
state = "aquatic", anc.state = "semi-terrestrial")
bt>
> Phylogenetic tree with 399 tips and 398 internal nodes.
>
> Tip labels:
> Samastacus_spinifrons, Virilastacus_araucanius, Virilastacus_retamali, Virilastacus_rucapihuelensis, Parastacus_nicoleti, Parastacus_pugnax, ...
>
> The tree includes a mapped, 2-state discrete character
> with states:
> aquatic, semi-terrestrial
>
> Rooted; includes branch lengths.cols2 <- setNames(c("#2F5597", "#956D45"), c("aquatic", "semi-terrestrial"))
## Burrowing tree
plot(bt, cols2, ftype = "off", lwd = 2, offset = 0.4, direction = "leftwards")
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
h <- max(nodeHeights(bt))
points(pp$xx[1:Ntip(bt)] - 0.02*h, pp$yy[1:Ntip(bt)], pch = 16,
col = cols2[burrowing[bt$tip.label]], cex = 0.6)
legend("topright", legend = c("aquatic burrowing", "semi-terrestrial burrowing"), pch = 21, pt.bg = cols2, pt.cex = 1.5, bty = "n", col = "transparent", cex = 1.5)
burrow <- readPNG("imgs/burrow3.png")
rasterImage(burrow, 70, 220, 260, 380)
layout(matrix(c(1, 1, 1, 2, 2, 2,
3, 3, 1, 2, 4, 4,
3, 3, 1, 2, 4, 4,
3, 3, 1, 2, 5, 5,
3, 3, 1, 2, 5, 5,
0, 0, 1, 2, 0, 0,
1, 1, 1, 2, 2, 2,
1, 1, 1, 2, 2, 2), nrow = 8, ncol = 6, byrow = TRUE))
## Body color tree
plot(tt, cols1, ftype = "off", lwd = 2, offset = 0.4)
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
h <- max(nodeHeights(tt))
points(pp$xx[1:Ntip(tt)] + 0.02*h, pp$yy[1:Ntip(tt)], pch = 16,
col = cols1[bodycol[tt$tip.label]], cex = 0.6)
legend("topleft", legend = c("conspicuous", "cryptic"), pch = 21, pt.bg = cols1,
pt.cex = 1.5, bty = "n", col = "transparent", cex = 1.5)
## Burrowing tree
plot(bt, cols2, ftype = "off", lwd = 2, offset = 0.4, direction = "leftwards")
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
h <- max(nodeHeights(bt))
points(pp$xx[1:Ntip(bt)] - 0.02*h, pp$yy[1:Ntip(bt)], pch = 16,
col = cols2[burrowing[bt$tip.label]], cex = 0.6)
legend("topright", legend = c("aquatic burrowing", "semi-terrestrial burrowing"), pch = 21, pt.bg = cols2, pt.cex = 1.5, bty = "n", col = "transparent", cex = 1.5)
## Adding pie charts
craydata_updated <- read.csv("col_data_updated.csv", row.names = 1)
unique(craydata_updated$body_color_zg)> [1] "blue" "cryptic" "blue_purple"
> [4] "red" "red_blue" "blue_red"
> [7] "cryptic_blue" "cryptic_blue " "blue_orange_red"
> [10] "cryptic_blue_orange_red" "cryptic_orange" "orange"
> [13] "orange_red" "orange_blue" "orange_cryptic"
> [16] "purple" "unpigmented"craytree <- read.tree("stern.2017.ml.new.tree")
## Changing levels of the factors "body_color"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "no"] <- "c"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "blue"] <- "b"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "orange"] <- "o"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "purple"] <- "p"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "red"] <- "r"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "blue_red"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "red_blue"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "cryptic_blue"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "cryptic_blue "] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "cryptic_blue_orange_red"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "cryptic_orange"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "orange_blue"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "orange_red"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "orange_cryptic"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "cryptic"] <- "c"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "blue_orange_red"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "blue_purple"] <- "g"
craydata_updated$body_color_zg[craydata_updated$body_color_zg == "unpigmented"] <- "c"
unique(craydata_updated$body_color_zg)> [1] "b" "c" "g" "r" "o" "p"craydata_updated$body_color_zg <- as.factor(craydata_updated$body_color_zg)
levels(craydata_updated$body_color_zg)> [1] "b" "c" "g" "o" "p" "r"craydata_updated$combined_burrow_category[craydata_updated$combined_burrow_category == "terrestrial_burrower"] <- "t"
craydata_updated$combined_burrow_category[craydata_updated$combined_burrow_category == "aquatic_burrower"] <- "a"
craydata_updated$combined_burrow_category <- as.factor(craydata_updated$combined_burrow_category)
levels(craydata_updated$combined_burrow_category)> [1] "a" "t"## Pruning the data and tree
chk <- name.check(craytree, craydata_updated)
summary(chk)> 67 taxa are present in the tree but not the data:
> Astacus_leptodactylus,
> Austropotamobius_italicus_carinthiacus,
> Austropotamobius_italicus_carsicus,
> Austropotamobius_italicus_italicus,
> Austropotamobius_torrentium,
> Cambarellus_blacki,
> ....
> 2 taxa are present in the data but not the tree:
> Austropotamobius_torrentium_torrentium,
> Lacunicambarus_miltus
>
> To see complete list of mis-matched taxa, print object.pruned.tree <- drop.tip(craytree, chk$tree_not_data)
pruned.data <- craydata_updated[!(rownames(craydata_updated) %in% chk$data_not_tree),, drop = FALSE]
name.check(pruned.tree, pruned.data)> [1] "OK"colsf <- setNames(c("skyblue", "black", "green4", "orange", "#6600CC", "red"), c("b", "c", "g", "o", "p", "r"))
par(mar = c(5, 2, 2, 0))
pie(table(pruned.data$body_color_zg), border = FALSE, col = colsf, labels = "", radius = 0.5)
text(0, 0.8, paste("Combined \n n = ", length(pruned.data$body_color_zg), sep = ""), font = 2)
legend("left", paste(round((table(pruned.data$body_color_zg)/length(pruned.data$body_color_zg)*100), 1), rep("%", length(unique(pruned.data$body_color_zg)), sep = "")), fill = colsf, bty = "n", cex = 0.8, border = FALSE)
par(mar = c(4, 2, 1.3, 0))
pie(table(pruned.data$body_color_zg[pruned.data$combined_burrow_category == "a"]), border = FALSE, labels = "", main = paste("Aquatic burrowing \n n = ", length(pruned.data$body_color_zg[pruned.data$combined_burrow_category == "a"]), sep = ""), col = colsf, radius = 0.98, cex.main = 0.8)
legend("left", paste(round((table(pruned.data$body_color_zg[pruned.data$combined_burrow_category == "a"])/length(pruned.data$body_color_zg[pruned.data$combined_burrow_category == "a"]))*100, 1), rep("%", length(unique(pruned.data$body_color_zg))), sep = ""), fill = colsf, bty = "n", cex = 0.8, border = FALSE)
pie(table(pruned.data$body_color_zg[pruned.data$combined_burrow_category == "t"]), border = FALSE, main = paste("Semi-terrestrial burrowing \n n = ", length(pruned.data$body_color_zg[pruned.data$combined_burrow_category == "t"]), sep = ""), col = colsf, labels = "", radius = 0.98, cex.main = 0.8)
legend("left", paste(round((table(pruned.data$body_color_zg[pruned.data$combined_burrow_category == "t"])/length(pruned.data$body_color_zg[pruned.data$combined_burrow_category == "t"]))*100, 1), rep("%", length(unique(pruned.data$body_color_zg))), sep = ""), fill = colsf, bty = "n", cex = 0.8, border = FALSE)
craydata.fam <- read.csv("col_data.csv", row.names = 1)
craydata.fam$species <- rownames(craydata.fam)
mlt <- melt(rel.tab, id.vars = c("species", "combined_burrow_category"), measure.vars = names(rel.tab[3:12]), direction = "long", value.name = "color", variable.name = "raterID")
mlt <- mlt[1:nrow(rel.tab), ]
rownames(mlt) <- mlt$species
craydata <- merge(mlt, craydata.fam[c(1, 7, 19)], by = "species")
rownames(craydata) <- craydata$species
craytree <- read.tree("stern.2017.ml.new.tree")
unique(craydata$color)> [1] "conspicuous" "cryptic"craydata$color <- as.factor(craydata$color)
levels(craydata$color)> [1] "conspicuous" "cryptic"craydata$family[craydata$family == "Cambaridea"] <- "Cambaridae"
craydata$family[craydata$family == "Parastacidea"] <- "Parastacidae"
craydata$family[craydata$family == "Astacidea"] <- "Astacidae"
craydata$family[craydata$family == "Cambaroididea"] <- "Cambaroididae"
craydata$family <- as.factor(craydata$family)
levels(craydata$family)> [1] "Astacidae" "Cambaridae" "Cambaroididae" "Parastacidae"## Pruning the data and tree
chk <- name.check(craytree, craydata)
summary(chk)> 67 taxa are present in the tree but not the data:
> Astacus_leptodactylus,
> Austropotamobius_italicus_carinthiacus,
> Austropotamobius_italicus_carsicus,
> Austropotamobius_italicus_italicus,
> Austropotamobius_torrentium,
> Cambarellus_blacki,
> ....
> 2 taxa are present in the data but not the tree:
> Austropotamobius_torrentium_torrentium,
> Lacunicambarus_miltus
>
> To see complete list of mis-matched taxa, print object.pruned.tree <- drop.tip(craytree, chk$tree_not_data)
pruned.data <- craydata[!(rownames(craydata) %in% chk$data_not_tree),, drop = FALSE]
name.check(pruned.tree, pruned.data)> [1] "OK"## Fan tree
pruned.data$body_color[pruned.data$body_color == "no"] <- "c"
pruned.data$body_color[pruned.data$body_color == "blue"] <- "b"
pruned.data$body_color[pruned.data$body_color == "orange"] <- "o"
pruned.data$body_color[pruned.data$body_color == "purple"] <- "p"
pruned.data$body_color[pruned.data$body_color == "red"] <- "r"
pruned.data$body_color[pruned.data$body_color == "blue_red"] <- "g"
pruned.data$body_color[pruned.data$body_color == "orange_red"] <- "g"
pruned.data$body_color[pruned.data$body_color == "blue_orange_red"] <- "g"
pruned.data$body_color[pruned.data$body_color == "blue_purple"] <- "g"
pruned.data$body_color[pruned.data$body_color == "unpigmented"] <- "c"
b <- pruned.data$species[pruned.data$body_color == "b"]
c <- pruned.data$species[pruned.data$body_color == "c"]
g <- pruned.data$species[pruned.data$body_color == "g"]
o <- pruned.data$species[pruned.data$body_color == "o"]
p <- pruned.data$species[pruned.data$body_color == "p"]
r <- pruned.data$species[pruned.data$body_color == "r"]
bodycol <- setNames(as.factor(pruned.data$body_color), rownames(pruned.data))
head(bodycol)> Astacoides_betsileoensis Astacoides_caldwelli Astacoides_crosnieri
> r c c
> Astacoides_madagascariensis Astacopsis_franklinii Astacopsis_gouldi
> c c c
> Levels: b c g o p rft <- paintBranches(pruned.tree, edge = sapply(b, match, pruned.tree$tip.label),
state = "b", anc.state = "c")
ft <- paintBranches(ft, edge = sapply(g, match, pruned.tree$tip.label),
state = "g")
ft <- paintBranches(ft, edge = sapply(o, match, pruned.tree$tip.label),
state = "o")
ft <- paintBranches(ft, edge = sapply(p, match, pruned.tree$tip.label),
state = "p")
ft <- paintBranches(ft, edge = sapply(r, match, pruned.tree$tip.label),
state = "r")
ft>
> Phylogenetic tree with 399 tips and 398 internal nodes.
>
> Tip labels:
> Samastacus_spinifrons, Virilastacus_araucanius, Virilastacus_retamali, Virilastacus_rucapihuelensis, Parastacus_nicoleti, Parastacus_pugnax, ...
>
> The tree includes a mapped, 6-state discrete character
> with states:
> b, c, g, o, p, r
>
> Rooted; includes branch lengths.colsf <- setNames(c("skyblue", "#956D45", "green4", "orange", "#6600CC", "red"), c("b", "c", "g", "o", "p", "r"))
plot(ft, colsf, type = "fan", ftype = "off", lwd = 2,
mar = c(2.5, 2, 2, 1), part = 0.98)
h <- max(nodeHeights(ft))
tick.spacing <- 20
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(h, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
points(pp$xx[1:Ntip(ft)], pp$yy[1:Ntip(ft)], pch = 16,
col = colsf[bodycol[ft$tip.label]], cex = 0.7)
text(scale, rep(-23, length(scale)), h - scale, cex = 0.6)
text(mean(scale), -38, "time (mya)", cex = 0.7)
legend("bottomleft", legend = c("blue", "cryptic", "polymorphic", "orange", "purple", "red"), pch = 21, pt.bg = colsf, pt.cex = 0.9, bty = "n", col = "transparent", cex = 0.8)
unique(pruned.data$family)> [1] Parastacidae Astacidae Cambaridae Cambaroididae
> Levels: Astacidae Cambaridae Cambaroididae ParastacidaeCambaridae <- pruned.data$species[pruned.data$family == "Cambaridae"]
Parastacidae <- pruned.data$species[pruned.data$family == "Parastacidae"]
Astacidae <- pruned.data$species[pruned.data$family == "Astacidae"]
Cambaroididae <- pruned.data$species[pruned.data$family == "Cambaroididae"]
node1 <- getMRCA(ft, Cambaridae)
node2 <- getMRCA(ft, Parastacidae)
node3 <- getMRCA(ft, Astacidae)
node4 <- getMRCA(ft, Cambaroididae)
nodes <- c(node1, node2, node3, node4)
labels <- c("Cambaridae", "Parastacidae", "Astacidae", "Cambaroididae")
for(i in 1:length(nodes))
arc.cladelabels(text = labels[i], node = nodes[i], ln.offset = 1.05, lab.offset = 1.1, mark.node = FALSE, lwd = 3, orientation = if(labels[i] %in% c("Astacidae", "Cambaroididae")) "horizontal" else "curved")
colsf2 <- setNames(c("skyblue", "black", "green4", "orange", "#6600CC", "red"), c("b", "c", "g", "o", "p", "r"))
plot(ft, colsf2, type = "fan", ftype = "off", lwd = 2,
mar = c(2.5, 2, 2, 1), part = 0.98)
h <- max(nodeHeights(ft))
tick.spacing <- 20
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(h, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
points(pp$xx[1:Ntip(ft)], pp$yy[1:Ntip(ft)], pch = 16,
col = colsf2[bodycol[ft$tip.label]], cex = 0.7)
text(scale, rep(-23, length(scale)), h - scale, cex = 0.6)
text(mean(scale), -38, "time (mya)", cex = 0.7)
legend("bottomleft", legend = c("blue", "cryptic", "polymorphic", "orange", "purple", "red"), pch = 21, pt.bg = colsf, pt.cex = 0.9, bty = "n", col = "transparent", cex = 0.8)
for(i in 1:length(nodes))
arc.cladelabels(text = labels[i], node = nodes[i], ln.offset = 1.05, lab.offset = 1.1, mark.node = FALSE, lwd = 3, orientation = if(labels[i] %in% c("Astacidae", "Cambaroididae")) "horizontal" else "curved")
## Adding different colors to polymorphic species
craydata.fam <- read.csv("col_data.csv", row.names = 1)
craydata.fam$species <- rownames(craydata.fam)
mlt <- melt(rel.tab, id.vars = c("species", "combined_burrow_category"), measure.vars = names(rel.tab[3:12]), direction = "long", value.name = "color", variable.name = "raterID")
mlt <- mlt[1:nrow(rel.tab), ]
rownames(mlt) <- mlt$species
craydata <- merge(mlt, craydata.fam[c(1, 7, 19)], by = "species")
rownames(craydata) <- craydata$species
craytree <- read.tree("stern.2017.ml.new.tree")
unique(craydata$color)> [1] "conspicuous" "cryptic"craydata$color <- as.factor(craydata$color)
levels(craydata$color)> [1] "conspicuous" "cryptic"craydata$family[craydata$family == "Cambaridea"] <- "Cambaridae"
craydata$family[craydata$family == "Parastacidea"] <- "Parastacidae"
craydata$family[craydata$family == "Astacidea"] <- "Astacidae"
craydata$family[craydata$family == "Cambaroididea"] <- "Cambaroididae"
craydata$family <- as.factor(craydata$family)
levels(craydata$family)> [1] "Astacidae" "Cambaridae" "Cambaroididae" "Parastacidae"poly <- craydata
## Changing levels of the factors "body_color"
poly$body_color[poly$body_color == "no"] <- "c"
poly$body_color[poly$body_color == "blue"] <- "b"
poly$body_color[poly$body_color == "orange"] <- "o"
poly$body_color[poly$body_color == "purple"] <- "p"
poly$body_color[poly$body_color == "red"] <- "r"
poly$body_color[poly$body_color == "blue_red"] <- "b+r"
poly$body_color[poly$body_color == "orange_red"] <- "o+r"
poly$body_color[poly$body_color == "blue_orange_red"] <- "b+o+r"
poly$body_color[poly$body_color == "blue_purple"] <- "b+p"
poly$body_color[poly$body_color == "unpigmented"] <- "c"
unique(poly$body_color)> [1] "r" "c" "o" "b+o+r" "b" "o+r" "b+p" "p" "b+r"poly$body_color <- as.factor(poly$body_color)
levels(poly$body_color)> [1] "b" "b+o+r" "b+p" "b+r" "c" "o" "o+r" "p" "r"## Pruning the data and tree
chk <- name.check(craytree, poly)
summary(chk)> 67 taxa are present in the tree but not the data:
> Astacus_leptodactylus,
> Austropotamobius_italicus_carinthiacus,
> Austropotamobius_italicus_carsicus,
> Austropotamobius_italicus_italicus,
> Austropotamobius_torrentium,
> Cambarellus_blacki,
> ....
> 2 taxa are present in the data but not the tree:
> Austropotamobius_torrentium_torrentium,
> Lacunicambarus_miltus
>
> To see complete list of mis-matched taxa, print object.poly.tree <- drop.tip(craytree, chk$tree_not_data)
poly.data <- poly[!(rownames(poly) %in% chk$data_not_tree),, drop = FALSE]
name.check(poly.tree, poly.data)> [1] "OK"polycol <- setNames(as.factor(poly.data$body_color), rownames(poly.data))
head(polycol)> Astacoides_betsileoensis Astacoides_caldwelli Astacoides_crosnieri
> r c c
> Astacoides_madagascariensis Astacopsis_franklinii Astacopsis_gouldi
> c c c
> Levels: b b+o+r b+p b+r c o o+r p rbr <- poly$species[poly$body_color == "b+r"]
or <- poly$species[poly$body_color == "o+r"]
bor <- poly$species[poly$body_color == "b+o+r"]
bp <- poly$species[poly$body_color == "b+p"]
pies <- matrix(0, Ntip(poly.tree), 5,
dimnames = list(poly.tree$tip.label, c("b", "c", "o", "p", "r")))
pies[b, "b"] <- 1
pies[c, "c"] <- 1
pies[br, c("b", "r")] <- 0.5
pies[or, c("o", "r")] <- 0.5
pies[bor, c("b", "o", "r")] <- 1/3
pies[bp, c("b", "p")] <- 0.5
pies[o, "o"] <- 1
pies[p, "p"] <- 1
pies[r, "r"] <- 1colsp <- setNames(c("skyblue", "black", "green4", "orange", "#6600CC", "red"), c("b", "c", "g", "o", "p", "r"))
plot(ft, colsp, type = "fan", ftype = "off", lwd = 2,
mar = c(2.5, 2, 2, 1), part = 0.98)
h <- max(nodeHeights(ft))
tick.spacing <- 20
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(h, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
text(scale, rep(-23, length(scale)), h - scale, cex = 0.6)
text(mean(scale), -38, "time (mya)", cex = 0.7)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
par(fg = "transparent")
colspie <- setNames(c("skyblue", "black", "orange", "#6600CC", "red"), c("b", "c", "o", "p", "r"))
tiplabels(pie = pies, piecol = colspie, cex = 0.22)
par(fg = "black")
legend("bottomleft", legend = c("blue", "cryptic", "orange", "purple", "red"), pch = 21, pt.bg = colspie, pt.cex = 0.9, bty = "n", col = "transparent", cex = 0.8)
for(i in 1:length(nodes))
arc.cladelabels(text = labels[i], node = nodes[i], ln.offset = 1.05, lab.offset = 1.1, mark.node = FALSE, lwd = 3, orientation = if(labels[i] %in% c("Astacidae", "Cambaroididae")) "horizontal" else "curved")
spp <- c("Euastacus_australasiensis", "Engaeus_cymus", "Cambarus_dubius", "Cambarus_harti", "Cherax_robustus",
"Distocambarus_carlsoni", "Faxonius_rusticus", "Procambarus_clarkii")
nodes2 <- sapply(spp, grep, x = ft$tip.label)
nodes2
labels2 <- LETTERS[1:length(spp)]
Australasiensis <- readPNG("imgs/Euastacus_australasiensis.png")
Engaeus <- readPNG("imgs/Engaeus_cymus.png")
Dubius <- readPNG("imgs/Cambarus_dubius.png")
Cambarus <- readPNG("imgs/Cambarus_harti.png")
Cherax <- readPNG("imgs/Cherax_robustus.png")
Distocambarus <- readPNG("imgs/Distocambarus_carlsoni.png")
Faxonius <- readPNG("imgs/Faxonius_rusticus.png")
Procambarus <- readPNG("imgs/Procambarus_clarkii.png")
plot(ft, plot = FALSE)
rasterImage(Australasiensis, -25, 40, 120, 100)
rasterImage(Engaeus, 120, 40, 245, 100)
rasterImage(Dubius, 245, 40, 370, 100)
rasterImage(Cambarus, 370, 40, 495, 100)
rasterImage(Cherax, -25, -18, 120, 40)
rasterImage(Distocambarus, 120, -18, 245, 40)
rasterImage(Faxonius, 245, -18, 370, 40)
rasterImage(Procambarus, 370, -18, 495, 40)
points(-6, 90, bg = "red", pch = 21, cex = 3)
text(-6, 90, "A")
points(135, 90, bg = "orange", pch = 21, cex = 3)
text(135, 90, "B")
points(260, 90, bg = "skyblue", pch = 21, cex = 3)
text(260, 90, "C")
points(385, 90, bg = "skyblue", pch = 21, cex = 3)
text(385, 90, "D")
points(-6, 30, bg = "purple", pch = 21, cex = 3)
text(-6, 30, "E")
points(135, 30, bg = "orange", pch = 21, cex = 3)
text(135, 30, "F")
points(260, 30, bg = "#956D45", pch = 21, cex = 3)
text(260, 30, "G")
points(385, 30, bg = "red", pch = 21, cex = 3)
text(385, 30, "H")
legend("topleft", legend = c("blue", "cryptic", "orange", "purple", "red"), pch = 21, pt.bg = colspie, pt.cex = 0.7, bty = "n", col = "transparent", cex = 0.7)
par(new = TRUE)
par(fg = "transparent")
plot(ft, "white", type = "fan", ftype = "off", mar = c(10, 4, 1, 4), part = 0.96)
par(fg = "black")
for(i in 1:length(nodes))
arc.cladelabels(text = labels[i], node = nodes[i], ln.offset = 1.065, lab.offset = 1.1, mark.node = FALSE, lwd = 1.5, col = "black", cex = 0.7, orientation = if(labels[i] %in% c("Astacidae", "Cambaroididae")) "horizontal" else "curved")
nodelabels(LETTERS[1:length(nodes2)], node = nodes2, frame = "circle", bg = c("red", "orange", "skyblue", "skyblue", "purple", "orange", "#956D45", "red"), cex = 1)
par(new = TRUE)
plot(ft, colsp, type = "fan", ftype = "off", mar = c(12, 6, 3, 6), part = 0.96)
h <- max(nodeHeights(ft))
tick.spacing <- 20
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(h, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
text(scale, rep(-28, length(scale)), h - scale, cex = 0.5, srt = 20)
text(mean(scale), -48, "time (mya)", cex = 0.5)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
par(fg = "transparent")
colspie <- setNames(c("skyblue", "black", "orange", "#6600CC", "red"), c("b", "c", "o", "p", "r"))
tiplabels(pie = pies, piecol = colspie, cex = 0.2)
par(new = TRUE)
par(fg = "transparent")
plot(ft, "white", type = "fan", ftype = "off", mar = c(11.8, 5.8, 2.8, 5.8), part = 0.96)
par(fg = "black")
for(i in 1:length(nodes2)){
xy <- add.arrow(ft, nodes2[i], col = "blue", arrl = 20, lwd = 1.5,
hedl = 5, offset = 2)
}
plot(ft, plot = FALSE)
rasterImage(Australasiensis, -25, 40, 120, 100)
rasterImage(Engaeus, 120, 40, 245, 100)
rasterImage(Dubius, 245, 40, 370, 100)
rasterImage(Cambarus, 370, 40, 495, 100)
rasterImage(Cherax, -25, -18, 120, 40)
rasterImage(Distocambarus, 120, -18, 245, 40)
rasterImage(Faxonius, 245, -18, 370, 40)
rasterImage(Procambarus, 370, -18, 495, 40)
points(-6, 90, bg = "red", pch = 21, cex = 3)
text(-6, 90, "A")
points(135, 90, bg = "orange", pch = 21, cex = 3)
text(135, 90, "B")
points(260, 90, bg = "skyblue", pch = 21, cex = 3)
text(260, 90, "C")
points(385, 90, bg = "skyblue", pch = 21, cex = 3)
text(385, 90, "D")
points(-6, 30, bg = "purple", pch = 21, cex = 3)
text(-6, 30, "E")
points(135, 30, bg = "orange", pch = 21, cex = 3)
text(135, 30, "F")
points(260, 30, bg = "#956D45", pch = 21, cex = 3)
text(260, 30, "G")
points(385, 30, bg = "red", pch = 21, cex = 3)
text(385, 30, "H")
legend("topleft", legend = c("blue", "cryptic", "orange", "purple", "red"), pch = 21, pt.bg = colspie,
pt.cex = 0.7, bty = "n", col = "transparent", cex = 0.7)
par(new = TRUE)
par(fg = "transparent")
plot(ft, "white", type = "fan", ftype = "off", mar = c(10, 4, 1, 4), part = 0.96)
par(fg = "black")
nodelabels(LETTERS[1:length(nodes2)], node = nodes2, frame = "circle", bg = c("red", "orange", "skyblue", "skyblue", "purple", "orange", "#956D45", "red"), cex = 1)
obj <- get("last_plot.phylo", envir = .PlotPhyloEnv)
par(new = TRUE)
plot(ft, colsp, type = "fan", ftype = "off", mar = c(11, 5, 2, 5), part = 0.96)
h <- max(nodeHeights(ft))
tick.spacing <- 20
min.tick <- 0
scale <- axis(1, pos = -5, at = seq(h, min.tick, by = -tick.spacing), cex.axis = 0.5, labels = FALSE)
pp <- get("last_plot.phylo", envir = .PlotPhyloEnv)
text(scale, rep(-28, length(scale)), h - scale, cex = 0.5, srt = 20)
text(mean(scale), -48, "time (mya)", cex = 0.5)
for(i in 1:length(scale)){
a1 <- 0
a2 <- 2*pi
draw.arc(0, 0, radius = scale[i], a1, a2, lwd = 1,
col = make.transparent("blue", 0.15))
}
par(fg = "transparent")
colspie <- setNames(c("skyblue", "black", "orange", "#6600CC", "red"), c("b", "c", "o", "p", "r"))
tiplabels(pie = pies, piecol = colspie, cex = 0.2)
par(fg = "black")