Isso é um R Markdown Notebook depositado no Bioforum. Quando executar o código, os resultados aparecerão abaixo dele!
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Cmd+Shift+Enter.
É imprescindível a instalação da versão R mais atual e, se estiver usando Windows, o compilador de C++ Rtools!
As coordenadas referem-se a um grupo de percevejos. Além disso, não esqueça da filogenia dessas espécies e o vetor referente às províncias biogeográficas atualizadas do Morrone (LÖWENBERG-NETO, 2015; MORRONE, 2022).
Observação: As funções atualizadas encontram-se em https://github.com/jozecaricardo/PanBioGeo
#Função que remove todos os objetos do ambiente:
rm(list = ls())
# instalação dos pacotes:
# install.packages('devtools')
library(devtools)
# nota: não há mais a necessidade de instalar o maptools
# install_github("liamrevell/phytools", force = TRUE)
# install.packages('geiger')
# install.packages('phytools')
# install.packages('viridis'); install.packages('mapdata')
install.packages("letsR"); install.packages('terra', repos='https://rspatial.r-universe.dev')
install.packages('fossil')
install.packages('vegan')
install.packages("spdep"); install.packages('dismo')
install.packages("TreeSearch"); install.packages('spatstat')
install.packages('ctv'); install.packages('Xming')
library(phytools)
library(terra)
library(geiger)
# Indicando o local onde se localiza as coordenadas, a árvore e o shapefile de interesse...
diret1 <- '~/Dropbox/papers_books/panbiogeography/' # formato MAC
diret1 <- 'D:/Dropbox/papers_books/panbiogeography/' # formato Windows
diret1 <- 'local de sua escolha/etc/'lycipta.coords <- read.table(file = paste0(diret1, 'Lycipta.csv', collapse = NULL), sep = ',', dec = '.', h = T)
# Visualização da estrutura dos dados:
# obs: Para utilização de banco de dados próprio, recomendamos organizar o arquivo com o seguinte cabeçalho e coordenadas em decimais: spp, lat, long (tudo minúsculo)! É muito importante deixar nessa ordem!
head(lycipta.coords, 10)# shapefile de trabalho:
diret2 = 'D:/Dropbox/papers_books/panbiogeography/shapes_pan/Lowenberg_Neto_2014_shapefile/'
# ou um endereço de sua escolha:
diret2 <- 'local de sua escolha/etc/' # importante terminar com a barra
morrone <- vect(paste0(diret2, 'NeotropicMap_Geo.shp'), crs = "+proj=longlat +datum=WGS84")
## PLOTANDO OS PONTOS
plot(morrone, axes = T, col = "gray")
# Os dados de Lycipta não estão na ordem que sugerimos. Logo, tivemos que inverter Long e Lat. Se os seus dados estiverem corretos, NÃO execute a linha de código seguinte.
points(lycipta.coords[,3:2], pch = 16, col = "red", cex = 1)
# Se seus dados estiverem corretos, execute a linha de código abaixo:
points(lycipta.coords[,2:3], pch = 16, col = "red", cex = 1)lycipta.tree <- read.tree(file = paste0(diret1, 'Lycipta.tre', collapse = NULL))
# X11(type = 'cairo') # ou windows() para Windows
#plot da árvore:
plotTree(lycipta.tree)\(\color{red}{\text{Só é necessário fazer isso se você tiver espécies com apenas um registro, bem como táxons que estejam na árvore mas não}}\)
\(\color{red}{\text{nos dados e vice-versa!}}\)
# carregar a função que edita a árvore e o arquivo com as coordenadas:
source(paste0(diret1, 'sing_data_frame.R'))
resul <- singleton.to.data.frame(data = lycipta.coords, phylogeny = lycipta.tree)
resul # apresenta duas coisas e, agora, as colunas lat e long mudaram de ordem. Isso deve acontecer!
#########################################################################
#########################################################################
# Para aqueles que não possuem árvore filogenética:
lat <- lycipta.coords[,2]
long <- lycipta.coords[,3]
names(lat) <- lycipta.coords[, 1]
names(long) <- lycipta.coords[, 1]
locais <- cbind(lat, long)
singletons <- NULL # apenas uma ocorrência
sppsORD <- sort(rownames(locais), decreasing = F)
for(i in 1:length(unique(rownames(locais)))){
singletons[i] <- length(grep(pattern = unique(sppsORD)[i], x = sppsORD))
}
# Atenção: se todos os valores do objeto singletons foram maiores que 1, significa que não existem espécies
# com apenas um ponto de ocorrência. Assim, execute o código 1A abaixo.
# código 1A:
if(all(singletons > 1) == TRUE){
spp <- data[, 1]
data_df <- data.frame(spp = spp, long = long, lat = lat)
}
# Entretanto, se tiver algum valor igual a um, execute o código 1B abaixo ao invés do 1A:
# código 1B:
singletons <- which(singletons == 1)
um_ponto <- unique(sppsORD)[singletons]
# now without singletons:
locaisNew <- locais
for(i in 1:length(um_ponto)){
linha <- which(rownames(locaisNew) == um_ponto[i])
locaisNew <- locaisNew[-linha,]
}
# data frame
lat <- locaisNew[, 1]
long <- locaisNew[, 2]
spp <- rownames(locaisNew)
data_df <- data.frame(spp = spp, long = long, lat = lat)# Carregar a função que gera os traços individuais
source(paste0(diret1, 'node_terminal.R'))
# ATENÇÃO ## a função abaixo irá gerar os MST's para todas as espécies presentes no arquivo inicial.
# Neste ponto, definimos a resolução espacial (com o argumento 'resol') que será mantida nas análises subsequentes (PAE_PCE)
# todos os táxons (com árvore)
mst_all_taxa <- terminal_node(coordin = resul$data_df, shape_file = morrone, sobrepo = F, caption = T, resol = c(10, 10), tree = resul$treeMod, seephylog = FALSE)
# todos os táxons (sem árvore)
mst_all_taxa <- terminal_node(coordin = data_df, shape_file = morrone, sobrepo = F, caption = T, resol = c(10, 10), seephylog = FALSE)
# a função gerará uma matriz de presenças e ausências
mst_all_taxa\(\color{blue}{\text{Podemos fazer isso apenas para nós terminais que quisermos...}}\)
# só nós terminais: E. l. triangulator
resul2A <- terminal_node(taxon = 'E_L_triangulator', coordin = resul$data_df, shape_file = morrone, sobrepo = F, caption = T, resol = c(10, 10), tree = resul$treeMod, seephylog = F)
head(resul2A)
# só nós terminais: E. l. triangulator + E. l. aceratos
resul2B <- terminal_node(taxon = c('E_L_triangulator', 'E_L_aceratos'), coordin = resul$data_df, shape_file = morrone, sobrepo = F, caption = T, resol = c(10, 10), tree = resul$treeMod, seephylog = F)
head(resul2B)\(\color{blue}{\text{Antes de quermos ver o MST de integrantes de um nó interno, temos que saber qual o número dos nós internos da árvore!}}\)
plotTree(resul$treeMod)
# ?labelnodes
labelnodes(1:(Ntip(resul$treeMod) + resul$treeMod$Nnode),
1:(Ntip(resul$treeMod) + resul$treeMod$Nnode), interactive = F, cex = 0.6, circle.exp = 0.4)\(\color{blue}{\text{Agora sim! Sabemos os números dos nós internos e, assim, vamos pedir para o R calcular os traços de todos os}}\)
\(\color{blue}{\text{integrantes de um grupo monofilético qualquer.}}\)
# internal nodes: node number 23
resul_node_23 <- terminal_node(coordin = resul$data_df, tree = resul$treeMod, shape_file = morrone, caption = T, resol = c(10, 10), nodes = 23, sobrepo = F, mintreeall = F)O MST do ancestral do grupo monofilético definido pelo nó 23 é o mesmo das espécies separadas integrantes desse nó?
# Altere o código para obter a resposta
resultado <- terminal_node(coordin = resul$data_df, tree = resul$treeMod, shape_file = morrone)# carregar a função que, a partir da matriz de presença e ausência, vai produzir traços generalizados com PAE-PCE:
source(paste0(diret1, 'pae_pce_pan.R'))
# Obs: a matriz resultante terá como base o tamanho de quadrícula definido no BLOCO 2
# producing the trees and a list with grid numbers and the characters (i.e., taxaset)
# RI > 0
# Please try to run at least more than one round, with N > 1...
pae1 <- pae_pce(preabsMat = mst_all_taxa, shapeFile = morrone, resolut = c(10, 10), gridView = TRUE, labelGrid = TRUE, legendSpecies = F, N = 10, sobrepo = F)
pae1 # tabela indicando a ocorrência das espécies nos grids# NOTA: rasterfiles (formato .tiff) serão produzidos, para cada iteração, no diretório "/out"
# raster de uma PAE-PCE:
paepce1 <- rast(file.choose()) # rodada 1 do PAE-PCE
paepce2 <- rast(file.choose()) # rodada 2 do PAE-PCE se existir etc.
# Associando cores às espécies
# cores <- viridis(n = length(unique(resul$data_df$spp)),
# alpha = 0.3, option = 'B')
# shape do Morrone:
plot(morrone, axes = T, lwd = 4) # biorregiões
# rasters dos traços generalizados:
plot(paepce1, add = T)
plot(paepce2, add = T) # se existir...
# espécies ocorrentes nos grids:
sp1 <- vect(file.choose(), crs = "+proj=longlat +datum=WGS84")
# ...
# de forma alternativa, podemos pegar os MSTs...
diret <- 'out'
for(i in unique(pae1[[1]]$spp)){
resul1_shape <- vect(paste0(diret, '/mst_',
i, '.shp'), crs = "+proj=longlat +datum=WGS84")
plot(resul1_shape, col = 'green', lwd = 3, lty = 2, add = T)
}
# podemos colocar apenas os pontos:
especies <- unique(resul$data_df[,1]) # adiciona nomes apenas uma vez, sem repetição
# associando cores às espécies
coresSpp <- setNames(object = viridis(n = length(especies),
option = 'B'), nm = especies)
# fazendo um loop...
for(i in 1:length(cores)){
subData <- subset(resul$data_df, spp == especies[i])
points(subData[, c(2, 3)], col = coresSpp[i], pch = 21)
}
# legenda:
legend(x = "bottomleft", legend = especies, pch = 19,
col = coresSpp, title = 'espécies',
title.col = 'red', pt.cex = 0.6, cex = 0.6)Escolha grupos-irmãos e verifique se sofreram vicariância, com a função terminal_node():
# Altere o código para obter a resposta
resul_node_23 <- terminal_node(coordin = resul$data_df, tree = resul$treeMod, shape_file = morrone, caption = T, resol = c(10, 10), nodes = 23, sobrepo = F, mintreeall = F)Escolha táxons terminais e verifique se sofreram vicariância, com a função terminal_node():
# Altere o código para obter a resposta
resul2A <- terminal_node(taxon = 'E_L_triangulator', coordin = resul$data_df, shape_file = morrone, sobrepo = F, caption = T, resol = c(10, 10), tree = resul$treeMod, seephylog = F)# chame a função is.allopatry...
source(paste0(diret1, 'is_allopatry.R'))
# a função produz um arquivo .txt no diretório /out, chamado "vicarius.txt": devemos SEMPRE usar o resultado mst_all_taxa (os MSTs de todas as espécies da árvore)!
is_allopatry(results = mst_all_taxa, phylogeny = resul$treeMod, coordin = resul$data_df)