forked from pwilliams0/Biogeography_and_global_diversity
-
Notifications
You must be signed in to change notification settings - Fork 0
/
2_PhyloB_Bird_24.R
67 lines (58 loc) · 2.47 KB
/
2_PhyloB_Bird_24.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#===========================================================#
# Calculate phylobetadiversity for birds, cell groups 2 & 4 #
#===========================================================#
library(tidyverse)
library(betapart)
library(picante)
library(ape)
# Split cells into 5 groups, every fifth cell
cell_nums <- unique(read.csv(file="Data/bird_cell_df.csv",
stringsAsFactors=FALSE)$cell_id)
group_1 <- cell_nums[c(TRUE, FALSE, FALSE, FALSE, FALSE)]
group_2 <- cell_nums[c(FALSE, TRUE, FALSE, FALSE, FALSE)]
group_3 <- cell_nums[c(FALSE, FALSE, TRUE, FALSE, FALSE)]
group_4 <- cell_nums[c(FALSE, FALSE, FALSE, TRUE, FALSE)]
group_5 <- cell_nums[c(FALSE, FALSE, FALSE, FALSE, TRUE)]
# Select set of cells to calculate phylobetadiversity
cell_list <- c(group_2, group_4)
# Create community matrix
bird_comm_mat <- read.csv(file="Data/bird_cell_df.csv",
stringsAsFactors=FALSE) %>%
dplyr::filter(cell_id %in% cell_list) %>%
group_by(cell_id, names_phylo) %>%
slice(1) %>% ungroup() %>%
mutate(present = 1) %>%
dplyr::select(cell_id, names_phylo, present) %>%
pivot_wider(names_from = names_phylo,
values_from = present) %>%
replace(.,is.na(.),0) %>%
mutate(cell_id = paste("X",cell_id,sep="")) %>%
column_to_rownames("cell_id") %>%
as.matrix()
bird_comm_mat <- bird_comm_mat[,order(colnames(bird_comm_mat))]
# Load phylogeny, prune for species in community matrix
bird_phylo <- picante::prune.sample(
bird_comm_mat,
ape::read.tree("Data/Raw/Stage2_MayrParSho_Hackett_mean phylogeny_ultrametric.tre"))
print("Done loading data")
# Run phylogenetic beta diversity (Sorenson)
start_time <- Sys.time()
bird_pb <- betapart::phylo.beta.pair(
bird_comm_mat,
bird_phylo,
index.family="sorensen")
end_time <- Sys.time()
print(end_time - start_time)
saveRDS(bird_pb, "Data/bird_pb_24.rds")
# sor is total phylobetdiveristy (turnover + nestedness)
# sim is phylobetdiveristy turnover
# sne is phylobetdiveristy nestedness
write.table(as.matrix(bird_pb$phylo.beta.sim),
file="Data/bird_pb_sim_24.txt", col.names=TRUE, row.name=TRUE, sep="\t", quote=FALSE)
#write.table(as.matrix(bird_pb$phylo.beta.sor),
# file="Data/bird_pb_sor_24.txt", col.names=TRUE, row.name=TRUE, sep="\t", quote=FALSE)
#write.table(as.matrix(bird_pb$phylo.beta.sne),
# file="Data/bird_pb_sne_24.txt", col.names=TRUE, row.name=TRUE, sep="\t", quote=FALSE)
print("Done with bird_24")
# [Runtime 11 hr 49 min]
# [253 GB memory used]