vegan包taxa2dist函数计算物种分类间距离及建树
vegan包taxa2dist函数计算物种分类间距离及建树
作者：修空调
审核：Listenlii
前言
最近在Listenlii的科研交流群提到了一个函数vegan::taxa2dist，可以直接根据物种分类单元信息计算距离，进而建树，这个功能我是第一次听说。
但是似乎该函数又有些bug（一个门的被搞到不同的支上去了），所以我决定盘一下源代码，学习一下。
代码
taxa2dist<-function (x, varstep = FALSE, check = TRUE, labels)
{
rich <- apply(x, 2, function(taxa) length(unique(taxa)))
S <- nrow(x)
if (check) {
keep <- rich < S & rich > 1
rich <- rich[keep]
x <- x[, keep, drop = FALSE]
}
i <- rev(order(rich))
x <- x[, i, drop = FALSE]
rich <- rich[i]
if (varstep) {
add <- -diff(c(nrow(x), rich, 1))
add <- add/c(S, rich)
add <- add/sum(add) * 100
}
else {
add <- rep(100/(ncol(x) + check), ncol(x) + check)
}
if (!is.null(names(add)))
names(add) <- c("Base", names(add)[-length(add)])
if (!check)
add <- c(0, add)
out <- matrix(add[1], nrow(x), nrow(x))
for (i in 1:ncol(x)) {
out <- out + add[i + 1] * outer(x[, i], x[, i], "!=")
}
out <- as.dist(out)
attr(out, "method") <- "taxa2dist"
attr(out, "steps") <- add
if (missing(labels)) {
attr(out, "Labels") <- rownames(x)
}
else {
if (length(labels) != nrow(x))
warning(gettextf("labels are wrong: needed %d, got %d",
nrow(x), length(labels)))
attr(out, "Labels") <- as.character(labels)
}
if (!check && any(out <= 0))
warning("you used 'check=FALSE' and some distances are zero: was this intended?")
out
}
名词
为了便于理解，解释一下本文提到的几个名词含义：
- 分类信息矩阵：由多个分类层级、分类单元组成的矩阵。分类单元为行，分类层级为列，如：
- 分类单元：可以理解为一个种、一个ASV或别的什么
- 分类层级：指界、门、纲、目、科、属、种。也可以是别的分类层级
参数
x: 分类信息矩阵，可以参考vegan包的实例数据dune.taxon。varstep: 是否针对不同的分类层级采用varstep法设定不同系数。系数的作用下面会讲。check:是否去除每个层级只出现了一次的分类单元和重复出现的分类单元。每个层级只出现了一次的分类单元意思是这个某个种的界门纲目科属在分类信息矩阵中都只出现了一次，如上表中的ASVfungi.labels: 人工设置输出矩阵的labels名称。默认采用分类信息矩阵的行名。
读懂代码
1. 代码并不长，其中最核心的一句是：
out <- out + add[i + 1] * outer(x[, i], x[, i], "!=")
其中outer(x[, i], x[, i], "!=")用来计算两个分类单元是否相等,相等距离为0(FALSE)，不相等为1(TRUE)（其实outer本来是用来计算向量外积的，具体可看参考文献1[1]）。
> demo
[1] "Magnoliidae" "Magnoliidae" "Magnoliidae" "Magnoliidae" "Bryidae" "Bryidae"
> outer(demo, demo, "!=")
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] FALSE FALSE FALSE FALSE TRUE TRUE
[2,] FALSE FALSE FALSE FALSE TRUE TRUE
[3,] FALSE FALSE FALSE FALSE TRUE TRUE
[4,] FALSE FALSE FALSE FALSE TRUE TRUE
[5,] TRUE TRUE TRUE TRUE FALSE FALSE
[6,] TRUE TRUE TRUE TRUE FALSE FALSE
源代码中使用了一个for循环上述计算过程（如下），是针对输入的分类信息矩阵逐个分类层级计算距离，然后相加：
out <- matrix(add[1], nrow(x), nrow(x))
for (i in 1:ncol(x)) {
out <- out + add[i + 1] * outer(x[, i], x[, i], "!=")
}
那么毫无疑问前面的add[i + 1]就是给这个逻辑矩阵加上一个系数，因此下一步我们看看这个系数是怎么加的。
2. add
上面提到了，是针对输入的分类信息矩阵逐个分类层级计算距离，所以显然add的长度是和提供的分类层级数的长度是一致的。
a. 当varstep为FALSE时，add会是一个由同一个数组成的数组，这个数具体多少不是那么重要，因为其不论多少，总之每个分类层级在计算距离上的权重是一致。
b. 当varstep为TRUE时, add的组成就不是一个数了。总的思想是在上一个级分类单元基础上，提供额外信息越多的分类单元占比越大。其具体计算过程如下：
(1) 计算每个分类层级的提供的信息数(即unique)，然后进行排序.
这步实际上是结合生物学背景，划分分类层级的顺序，原理也很简单，例如可能存在5门18属，但是不可能出现18门5属。
> rich <- apply(x, 2, function(taxa) length(unique(taxa)))
> rich
Genus Family Order Superorder Subclass
27 15 10 6 2
> i <- rev(order(rich))
> i
[1] 1 2 3 4 5
> x <- x[, i, drop = FALSE]
> x
Genus Family Order Superorder Subclass
Achimill Achillea Asteraceae Asterales Asteranae Magnoliidae
Agrostol Agrostis Poaceae Poales Lilianae Magnoliidae
Airaprae Aira Poaceae Poales Lilianae Magnoliidae
Alopgeni Alopecurus Poaceae Poales Lilianae Magnoliidae
Anthodor Anthoxanthum Poaceae Poales Lilianae Magnoliidae
Bellpere Bellis Asteraceae Asterales Asteranae Magnoliidae
Bromhord Bromus Poaceae Poales Lilianae Magnoliidae
Chenalbu Chenopodium Amaranthaceae Caryophyllales Caryophyllanae Magnoliidae
Cirsarve Cirsium Asteraceae Asterales Asteranae Magnoliidae
Comapalu Comarum Rosaceae Rosales Rosanae Magnoliidae
Eleopalu Eleocharis Cyperaceae Poales Lilianae Magnoliidae
Elymrepe Elymus Poaceae Poales Lilianae Magnoliidae
Empenigr Empetrum Ericaceae Ericales Asteranae Magnoliidae
Hyporadi Hypochaeris Asteraceae Asterales Asteranae Magnoliidae
Juncarti Juncus Juncaceae Poales Lilianae Magnoliidae
Juncbufo Juncus Juncaceae Poales Lilianae Magnoliidae
Lolipere Lolium Poaceae Poales Lilianae Magnoliidae
Planlanc Plantago Plantaginaceae Lamiales Asteranae Magnoliidae
Poaprat Poa Poaceae Poales Lilianae Magnoliidae
Poatriv Poa Poaceae Poales Lilianae Magnoliidae
Ranuflam Ranunculus Ranunculaceae Ranunculales Ranunculanae Magnoliidae
Rumeacet Rumex Polygonaceae Caryophyllales Caryophyllanae Magnoliidae
Sagiproc Sagina Caryophyllaceae Caryophyllales Caryophyllanae Magnoliidae
Salirepe Salix Salicaceae Malpighiales Rosanae Magnoliidae
Scorautu Scorzoneroides Asteraceae Asterales Asteranae Magnoliidae
Trifprat Trifolium Fabaceae Fabales Rosanae Magnoliidae
Trifrepe Trifolium Fabaceae Fabales Rosanae Magnoliidae
Vicilath Vicia Fabaceae Fabales Rosanae Magnoliidae
Bracruta Brachythecium Brachytheciaceae Hypnales Bryanae Bryidae
Callcusp Calliergonella Hypnaceae Hypnales Bryanae Bryidae
(2) 计算每层分类单元递增的信息情况
#每个分类层级的信息总数，注意最前和最后，开头的30是总物种数，结尾的1可以看做所有的物种都是从1个祖先分化的
> c(nrow(x), rich, 1)
Genus Family Order Superorder Subclass
30 27 15 10 6 2 1
> add <- -diff(c(nrow(x), rich, 1))
> add
Genus Family Order Superorder Subclass
3 12 5 4 4 1
(3)再用逐层递增的信息数除以该层物种总数
add <- add/c(S, rich)
> add
Genus Family Order Superorder Subclass
0.1000000(3/30) 0.4444444(12/27) 0.3333333(5/15) 0.4000000(4/10) 0.6666667(4/6) 0.5000000(1/2)
最后除以总信息数，获得系数
> add <- add/sum(add) * 100
> add
Genus Family Order Superorder Subclass
4.090909 18.181818 13.636364 16.363636 27.272727 20.454545
为什么会分类错误？
domain | phylum | class | order | family | genus | |
|---|---|---|---|---|---|---|
ASV694 | Bacteria | RCP2-54 | unclassified | unclassified | unclassified | unclassified |
ASV389 | Bacteria | Proteobacteria | Gammaproteobacteria | Burkholderiales | Rhodocyclaceae | unclassified |
ASV1571 | Bacteria | Proteobacteria | Alphaproteobacteria | Caulobacterales | Caulobacteraceae | Phenylobacterium |
显然上表中ASV389应该与ASV1571分为一类，但是使用taxa2dist计算时却是ASV389与ASV694聚为了一类。
先不考虑varstep，ASV694与ASV389的距离为4，（domain和genus是一样的）；
ASV1571和ASV389的距离也为4（domain和phylum是一样的）。
那么显然当考虑varstep时genus和phylum的系数就是关键了。
而原始数据中genus这级的“unclassified”过多，R计算时认为phylum提供的信息比genus多，造成出现了phylum的系数比genus大这种情况，从而导致了ASV389与ASV694距离更近这样的错误的出现。
解决方法：
只要保证genus这一级提供的信息比phylum多就可以解决这个问题。
例如将“unclassified”改为"unclassified_Rhodocyclaceae"和"unclassified_RCP2-54"。
致谢
特别感谢华南农大的连腾祥老师提出这个问题！
参考文献
[1] https://www.jianshu.com/p/961e0bf6ed68