picante: R tools for integrating phylogenies and ecology

生態学と系統学の統合的ツール

CRAN: http://cran.r-project.org/web/packages/picante/index.html

> library(picante)

Loading required package: vegan
Loading required package: permute

Attaching package: 'permute'

The following object is masked from 'package:igraph':

    permute

The following object is masked from 'package:devtools':

    check

Loading required package: lattice
This is vegan 2.3-2

Attaching package: 'vegan'

The following object is masked from 'package:VGAM':

    calibrate

The following object is masked from 'package:igraph':

    diversity

Loading required package: nlme

> data("IvesGodfray")
> data("phylocom")

バージョン: 1.6.2

関数名	概略
`IvesGodfray`	Host-parasitoid food web data
`Kcalc`	K statistic of phylogenetic signal
`color.plot.phylo`	Color tip labels based on trait
`comdist`	Calculates inter-community mean pairwise distance
`comdistnt`	Calculates inter-community mean nearest taxon distance
`comm.phylo.cor`	Correlations between species co-occurrence and phylogenetic distances
`comm.phylo.qr`	Quantile regression slopes between species co-occurrence and phylogenetic distances
`cor.table`	Table of correlations and P-values
`df2vec`	Picante utility functions
`evol.distinct`	Species' evolutionary distinctiveness
`expected.pd`	Expected PD, PD Variance, and Edge Abundance Distribution of a phylogeny
`match.phylo.data`	Match taxa in phylogeny and data
`matrix2sample`	Convert community data matrix to Phylocom sample
`mntd`	Mean nearest taxon distance
`mpd`	Mean pairwise distance
`multiPhylosignal`	Calculates phylogenetic signal for data.frame of traits
`pblm`	Phylogenetic Bipartite Linear Model
`pcd`	Phylogenetic Community Dissimilarity
`pd`	Calculate Faith's Phylogenetic Diversity
`phyEstimate`	Phylogenetic estimation of traits for unobserved taxa
`phylocom`	Phylocom default data
`phylosignal`	Measure phylogenetic signal
`phylosor`	Phylogenetic index of beta-diversity PhyloSor
`phylosor.rnd`	Null PhyloSor values of phylogenetic beta-diversity
`phylostruct`	Permutations to Test for Phylogenetic Signal in Community Composition
`picante-package`	picante: R tools for integrating phylogenies and ecology
`prune.sample`	Prune tree to match community data or trait data
`psd`	Phylogenetic Species Diversity Metrics
`randomizeMatrix`	Null models for community data matrix randomization
`raoD`	Rao's quadratic entropy
`readsample`	Read Phylocom sample
`sample2matrix`	Convert Phylocom sample to community data matrix
`ses.mntd`	Standardized effect size of MNTD
`ses.mpd`	Standardized effect size of MPD
`ses.pd`	Standardized effect size of PD
`specaccum.psr`	Phylogenetic Species Richness Sample-Based Rarefaction Curve
`species.dist`	Species co-occurrence distances
`sppregs`	Regressions to Separate Phylogenetic Attraction and Repulsion
`tax.distinctiveness`	Taxonomic distinctiveness sensu Vane-Wright or May
`traitgram`	Draw phylogeny with nodes at trait positions
`unifrac`	Unweighted UniFrac distance between communities
`writesample`	Write a Phylocom community sample file
`writetraits`	Write a Phylocom traits formatted file

IvesGodfray

> data("IvesGodfray")
> IvesGodfray %>% {
+   print(class(.))
+   str(., max.levels = 1)
+ }

