【问题标题】:Bray-Curtis Pairwise Analysis in RR中的Bray-Curtis成对分析
【发布时间】:2018-07-05 22:37:56
【问题描述】:

我正在尝试使用 R 中的 Vegan 包计算和可视化配对/合并站点中社区之间的 Bray-Curtis 差异。

下面是一个简化的示例数据框:

Site = c("A", "B", "C", "D", "E", "F", "G", "H", "I", "J") 

PoolNumber = c(1, 3, 4, 2, 4, 1, 2, 3, 4, 4) 

Sp1 = c(3, 10, 7, 0, 12, 9, 4, 0, 4, 3) 

Sp2 = c(2, 1, 17, 1, 2, 9, 3, 1, 6, 7)

Sp3 = c(5, 12, 6, 10, 2, 4, 0, 1, 3, 3)

Sp4 = c(9, 6, 4, 8, 13, 5, 2, 20, 13, 3)

df = data.frame(Site, PoolNumber, Sp1, Sp2, Sp3, Sp4)

“Site”是一个变量,表示每个样本的采集位置 “Sp”列表示每个站点的物种丰度值。 我想比较具有相同“PoolNumber”的网站对,并为每次比较获取一个相异值。

大多数示例建议我应该创建一个仅包含“Sp”列的矩阵并使用以下代码:

matrix <- df[,3:6]

braycurtis = vegdist(matrix, "bray")

hist(braycurtis)

但是,如果我消除带有“PoolNumber”和“Site”的列,我不确定如何告诉 R 要比较哪些行。这是否涉及按“PoolNumber”组织,将其用作行名,然后编写一个循环来比较每 2 行? 我也发现输出难以解释。较低的 Bray-Curtis 值表示更多相似的社区(接近 0 的值),而较高的值(接近 1)表示更多不同的社区,但是有没有办法判断方向性,哪一个更多样化?

我是 R 初学者,因此对于任何滥用术语/格式的行为,我深表歉意。欢迎提出所有建议。

谢谢

【问题讨论】:

  • @JariOksanen,谢谢你的回答!我最终做了一些与你的建议非常相似的事情,对每一对进行子集化,然后使用 vegdist 函数。我发现使用 heatmap(as.matrix(braycurtis)) 可视化整个数据集中的差异很有用,并查看了池中每个站点的 Shannon Weaver 多样性,以更好地了解差异。 Vegan 对这些分析非常有用。感谢您创建和维护它!
  • 我知道这个评论帖子不久前就结束了,但我很想知道你会如何解释热图?谢谢!

标签: r vegan


【解决方案1】:

你的意思是你想得到一个不相似的子集,等于PoolNumbervegdist 函数将为您提供所有不同之处,您可以从中选择您的配对。当您首先将相异性转换为对称矩阵,然后从该对称矩阵中选择您的子集时,这是最简单的:

braycurtis <- vegdist(df[,3:6])
as.matrix(braycurtis)[df$PoolNumber==4,df$PoolNumber==4]
as.dist(as.matrix(braycurtis)[df$PoolNumber==4,df$PoolNumber==4])

如果您只想获得平均值,vegan::meandist 函数将为您提供:

meandist(braycurtis, df$PoolNumber)

这里的对角线值将是PoolNumber 内的平均差异和不同PoolNumbers 之间的非对角线平均差异。看vegan::meandist的代码就知道是怎么做到的了。

Bray-Curtis 相异(与所有正常相异一样)是一种对称度量,它不知道多样性的概念。您可以评估每个站点的多样化程度,但您需要先告诉我们“多样化”是什么意思(多样性还是其他?)。然后,您只需在计算中使用这些值。

如果您只想查看项目(物种)的数量,以下函数将为您提供下三角的差异(上三角的值将与一个符号的开关相同):

designdist(df[,3:6], "A-B", "binary")

或者,您可以使用逐行统计信息并查看它们的差异。这是一个 Shannon-Weaver 多样性指数的例子:

H <- diversity(df[,3:6])
outer(H, H, "-")

要获取子集,与使用 Bray-Curtis 索引的工作类似。

【讨论】:

    猜你喜欢
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 2021-09-18
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    • 1970-01-01
    相关资源
    最近更新 更多