phyloregion
is available from the Comprehensive R Archive Network, so you can use the following line of code to install and run it:
install.packages("phyloregion")
Alternatively, you can install the development version of phyloregion
hosted on GitHub. To do this, you will need to install the devtools
package. In R, type:
if (!requireNamespace("remotes", quietly = TRUE))
install.packages("remotes")
remotes::install_github("darunabas/phyloregion")
When installed, load the package in R:
library(phyloregion)
phyloregion
The workflow of the phyloregion
package demonstrates steps from preparation of different types of data to visualizing the results of biogeographical regionalization, together with tips on selecting the optimal method for achieving the best output, depending on the types of data used and research questions.
In R, phylogenetic relationships among species / taxa are often represented as a phylo object implemented in the ape
package(Paradis and Schliep 2018). Phylogenies (often in the Newick or Nexus formats) can be imported into R with the read.tree
or read.nexus
functions of the ape
package(Paradis and Schliep 2018).
library(ape)
library(Matrix)
library(sp)
data(africa)
<- africa$comm
sparse_comm
<- africa$phylo
tree <- keep.tip(tree, intersect(tree$tip.label, colnames(sparse_comm)))
tree par(mar=c(2,2,2,2))
plot(tree, show.tip.label=FALSE)
The phyloregion
package has functions for manipulating three kinds of distribution data: point records, polygons and raster layers. An overview can be easily obtained with the functions points2comm
, polys2comm
and raster2comm
for point records, polygons, or raster layers, respectively. Depending on the data source, all three functions ultimately provide convenient interfaces to convert the distribution data to a community matrix at varying spatial grains and extents for downstream analyses.
We will play around with these functions in turn.
points2comm
Here, we will generate random points in geographic space, similar to data obtained from museum records, GBIF, iDigBio, or CIESIN which typically have columns of geographic coordinates for each observation.
<- readRDS(system.file("nigeria/nigeria.rds", package = "phyloregion"))
s
set.seed(1)
<- data.frame(sp::spsample(s, 10000, type = "nonaligned"))
m names(m) <- c("lon", "lat")
<- paste0("sp", sample(1:1000))
species $taxon <- sample(species, size = nrow(m), replace = TRUE)
m
<- points2comm(dat = m, mask = s, res = 0.5, lon = "lon", lat = "lat",
pt species = "taxon")
head(pt[[1]][1:5, 1:5])
## 5 x 5 sparse Matrix of class "dgCMatrix"
## sp1 sp10 sp100 sp1000 sp101
## v100 . . . . 1
## v101 . . . . .
## v102 . . . . .
## v103 . 1 . . .
## v104 . . . . .
polys2comm
This function converts polygons to a community matrix at varying spatial grains and extents for downstream analyses. Polygons can be derived from the IUCN Redlist spatial database (https: //www.iucnredlist.org/resources/spatial-data-download), published monographs or field guides validated by taxonomic experts. To illustrate this function, we will use the function random_species
to generate random polygons for five random species over the landscape of Nigeria as follows:
<- readRDS(system.file("nigeria/nigeria.rds", package="phyloregion"))
s <- random_species(100, species=5, shp=s)
sp <- polys2comm(dat = sp, species = "species", trace=0) pol
##
|
| | 0%
|
|======== | 20%
|
|================ | 40%
|
|======================== | 60%
|
|================================ | 80%
|
|========================================| 100%
head(pol[[1]][1:5, 1:5])
## 5 x 5 sparse Matrix of class "dgCMatrix"
## species1 species2 species3 species4 species5
## v27905 . . 1 1 .
## v27906 1 . 1 1 .
## v27907 1 . 1 1 .
## v27908 1 . 1 1 .
## v27909 1 . 1 1 .
raster2comm
This third function, converts raster layers (often derived from species distribution modeling, such as aquamaps(Kaschner et al. 2008)) to a community matrix.
<- system.file("NGAplants", package="phyloregion")
fdir <- file.path(fdir, dir(fdir))
files <- raster2comm(files) ras
##
|
| | 0%
|
|===== | 12%
|
|========== | 25%
|
|=============== | 38%
|
|==================== | 50%
|
|========================= | 62%
|
|============================== | 75%
|
|=================================== | 88%
|
|========================================| 100%
head(ras[[1]])
## 6 x 8 sparse Matrix of class "dgCMatrix"
## Chytranthus_gilletii Commelina_ramulosa Cymbopogon_caesius
## v1 . . .
## v101 . 1 .
## v102 . . .
## v103 . . .
## v104 . . .
## v105 . . .
## Dalechampia_ipomoeifolia Grewia_barombiensis Indigofera_berhautiana
## v1 . . .
## v101 . . 1
## v102 . . 1
## v103 . . 1
## v104 . . 1
## v105 1 . 1
## Indigofera_erecta Indigofera_oubanguiensis
## v1 1 .
## v101 . .
## v102 . .
## v103 . .
## v104 . .
## v105 . .
The object ras
above also returns two objects: a community data frame and a shapefile of grid cells with the numbers of species per cell and can be plotted as a heatmap using our nice plot_swatch
function as follows:
<- readRDS(system.file("nigeria/SR_Naija.rds", package = "phyloregion"))
s par(mar=rep(0,4))
plot_swatch(s, values = s$SR, leg=1, border=NA,
col = hcl.colors(n=20, palette = "Blue-Red 3", rev=FALSE))
## Warning in wkt(obj): CRS object has no comment
Community data are commonly stored in a matrix with the sites as rows and species / operational taxonomic units (OTUs) as columns. The elements of the matrix are numeric values indicating the abundance/observations or presence/absence (0/1) of OTUs in different sites. In practice, such a matrix can contain many zero values because species are known to generally have unimodal distributions along environmental gradients (Ter Braak and Prentice 2004), and storing and analyzing every single element of that matrix can be computationally challenging and expensive.
phyloregion
differs from other R packages (e.g. vegan (Oksanen et al. 2019), picante (Kembel et al. 2010) or betapart(Baselga and Orme 2012)) in that the data are not stored in a (dense) matrix
or data.frame
but as a sparse matrix making use of the infrastructure provided by the Matrix package (Bates and Maechler 2019). A sparse matrix is a matrix with a high proportion of zero entries(Duff 1977), of which only the non-zero entries are stored and used for downstream analysis.
A sparse matrix representation has two advantages. First the community matrix can be stored in a much memory efficient manner, allowing analysis of larger datasets. Second, for very large datasets spanning thousands of taxa and spatial scales, computations with a sparse matrix are often much faster.
The phyloregion
package contains functions to conveniently change between data formats.
library(Matrix)
data(africa)
<- africa$comm
sparse_comm <- as.matrix(sparse_comm)
dense_comm object.size(dense_comm)
## 4216952 bytes
object.size(sparse_comm)
## 885952 bytes
Here, the data set in the dense matrix representation consumes roughly five times more memory than the sparse representation.
We demonstrate the utility of phyloregion
in mapping standard conservation metrics of species richness, weighted endemism (weighted_endemism
) and threat (map_traits
) as well as fast computations of phylodiversity measures such as phylogenetic diversity (PD
), phylogenetic endemism (phylo_endemism
), and evolutionary distinctiveness and global endangerment (EDGE
). The major advantage of these functions compared to available tools e.g. biodiverse (Laffan, Lubarsky, and Rosauer 2010), is the ability to utilize sparse matrix that speeds up the analyses without exhausting computer memories, making it ideal for handling any data from small local scales to large regional and global scales.
weighted_endemism
Weighted endemism is species richness inversely weighted by species ranges(Crisp et al. 2001),(Laffan and Crisp 2003),(Barnabas H. Daru et al. 2020).
library(raster)
##
## Attaching package: 'raster'
## The following objects are masked from 'package:ape':
##
## rotate, zoom
data(africa)
<- weighted_endemism(africa$comm)
Endm head(Endm)
## v3635 v3636 v3637 v3638 v3639 v3640
## 1.770041 2.637894 1.825862 1.270093 1.043782 0.259324
<- merge(africa$polys, data.frame(grids=names(Endm), WE=Endm), by="grids")
m <- m[!is.na(m@data$WE),]
m
par(mar=rep(0,4))
plot_swatch(m, values = m$WE, leg = 3, border = NA,
col = hcl.colors(n=20, palette = "Blue-Red 3", rev=FALSE))
## Warning in wkt(obj): CRS object has no comment
PD
– phylogenetic diversityPhylogenetic diversity (PD
) represents the length of evolutionary pathways that connects a given set of taxa on a rooted phylogenetic tree (Faith 1992). This metric is often characterised in units of time (millions of years, for dated phylogenies). We will map PD for plants of southern Africa.
data(africa)
<- africa$comm
comm <- africa$phylo
tree <- africa$polys
poly
<- PD(comm, tree)
mypd head(mypd)
## v3635 v3636 v3637 v3638 v3639 v3640
## 4226.216 5372.009 4377.735 3783.992 3260.111 1032.685
<- merge(poly, data.frame(grids=names(mypd), pd=mypd), by="grids")
M <- M[!is.na(M@data$pd),]
M head(M)
## grids pd
## 1 v3635 4226.216
## 2 v3636 5372.009
## 3 v3637 4377.735
## 4 v3638 3783.992
## 5 v3639 3260.111
## 6 v3640 1032.685
par(mar=rep(0,4))
plot_swatch(M, values = M$pd, border=NA, leg=3,
col = hcl.colors(n=20, palette = "Blue-Red 3", rev=FALSE))
## Warning in wkt(obj): CRS object has no comment
phylo_endemism
– phylogenetic endemismPhylogenetic endemism is not influenced by variations in taxonomic opinion because it measures endemism based on the relatedness of species before weighting it by their range sizes(Rosauer et al. 2009),(Barnabas H. Daru et al. 2020).
library(raster)
data(africa)
<- africa$comm
comm <- africa$phylo
tree <- africa$polys
poly
<- phylo_endemism(comm, tree)
pe head(pe)
## v3635 v3636 v3637 v3638 v3639 v3640
## 32.536530 45.262625 35.004944 27.603721 23.183947 6.439589
<- merge(poly, data.frame(grids=names(pe), pe=pe), by="grids")
mx <- mx[!is.na(mx@data$pe),]
mx head(mx)
## grids pe
## 1 v3635 32.536530
## 2 v3636 45.262625
## 3 v3637 35.004944
## 4 v3638 27.603721
## 5 v3639 23.183947
## 6 v3640 6.439589
par(mar=rep(0,4))
plot_swatch(mx, values = mx$pe, border=NA, leg=3,
col = hcl.colors(n=20, palette = "Blue-Red 3", rev=FALSE))
## Warning in wkt(obj): CRS object has no comment
EDGE
– Evolutionary Distinctiveness and Global EndangermentThis function calculates EDGE by combining evolutionary distinctiveness (ED; i.e., phylogenetic isolation of a species) with global endangerment (GE) status as defined by the International Union for Conservation of Nature (IUCN).
data(africa)
<- africa$comm
comm <- africa$IUCN
threat <- africa$phylo
tree <- africa$polys
poly
<- EDGE(threat, tree, Redlist = "IUCN", species="Species")
x head(x)
## Abutilon_angulatum_OM1934 Abutilon_sonneratianum_LTM034
## 2.903551 2.903551
## Acalypha_glabrata_glabrata_OM441 Acalypha_glabrata_pilosa_OM1979
## 2.480505 2.480505
## Acalypha_sonderiana_OM2163 Acokanthera_oblongifolia_OM2240
## 2.914481 2.211561
<- map_trait(comm, x, FUN = sd, shp=poly)
y
par(mar=rep(0,4))
plot_swatch(y, y$traits, border=NA, leg=3,
col = hcl.colors(n=20, palette = "Blue-Red 3", rev=FALSE))
## Warning in wkt(obj): CRS object has no comment
The three commonly used methods for quantifying -diversity, the variation in species composition among sites, – Simpson, Sorenson and Jaccard(Laffan et al. 2016). The phyloregion
’s functions beta_diss
and phylobeta
compute efficiently pairwise dissimilarities matrices for large sparse community matrices and phylogenetic trees for taxonomic and phylogenetic turnover, respectively. The results are stored as distance objects for subsequent analyses.
phyloregion
offers a fast means of computing phylogenetic beta diversity, the turnover of branch lengths among sites, making use of and improving on the infrastructure provided by the betapart
package(Baselga and Orme 2012) allowing a sparse community matrix as input.
data(africa)
<- africa$comm
sparse_comm
<- africa$phylo
tree <- keep.tip(tree, intersect(tree$tip.label, colnames(sparse_comm)))
tree <- phylobeta(sparse_comm, tree) pb
<- phyloregion(pb[[1]], shp=africa$polys) y
plot_NMDS(y, cex=3)
text_NMDS(y)
par(mar=rep(0,4))
plot(y, palette="NMDS")
## Warning in wkt(obj): CRS object has no comment
sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] raster_3.4-5 sp_1.4-5 Matrix_1.3-2 ape_5.4-1
## [5] phyloregion_1.0.6
##
## loaded via a namespace (and not attached):
## [1] phangorn_2.5.5 xfun_0.21 bslib_0.2.4
## [4] splines_4.0.4 lattice_0.20-41 clustMixType_0.2-9
## [7] colorspace_2.0-0 vctrs_0.3.6 geometry_0.4.5
## [10] htmltools_0.5.1.1 yaml_2.2.1 mgcv_1.8-34
## [13] utf8_1.1.4 rlang_0.4.10 jquerylib_0.1.3
## [16] pillar_1.5.0 RColorBrewer_1.1-2 foreach_1.5.1
## [19] lifecycle_1.0.0 stringr_1.4.0 rgeos_0.5-5
## [22] codetools_0.2-18 evaluate_0.14 knitr_1.31
## [25] magic_1.5-9 permute_0.9-5 doParallel_1.0.16
## [28] rcdd_1.2-2 parallel_4.0.4 fansi_0.4.2
## [31] highr_0.8 itertools_0.1-3 Rcpp_1.0.6
## [34] vegan_2.5-7 betapart_1.5.2 jsonlite_1.7.2
## [37] abind_1.4-5 picante_1.8.2 fastmatch_1.1-0
## [40] digest_0.6.27 stringi_1.5.3 dismo_1.3-3
## [43] grid_4.0.4 rgdal_1.5-23 quadprog_1.5-8
## [46] tools_4.0.4 magrittr_2.0.1 sass_0.3.1
## [49] tibble_3.1.0 randomForest_4.6-14 cluster_2.1.1
## [52] crayon_1.4.1 pkgconfig_2.0.3 MASS_7.3-53.1
## [55] ellipsis_0.3.1 rmarkdown_2.7 iterators_1.0.13
## [58] R6_2.5.0 igraph_1.2.6 nlme_3.1-152
## [61] compiler_4.0.4