[1] "list"
List of 3
 $ host        :'data.frame':    12 obs. of  12 variables:
  ..$ X1 : int [1:12] 6 5 0 2 3 2 4 1 1 1 ...
  ..$ X2 : int [1:12] 5 6 0 2 3 2 4 1 1 1 ...
  ..$ X3 : int [1:12] 0 0 1 0 0 0 0 0 0 0 ...
  ..$ X4 : int [1:12] 2 2 0 4 2 3 2 1 1 1 ...
  ..$ X5 : int [1:12] 3 3 0 2 4 2 3 1 1 1 ...
  ..$ X6 : int [1:12] 2 2 0 3 2 4 2 1 1 1 ...
  ..$ X7 : int [1:12] 4 4 0 2 3 2 6 1 1 1 ...
  ..$ X8 : int [1:12] 1 1 0 1 1 1 1 4 3 3 ...
  ..$ X9 : int [1:12] 1 1 0 1 1 1 1 3 5 4 ...
  ..$ X10: int [1:12] 1 1 0 1 1 1 1 3 4 5 ...
  ..$ X11: int [1:12] 4 4 0 2 3 2 5 1 1 1 ...
  ..$ X12: int [1:12] 1 1 0 1 1 1 1 2 2 2 ...
 $ parasitoid  :'data.frame':    27 obs. of  27 variables:
  ..$ X5 : int [1:27] 4 3 2 2 1 1 0 0 0 0 ...
  ..$ X6 : int [1:27] 3 4 2 2 1 1 0 0 0 0 ...
  ..$ X12: int [1:27] 2 2 4 3 1 1 0 0 0 0 ...
  ..$ X25: int [1:27] 2 2 3 4 1 1 0 0 0 0 ...
  ..$ X26: int [1:27] 1 1 1 1 3 2 0 0 0 0 ...
  ..$ X27: int [1:27] 1 1 1 1 2 3 0 0 0 0 ...
  ..$ X1 : int [1:27] 0 0 0 0 0 0 7 6 5 5 ...
  ..$ X3 : int [1:27] 0 0 0 0 0 0 6 7 5 5 ...
  ..$ X2 : int [1:27] 0 0 0 0 0 0 5 5 7 6 ...
  ..$ X4 : int [1:27] 0 0 0 0 0 0 5 5 6 7 ...
  ..$ X7 : int [1:27] 0 0 0 0 0 0 4 4 4 4 ...
  ..$ X8 : int [1:27] 0 0 0 0 0 0 4 4 4 4 ...
  ..$ X9 : int [1:27] 0 0 0 0 0 0 4 4 4 4 ...
  ..$ X17: int [1:27] 0 0 0 0 0 0 3 3 3 3 ...
  ..$ X18: int [1:27] 0 0 0 0 0 0 3 3 3 3 ...
  ..$ X23: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
  ..$ X24: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
  ..$ X22: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
  ..$ X19: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
  ..$ X20: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
  ..$ X21: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
  ..$ X10: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
  ..$ X11: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
  ..$ X13: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
  ..$ X14: int [1:27] 0 0 0 0 0 0 1 1 1 1 ...
  ..$ X15: int [1:27] 0 0 0 0 0 0 1 1 1 1 ...
  ..$ X16: int [1:27] 0 0 0 0 0 0 2 2 2 2 ...
 $ interactions:'data.frame':    12 obs. of  28 variables:
  ..$ X5     : num [1:12] 0 0 0 0 0 ...
  ..$ X6     : num [1:12] 0 0 0 0 0 ...
  ..$ X12    : num [1:12] 0 0 1.51 7.19 2.75 ...
  ..$ X25    : num [1:12] 25.47 1.69 0 0 0 ...
  ..$ X26    : num [1:12] 0 0 0.0483 0.0136 0.0209 ...
  ..$ X27    : num [1:12] 3.868 0.182 0 0 0 ...
  ..$ X1     : num [1:12] 1.3471 0.1198 0.0937 0.2327 0.1262 ...
  ..$ X3     : num [1:12] 0.4199 0.0304 0 0.0407 0 ...
  ..$ X2     : num [1:12] 0 0 0 0.0136 0 ...
  ..$ X4     : num [1:12] 0 0.0304 0.1731 1.115 0.2095 ...
  ..$ X7     : num [1:12] 0.7729 0.1773 0.0923 0.2893 0.0133 ...
  ..$ X8     : num [1:12] 1.515 0.136 0.736 1.212 0.149 ...
  ..$ X9     : num [1:12] 0 0 0 0.0407 0 ...
  ..$ X17    : num [1:12] 0 0 0 0 0 ...
  ..$ X18    : num [1:12] 0.2838 0 0.0866 0.0136 0 ...
  ..$ X23    : num [1:12] 0.0697 0 0.0758 0.2972 0.2809 ...
  ..$ X24    : num [1:12] 0.3533 0 0.029 0.2031 0.0643 ...
  ..$ X22    : num [1:12] 0.147 0 0 0 0 ...
  ..$ X19    : num [1:12] 0 0 0 0.0362 0 ...
  ..$ X20    : num [1:12] 0.1009 0 0 0.0497 0 ...
  ..$ X21    : num [1:12] 0.0625 0 0 0.0362 0 ...
  ..$ X10    : num [1:12] 1.198 0.0911 0.0816 0.3546 0.3097 ...
  ..$ X11    : num [1:12] 0 0 0 0 0 0 0 0 0 0 ...
  ..$ X13    : num [1:12] 1.5063 0.1519 0.0136 0.0893 0.0653 ...
  ..$ X14    : num [1:12] 0.0313 0 0 0 0 ...
  ..$ X15    : num [1:12] 0.0536 0 0 0 0 ...
  ..$ X16    : num [1:12] 0 0 0 0 0 ...
  ..$ density: num [1:12] 74.83 6.81 11.18 27.95 7.79 ...

