根据物种出生-死亡表计算总的系统发育分支长度之和

2024-05-09 • 问答

背景

这个问题有点复杂，因此我将首先介绍背景。为了生成物种出生死亡表（ L 表）的示例，我建议使用dd_sim()包中的DDD函数。

library(DDD)
library(tidyverse)
library(picante)

result <- dd_sim(c(0.2,0.1,20),10) 
# with birth rate 0.2,death rate 0.1,carrying capacity 20 and overall 10 million years.
L <- result$L

L

            [,1] [,2] [,3]       [,4]
 [1,] 10.0000000    0   -1 -1.0000000
 [2,] 10.0000000   -1    2  2.7058965
 [3,]  8.5908774    2    3  6.6301616
 [4,]  8.4786474    3    4  3.3866813
 [5,]  8.4455262   -1   -5 -1.0000000
 [6,]  8.3431071    4    6  3.5624756
 [7,]  5.3784683    2    7  0.6975934
 [8,]  3.8950593    6    8 -1.0000000
 [9,]  1.5032100   -5   -9 -1.0000000
[10,]  0.8393589    7   10 -1.0000000
[11,]  0.6118985   -5  -11 -1.0000000

L 表有4列：

第一列是物种在Mya出生的时间
第二列是物种父本的标签；正值和负值指示该物种属于左冠系还是右冠系
第三列是子代物种本身的标签；正值和负值指示该物种属于左冠系还是右冠系

我做什么

有了这个L表，我现在有了一个现存的社区。我想计算其系统发育多样性（也称为Faith指数）

使用DDD和picante函数可以做到：

# convert L table into community data matrix
comm = as.data.frame(L) %>% dplyr::select(V3:V4) %>%
dplyr::rename(id = V3,pa = V4) %>%
dplyr::mutate(id = paste0("t",abs(id))) %>%
dplyr::mutate(pa = dplyr::if_else(pa == -1,1,0)) %>%
dplyr::mutate(plot = 0) %>%
dplyr::select(plot,pa,id) %>%
picante::sample2matrix()

# convert L table into phylogeny
phy = DDD::L2phylo(L,dropextinct = T)

# calculate Faith's index using pd() function
Faith = picante::pd(comm,phy)

问题

尽管我实现了目标，但该过程似乎是多余且耗时的。我必须来回转换原始的L表，因为我必须使用现有函数。

根据定义，信念指数基本上是群落总系统发育分支长度的总和，所以我的问题是：

是否可以直接从L表计算Faith的索引？

谢谢前进！

您可以简单地使用phy$edge.length生成的phylo对象的DDD::L2phylo组件：

## Measuring the sum of the branch lengths from `phy`
sum_br_length <- sum(phy$edge.length)
sum_br_length == Faith$PD
# [1] TRUE

## Measuring the sum of the branch length from `L`
sum_br_length <- sum(DDD::L2phylo(L,dropextinct = TRUE)$edge.length)
sum_br_length == Faith$PD
# [1] TRUE

一些微基准测试很有趣：

library(microbenchmark)
## Function 1
fun1 <- function(L) {
    comm = as.data.frame(L) %>% dplyr::select(V3:V4) %>%
    dplyr::rename(id = V3,pa = V4) %>%
    dplyr::mutate(id = paste0("t",abs(id))) %>%
    dplyr::mutate(pa = dplyr::if_else(pa == -1,1,0)) %>%
    dplyr::mutate(plot = 0) %>%
    dplyr::select(plot,pa,id) %>%
    picante::sample2matrix()

    # convert L table into phylogeny
    phy = DDD::L2phylo(L,dropextinct = T)

    # calculate Faith's index using pd() function
    Faith = picante::pd(comm,phy)
    return(Faith$PD)
}
## Function 2
fun2 <- function(L) {
    phy <- DDD::L2phylo(L,dropextinct = T)
    return(sum(phy$edge.length))
}
## Function 3
fun3 <- function(L) {
    return(sum(DDD::L2phylo(L,dropextinct = TRUE)$edge.length))
}

## Do all of them give the same results
fun1(L) == Faith$PD
# [1] TRUE
fun2(L) == Faith$PD
# [1] TRUE
fun3(L) == Faith$PD
# [1] TRUE

## Which function fastest?
microbenchmark(fun1(L),fun2(L),fun3(L))
# Unit: milliseconds
#     expr      min       lq     mean   median       uq       max neval
#  fun1(L) 6.486462 6.900641 8.273386 7.445334 8.667535 16.888429   100
#  fun2(L) 1.627854 1.683204 2.215531 1.771219 2.229408  9.522366   100
#  fun3(L) 1.630635 1.663181 2.229206 1.859733 2.448196  7.573001   100

我检查了pd::sample2matrix以了解其内部功能。 tapply调用和以下行似乎是唯一必要的部分。

library(DDD)
library(tidyverse)
library(picante)
#> Loading required package: ape
#> Loading required package: vegan
#> Loading required package: permute
#> Loading required package: lattice
#> This is vegan 2.5-6
#> Loading required package: nlme
#> 
#> Attaching package: 'nlme'
#> The following object is masked from 'package:dplyr':
#> 
#>     collapse
set.seed(100)
result <- dd_sim(c(0.2,0.1,20),10) 
# with birth rate 0.2,death rate 0.1,carrying capacity 20 and overall 10 million years.
L <- result$L

# convert L table into community data matrix
comm_original = as.data.frame(L) %>% dplyr::select(V3:V4) %>%
  dplyr::rename(id = V3,pa = V4) %>%
  dplyr::mutate(id = paste0("t",abs(id))) %>%
  dplyr::mutate(pa = dplyr::if_else(pa == -1,0)) %>%
  dplyr::mutate(plot = 0) %>%
  dplyr::select(plot,id) %>%
  picante::sample2matrix()


# Instead of using dplyr,we'll do some base R operations
# on L. The code doesn't look as nice,but it should be
# significantly faster.
pa <- ifelse(L[,4] == -1,0)
plot <- rep(0,length(pa))
id <- paste0("t",abs(L[,3]))
comm_new <- tapply(pa,list(plot,id),sum)
comm_new[is.na(comm_new)] <- 0
# convert L table into phylogeny
phy = DDD::L2phylo(L,dropextinct = T)

# calculate Faith's index using pd() function
picante::pd(comm_original,phy)
#>         PD SR
#> 0 29.82483  6

picante::pd(comm_new,phy)
#>         PD SR
#> 0 29.82483  6
Created on 2019-11-17 by the reprex package (v0.3.0)

编辑：original()是您最初构建comm的方式，new()是上面给出的方式。如果您将其交换，看来您可以期望将速度提高2倍。我知道这并不是很大的收益，具体取决于工作负载的大小，但总比没有好。

expression      min   median `itr/sec` mem_alloc `gc/sec` n_itr  n_gc total_time result               memory          time    gc            
  <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl> <int> <dbl>   <bch:tm> <list>               <list>          <list>  <list>        
1 original()   9.76ms  10.24ms      96.1     552KB     2.04    47     1      489ms <df[,1125] [1 x 1,1~ <df[,3] [107 x~ <bch:t~ <tibble [48 x~
2 new()        4.57ms   4.84ms     201.      464KB     2.07    97     1      483ms <dbl[,~ <df[,3] [63 x ~ <bch:t~ <tibble [98 x~

根据物种出生-死亡表计算总的系统发育分支长度之和

背景

我做什么

问题

jiejieyanyanhaohao 回答：根据物种出生-死亡表计算总的系统发育分支长度之和

大家都在问