-
Notifications
You must be signed in to change notification settings - Fork 0
/
RevBayes_BSDR.Rmd
128 lines (104 loc) · 5.95 KB
/
RevBayes_BSDR.Rmd
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
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
---
output: github_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE,
eval = F, warning=F, message=F,
fig.path = "RevBayes_BSDR_files/")
```
# RevBayes_BSDR
**Author**: Léo-Paul Dagallier
**Last update**: `r format(Sys.Date())`
***
**Note**: the RevBayes logs have been compressed. To reproduce this notebook, you will first need to decompress the file `output/RevBayes/monodoreae3_BDSE_RevBayes_logs.zip`.
# RevBayes
## Input data
Prepare the path variables:
- in bash:
```{bash, eval = F}
path_to_output="outputs/";
cd $path_to_output
mkdir RevBayes
cd RevBayes
```
- in R:
```{r, eval=TRUE}
path_to_tree = c("data/name_MCC_monodoreae3_monod_pruned.newick")
path_to_output = c("outputs/RevBayes/")
data_suffix = c("monodoreae3")
```
## Branch-Specific Diversification Rate Estimation
See: https://revbayes.github.io/tutorials/divrate/branch_specific.html.
### With both speciation and extinction rate shifts
Prepare several analysis:
6 vs 10 rate categories
1 rate shift, 4 rate shifts, 10 rate shifts
#### Run the analysis
```{bash, eval = F}
mpirun -np 2 rb-mpi mcmc_BDSE_6RateCat_1shift_monodoreae3.Rev
mpirun -np 2 rb-mpi mcmc_BDSE_6RateCat_4shifts_monodoreae3.Rev
mpirun -np 2 rb-mpi mcmc_BDSE_6RateCat_10shifts_monodoreae3.Rev
mpirun -np 2 rb-mpi mcmc_BDSE_10RateCat_1shift_monodoreae3.Rev
mpirun -np 2 rb-mpi mcmc_BDSE_10RateCat_4shift_monodoreae3.Rev
mpirun -np 2 rb-mpi mcmc_BDSE_10RateCat_10shift_monodoreae3.Rev
#2750
mpirun -np 2 rb-mpi mcmc_BDSE_6RateCat_4shifts_fixed_monodoreae3.Rev
mpirun -np 2 rb-mpi mcmc_BDSE_6RateCat_4shifts_free_monodoreae3.Rev
mpirun -np 2 rb-mpi mcmc_BDSE_6RateCat_4shifts_monodoreae3.Rev
```
#### Plot the results
```{r plot-prepare, echo=TRUE, warning=F, message=FALSE}
library(RevGadgets)
library(phytools)
library(tibble)
library(ggtree)
library(treeio)
library(ggplot2)
library(RColorBrewer)
source(file = "R/plot_branch_rates_tree2.R")
# load the files:
my_tree_file = path_to_tree
my_branch_rates_file = paste0(path_to_output, data_suffix, "_BDSE_6RateCat_1shift_rates.log")
my_branch_rates_files = list.files(path_to_output)[grep("rates.log", list.files(path_to_output))]
# set the colors:
Colors <- colorRampPalette(rev(c('darkred',brewer.pal(n = 8, name = "Spectral"),'darkblue')))(100)
```
Branch-specific **speciation** rates:
```{r plot-speciation, echo=TRUE, warning=F}
plot_branch_rates_tree2(tree_file=my_tree_file, branch_rates_file=my_branch_rates_file, parameter_name = "lambda", trans = "identity", colors = Colors) + geom_tiplab(color = "black", size = 2) + scale_color_gradientn("Speciation rate", colors = Colors, trans = "identity") + theme(legend.position=c(0.2,0.85))+ scale_x_continuous(limits = c(0,35))
ggsave(paste0(path_to_output, "RevBayes_", data_suffix, "_BDSE_6RateCat_Speciation_rate.pdf"), width=15, height=15, units="cm")
```
For all the analyses:
```{r plot-speciation-all, echo=TRUE, warning=F, eval=F}
for (f in my_branch_rates_files){
plot_branch_rates_tree2(tree_file=my_tree_file, branch_rates_file=paste0(path_to_output, f), parameter_name = "lambda", trans = "identity", colors = Colors) + geom_tiplab(color = "black", size = 2) + scale_color_gradientn("Speciation rate", colors = Colors, trans = "identity") + theme(legend.position=c(0.2,0.85))+ scale_x_continuous(limits = c(0,35))
ggsave(paste0(path_to_output, "RevBayes_BSDR_Speciation_rate_", gsub(x = f, pattern = ".log", replacement = ""), ".pdf"), width=15, height=15, units="cm")
}
```
Branch-specific **extinction** rates:
```{r plot-extinction, echo=TRUE, warning=F}
plot_branch_rates_tree2(tree_file=my_tree_file, branch_rates_file=my_branch_rates_file, parameter_name = "mu", trans = "identity", colors = Colors) + geom_tiplab(color = "black", size = 2) + scale_color_gradientn("Extinction rate", colors = Colors, trans = "identity") + theme(legend.position=c(0.2,0.85))+ scale_x_continuous(limits = c(0,35))
ggsave(paste0(path_to_output, "RevBayes_", data_suffix, "_BDSE_6RateCat_Extinction_rate.pdf"), width=15, height=15, units="cm")
```
For all the analyses:
```{r plot-extinction-all, echo=TRUE, warning=F, eval=F}
for (f in my_branch_rates_files){
plot_branch_rates_tree2(tree_file=my_tree_file, branch_rates_file=paste0(path_to_output,f), parameter_name = "mu", trans = "identity", colors = Colors) + geom_tiplab(color = "black", size = 2) + scale_color_gradientn("Extinction rate", colors = Colors, trans = "identity") + theme(legend.position=c(0.2,0.85))+ scale_x_continuous(limits = c(0,35))
ggsave(paste0(path_to_output, "RevBayes_BSDR_Extinction_rate_", gsub(x = f, pattern = ".log", replacement = ""), ".pdf"), width=15, height=15, units="cm")
}
```
Branch-specific **net diversification** rates:
```{r plot-netdiv, echo=TRUE, warning=F}
plot_branch_rates_tree2(tree_file=my_tree_file, branch_rates_file=my_branch_rates_file, parameter_name = "net_div", trans = "identity", colors = Colors) + geom_tiplab(color = "black", size = 2) + scale_color_gradientn("Net Diversification rate", colors = Colors, trans = "identity") + theme(legend.position=c(0.2,0.85))+ scale_x_continuous(limits = c(0,35))
ggsave(paste0(path_to_output, "RevBayes_", data_suffix, "_BDSE_6RateCat_Net_Diversification_rate.pdf"), width=15, height=15, units="cm")
```
For all the analyses:
```{r plot-netdiv-all, echo=TRUE, warning=F, eval=F}
for (f in my_branch_rates_files){
plot_branch_rates_tree2(tree_file=my_tree_file, branch_rates_file=paste0(path_to_output,f), parameter_name = "net_div", trans = "identity", colors = Colors) + geom_tiplab(color = "black", size = 2) + scale_color_gradientn("Net diversification rate", colors = Colors, trans = "identity") + theme(legend.position=c(0.2,0.85))+ scale_x_continuous(limits = c(0,35))
ggsave(paste0(path_to_output, "RevBayes_BSDR_Net_Diversification_rate_", gsub(x = f, pattern = ".log", replacement = ""), ".pdf"), width=15, height=15, units="cm")
}
```