Genus Lygodium
Exploring Lygodium Phylogeny and Ancestral Traits
Chloroplast trnL intron and trnL-F intergenic spacer
John: jmh09r@my.fsu.edu
date()## [1] "Wed Jan 11 12:24:40 2017"Options
library(knitr)
opts_knit$set(verbose = FALSE)Set Directory
setwd("D:/Lygodium/PhyloMicro")Load Library
Search and retrieve from NCBI
library(rentrez)NCBI Pages for Primary Data Source
Population Set
113911638: https://www.ncbi.nlm.nih.gov/popset/113911638 #Madeira
Combine Population and Individual Sequences
ACC = entrez_fetch(db = "nucleotide",
id = c("AF448934", "KT423427",
"AF448932", "KT423393",
"KT423514", "KT423515"),
rettype="fasta")
PopSet = entrez_fetch(db="popset",
id = c("113911638"),
rettype="fasta")
#trnLF = PopSet
trnLF = paste(PopSet, ACC, sep=" ")Save to Text File
write(trnLF, "./trnLF.fasta")Load and Align Sequences
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ape)
library(phangorn)
tf = tempfile()
write(trnLF, tf)
trnLF0 = read.dna(tf, format = "fasta")
trnLF_aligned = muscle(trnLF0, exec = "C:/Windows/muscle")Clean-up Labels and Quick Plot
BaseTree = nj(dist.dna(trnLF_aligned))
BaseTree = root(BaseTree, outgroup = 40)
Lab.df = as.data.frame(BaseTree$tip.label)
names(Lab.df) = "Lab"
Lab.df = Lab.df %>%
mutate(ACC = stringr::word(Lab.df$Lab, 1),
Gen = stringr::word(Lab.df$Lab, 2),
Sp = stringr::word(Lab.df$Lab, 3),
Epi = paste(Gen, Sp, ACC, sep = " "))
uglywords = ".1"
Lab.df$Labels = tm::removeWords(Lab.df$Epi, uglywords)
Lab.df$Labels = gsub("Lygodium", "L.", Lab.df$Epi)
Lab.df$Labels = gsub("Osmundastrum", "O.", Lab.df$Labels)
BaseTree$tip.label = Lab.df$Labels
plot(BaseTree, cex = 0.5)
Pairwise Distance
Unweighted Pair Group Method with Arithmetic Mean hierarchical clustering (UPGMA).
trnLF1 = as.phyDat(trnLF_aligned)
dm = dist.ml(trnLF1)
treeUPGMA = upgma(dm)
parsimony(treeUPGMA, trnLF1)## [1] 535treeUPGMA$tip.label = Lab.df$Labels
plot(treeUPGMA, cex = 0.5)
Improving Parsimony
trnLF1 = as.phyDat(trnLF_aligned)
dm = dist.ml(trnLF1)
treeUPGMA = upgma(dm)
parsimony(treeUPGMA, trnLF1)## [1] 535treePars = optim.parsimony(treeUPGMA, trnLF1)## Final p-score 521 after 8 nni operationstreeRatchet = pratchet(trnLF1, trace = 0)
parsimony(c(treePars, treeRatchet), trnLF1)## [1] 521 521treePars2 = treePars
treePars2$tip.label = Lab.df$Labels
plot(root(treePars2, 40), "fan", cex = 0.5)
Explore Models of Evolution
Nucleotide substitution models
mt = modelTest(trnLF1)## 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"mt2 = mt #copy
arrange(mt2, BIC) #sort by BIC## Model df logLik AIC AICw AICc AICcw
## 1 GTR+G 90 -4194.237 8568.474 7.278480e-01 8583.754 7.609055e-01
## 2 GTR+G+I 91 -4194.237 8570.474 2.677210e-01 8586.108 2.344621e-01
## 3 GTR+I 90 -4199.339 8578.677 4.430748e-03 8593.957 4.631984e-03
## 4 HKY+G 86 -4214.591 8601.183 5.746349e-08 8615.090 1.193381e-07
## 5 HKY+G+I 87 -4214.592 8603.183 2.113650e-08 8617.427 3.709504e-08
## 6 HKY+I 86 -4219.357 8610.714 4.894375e-10 8624.621 1.016446e-09
## 7 GTR 89 -4210.525 8599.050 1.669339e-07 8613.980 2.078650e-07
## 8 HKY 85 -4230.930 8631.859 1.253613e-14 8645.434 3.074062e-14
## 9 K80+G 83 -4256.935 8679.870 4.705766e-25 8692.793 1.598407e-24
## 10 SYM+G 87 -4243.743 8661.485 4.622382e-21 8675.729 8.112386e-21
## 11 K80+G+I 84 -4256.935 8681.871 1.730889e-25 8695.117 5.000815e-25
## 12 SYM+G+I 88 -4243.743 8663.486 1.700254e-21 8678.070 2.516212e-21
## 13 K80+I 83 -4262.531 8691.061 1.748157e-27 8703.984 5.937961e-27
## 14 SYM+I 87 -4249.594 8673.188 1.329330e-23 8687.432 2.333004e-23
## 15 K80 82 -4275.337 8714.674 1.303397e-32 8727.278 5.193795e-32
## 16 SYM 86 -4262.727 8697.454 7.150192e-29 8711.361 1.484926e-28
## 17 F81+G 85 -4326.783 8823.566 2.948388e-56 8837.140 7.229925e-56
## 18 F81+G+I 86 -4326.783 8825.566 1.084476e-56 8839.473 2.252201e-56
## 19 F81+I 85 -4331.058 8832.115 4.102560e-58 8845.690 1.006014e-57
## 20 F81 84 -4341.295 8850.591 3.992273e-62 8863.837 1.153431e-61
## 21 JC+G 82 -4362.695 8889.390 1.499694e-70 8901.994 5.976001e-70
## 22 JC+G+I 83 -4362.695 8891.390 5.516181e-71 8904.314 1.873680e-70
## 23 JC+I 82 -4367.289 8898.578 1.516582e-72 8911.182 6.043297e-72
## 24 JC 81 -4378.576 8919.153 5.166214e-77 8931.441 2.409903e-76
## BIC
## 1 9023.762
## 2 9030.821
## 3 9033.965
## 4 9036.236
## 5 9043.295
## 6 9045.767
## 7 9049.280
## 8 9061.853
## 9 9099.747
## 10 9101.597
## 11 9106.806
## 12 9108.656
## 13 9110.938
## 14 9113.300
## 15 9129.492
## 16 9132.508
## 17 9253.560
## 18 9260.619
## 19 9262.110
## 20 9275.526
## 21 9304.208
## 22 9311.267
## 23 9313.396
## 24 9328.912(fit = eval(get("GTR+G", env), env)) #select model, save as "fit"##
## loglikelihood: -4194.237
##
## unconstrained loglikelihood: -5570.957
## Discrete gamma model
## Number of rate categories: 4
## Shape parameter: 1.311994
##
## Rate matrix:
## a c g t
## a 0.0000000 1.1146883 2.5598350 0.2884088
## c 1.1146883 0.0000000 0.5322027 3.5141055
## g 2.5598350 0.5322027 0.0000000 1.0000000
## t 0.2884088 3.5141055 1.0000000 0.0000000
##
## Base frequencies:
## 0.3255293 0.1694899 0.2136642 0.2913166Optimze Selected Model
(GTR+G)
Maximum Likelihood Tree
# optimise model parameters without fitting edges
fit1 = optim.pml(fit,
optNni = FALSE, optBf = TRUE,
optQ = TRUE, optInv = TRUE,
optGamma = TRUE, optEdge = FALSE,
optRate = TRUE,
control = pml.control(epsilon = 1e-08,
maxit = 10, trace = 0))
# fix substitution model and fit tree
fit2 = optim.pml(fit1,
optNni = TRUE, optBf = FALSE,
optQ = FALSE, optInv = FALSE,
optGamma = FALSE, optEdge = TRUE,
control = pml.control(epsilon = 1e-08,
maxit = 10, trace = 0))## Warning: I unrooted the tree# fine tune
fit3 = optim.pml(fit2,
optNni = TRUE, optBf = TRUE,
optQ = TRUE, optInv = TRUE,
optGamma = TRUE, optEdge = TRUE,
optRate = FALSE,
control = pml.control(epsilon = 1e-08,
maxit = 10, trace = 0))
fit4 = fit3$tree
fit4$tip.label## [1] "DQ845205.1 Lygodium microphyllum voucher L42A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [2] "DQ845206.1 Lygodium microphyllum voucher L44A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [3] "DQ845209.1 Lygodium microphyllum voucher L26A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [4] "DQ845210.1 Lygodium microphyllum voucher L19C tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [5] "DQ845211.1 Lygodium microphyllum voucher L16A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [6] "DQ845212.1 Lygodium microphyllum voucher L10A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [7] "DQ845213.1 Lygodium microphyllum voucher L13A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [8] "DQ845214.1 Lygodium articulatum voucher L97A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [9] "DQ845215.1 Lygodium articulatum voucher L98A tRNA-Leu (trnL) gene and intergenic spacer, partial sequence; chloroplast"
## [10] "DQ845216.1 Lygodium reticulatum voucher L141A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [11] "DQ845217.1 Lygodium reticulatum voucher L140A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [12] "DQ845219.1 Lygodium cubense voucher L160B tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [13] "DQ845220.1 Lygodium oligostachyum voucher L95A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [14] "DQ845222.1 Lygodium smithianum voucher L8A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [15] "DQ845223.1 Lygodium volubile voucher L70A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [16] "DQ845225.1 Lygodium lanceolatum voucher L135A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [17] "DQ845226.1 Lygodium venustum voucher L90A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [18] "DQ845228.1 Lygodium circinatum voucher L130A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [19] "DQ845229.1 Lygodium palmatum voucher L150A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [20] "DQ845231.1 Lygodium japonicum voucher L57B tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [21] "DQ845232.1 Lygodium japonicum voucher L51B tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [22] "DQ845235.1 Lygodium flexuosum voucher L180A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [23] "DQ845238.1 Lygodium flexuosum voucher L82AS tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [24] "DQ845239.1 Lygodium flexuosum voucher L82C tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [25] "DQ845240.1 Lygodium polystachyum voucher L81A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [26] "AF448934.1 Lygodium microphyllum chloroplast trnL-trnF intergenic spacer, partial sequence"
## [27] "KT423427.1 Schizaea elegans isolate PL5550 trnL-trnF intergenic spacer region, partial sequence; chloroplast"
## [28] "AF448932.1 Schizaea elegans chloroplast trnL-trnF intergenic spacer, partial sequence"
## [29] "KT423393.1 Schizaea pectinata isolate BR245 trnL-trnF intergenic spacer region, partial sequence; chloroplast"
## [30] "KT423514.1 Lygodium palmatum isolate SHsn trnL-trnF intergenic spacer region, partial sequence; chloroplast"
## [31] "KT423515.1 Lygodium reticulatum isolate SM8578 trnL-trnF intergenic spacer region, partial sequence; chloroplast"
## [32] "DQ845207.1 Lygodium microphyllum voucher L101B tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [33] "DQ845208.1 Lygodium microphyllum voucher L105A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [34] "DQ845218.1 Lygodium reticulatum voucher L143A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [35] "DQ845221.1 Lygodium oligostachyum voucher L96A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [36] "DQ845224.1 Lygodium volubile voucher L71A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [37] "DQ845227.1 Lygodium venustum voucher L91A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [38] "DQ845230.1 Lygodium palmatum voucher L150C tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [39] "DQ845233.1 Lygodium japonicum voucher L56A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [40] "DQ845234.1 Lygodium japonicum voucher L63A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [41] "DQ845236.1 Lygodium flexuosum voucher L22A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"
## [42] "DQ845237.1 Lygodium flexuosum voucher L80A tRNA-Leu (trnL) gene and trnL-trnF intergenic spacer, partial sequence; chloroplast"fit4 = root(fit4, 29, resolve.root = TRUE) # reroots the tree
Lab.df = as.data.frame(fit4$tip.label)
names(Lab.df) = "Lab"
Lab.df = Lab.df %>%
mutate(ACC = stringr::word(Lab.df$Lab, 1),
Gen = stringr::word(Lab.df$Lab, 2),
Sp = stringr::word(Lab.df$Lab, 3),
Epi = paste(Gen, Sp, ACC, sep = " "))
uglywords = ".1"
Lab.df$Labels = tm::removeWords(Lab.df$Epi, uglywords)
Lab.df$Labels = gsub("Lygodium", "L.", Lab.df$Epi)
Lab.df$Labels = gsub("Schizaea", "S.", Lab.df$Labels)
fit4$tip.label = Lab.df$Labels
plot(fit4,
cex = 0.5,
edge.color = "black",
tip.color = "black",
use.edge.length = TRUE,
label.offset = 0.005,
show.node.label = TRUE)
add.scale.bar()
title("Optimized Parsimony \n(GTR+G)")
Bootstrap Optimized Model
set.seed(1976)
bs = bootstrap.pml(fit3,
bs = 10,
optNni = TRUE,
control = pml.control(trace = 0))
Lab.df = as.data.frame(attr(bs, "TipLabel"))
names(Lab.df) = "Lab"
Lab.df = Lab.df %>%
mutate(ACC = stringr::word(Lab.df$Lab, 1),
Gen = stringr::word(Lab.df$Lab, 2),
Sp = stringr::word(Lab.df$Lab, 3),
Epi = paste(Gen, Sp, ACC, sep = " "))
uglywords = ".1"
Lab.df$Labels = tm::removeWords(Lab.df$Epi, uglywords)
Lab.df$Labels = gsub("Lygodium", "L.", Lab.df$Epi)
Lab.df$Labels = gsub("Schizaea", "S.", Lab.df$Labels)
bs2 = bs #Copy
attr(bs2, "TipLabel") = as.character(Lab.df$Labels)Plot Results
Optimized and bootstrapped.
plotBS(midpoint(fit4), bs2,
type="p", cex=0.4,
bs.adj = c(1.25, 1.25),
bs.col = "black")
add.scale.bar()
title("Maximum Parsimony")
Consensus Network
cnet = consensusNet(bs2, p=0.2)
set.seed(123)
plot(cnet, "2D",
show.tip.label = FALSE,
show.edge.label = FALSE,
cex = 1, tip.col = "red")
Temperature Tolerance
Reduced Tree
Selecting a representative member of each Lygodium species and mximizing parsimony. For a fuller description of the below steps, see: https://rpubs.com/JMHumphreys/PhyllotisTree
RedACC = entrez_fetch(db = "nucleotide",
id = c("DQ845205", "DQ845214", "DQ845216", "DQ845219", "DQ845220",
"DQ845222", "DQ845223", "DQ845225", "DQ845226", "DQ845228",
"DQ845229", "DQ845231", "DQ845235", "DQ845240"),
rettype="fasta")
write(RedACC, "./trnLF.fasta")
tf = tempfile()
write(RedACC, tf)
RedACC0 = read.dna(tf, format = "fasta")
RedACC_aligned = muscle(RedACC0, exec = "C:/Windows/muscle")
TrnLF = as.phyDat(RedACC_aligned)
dm = dist.ml(TrnLF)
RedTree = upgma(dm)Update Labels and Plot
Lab.df = as.data.frame(RedTree$tip.label)
names(Lab.df) = "Lab"
Lab.df = Lab.df %>%
mutate(ACC = stringr::word(RedTree$tip.label, 1),
Gen = stringr::word(RedTree$tip.label, 2),
Sp = stringr::word(RedTree$tip.label, 3),
Epi = paste(Gen, Sp, ACC, sep = " "))
uglywords = ".1"
Lab.df$Labels = tm::removeWords(Lab.df$Epi, uglywords)
Lab.df$Labels = gsub("Lygodium", "L.", Lab.df$Epi)
RedTree$tip.label = Lab.df$Labels
plot.new()
plot(RedTree)
Habitat Tree
Constructing a “Habitat Tree” using environmental variables (mostly climatic).
Load Packages:
suppressMessages(library(rgdal))
suppressMessages(library(reshape2))
suppressMessages(library(raster))
suppressMessages(library(fields))
suppressMessages(library(rasterVis))
suppressMessages(library(rgeos))
suppressMessages(library(Thermimage))
suppressMessages(library(maptools))
suppressMessages(library(mapproj))
suppressMessages(library(spdep))
suppressMessages(library(ggplot2))
suppressMessages(library(GISTools))
suppressMessages(library(lattice))
suppressMessages(library(gridExtra))
suppressMessages(library(dismo))
suppressMessages(library(dplyr))Get Point Data
Lyg.df = read.csv("D:/Lygodium/GBIF_data/MDat/REvCoord/Comb122816.csv",
header = TRUE, sep=",")
#Observation Summary
Count1 = as.data.frame(
Lyg.df %>%
group_by(Species) %>%
summarise(no_rows = length(Species)))
Count1## Species no_rows
## 1 Lygodium articulatum 2288
## 2 Lygodium auriculatum 27
## 3 Lygodium circinnatum 692
## 4 Lygodium cubense 18
## 5 Lygodium flexuosum 1088
## 6 Lygodium japonicum 11400
## 7 Lygodium lanceolatum 111
## 8 Lygodium microphyllum 25020
## 9 Lygodium oligostachyum 26
## 10 Lygodium palmatum 120
## 11 Lygodium polystachyum 73
## 12 Lygodium reticulatum 288
## 13 Lygodium smithianum 160
## 14 Lygodium venustum 1757
## 15 Lygodium volubile 908
## 16 Schizaea elegans 641
## 17 Schizaea pectinata 238Get Country Boundaries
World = readOGR(dsn = "D:/Lygodium/GBIF_data/MDat/World",
layer = "World",
stringsAsFactors = FALSE)## OGR data source with driver: ESRI Shapefile
## Source: "D:/Lygodium/GBIF_data/MDat/World", layer: "World"
## with 7 features
## It has 4 fieldsWorld = spTransform(World, "+proj=merc +a=6378137 +b=6378137 +lat_ts=0 +lon_0=0
+x_0=0 +y_0=0 +k=1 +units=km +nadgrids=@null +no_defs")
ext = as.vector(c(-180.00000, 180.00000,
-55.90222, 83.60220))
Countries = map("world", fill = TRUE,
xlim = ext[1:2],
ylim = ext[3:4],
plot = FALSE)
IDs = sapply(strsplit(Countries$names, ":"), function(x) x[1])
LL84 = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
CountriesP = map2SpatialPolygons(Countries, IDs = IDs,
proj4string = CRS(projection(LL84 )))
pid = sapply(slot(CountriesP, "polygons"),
function(x) slot(x, "ID"))
p.df = data.frame( ID=1:length(CountriesP),
row.names = pid)
Countries = SpatialPolygonsDataFrame(CountriesP, p.df)Remove Observations over the Ocean
LL = cbind(Lyg.df$Long, Lyg.df$Lat)
Phy.pnt = SpatialPointsDataFrame(LL, Lyg.df)
proj4string(Phy.pnt) = LL84
N = length(Phy.pnt[,1])
N ## [1] 44855plot(Countries, col="tan", main = "Observations")
plot(Phy.pnt, cex=0.3, pch=19, col=factor(Phy.pnt$Species), add=T)
Population Density
Population Density will be a proxy for sampling bias.
http://sedac.ciesin.columbia.edu/data/set/gpw-v3-population-density
Pop = raster("D:/Lygodium/gl_gpwv3_pdens_00_wrk_25/gldens00/glds00g/w001001.adf")Elevation
SRTM (NASA)
Elev = raster("D:/BulkData_dl/Global_LC/srtmv4_30s/w001001.adf")Land Cover
As well as a few elevational indicies.
IMI = “Integrated Mositure Index”
HLI = “Heat Loading Index”
CTI = “Compound Topographic Index”
These were constructed as part of a Lygodium-related modeling project.
LCover = raster("D:/BulkData_dl/Global_LC/Globcover2009_V2.3_Global_/GLOBCOVER_L4_200901_200912_V2.3.tif")
IMI = raster("D:/BulkData_dl/Global_LC/SRTMindx/imi.tif")
HLI = raster("D:/BulkData_dl/Global_LC/SRTMindx/hli.tif")
CTI = raster("D:/BulkData_dl/Global_LC/SRTMindx/cti.tif")Soils
Soil PH, Oraganic Content, Cation Exchange Capacity, and Aluminum content.
PH = raster("D:/BulkData_dl/Global_LC/SRTMindx/Soils/PH_aq.tif")
OC = raster("D:/BulkData_dl/Global_LC/SRTMindx/Soils/ORGC.tif")
CE = raster("D:/BulkData_dl/Global_LC/SRTMindx/Soils/CACC.tif")
AL = raster("D:/BulkData_dl/Global_LC/SRTMindx/Soils/ALSA.tif")
Soil.stk = stack(PH, OC, CE, AL)
names(Soil.stk) = c("PH", "OC", "CE", "AL")Percent Canopy
pC = raster("D:/BulkData_dl/Global_LC/SRTMindx/pCanopy.tif")Background Points
r2 = raster(res = 0.75,
crs = proj4string(Phy.pnt))
extent(r2) = extent(Countries)
Domain.r = rasterize(Countries, r2,
field = 0,
background = NA)
Grid.pnts = rasterToPoints(Domain.r, sp = TRUE)
Grid.pnts@data = Grid.pnts@data %>%
mutate(Long = Grid.pnts@coords[,1],
Lat = Grid.pnts@coords[,2],
OBS = 0,
Species = "Quad",
Source = "Quad")
Grid.pnts@data = Grid.pnts@data[2:6]Combine Observation and Background Points
Phy.pnt$OBS = 1
PhyPA = spRbind(Phy.pnt, Grid.pnts)
PhyPA@data = PhyPA@data %>%
mutate(Long = as.numeric(PhyPA@coords[,1]),
Lat = as.numeric(PhyPA@coords[,2]))Extract Potential Fixed Effects
setwd("D:/Phyllotis/WClim")
Bio = getData('worldclim',
var='bio',
res=2.5)
setwd("D:/Lygodium")
BioClim = as.data.frame(
extract(Bio, spTransform(PhyPA, CRS(proj4string(Bio))), method = "simple"))
Soil.df = as.data.frame(
extract(Soil.stk, spTransform(PhyPA, CRS(proj4string(Soil.stk))), method = "simple"))
Temp.df = cbind(BioClim, Soil.df) %>%
mutate(Pop = extract(Pop, spTransform(PhyPA, CRS(proj4string(Bio))), method = "simple"),
Elev = extract(Elev, spTransform(PhyPA, CRS(proj4string(Elev))), method = "simple"),
IMI = extract(IMI, spTransform(PhyPA, CRS(proj4string(IMI))), method = "simple"),
HLI = extract(HLI, spTransform(PhyPA, CRS(proj4string(HLI))), method = "simple"),
CTI = extract(CTI, spTransform(PhyPA, CRS(proj4string(CTI))), method = "simple"),
pC = extract(pC, spTransform(PhyPA, CRS(proj4string(pC))), method = "simple"))
Temp.df = as.data.frame(apply(Temp.df, 2, function(x) {
x[which(is.na(x))] = mean(x,na.rm=TRUE);x}))
PhyPA@data = cbind(PhyPA@data, Temp.df)Removing a Non-Lygodium Points
Env.df = PhyPA@data
##Try
nEnv.df = Env.df %>%
filter(Species != "Quad" &
Species != "Schizaea elegans" &
Species != "Schizaea pectinata")
sEnv.df = nEnv.df[6:34]
sEnv.df$Labs = nEnv.df$SpeciesSample Cluster and Re-Code as Factor
Environmental variables are dicretized using clustering algorhythms and then coverted to a matrix and analysed as “Habitat DNA” in phylogentic approach. Some information is lost through the clustering procedure, but, the dataset is too large otherwise.
library(FactoMineR)
by_Spp = sEnv.df %>% group_by(Labs)
set.seed(1976)
DF = as.data.frame(sample_n(by_Spp, 500, replace = TRUE))
DF2 = DF[1:29]
data.cat = DF2
for (i in 1:ncol(DF2)){
vari = data.cat[,i,drop=FALSE]
res.hcpc = HCPC(vari, nb.clust=-1, graph=FALSE)
maxi = unlist(by(res.hcpc$data.clust[,1], res.hcpc$data.clust[,2], max))
breaks = c(min(vari), maxi)
aaQuali = cut(vari[,1], breaks, include.lowest=TRUE)
nL = length(levels(aaQuali))
LevN = as.character(levels(aaQuali))
RepF = LETTERS[1:nL]
xRet = factor(aaQuali, levels = LevN, labels = RepF)
data.cat[,i] = xRet
}
data.cat$Spp = DF$Labs
melt.df = melt(data.cat, "Spp")## Warning: attributes are not identical across measure variables; they will
## be droppedCast.df = dcast(melt.df, Spp + variable ~ value, length)
for (i in 1:15) {
Temp3 = Cast.df
SppName = levels(Temp3$Spp)[i]
Temp3 = Temp3 %>%
filter(Spp == paste(SppName))
Temp3 = Temp3[,(3:8)]
GeneSeq = colnames(Temp3)[max.col(Temp3, ties.method="first")]
if (i == 1) {
Data.df = GeneSeq
}
else {Data.df = rbind(Data.df, GeneSeq) }
}
row.names(Data.df) = levels(DF$Labs)[1:15]
dim(Data.df) #15 Species, each with 29 "Habitat Sequences"## [1] 15 29head(Data.df)## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## Lygodium articulatum "A" "B" "B" "C" "A" "B" "C" "A" "B" "A"
## Lygodium auriculatum "C" "A" "E" "A" "C" "F" "A" "C" "C" "F"
## Lygodium circinnatum "C" "A" "F" "A" "C" "F" "A" "C" "C" "E"
## Lygodium cubense "C" "B" "D" "B" "D" "E" "B" "C" "C" "E"
## Lygodium flexuosum "C" "B" "D" "B" "D" "E" "B" "C" "C" "F"
## Lygodium japonicum "B" "D" "B" "D" "D" "B" "D" "C" "B" "E"
## [,11] [,12] [,13] [,14] [,15] [,16] [,17] [,18] [,19]
## Lygodium articulatum "B" "B" "A" "D" "A" "B" "C" "B" "C"
## Lygodium auriculatum "E" "B" "B" "B" "C" "C" "C" "B" "B"
## Lygodium circinnatum "E" "B" "B" "D" "A" "B" "C" "B" "B"
## Lygodium cubense "D" "A" "A" "A" "C" "A" "A" "B" "A"
## Lygodium flexuosum "E" "A" "B" "A" "D" "B" "A" "A" "A"
## Lygodium japonicum "B" "A" "A" "C" "B" "A" "B" "B" "B"
## [,20] [,21] [,22] [,23] [,24] [,25] [,26] [,27] [,28]
## Lygodium articulatum "B" "A" "C" "C" "A" "B" "B" "B" "A"
## Lygodium auriculatum "B" "A" "B" "C" "A" "A" "B" "A" "B"
## Lygodium circinnatum "B" "A" "A" "B" "A" "A" "B" "A" "B"
## Lygodium cubense "D" "A" "C" "A" "A" "A" "B" "A" "B"
## Lygodium flexuosum "B" "A" "A" "A" "A" "A" "B" "A" "B"
## Lygodium japonicum "B" "A" "A" "A" "A" "A" "B" "A" "B"
## [,29]
## Lygodium articulatum "D"
## Lygodium auriculatum "A"
## Lygodium circinnatum "D"
## Lygodium cubense "B"
## Lygodium flexuosum "A"
## Lygodium japonicum "C"Convert to DNA Matrix and Plot
GeogT.df = Data.df
GeogT = phyDat(as.matrix(GeogT.df),
type="USER",
levels= LETTERS[1:6])
dm = dist.ml(GeogT)
GeogTree = upgma(dm)
parsimony(GeogTree, GeogT)## [1] 127GeogTreeP = optim.parsimony(GeogTree, GeogT)## Final p-score 124 after 1 nni operationsGeogTreeR = pratchet(GeogT, trace = 0)
parsimony(c(GeogTreeP, GeogTreeR), GeogT)## [1] 124 119GeogTreeR$tip.label = gsub("Lygodium", "L.", GeogTreeR$tip.label)
plot(GeogTreeR, "unrooted", main = "Habitat Tree \n(unrooted)", cex = 0.5)
Remove L. auriculatum
Later the “Habitat Tree” is compared to the genetic species tree, which lacks sequences for L. auriculatum.
GeogT2.df = Data.df[-2,]
GeogT2 = phyDat(as.matrix(GeogT2.df),
type="USER",
levels= LETTERS[1:6])
dm = dist.ml(GeogT2)
GeogTree2 = upgma(dm)
parsimony(GeogTree2, GeogT2)## [1] 120GeogTree2P = optim.parsimony(GeogTree2, GeogT2)## Final p-score 114 after 4 nni operationsGeogTreeR2 = pratchet(GeogT2, trace = 0)
parsimony(c(GeogTree2P, GeogTreeR2), GeogT2)## [1] 114 113GeogTreeR2$tip.label = gsub("Lygodium", "L.", GeogTreeR2$tip.label)
GeogTree2$tip.label = gsub("Lygodium", "L.", GeogTree2$tip.label)
GeogTreeR2 = root(GeogTreeR2,
outgroup = "L. articulatum")
plot(GeogTreeR2, "phylogram", cex = 0.5)
Compare Trees and Untangle
Compare Habitat Tree to Species Tree. First, set up the trees ensuring lables correspond and convert to dendrogram. Next, plot the base entanglement.
library(dendextend)##
## ---------------------
## Welcome to dendextend version 1.3.0
## Type browseVignettes(package = 'dendextend') for the package vignette.
## The github page is: https://github.com/talgalili/dendextend/
##
## Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
## Or contact: <tal.galili@gmail.com>
##
## To suppress this message use: suppressPackageStartupMessages(library(dendextend))
## ---------------------##
## Attaching package: 'dendextend'## The following object is masked from 'package:raster':
##
## rotate## The following objects are masked from 'package:ape':
##
## ladderize, rotate## The following object is masked from 'package:stats':
##
## cutreeTanRedTree = RedTree
TanRedTree$tip.label = c("L. microphyllum", "L. articulatum", "L. reticulatum", "L. cubense",
"L. oligostachyum", "L. smithianum", "L. volubile", "L. lanceolatum",
"L. venustum", "L. circinatum", "L. palmatum", "L. japonicum",
"L. flexuosum", "L. polystachyum")
GeogTree2$tip.label = c("L. articulatum", "L. circinatum",
"L. cubense", "L. flexuosum", "L. japonicum", "L. lanceolatum",
"L. microphyllum", "L. oligostachyum", "L. palmatum",
"L. polystachyum", "L. reticulatum", "L. smithianum",
"L. venustum", "L. volubile")
dend = as.dendrogram(TanRedTree)
dend2 = as.dendrogram(GeogTree2)
dend_diff(dend, dend2)
dends = dendlist(dend, dend2)
dends %>% entanglement # lower is better## [1] 0.6093824Reduce Entanglement
A series of methods are applied to locate the arrangment that best simplifies the entanglment. The goal is to change the layout to achieve the lowest entanglemnt score.
dends %>% dendextend::untangle(method = "DendSer") %>% entanglement # lower is better## [1] 0.6093824dends %>% dendextend::untangle(method = "step1side") %>% entanglement # lower is better## [1] 0.2267948dends %>% dendextend::untangle(method = "step1side") %>%
tanglegram(common_subtrees_color_branches = TRUE)
x = dends
x %>% plot(main = paste("entanglement =", round(entanglement(x), 2)))
x = dends %>% dendextend::untangle(method = "ladderize")
x %>% plot(main = paste("entanglement =", round(entanglement(x), 2)))
set.seed(1976)
x = dends %>% dendextend::untangle(method = "random", R = 10)
x %>% plot(main = paste("entanglement =", round(entanglement(x), 2)))
x = dends %>% dendextend::untangle(method = "step2side")
x %>% plot(main = paste("entanglement =", round(entanglement(x), 2)))
Plot Best Entanglegram
x %>% set("rank_branches") %>%
tanglegram(main_left = "Species Tree",
main_right = "Habitat Tree",
columns_width = c(3,2,3),
dLeaf_left = -0.1,
dLeaf_right = 0.1,
margin_outer = 4,
margin_inner = 7,
color_lines = c(1:14),
common_subtrees_color_branches = FALSE)
Ancestral Reconstruction
Comapring the contnuous environmental variables to the reduced species tree.
Temperature Tolerance
The below Minimum Thermal Thresholds were estimated in conjunction with a different distribution modeling project. Here we add the estimated temperature tolerances to a dataset with correponding species.
Env.df = as.data.frame(
aggregate(PhyPA@data[, 6:34], list(PhyPA$Species), median))
Env.df = Env.df[-c(2, 16, 17, 18),]
nEnv.df = Env.df[2:30]
LNames = c("L. articulatum", "L. circinatum", "L. cubense", "L. flexuosum",
"L. japonicum", "L. lanceolatum", "L. microphyllum", "L. oligostachyum", "L. palmatum",
"L. polystachyum", "L. reticulatum", "L. smithianum", "L. venustum", "L. volubile")
###Min Temperature Threshold from SDMs
Tvals = c(-2.25, -1.05, 0.85, -3.15,
-4.65, -0.75, -0.35, 3.55, -6.5,
-0.85, -0.05, -0.95, 0.85, -0.95)
Env.df$Spp = LNames
Env.df$Temp = Tvals
TreeLabs = as.data.frame(RedTree$tip.label)
names(TreeLabs) = "Label"
TreeLabs$Order = 1:nrow(TreeLabs)
TreeLabs = arrange(TreeLabs, Label)
Env.df = cbind(TreeLabs, Env.df)
Env.df = arrange(Env.df, Order)Ancestral Trait (Model Min Temp)
library(phytools)##
## Attaching package: 'phytools'## The following object is masked from 'package:dendextend':
##
## untangle## The following object is masked from 'package:Matrix':
##
## expmTipValue = as.vector(Env.df$Temp)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue,
vars=TRUE, CI=TRUE)
TemAT## Ancestral character estimates using fastAnc:
## 15 16 17 18 19 20 21
## -2.435067 -2.493176 -1.822076 -0.796332 -1.563451 -1.503770 -0.924381
## 22 23 24 25 26 27
## -1.206525 -0.572228 0.141048 -1.605125 0.209045 0.757104
##
## Variances on ancestral states:
## 15 16 17 18 19 20 21 22
## 9.809172 5.069810 3.779362 3.284645 2.994213 2.197042 1.254685 1.079077
## 23 24 25 26 27
## 1.160673 0.803029 1.045921 0.727699 0.642390
##
## Lower & upper 95% CIs:
## lower upper
## 15 -8.573708 3.703574
## 16 -6.906358 1.920007
## 17 -5.632430 1.988278
## 18 -4.348556 2.755892
## 19 -4.954995 1.828093
## 20 -4.408964 1.401425
## 21 -3.119831 1.271068
## 22 -3.242547 0.829496
## 23 -2.683824 1.539369
## 24 -1.615345 1.897441
## 25 -3.609623 0.399372
## 26 -1.462937 1.881028
## 27 -0.813821 2.328028obj = contMap(RedTree, TipValue, plot=FALSE)
dev.off()## null device
## 1plot(obj,
sub = "Min Temp",
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))Annual Mean Temperature
TipValue = as.vector(Env.df$bio1/10)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
obj = setMap(obj, invert=TRUE)
plot(obj,
title = "Mean Annual Temperature",
type="phylogram",
legend= 0.7*max(nodeHeights(RedTree)))
title(main=expression("Annual Minimum Temperature"*~degree*C))
Mean Diurnal Range
(Mean of monthly (max temp - min temp))
TipValue = as.vector(Env.df$bio2)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Isothermality
(BIO2/BIO7) (*100)
TipValue = as.vector(Env.df$bio3)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
emperature Seasonality
(standard deviation *100)
TipValue = as.vector(Env.df$bio4)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Max Temperature of Warmest Month
TipValue = as.vector(Env.df$bio5)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Min Temperature of Coldest Month
TipValue = as.vector(Env.df$bio6/10)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Temperature Annual Range
(BIO5-BIO6)
TipValue = as.vector(Env.df$bio7)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Mean Temperature of Wettest Quarter
TipValue = as.vector(Env.df$bio8)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Mean Temperature of Driest Quarter
TipValue = as.vector(Env.df$bio9)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Mean Temperature of Warmest Quarter
TipValue = as.vector(Env.df$bio10)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Mean Temperature of Coldest Quarter
TipValue = as.vector(Env.df$bio11)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Annual Precipitation
TipValue = as.vector(Env.df$bio12)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Precipitation of Wettest Month
TipValue = as.vector(Env.df$bio13)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Precipitation of Driest Month
TipValue = as.vector(Env.df$bio14)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Precipitation Seasonality
(Coefficient of Variation)
TipValue = as.vector(Env.df$bio15)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Precipitation of Wettest Quarter
TipValue = as.vector(Env.df$bio16)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Precipitation of Driest Quarter
TipValue = as.vector(Env.df$bio17)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Precipitation of Warmest Quarter
TipValue = as.vector(Env.df$bio18)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Precipitation of Coldest Quarter
TipValue = as.vector(Env.df$bio19)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Soil PH)
TipValue = as.vector(Env.df$PH)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Soil Organic Content (g/m^3)
TipValue = as.vector(Env.df$OC)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Soil Conductivity
TipValue = as.vector(Env.df$CE)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Soil Aluminum
TipValue = as.vector(Env.df$AL)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Elevation
TipValue = as.vector(Env.df$Elev)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Heat Loading Index
TipValue = as.vector(Env.df$HLI)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Compound Topographic Index
TipValue = as.vector(Env.df$CTI)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Percent Canopy Cover
TipValue = as.vector(Env.df$pC)
names(TipValue) = Env.df$Label
TemAT = fastAnc(RedTree, TipValue, vars=FALSE,CI=FALSE)
obj = contMap(RedTree, TipValue, plot=FALSE)
plot(obj,
type="phylogram",
legend=0.7*max(nodeHeights(RedTree)))
Session Info
sessionInfo()## R version 3.3.2 (2016-10-31)
## 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] grid stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] phytools_0.5-64 dendextend_1.3.0 FactoMineR_1.34
## [4] dismo_1.1-1 gridExtra_2.2.1 GISTools_0.7-4
## [7] MASS_7.3-45 ggplot2_2.2.0 spdep_0.6-8
## [10] Matrix_1.2-7.1 mapproj_1.2-4 maptools_0.8-40
## [13] Thermimage_2.2.3 rgeos_0.3-21 rasterVis_0.41
## [16] latticeExtra_0.6-28 RColorBrewer_1.1-2 lattice_0.20-34
## [19] fields_8.4-1 maps_3.1.1 spam_1.4-0
## [22] raster_2.5-8 reshape2_1.4.2 rgdal_1.2-4
## [25] sp_1.2-3 phangorn_2.1.1 ape_4.0
## [28] dplyr_0.5.0 rentrez_1.0.4 knitr_1.15.1
##
## loaded via a namespace (and not attached):
## [1] nlme_3.1-128 gmodels_2.16.2
## [3] httr_1.2.1 rprojroot_1.1
## [5] prabclus_2.2-6 numDeriv_2016.8-1
## [7] tools_3.3.2 backports_1.0.4
## [9] R6_2.2.0 DBI_0.5-1
## [11] lazyeval_0.2.0 colorspace_1.3-2
## [13] trimcluster_0.1-2 nnet_7.3-12
## [15] mnormt_1.5-5 curl_2.3
## [17] animation_2.4 expm_0.999-0
## [19] flashClust_1.01-2 NLP_0.1-9
## [21] slam_0.1-40 diptest_0.75-7
## [23] scales_0.4.1 DEoptimR_1.0-8
## [25] hexbin_1.27.1 mvtnorm_1.0-5
## [27] tm_0.6-2 robustbase_0.92-7
## [29] quadprog_1.5-5 stringr_1.1.0
## [31] digest_0.6.10 foreign_0.8-67
## [33] rmarkdown_1.2 htmltools_0.3.5
## [35] plotrix_3.6-3 combinat_0.0-8
## [37] zoo_1.7-13 jsonlite_1.1
## [39] mclust_5.2 gtools_3.5.0
## [41] magrittr_1.5 modeltools_0.2-21
## [43] leaps_2.9 Rcpp_0.12.8
## [45] munsell_0.4.3 scatterplot3d_0.3-37
## [47] stringi_1.1.2 whisker_0.3-2
## [49] clusterGeneration_1.3.4 flexmix_2.3-13
## [51] plyr_1.8.4 parallel_3.3.2
## [53] gdata_2.17.0 deldir_0.1-12
## [55] splines_3.3.2 igraph_1.0.1
## [57] boot_1.3-18 fpc_2.1-10
## [59] stats4_3.3.2 fastmatch_1.0-4
## [61] LearnBayes_2.15 XML_3.98-1.5
## [63] evaluate_0.10 msm_1.6.4
## [65] gtable_0.2.0 kernlab_0.9-25
## [67] assertthat_0.1 coda_0.19-1
## [69] survival_2.40-1 class_7.3-14
## [71] viridisLite_0.1.3 tibble_1.2
## [73] cluster_2.0.5