Kcalc

Blombergの系統シグナル$K$の算出。$K$のみ。有意水準Pは求めない（より細かな解析はphylosignal()）

> set.seed(71)
> randtree <- ape::rcoal(4)
> randtraits <- ape::rTraitCont(randtree)
> Kcalc(randtraits, randtree)

          [,1]
[1,] 0.2697736

cor.table

形質データの標準値と対比値それぞれの相関を求める

Arguments

x
cor.method... pearson or spearman
cor.type... standard or contrast

> # ape::pic() の仮想形質データ
> tmp_df <- data.frame(
+   X = c(4.09434, 3.61092, 2.37024, 2.02815, -1.46968),
+   Y = c(4.74493, 3.33220, 3.36730, 2.89037, 2.30259))
> cor.table(x          = tmp_df, 
+           cor.method = "pearson",
+           cor.type   = "contrast")

$r
          X         Y
X 1.0000000 0.8479677
Y 0.8479677 1.0000000

$df
[1] 4

$P
           X          Y
X 0.00000000 0.03291371
Y 0.03291371 0.00000000

ape::pic()と異なり、対比値を求める際に枝の長さを考慮していない（単純な形質値の差を使用？）

multiPhylosignal

系統シグナルの算出

pcd

系統的な群集非類似度

> phylocom %$% pcd(comm = sample, tree = phylo)

$PCD
           clump1   clump2a   clump2b    clump4      even
clump2a 1.2127013                                        
clump2b 1.1380509 1.1117960                              
clump4  1.6254549 0.9951988 0.9337467                    
even    1.5154868 1.3749415 1.2782336 0.7005061          
random  1.5133958 1.3614215 1.2448956 1.2047632 0.9288258

$PCDc
           clump1   clump2a   clump2b    clump4      even
clump2a 0.7352941                                        
clump2b 0.7352941 0.7352941                              
clump4  1.1029412 0.7352941 0.7352941                    
even    1.1029412 1.1029412 1.1029412 0.7352941          
random  1.1029412 1.1029412 1.1029412 1.1029412 1.1029412

$PCDp
           clump1   clump2a   clump2b    clump4      even
clump2a 1.6492737                                        
clump2b 1.5477492 1.5120425                              
clump4  1.4737458 1.3534704 1.2698955                    
even    1.3740414 1.2466136 1.1589318 0.9526883          
random  1.3721456 1.2343555 1.1287053 1.0923186 0.8421354

$PSVmncd
[1] 0.7482580 0.6806116 0.6265613 0.5681849 0.5169596 0.4756883 0.4334882
[8] 0.3938833

$PSVpool
[1] 0.8233333

pd

Faithの系統的多様度を算出

> phylocom %$% pd(samp = sample, tree = phylo,
+                 include.root = TRUE)

        PD SR
clump1  16  8
clump2a 17  8
clump2b 18  8
clump4  22  8
even    30  8
random  27  8

phylocom

系統、群集および形質値を含んだ標準データセット

Format

3つの要素を含んだリスト

phylo... 系統樹データ（phyloクラス）
sample... 群集ごとの種データ（matrixクラス）
traits... sampleに含まれる種の４形質のデータ（data.frameクラス）

> data("phylocom")
> phylocom %>% {
+   class(.) %>% print(.)
+   str(.)
+ }

[1] "list"
List of 3
 $ phylo :List of 5
  ..$ edge       : num [1:62, 1:2] 33 34 35 36 37 37 36 38 38 35 ...
  ..$ tip.label  : chr [1:32] "sp1" "sp2" "sp3" "sp4" ...
  ..$ Nnode      : int 31
  ..$ node.label : chr [1:31] "A" "B" "C" "D" ...
  ..$ edge.length: num [1:62] 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "class")= chr "phylo"
 $ sample: num [1:6, 1:25] 1 1 1 1 1 0 0 2 0 1 ...
  ..- attr(*, "dimnames")=List of 2
  .. ..$ : chr [1:6] "clump1" "clump2a" "clump2b" "clump4" ...
  .. ..$ : chr [1:25] "sp1" "sp10" "sp11" "sp12" ...
 $ traits:'data.frame':    32 obs. of  4 variables:
  ..$ traitA: int [1:32] 1 1 2 2 2 2 2 2 1 1 ...
  ..$ traitB: int [1:32] 1 1 1 1 2 2 2 2 3 3 ...
  ..$ traitC: int [1:32] 1 2 3 4 1 2 3 4 1 2 ...
  ..$ traitD: int [1:32] 0 0 0 0 0 0 0 0 1 1 ...

