做宏基因组或者扩增子分析的朋友到了某个阶段大概率都会碰到同一个需求样本多了OTU表也有了想看看物种和物种之间到底有哪些联动关系哪些物种总是“同进同退”哪些又是“此消彼长”。我最早做这个用的是Rigraph硬写后来又试过Cytoscape插件都折腾得够呛。后来接触到SMETANA这个工具才发现这种“丰度表 - 网络图”的分析其实可以一条命令跑完。这篇东西就围绕SMETANA的安装和使用展开把我从装环境、整理数据、跑通流程、看懂结果到把结果整理成投稿图的全过程记录下来。1. SMETANA是什么它到底帮你解决了什么问题1.1 从一张丰度表到一张物种网络图先说说SMETANA解决的核心问题。你手头有一张特征-样本丰度矩阵行是OTU/ASV或者物种列是样本值通常是测序得到的reads数、相对丰度等。我们想知道的是哪些物种在样本之间呈现一致的丰度变化趋势哪些物种的变化方向恰好相反这种关系一旦被系统性地算出来就可以构成一张网络——节点是物种边是物种间的共丰度关系正相关和负相关还能用不同颜色或线型区分。这个网络不是画着好看的。它背后承载的生态学问题是物种之间的互作可能是什么样哪些物种可能是群落里的关键类群群落结构对外界环境的响应是否集中在某一群物种上SMETANA干的事情就是从输入矩阵出发把相关性计算、显著性检验、网络构建、模块识别这些步骤串成一条流水线最终输出边列表、节点属性表和模块划分结果。它的价值在于把“统计计算”和“网络构建”这两层做了整合你不用再手动拿着相关矩阵到别的软件里二次加工。1.2 和SparCC、CoNet这类工具的差别很多人听到“共丰度网络”第一反应是“我用SparCC或者CoNet不就行了”我一开始也是这么想的但实际用下来发现定位差别挺明显。SparCC是一个相关性计算方法它对丰度数据属于“成分数据”——也就是各物种丰度总和固定导致的伪相关——做了专门校正计算完给的是一个相关矩阵。矩阵后续怎么变成网络、怎么画图、怎么找模块它不管。CoNet是Cytoscape里的插件交互式玩很舒服适合临时看数据但要是你想在服务器上批量处理几十份样本、把网络分析嵌进标准流程里它就不太顺。SMETANA更像一套命令行流水线。它把从矩阵到网络的中间环节都接管了只留少数几个关键参数给你。这种设计的取舍很明确牺牲一部分灵活性换来流程的标准化和可重复性。所以我的建议是探索性分析可以用CoNet这种交互工具正式项目、需要流程化重现的时候用SMETANA这种命令行工具更省心。而且需要说明一点不同算法算出来的网络会有差异这是正常的选定一种工具后尽量别混着比否则结果解释会很混乱。2. 环境准备与安装两条安装路径和几个查错顺序2.1 安装前的环境检查SMETANA是Python写的工具安装前先确认三件事第一Python版本。我用的时候是3.9整体稳定建议至少3.8以上。第二有没有独立的虚拟环境。这点非常关键SMETANA依赖numpy、pandas、scipy、networkx、matplotlib这一类库很多做生物信息的人机器上已经装了一堆包版本很乱直接在base环境里装很容易出现依赖冲突。第三确认pip源可用如果是在内网服务器装最好先配置好可用的pip镜像源。创建独立环境的命令很简单conda create -n smetana_env python3.9 conda activate smetana_env2.2 路径一pip直接安装如果PyPI上有对应包最省事的方式就是pip安装pip install smetana装完之后立刻验证一下smetana --help能弹出参数帮助基本就成功一大半。如果这一步报错No module named smetana要么是安装过程没结束要么是当前激活的Python环境跟pip安装的环境不是同一个。检查which python和which pip是不是指向同一个环境是我每次排查这种问题首先做的事。2.3 路径二GitHub源码安装有时候因为网络或版本原因pip源里搜不到那就用源码安装。流程是把GitHub上的官方仓库clone下来在项目目录里执行安装命令git clone https://github.com/EnvGen/smetana.git cd smetana python setup.py install我后来更习惯用可编辑模式pip install -e .这种模式的好处是如果某个依赖缺失你会立刻看到报错而且修改源码后不用重新安装。缺点是如果你不熟悉Python依赖需要自己逐个把缺的包补齐。常见到的报错无非是某个ModuleNotFoundError缺什么就补什么没有太多技术含量。2.4 安装后自检用--help做冒烟测试安装完成不等于万事大吉。我有一次在服务器上装完命令行能呼出帮助但一跑真实数据就在某一步崩溃。后来发现是安装时缺了一个图像相关的底层库导致出图步骤报错。从那以后我学乖了安装完先做一个“冒烟测试”拿一份最小的人造数据把完整流程跑一遍确认每个输出文件都生成了再放真实数据进去。这一步最多花你五分钟却能避免你把一个坏环境当成数据问题排查一下午。3. 输入数据整理格式要求、BIOM转换和过滤策略3.1 SMETANA的标准输入格式SMETANA的输入核心是一张特征-样本丰度矩阵通常以CSV或TSV形式提供。第一列是特征IDOTU、ASV或物种名第一行是样本ID中间的单元格全部是数值。听起来很简单但实际使用中我见过各种“出格”的输入有带注释行的有列名带空格的有第一列没有表头单元格的这些都会在解析阶段引发各种问题。我的建议是输入文件尽量保持“朴素”第一行单纯是样本ID第一列单纯是特征ID不要有多余的行。一份合格的测试数据长这样OTU_ID,S1,S2,S3,S4,S5 OTU1,120,80,200,150,60 OTU2,30,200,20,80,150 OTU3,0,50,0,20,0 OTU4,90,20,130,0,1803.2 从BIOM和QIIME2导出到可用CSV很多人用的上游流程是QIIME2得到的特征表是BIOM格式。BIOM不能直接被SMETANA读需要先转换。QIIME2的导出命令是qiime tools export --input-path feature-table.biom --output-path exported_table导出目录里会有一个BIOM文件再用biom工具转成TSVbiom convert -i exported_table/feature-table.biom -o otu_table.tsv --to-tsv如果你不想额外装biom包用pandas加biom包也可以import pandas as pd from biom import load_table table load_table(feature-table.biom) df table.to_dataframe(denseTrue) df.to_csv(otu_table.csv)这里我特别提醒一个细节BIOM文件里特征表默认可能是“特征为行、样本为列”的布局转换后记得确认一下方向。SMETANA大多数版本按特征在行、样本在列处理方向反了会把特征ID当成样本ID后面的输出会莫名其妙。3.3 用pandas做低丰度、低出现率过滤输入数据不是越全越好。如果某个OTU只在极少数样本里出现或者总丰度低到可以忽略它对相关矩阵计算不仅没有贡献还会引入噪声导致很多虚假相关或计算错误。我一般用pandas做两层过滤先按“出现率”过滤再按“总丰度”过滤。import pandas as pd df pd.read_csv(otu_table.csv, index_col0) # 至少在20%的样本中出现 min_prevalence int(df.shape[1] * 0.2) df df[df.astype(bool).sum(axis1) min_prevalence] # 总丰度不低于50 df df[df.sum(axis1) 50] # 另外去掉全零行避免零方差问题 df df[df.sum(axis1) 0] df.to_csv(otu_table.filtered.csv)过滤的阈值没有绝对标准要看你的样本量和测序深度。样本多就适当把出现率标准提高一点样本少就要放宽否则剩下的特征太少网络没有意义。我个人的体会是过滤之后特征数量控制在几百到两千以内比较合适既要保证网络有结构又不至于让计算量失控。4. 实操全流程测试数据、核心参数和输出文件4.1 为什么不建议直接拿真实数据试跑工具刚装好大家第一反应都是拿自己的真实数据立刻跑。我劝你先忍住。真实数据里往往藏着各种坑——空值、异常值、列名不规范、某些OTU丰度分布极差。如果直接跑报错了你根本分不清是安装问题、格式问题还是数据本身的问题。先手动造一份十几行的小表跑一遍全流程你就能对SMETANA每个阶段做什么有个底。之后再上真实数据遇到问题就知道该往哪个方向排查。再造一次测试数据就用前面那份CSV就行保存成test_table.csv。4.2 一条最简命令的完整解释进入测试数据所在目录后我通常这样跑smetana -i test_table.csv -o smetana_out -p 0.05 --threads 4这几个参数是我用的所有版本里都存在的核心项-i输入丰度表路径这是必填项。-o输出目录如果目录不存在会自动新建。-p显著性阈值默认是0.05用来过滤掉那些不够显著的边。--threads并行计算线程数按机器CPU核心数设现阶段实验性使用4个线程就够了。需要说明的是不同版本对参数命名可能会调整最准确的参数列表一定要看你装的那个版本的smetana --help输出。这个习惯要养成不要死记我文章里的命令。跑完后终端会打出一些日志信息显示每个阶段的进度。如果只是跑一份小数据可能几秒就结束了。真正的大数据最耗时的是两两相关性的计算两千个特征就要跑一阵子。4.3 输出文件的正确打开方式跑完之后打开输出目录你可能会看到类似这样的文件边列表文件比如network.tsv或edges.csv、节点属性文件比如node_metrics.csv、模块文件比如modules.csv以及一个快速预览用的网络图文件。每个文件的命名和格式在不同版本略有差别但核心内容是一致的。边列表大概长这样node1,node2,rho,p_value OTU1,OTU2,0.88,0.001 OTU4,OTU5,-0.75,0.012每一行代表一条经过筛选后保留的边。看到这里恭喜你工具已经完全跑通了。接下来要做的才是真正需要生物信息学和生态学知识的部分怎么解释这些输出。5. 结果解读与可视化从边列表到生态故事5.1 边列表字段与相关性正负的含义边列表是整个网络分析的核心。每一行两个节点代表一条边rho列是相关系数p_value列则告诉你这条相关性是否显著。正相关意味着两个物种在样本间变化趋势趋同可能是互利共生、共享生态位或者对同一环境因素有相似响应负相关则意味着此消彼长可能是竞争资源、偏害作用或者对某环境因子有相反偏好。但我要泼一盆冷水网络分析得到的相关性只是假说生成工具不是因果证据。你发现A和B正相关不能直接写“A和B存在互作”只能说“A和B的丰度存在显著共变化推测可能存在生态联系”。这种表述在文章里才是稳妥的审稿人也挑不出毛病。5.2 节点属性表里哪些指标值得关注SMETANA输出的节点属性表一般包含多个拓扑指标其中三个最值得关注degree度节点直接相连的邻居数量。度高的节点常被称为hub在网络里地位显眼。如果hub物种丰度发生变化整个网络的结构都可能受影响。betweenness centrality介数中心性衡量节点承担“桥梁”角色的程度。介数高的节点连接不同群落模块是物种互作链上容易成为“堵点”的位置。clustering coefficient聚类系数反映节点的邻居之间互相连接的紧密程度。聚类系数高的区域通常被认为是比较稳定的功能群。实际分析里我不建议把全部指标堆到论文里而是结合研究问题选择一两个。比如你想找“核心物种”就重点看degree你想找“模块间的桥接物种”就看betweenness。把这些高指标节点挑出来后再回到原始丰度表看它们的分布特征往往能给出更有价值的线索。5.3 模块划分怎么解释模块是网络中连接紧密的子图模块内的节点彼此之间的相关性比较强模块与模块之间的联系则相对稀疏。生态学上一个模块往往被解释为一群栖息偏好相近或互作密切的物种组合。SMETANA输出的模块文件会告诉你每个特征归属于哪个模块。解读模块信息时我的做法是把每个模块当成一个整体计算它在各个样本里的平均丰度再去和环境因子做相关性分析。如果发现某个模块在高污染样本里丰度显著偏高那这个模块里的物种就值得做进一步的分类学注释和功能预测它们很可能代表了一组对特定环境条件有相似响应的类群。这种“先找模块再找模块-环境关联”的分析路径比单独盯着一两个OTU要稳得多。5.4 用Cytoscape快速绘制投稿级网络图SMETANA自带的快速预览图适合自查投稿最好还是用Cytoscape重画。流程不复杂打开Cytoscape选择File - Import - Table from File导入边列表文件。确认第一行是表头Source/Target列对应两个节点。导入后使用Layout菜单里的Prefuse Force Directed Layout或yFiles Organic Layout重新排布让网络结构更清楚。在Style面板里把边的颜色映射到rho值正相关一个颜色负相关另一个颜色边的粗细映射到相关系数绝对值。把节点大小映射到degree或betweenness突出核心节点。这一步操作本身不涉及分析但对结果表达影响很大。一张颜色区分正负相关、节点大小体现重要性的网络图比通篇纯文本描述强得多。6. 实际运行中容易踩的坑和排查思路6.1 UnicodeDecodeError文件编码问题症状输入CSV在命令行运行时抛UnicodeDecodeError。最常见场景是文件在Windows上用Excel编辑过保存编码变成了GBK或带BOM的格式。排查时先用文本编辑器或file命令看编码然后统一转成UTF-8。# 查看编码 file otu_table.csv # 转换编码 iconv -f GBK -t UTF-8 otu_table.csv otu_table_utf8.csv这个坑看着小实际发生率不低。我的习惯是所有进入分析流程的CSV文件先统一用UTF-8无BOM格式保存能省掉后面一大堆莫名其妙的报错。6.2 零方差行导致的相关计算崩溃症状计算相关矩阵时程序报错提示输入包含常量值或结果出现大量NaN。原因通常是输入数据里有全零行或者某个特征在所有样本里都是一个固定值。相关性计算遇到零方差的数据时数学上是未定义的程序自然跑不下去。排查过程用pandas检查每一行的标准差找出恒为0的行。df[df.std(axis1) 0]解决方式很简单过滤掉这些行。这一步最好在进入SMETANA之前就做完我在前面第3节给的过滤代码里已经包含了。记住零方差行不是“数据”是“错误输入”留着只会让程序崩溃。6.3 依赖版本冲突症状源码安装后运行import smetana或者pip install -e .时报错提示某个numpy或pandas的API找不到。这通常是环境中某个依赖版本过新/过旧导致的。比如太新的numpy移除了旧API太新的pandas改变了某些行为。排查思路先看报错栈里提到的是哪个库再查看当前版本pip list | grep -E numpy|pandas|scipy|networkx解决在虚拟环境里重新装一个兼容版本组合。不同SMETANA版本要求的依赖组合不一样强烈建议装一个全新环境不要图省事直接覆写旧环境。这也是为什么我一直强调虚拟环境的原因——就算装坏了直接删掉重建就行不影响其他项目。6.4 上千OTU时的内存问题症状输入特征超过几千个时计算过程中内存占用暴涨最后被系统OOM杀掉。根本原因是两两相关计算的空间复杂度是O(n^2)特征数翻一倍内存占用可能涨四倍。解决思路第一把特征数通过过滤降到合理范围第二如果必须要保留大矩阵就减少线程数有些计算库在线程多时内存开销更大第三实在不行在服务器上分析别用自己笔记本硬扛。按我的经验几千个特征是比较好跑的区间上万特征就需要认真考虑计算资源了。6.5 结果全部是NaN或边列表为空症状程序正常跑完但边列表是空的或者节点属性表里大量NaN。发生这种问题第一反应看输入数据本身是不是稀疏度太高是不是绝大多数样本的绝大多数特征是0相关性计算的统计功效不够会导致所有相关性检验都不显著最后一条边都不保留。排查时我会先做一个最直接检查统计输入矩阵的非零比例。如果非零比例低于5%说明数据太稀疏建议加大数据过滤力度或者考虑先做数据聚合比如把OTU按属水平合并把稀疏矩阵变得更稠密。这一步比在那里调p值阈值有意义得多。7. 进阶把SMETANA放进你的自动化分析流程7.1 在Jupyter里通过Python API调用命令行适合一键跑但如果你想在Notebook里做数据过滤、网络分析和可视化联动可以尝试调用它的Python API。一般逻辑是import smetana后调用相应函数具体的函数名和参数在不同版本会有变化所以我不在这里写死。在Notebook里的优势是能立刻看到中间结果比如过滤后剩多少特征、边列表的前几行长什么样对调试和教学都友好。不过要提醒一点如果最终要在服务器上跑正式分析我还是建议回到命令行。命令行日志清晰出错时定位快也方便写进流程文件。Notebook更适合小规模探索。7.2 参数记录与重跑的可复现性做网络分析最怕的是“跑完忘了参数”。过滤阈值、显著性阈值、相关性方法任何一个参数变了结果都会变。文章里如果说“用的是SMETANA默认参数”其实是没有意义的因为默认值也会随版本变化。我的做法是每次分析都记录一份分析配置哪怕就写成一行注释# SMETANA run 2025-XX-XX # input: otu_table.filtered.csv (344 features) # filter: prevalence20%, total50 # smetana -i otu_table.filtered.csv -o output -p 0.05 --threads 4 smetana -i otu_table.filtered.csv -o output -p 0.05 --threads 4有条件的话把这个命令写进Snakemake或Nextflow流程这样所有参数都固化在文件里以后换数据重跑或者跑完忘了细节回看配置文件就一目了然。这种习惯在项目周期长、分析任务多的场景下特别有用能省掉很多不必要的返工和扯皮。按我自己的经验SMETANA这类工具最让人舒服的地方在于它把“零散的中间步骤”收拢成了一条流水线安装和数据格式熟练之后从丰度表到网络图通常就是一条命令的事。真正需要花心思的反而是输入数据的清洗、网络结果的生态学解释这两头。建议刚开始接触的同学务必花半小时看一遍自身版本smetana --help的输出再找一份几MB级别的小数据完整跑一遍把每个输出文件打开看一遍你就能建立很直观的“数据到结果”的映射感。后面再换复杂数据心里就有底了。