1. 项目概述从群落数据到距离矩阵如果你刚开始接触生态学、微生物组学或者任何涉及群落物种组成分析的研究那么“距离”这个概念很快就会成为你绕不开的核心。我们手里通常有一张表格行是样本比如不同土壤、不同水体、不同病人的肠道列是物种比如细菌OTU、植物种类、基因家族表格里填的是每个物种在每个样本中的丰度。直接比较两个样本看它们像不像总不能一个个数物种吧这时候我们就需要一个量化的“尺子”来测量样本间的差异这把尺子就是距离或相异性。vegdist()函数来自R语言中生态学数据分析的“瑞士军刀”包——vegan就是专门用来制造这把尺子的工具。它不生产数据它只是距离的计算工。简单来说你给它一个物种丰度矩阵它就能给你算出一个样本两两之间的距离矩阵。这个距离矩阵是后续进行排序分析如NMDS、PCoA、聚类分析、差异检验如PERMANOVA的绝对基石。我刚开始用的时候以为这不过就是个算距离的函数随便选个方法就行。结果在分析微生物组数据时用错了距离算法导致后续的PERMANOVA结果完全解释不通白白浪费了两周时间排查数据问题。所以今天我就结合这些年踩过的坑和积累的经验把vegdist()里里外外讲透让你不仅能“会用”更能“懂用”知道在什么场景下该选哪把“尺子”。2. 核心概念与函数参数全解在深入代码之前我们必须统一思想什么是“距离”在群落生态学中我们通常计算的是相异性值越大表示两个样本越不相似。有些指数如Bray-Curtis取值范围是0到1有些如欧氏距离则没有上限。vegdist()函数的核心任务就是根据你指定的算法将样本对的物种组成向量转化为一个代表差异的数值。2.1 函数基本语法与参数vegdist()函数的基本调用格式如下vegdist(x, method bray, binary FALSE, diag FALSE, upper FALSE, na.rm FALSE, ...)别看参数不多每一个都至关重要选错了直接影响结果。x: 这是输入数据通常是一个数值矩阵或数据框。行是样本列是物种/变量。这是最容易出错的地方务必确保你的数据矩阵是这个方向。method: 这是核心中的核心指定计算距离的方法。vegan包内置了丰富的选项也是我们重点讲解的对象。默认是brayBray-Curtis相异性。binary: 逻辑值TRUE/FALSE。如果设为TRUE会在计算距离前将丰度数据转换为“有1/无0”的二元数据。这相当于先做了一次“存在与否”的转换再计算距离。对于某些关注物种有无而非丰度的研究如分布地理学很有用。diag: 逻辑值。是否在输出的距离矩阵中打印对角线上的值样本到自身的距离通常为0。upper: 逻辑值。是否以上三角矩阵的形式打印输出只显示矩阵右上角部分。na.rm: 逻辑值。是否在计算时移除缺失值NA。对于群落数据NA可能意味着未检测或真实缺失需要根据研究背景谨慎处理。注意diag和upper参数只影响距离矩阵的显示方式不影响其作为dist对象的内在结构和后续分析。在R中dist对象是一种高效存储对称距离矩阵下三角部分的数据类型。2.2 关键方法method深度解析选择哪种距离方法取决于你的数据特性和科学问题。下面我把最常用的几种方法掰开揉碎了讲。2.2.1 Bray-Curtis 相异性 (method bray)这是生态学中最经典、最常用的距离之一尤其适用于群落丰度数据。公式思想计算两个样本共有物种的丰度绝对值差之和然后除以两个样本的总丰度之和。公式为BC sum|A_i - B_i| / sum(A_i B_i)。其中A_i和B_i是物种i在两个样本中的丰度。特点与适用场景关注群落组成它对物种组成的变化敏感既考虑物种有无也考虑丰度差异。不受零值过度影响大量两个样本都没有的物种双零值不会增加它们的相似性这是与某些相关系数距离的关键区别。这在生态学上是合理的两个沙漠都没有某树种并不能说明它们相似。对丰度敏感一个物种在样本A中是100在样本B中是10与在A中是10在B中是1所产生的差异贡献是不同的。范围在[0,1]0表示两个样本完全相同1表示完全不同没有共有物种。实操心得对于绝大多数微生物16S rRNA基因测序得到的OTU/ASV丰度表、宏基因组物种组成表Bray-Curtis通常是首选的基准距离。它非常稳健我建议在初次分析时都从它开始。2.2.2 Jaccard 相异性 (method jaccard)这是一个基于物种有无二元数据的相异性指数。公式思想J (b c) / (a b c)。其中a是两个样本共有的物种数b是样本1有而样本2无的物种数c是样本2有而样本1无的物种数。特点与适用场景忽略丰度只关心有无它把所有的丰度信息都丢掉了只关注物种是否存在。这对于某些类型的DNA指纹数据如DGGE、T-RFLP或关注物种分布格局的研究特别有用。双零值无影响和Bray-Curtis一样双零值不贡献相似性。范围在[0,1]。与Bray-Curtis的关系当你设置binary TRUE并使用method bray时你计算的就是Jaccard距离。因为Bray-Curtis公式在二元数据下会简化为Jaccard。这是一个需要记住的等价关系。实操心得当你怀疑样本间的差异主要源于物种的“出现-消失”比如强环境过滤导致某些物种完全不存在而不是丰度的增减时可以用Jaccard距离来验证。与Bray-Curtis的结果对比如果两者格局差异很大可能说明你数据中稀有物种低丰度但广泛存在的影响很显著。2.2.3 欧氏距离 (method euclidean)这是最广为人知的几何距离但在群落数据分析中需要格外小心。公式就是多维空间中点与点的直线距离sqrt(sum((A_i - B_i)^2))。特点与适用场景对丰度绝对值敏感一个在所有样本中丰度都很高的物种其微小的相对变化就能对欧氏距离产生巨大贡献。这可能导致结果被少数高丰度物种主导。受双零值影响两个样本都没有的物种差值为0不增加距离。这看起来合理但在高维稀疏的群落数据中这会导致拥有大量共有“零值”的样本被拉近可能产生误导。没有上限。在群落数据中的问题原始丰度数据直接计算欧氏距离通常不是好主意。它没有标准化过程样本总测序深度文库大小的差异会严重影响距离。样本A总读段100万样本B总读段10万即使组成比例相似欧氏距离也会很大。如何正确使用欧氏距离通常用于经过转化的数据。例如对丰度数据进行Hellinger转化decostand(x, hellinger)后再计算欧氏距离是一个非常有效且数学性质良好的方法常用于基于冗余分析RDA的模型。所以不要轻易对原始OTU表用method euclidean。2.2.4 UniFrac 距离这是一个基于系统发育信息的距离特别适用于微生物组数据。它分为未加权UniFrac只考虑分支有无和加权UniFrac同时考虑分支长度和物种丰度。vegdist()函数本身不直接计算UniFrac距离但vegan包通过distance()函数或专门的phyloseq包可以方便地调用。由于其重要性这里简要提及其思想核心思想比较两个样本的微生物群落时不仅看物种是否相同还看它们系统发育上的远近。丢失一个独有物种和丢失一个在多个样本中常见的物种的近亲意义是不同的。适用场景当你拥有物种的系统发育树如16S数据通过QIIME2、mothur等流程生成时强烈建议使用UniFrac距离。它能揭示基于进化关系的群落差异。为了更直观地对比我将常用方法总结如下表方法 (method)核心关注点对双零值的处理对丰度的敏感性典型适用场景注意事项bray(Bray-Curtis)物种组成与丰度忽略不增加相似性敏感绝大多数群落丰度数据微生物组、动植物群落默认选择稳健通用jaccard物种有无存在/缺失忽略不敏感物种分布研究、二元化数据相当于binaryTRUE时的brayeuclidean(欧氏距离)多维空间几何距离视为相同差为0极度敏感对绝对值不推荐直接用于原始群落数据需先进行数据标准化/转化如Hellingermanhattan(曼哈顿距离)丰度绝对差异之和视为相同差为0敏感对绝对值可作为某些分析的替代但不如Bray-Curtis常用受总丰度影响大kulczynski对稀有物种更宽容的Bray-Curtis变体忽略敏感群落中存在大量低丰度物种时比Bray-Curtis更稳定unifrac(需通过其他函数)系统发育分支的共享情况由系统发育树定义加权版本敏感拥有系统发育树的微生物组数据能揭示进化尺度的差异3. 完整实操流程从数据到距离矩阵理论说再多不如亲手跑一遍。我们假设你手头有一个名为otu_table.csv的OTU丰度表现在我们来完成从数据导入、检查、预处理到计算距离的全过程。3.1 数据准备与检查首先加载必要的包并读入数据。# 安装并加载vegan包 # install.packages(vegan) # 如果未安装需先运行此命令 library(vegan) # 读入数据。假设你的数据是CSV格式第一列是OTU ID第一行是样本名 # 注意read.csv默认会把行名放在第一列我们需要正确处理 otu_raw - read.csv(otu_table.csv, row.names 1, check.names FALSE) # 查看数据前6行和前6列了解数据结构 head(otu_raw[, 1:6]) dim(otu_raw) # 查看数据维度行数(物种数) x 列数(样本数)关键检查点1数据方向。dim()输出的结果通常应该是[物种数目, 样本数目]。vegdist()要求样本在行物种在列。如果你的数据是转置的样本在列需要先转置otu_raw - t(otu_raw)。关键检查点2数据格式。用str(otu_raw)或class(otu_raw)查看确保它是一个data.frame或matrix且内部的数值都是numeric整数或小数。如果有非数值列如分类信息需要先拆分出去。关键检查点3缺失值与零。用sum(is.na(otu_raw))检查是否有NA。群落数据中NA可能代表未检测通常需要根据情况处理如视为0或移除。用sum(otu_raw 0) / length(otu_raw)可以计算数据的稀疏度零的比例微生物组数据通常非常稀疏70%的零。3.2 数据预处理标准化对于群落数据直接计算距离前往往需要进行标准化以消除样本间总测序深度文库大小不同带来的影响。最常用的方法是总和标准化Total Sum Scaling即将每个样本的计数转换为相对丰度百分比。# 方法1使用vegan包的decostand函数进行总和标准化 otu_relab - decostand(otu_raw, method total) # 检查每个样本的总和现在应该是1或100如果乘以100的话 colSums(otu_relab)[1:5] # 方法2手动计算原理相同 # otu_relab - apply(otu_raw, 2, function(x) x / sum(x)) # 注意apply的第二个参数2表示按列样本计算。如果你的数据样本在行则应为1。重要提示vegdist()函数在计算某些距离如bray,jaccard时内部会进行与算法相关的标准化处理。例如Bray-Curtis计算时每个样本的丰度在公式分母中已经被总和考虑了。因此对于Bray-Curtis你既可以使用原始计数也可以使用相对丰度两者计算出的距离矩阵在数值上可能不同但样本间的相对距离关系排序通常高度一致。我个人的习惯是对于Bray-Curtis我倾向于使用原始计数让函数内部处理而对于计划使用欧氏距离的数据则必须先进行标准化或Hellinger转化。3.3 计算距离矩阵并解读现在我们来计算最常用的Bray-Curtis距离矩阵。# 使用相对丰度数据计算Bray-Curtis距离 dist_bray - vegdist(t(otu_relab), method bray) # 注意如果otu_relab是物种为行样本为列需要转置(t)它因为vegdist要求样本在行。 # 如果上一步你的数据已经是样本在行则不需要t()。 # 查看距离对象 dist_bray # 输出是一个‘dist’对象只显示了下三角部分节省空间。 # 你可以看到距离的大致范围确认是否在[0,1]之间。 # 查看距离矩阵的维度样本数 attr(dist_bray, Size) # 或者用nrow(as.matrix(dist_bray)) # 将dist对象转换为矩阵以便查看具体值 dist_matrix - as.matrix(dist_bray) # 查看前4个样本之间的距离 dist_matrix[1:4, 1:4]这个dist_bray对象就是后续所有分析的起点。你可以把它输入到metaMDS()函数做NMDS排序输入到hclust()函数做层次聚类或者输入到adonis2()函数做PERMANOVA分析。3.4 不同距离方法的对比计算为了理解不同方法带来的差异我们可以同时计算几种距离并简单比较。# 计算几种常用距离 dist_jaccard - vegdist(t(otu_relab), method jaccard) # 注意这里用的还是丰度数据但jaccard方法会将其视为二元数据不 # 等一下这里有坑对于丰度数据直接使用method“jaccard”vegan实际上使用的是基于丰度的Jaccard变体如“jaccard”对应的是“binomial”。如果想用经典的二元Jaccard应该 # 正确计算经典二元Jaccard距离先转换为有无 otu_pa - decostand(otu_raw, method pa) # “pa”即 presence-absence (0/1) dist_jaccard_binary - vegdist(t(otu_pa), method jaccard) # 或者用 vegdist(t(otu_raw), methodjaccard, binaryTRUE) # 计算欧氏距离在相对丰度数据上 dist_euclidean - vegdist(t(otu_relab), method euclidean) # 计算曼哈顿距离 dist_manhattan - vegdist(t(otu_relab), method manhattan) # 简单比较查看第一个样本与其他样本在不同距离下的值 sample1_distances - data.frame( Sample rownames(dist_matrix)[2:attr(dist_bray, Size)], Bray_Curtis dist_matrix[2:attr(dist_bray, Size), 1], Jaccard_Binary as.matrix(dist_jaccard_binary)[2:attr(dist_bray, Size), 1], Euclidean as.matrix(dist_euclidean)[2:attr(dist_bray, Size), 1] ) head(sample1_distances)通过这个对比你可以直观感受不同度量标准下样本间“差异”的绝对大小和排序是否一致。通常Bray-Curtis和二元Jaccard的结果可能差异较大这反映了丰度信息的重要性。4. 高级应用与结果可视化得到距离矩阵不是终点而是起点。这里介绍两个最直接的应用可视化与统计检验。4.1 距离矩阵的可视化热图与聚类树热图是展示距离矩阵最直观的方式。# 加载绘图需要的包 library(pheatmap) # 或者用ggplot2扩展包但pheatmap最简单 # 使用pheatmap绘制距离矩阵热图 pheatmap(as.matrix(dist_bray), cluster_rows TRUE, # 对行样本聚类 cluster_cols TRUE, # 对列样本聚类通常与行一致 clustering_distance_rows dist_bray, # 聚类使用的距离 clustering_distance_cols dist_bray, clustering_method average, # 聚类方法可选ward.D, complete, average等 main Bray-Curtis Dissimilarity Heatmap, color colorRampPalette(c(navy, white, firebrick3))(100) # 自定义颜色梯度 )从热图中你可以快速看出哪些样本彼此更相似颜色偏蓝/浅哪些差异更大颜色偏红。聚类树状图展示了样本的层次分组关系。4.2 基于距离的统计检验PERMANOVA入门如果你想检验不同分组如处理组 vs 对照组的群落结构是否有显著差异PERMANOVA通过vegan包的adonis2函数实现是最常用的方法。它的原理是基于距离矩阵进行方差分析。# 假设你有一个分组信息的数据框metadata其中有一列“Group”表示样本的分组 # metadata的行名需要与距离矩阵的样本名一致 # 确保样本顺序一致 rownames(metadata) - metadata$SampleID # 假设SampleID是样本名列 common_samples - intersect(rownames(metadata), rownames(dist_matrix)) metadata - metadata[common_samples, ] dist_matrix_sub - as.matrix(dist_bray)[common_samples, common_samples] # 运行PERMANOVA # 注意adonis2要求输入数据是dist对象而不是矩阵 permanova_result - adonis2(dist_bray ~ Group, data metadata, permutations 999) # 公式 dist ~ Group 表示检验Group分组对距离的解释程度 # permutations 999 表示使用999次置换检验来计算p值 # 查看结果 print(permanova_result)结果中你会关注R2值类似于回归中的R-squared表示分组变量能解释的距离方差的比例和Pr(F)值p值。一个显著的p值如0.05表明不同分组间的群落结构存在统计学差异。重要警告PERMANOVA的一个关键前提是组内离散度同质性类似于方差齐性。如果不同分组的样本在其多维空间内的分散程度离散度差异很大PERMANOVA的结果可能不可靠容易产生假阳性。在报告PERMANOVA结果前务必用betadisper()函数检验离散度同质性。# 检验离散度同质性 dispersion - betadisper(dist_bray, group metadata$Group) anova(dispersion) # 查看离散度差异的ANOVA检验结果 permutest(dispersion, permutations 999) # 置换检验更稳健 plot(dispersion) # 可视化各组离散度主坐标分析如果离散度检验结果显著p0.05说明组间离散度不同此时PERMANOVA的结果需要谨慎解释或者考虑使用对离散度差异不敏感的其他方法。5. 常见陷阱、问题排查与经验总结即使理解了原理实操中依然会遇到各种问题。下面是我总结的几个高频“坑点”。5.1 错误Error in rowSums(x, na.rm TRUE) : x必需是数值问题描述运行vegdist()时出现此错误。原因排查数据非数值最常见。你的数据框中可能混入了字符型列如分类学信息“k__Bacteria;p__Firmicutes”。用str(your_data)检查每一列的数据类型。数据方向错误vegdist()对每行求和如果某一行全是非数值比如物种名就会报错。解决方案# 确保只将数值部分通常是OTU丰度表传递给vegdist # 假设你的数据框df前7列是分类信息从第8列开始是样本丰度 otu_numeric - df[, 8:ncol(df)] # 或者如果分类信息是行名直接使用整个数据框但需确保全是数字 # 转换所有列为数值型如果确信可以 otu_numeric - apply(otu_numeric, 2, as.numeric) dist - vegdist(otu_numeric, methodbray)5.2 错误距离值异常全为0、NaN或非常大全为0或NaN可能原因1数据所有值都相同或几乎全为0。检查数据summary(as.vector(as.matrix(your_data)))。可能原因2使用了binaryTRUE但数据已经是0/1且样本间物种组成完全相同概率极低。可能原因3数据中有大量NA且na.rmFALSE默认。计算时遇到NA会导致结果为NA。使用sum(is.na(your_data))检查并用na.rmTRUE或事先处理NA如用0填充但需有生物学依据。距离值非常大如欧氏距离根本原因未进行标准化样本总丰度差异巨大。务必先进行标准化如decostand(x, total)或使用对总丰度不敏感的距离如Bray-Curtis。5.3 问题应该选择哪种距离方法这是最常被问到的问题。我的决策流程通常是默认起点对于绝大多数群落丰度数据Bray-Curtis(methodbray) 是第一选择。它平衡了稳健性和解释性。关注物种有无如果你的科学问题更关注物种的分布、存在与否例如研究物种的地理分布界限或者你的数据本质就是二元化的如PCR产物电泳的有无使用二元Jaccard(binaryTRUE或对0/1数据用methodjaccard)。拥有系统发育树对于微生物组数据如果拥有可靠的系统发育树一定要使用UniFrac距离通过phyloseq::distance()或GUniFrac包。它能提供纯组成分析无法揭示的进化维度信息。用于线性模型如果你计划使用基于欧氏距离的线性模型方法如RDA那么对数据做Hellinger转化后计算欧氏距离是一个经典且数学性质良好的选择 (dist(decostand(x, hellinger))注意这里用stats::dist因为vegdist的欧氏距离与dist相同)。敏感性分析在关键研究中可以尝试2-3种不同的距离算法如Bray-Curtis, Jaccard, UniFrac看看主要结论如分组是否分离是否一致。如果结论一致则结果非常稳健。如果不一致则需要深入思考哪种距离更贴合你的生物学问题。5.4 性能与大数据处理当你的样本量很大如1000时计算距离矩阵和后续的置换检验如PERMANOVA会非常耗时。计算加速vegdist()函数本身是用C代码编写的效率已经很高。对于超大规模数据可以考虑使用parallel包进行并行计算如果算法支持。使用专门为大数据设计的包如bigstatsr或BiocParallel但可能需要自定义距离函数。在云计算平台如RStudio Server on AWS上使用更高配置的实例。降低维度在计算距离前可以考虑先过滤掉极低丰度或低出现率的物种如在所有样本中总丰度10或出现样本数5%的OTU这能显著减少数据维度且对整体群落模式影响通常很小。# 示例过滤掉总丰度小于20的OTU otu_filtered - otu_raw[rowSums(otu_raw) 20, ]5.5 我的个人经验与最终建议标准化是习惯不是教条我养成的习惯是对于任何新的群落数据集在计算距离前都会先用decostand(x, total)看一眼相对丰度。这能帮我理解数据的规模。即使用Bray-Curtis我也会对比一下使用原始计数和相对丰度结果的距离矩阵相关性mantel检验确保结论不受影响。距离矩阵的对象类型记住vegdist()返回的是一个dist对象不是矩阵。很多函数如hclust(),adonis2()直接接受dist对象。当你需要提取特定样本对的距离或进行矩阵运算时再用as.matrix()转换。保存中间结果距离矩阵的计算特别是对于大型数据或像UniFrac这样的复杂距离可能很耗时。计算完成后用saveRDS(dist_bray, file my_dist_matrix.rds)将其保存到本地。下次分析时用readRDS()加载可以节省大量时间。理解“零”的含义在群落数据中零可能代表“真不存在”、“存在但未检测到”或“测序深度不足”。对待零值的方式如在Bray-Curtis中忽略双零是生态学指数设计的智慧也是与普通统计距离的根本区别。始终带着“生态学意义”去选择和理解距离。最后没有“唯一正确”的距离。vegdist()提供了一系列工具最好的选择源于你对数据的理解、科学问题的界定以及方法本身的前提假设。从Bray-Curtis开始结合具体问题尝试其他方法并学会用mantel()函数比较不同距离矩阵的相关性你会逐渐培养出对群落距离的直觉让这把“尺子”真正为你所用量出数据背后真实的生物学故事。