> phylocom$sample %>% dplyr::glimpse()

Observations: 6
Variables: 25
$ sp1  (dbl) 1, 1, 1, 1, 1, 0
$ sp10 (dbl) 0, 2, 0, 1, 0, 0
$ sp11 (dbl) 0, 2, 0, 0, 0, 0
$ sp12 (dbl) 0, 2, 0, 0, 0, 1
$ sp13 (dbl) 0, 0, 0, 0, 1, 0
$ sp14 (dbl) 0, 0, 0, 0, 0, 4
$ sp15 (dbl) 0, 0, 0, 0, 0, 2
$ sp17 (dbl) 0, 0, 2, 2, 1, 3
$ sp18 (dbl) 0, 0, 2, 2, 0, 0
$ sp19 (dbl) 0, 0, 2, 0, 0, 0
$ sp2  (dbl) 1, 1, 1, 1, 0, 1
$ sp20 (dbl) 0, 0, 2, 0, 0, 0
$ sp21 (dbl) 0, 0, 0, 0, 1, 0
$ sp22 (dbl) 0, 0, 0, 0, 0, 1
$ sp24 (dbl) 0, 0, 0, 0, 0, 2
$ sp25 (dbl) 0, 0, 0, 2, 1, 0
$ sp26 (dbl) 0, 0, 0, 2, 0, 0
$ sp29 (dbl) 0, 0, 0, 0, 1, 0
$ sp3  (dbl) 1, 1, 1, 0, 0, 0
$ sp4  (dbl) 1, 1, 1, 0, 0, 0
$ sp5  (dbl) 1, 0, 0, 0, 1, 2
$ sp6  (dbl) 1, 0, 0, 0, 0, 0
$ sp7  (dbl) 1, 0, 0, 0, 0, 0
$ sp8  (dbl) 1, 0, 0, 0, 0, 0
$ sp9  (dbl) 0, 2, 0, 1, 1, 0

phylosignal

Blombergの系統シグナル$K$を求める。$K$を計算するだけならKcalc()で良い

Arguments

x
phy
reps... ランダマイゼーション繰り返し回数
chackdata
...

> set.seed(100)
> randtree <- rcoal(20)
> randtraits <- rTraitCont(randtree)
> phylosignal(randtraits[randtree$tip.label], randtree)

          K PIC.variance.obs PIC.variance.rnd.mean PIC.variance.P
1 0.3538194        0.0134462            0.05240829          0.006
  PIC.variance.Z
1      -1.476381

specaccum.psr

ref) vegan::specaccum()

Arguments

samp... 群集行列データ
tree... 系統樹オブジェクト
permutations
methods
...

> set.seed(71)
> phylocom %$% vegan::specaccum(sample, 
+                               permutations = 200, 
+                               method       = "random") %>% 
+   {
+     print(.)
+     accum.sr <<- .
+   }
> set.seed(71)
> phylocom %$% specaccum.psr(sample, 
+                            tree         = phylo, 
+                            permutations = 100, 
+                            method       = "random") %>% 
+   {
+     print(.)
+     accum.psr <<- .
+   }
> 
> accum.sr %>% {
+   plot(., col = "blue")
+   points(.$sites, accum.sr$richness, pch = 19, col = "blue")
+ }
> accum.psr %>% {
+   plot(., add = TRUE, col = "red")
+   points(.$sites, accum.psr$richness, pch = 19, col = "red")
+ }

traitgram

系統樹のノードにおける形質値のばらつきをプロットする

Arguments

x
phy
xaxt
underscore
show.names... logic. tip名を表示するか
show.xaxis.values
method... 系統樹上の中間形質値の算出方法。ML or pic

> set.seed(100)
> randtree <- rcoal(20)
> randtree$edge.length
> (randtraits <- rTraitCont(randtree))
> traitgram(x = randtraits, phy = randtree, method = "ML")
> # traitgram(randtraits, randtree, method = "pic")