引言
这年头做微生物组研究,离了扩增子测序几乎是寸步难行。16S和ITS这两个经典的标记基因,在过去十几年里撑起了肠道、土壤、水体、植物根际等无数微生态研究方向的基本盘。但恰恰是这个“基本盘”,这几年正在经历一轮非常明显的技术换挡:一边是二代测序凭借高深度、低成本继续守着自己的领地,另一边是三代入读长测序带着全长16S/ITS的超高分辨率和注释精度冲进来抢地盘。2026年4月份的这场“微生物组-扩增子二、三代16S/ITS分析和可视化技术研讨会”,把两代测序的数据处理分析、下游可视化和实际应用场景放到一起讲,说穿了就是对这套正在剧烈变化的方法学体系做一次系统性的梳理。
这篇文章我打算结合研讨会的核心议题和自己实际做项目的经验,把扩增子分析从实验设计、数据质控、降噪聚类、物种注释到可视化的完整链路拆开讲一遍。两代测序之间的差异在哪里、分析流程上哪些环节必须分开处理、哪些图表能讲清楚故事、哪些坑是新手几乎必踩的,都会逐一交代清楚。无论你是刚进实验室的研究生,还是已经被ASV表和多样性指数折磨过几轮的从业者,这篇文章应该能帮你在自己的项目里省下不少试错时间。
1. 扩增子研究的设计思路与核心逻辑
1.1 16S与ITS:两个标记基因,两套分析逻辑
很多刚开始接触微生物组的人,容易把16S和ITS混为一谈。实际上,这两个标记基因的应用场景、引物选择、数据库以及后续分析的重点方向都有不小的差异。
16S rRNA基因编码的是原核生物核糖体小亚基RNA,全长约1542bp,包含九个高变区(V1-V9)。它的优势在于保守区和高变区交错排列——保守区可以用来设计通用引物,高变区的序列差异则提供了物种区分的基础。二代平台通常选择V3-V4或V4区段,读长在300bp左右就能获得不错的区分度;三代平台得益于长读长,可以覆盖完整的16S全长,直接拿到接近于种水平的分类分辨率。
ITS(Internal Transcribed Spacer)则是真核生物核糖体DNA上位于18S和28S之间的内转录间隔区,细分为ITS1和ITS2两个亚区。它的特点是进化速率比16S快得多,序列变异丰富,尤其适合真菌的物种鉴定。真菌ITS数据库(如UNITE)在分类学注释上的标准化程度也比较高,和16S数据库(如Greengenes、SILVA、RDP)的体系有明显区别。
实操中一个经常被忽略的细节是:16S分析中很多通用工具链的默认参数是针对细菌设计的,如果做的是古菌或者线粒体、叶绿体污染比较严重的样本(比如植物根系、叶片表面),必须在质控阶段或者注释之后做专门的过滤。ITS分析则要额外留意真菌ITS区的长度变异问题,ITS1和ITS2的长度在不同的分类群之间差异很大,这会直接影响聚类和降噪的参数选择。
1.2 二代与三代测序:读长、通量与成本三方平衡
二代测序在扩增子领域的主流方案是Illumina MiSeq和NovaSeq,PE250/PE300模式产生的双端读长可以覆盖V3-V4和ITS区间。它的核心优势是单样本测序深度高、价格相对便宜,一个run可以pooling大量样本,适合大规模人群队列或田间试验的筛查研究。
三代测序则以PacBio的HiFi测序和Nanopore的长读长测序为代表。PacBio Sequel平台通过CCS模式把单分子测序的多次遍历聚合成高准确率的HiFi reads,读长通常在1.5-2.5kb之间,刚好完整覆盖16S或ITS全长。Nanopore的优势则是测序仪便携、实时产出数据,缺点是单碱基准确率相比HiFi还是低一些,但当前的最新版本在化学和碱基识别算法改进后已经有了明显提升。
我的经验是:如果项目目标是小样本量、需要精确到种水平的群落组成,或者研究对象的16S和ITS基因存在显著的拷贝数异质性,三代全长测序的优势非常明显。如果是几百上千例的大规模比较研究,预算有限,二代的V3-V4方案依然是性价比之王。两者之间的关系不是替代,而是互补。研讨会把两代测序放在同一个背景下讲,本质上是希望从业者在设计实验时不要被单一技术路线锁死。
1.3 实验设计中的几个前置决定
分析做得再漂亮,实验设计阶段图省事,后续任何质控都救不回来。扩增子的实验设计有几个前置决定会影响全套分析流程:
第一是引物选择。细菌16S通用引物中,515F/806R(针对V4)、341F/785R(针对V3-V4)都是经典选择。真菌ITS分析中,ITS1F/ITS2和ITS3/ITS4分别偏向ITS1和ITS2区域。引物的选择必须和目标数据库配套,不然注释阶段会出现大量unclassified。第二是样本量的确定。扩增子测序的生物学重复建议至少6个,如果群落组间差异较小,最好10个以上;测序深度方面,细菌16S V3-V4建议每个样本至少3-5万条有效reads,真菌ITS由于基因组中rRNA重复单元数量差异大,测序深度要求可以适当放宽但也不要低于2万。第三是建库策略。二代的双Index方案、三代的Barcode方案都要注意避免索引错配,样本处理过程中的污染控制对低生物量样本尤其关键。
这几个决定在做研讨会内容拆解时被反复强调,不是没有原因的。太多项目最后卡在分析出的结果解释不了,往前追溯都是实验设计埋下的隐患。
2. 核心分析环节解析与实操要点
2.1 原始数据质控:两代测序的差异化处理
扩增子分析的第一步永远是质控,但二代和三代数据的质控思路差异很大。
二代双端数据在导入QIIME 2之前,第一步是用cutadapt或QIIME 2内置插件去除引物和接头序列。这一步很多人会跳过或者图省事,觉得反正后面有质量控制。但实际上引物没去除干净,会直接影响下游的denoising,因为DADA2这类算法会把引物区域错误的碱基当作真实变异。参数上,cutadapt的--error-rate一般设置在0.1左右,如果引物长度较短或者有简并碱基,可以适当放宽。PEAR或Vsearch的paired-end合并阶段,最小overlap建议30bp,允许的错配率控制在5%以内。
三代数据质控的核心是HiFi read的生成环节。PacBio平台需要先通过CCS软件把subreads聚合成HiFi reads,常用的参数是--hifi-kinetics -j 64 --min-rq 0.99,意思是只保留准确率不低于99%的reads。经过这个步骤后,再进行引物和Barcode的去除。Nanopore数据则用porechop或chopper做adaptor trimming和长度过滤,之后还要根据序列质量进行二次筛选。三代数据的过滤阈值一般建议长度范围设为1200-1800bp(16S全长)或200-1200bp(ITS,因为不同真菌ITS区域长度波动大),这一条是经验值,不同引物和不同生态类型可以微调。
2.2 降噪与聚类:ASV还是OTU
这个选择是扩增子分析领域过去几年最核心的方法论争论。传统的OTU聚类基于97%序列相似性阈值,优点是计算开销小、跨研究可比性强,缺点是把真实的生物学变异和测序误差混在一起,无法区分菌株水平的细微差异。DADA2为代表的denoising算法则通过误差模型校正测序错误,输出单核苷酸分辨率的ASV表,分辨率更高,但也会因为测序错误或罕见变异产生大量低丰度ASV,导致样本间菌群组成表观上看起来更加分散。
实际项目中,我建议的做法是两条腿走路:常规分析以ASV为主要单位,同时额外产出一份按97%聚类生成的OTU表。下游分析中alpha多样性和beta多样性可以分别跑一遍,对照结果。如果两种策略得出的群落差异结论不一致,通常是样本量不足或组间差异太弱,需要在实验设计上补样本,而不是继续在参数里找答案。另一方面,三代全长数据分析更推荐直接走ASV路线。因为全长序列的分辨率本来就高于V3-V4区段,再按97%聚类反而会损失大量物种信息。
QIIME 2的DADA2插件是当前做denoising最常见的选择,几个关键参数值得反复调:--p-trunc-len需要根据测序质量曲线设置,一般V3-V4可以设为260-280bp;--p-max-ee建议2.0以内,可以严格一点到1.0;--p-trim-left设置为引物长度即可。三代数据用q2-dada2也可以处理,但需要先用q2-cutadapt完成引物去除,并在--p-trunc-len这里设置一个较大的值(比如1500)表示默认不截断。对于大量三代数据,QIIME 2自带的q2-quality-filter配合q2-deblur也是一个备选方案,但对输入质量要求较高。
2.3 物种注释与数据库选择
物种注释阶段最影响结果解读的变量是数据库。16S分析的常用数据库有Greengenes 13.8、SILVA 138和RDP 11.5。Greengenes的注释体系在早期研究中使用广泛,兼容性最好,但近几年更新停滞;SILVA的覆盖范围更广、注释分类层级更完整,是目前的主流选择;RDP有独立的RDP Classifier算法,适合做快速分类。真菌ITS注释基本默认使用UNITE数据库,目前推荐使用UNITE 9.0版,它提供了动态聚类(dynamic)版本的分类单元,对未知真菌的处理更精细。
NCBI的新版16S数据库(如RefSeq Targeted Loci)准确率最高,但覆盖率相对较低,适合做精细鉴定复核而非全样本注释。QIIME 2中使用q2-feature-classifier进行分类注释,先训练classifier再对特征序列进行预测,分类置信度阈值建议设置在0.7-0.8之间。注意classifier的训练集必须和引物区域匹配,用V3-V4引物扩增的序列却用V4训练的classifier去注释,结果必然出现偏差。
2.4 多样性分析:从指数到比较的可视化拼接
Alpha多样性(组内多样性)常用指标包括Observed Features、Chao1、Shannon、Simpson、Faith's PD等。Observed Features对测序深度敏感,Chao1是对物种丰富度的估计值,Shannon和Simpson综合反映丰富度和均匀度。QIIME 2中执行core-metrics-phylogenetic时会一次性输出多个核心指标,包括Beta多样性距离矩阵(Bray-Curtis、Jaccard、Unweighted UniFrac、Weighted UniFrac)。选择哪个Beta多样性指标要结合研究问题:UniFrac系列考虑了物种间的系统发育距离,适合宏进化层面的群落比较;Bray-Curtis则从丰度差异入手,更贴近群落组成的直接变化。
实际报告里我会同时展示Unweighted和Weighted UniFrac的结果,前者反映稀有物种的组成差异,后者受优势物种影响更大。两个指标结果方向一致说明群落差异稳健;如果不一致,则需要在讨论中进行解释,不能回避。差异显著性检验中,PERMANOVA(Adonis)是绝对主流,QIIME 2里--p-method permanova即可完成;但PERMANOVA对组间离散度敏感,必要时补充执行betadisper检验离散度是否组间存在差异。
3. 可视化实战:把微生物组数据讲明白
3.1 可视化工具链的选择
扩增子分析的下游可视化,我个人的主力方案是R语言生态系统,主要是phyloseq和ggplot2这套组合。phyloseq可以非常方便地导入QIIME 2的产物——通过qza_to_phyloseq函数把特征表、元数据、分类注释和系统发育树一次性整合成phyloseq对象,之后的多样性计算和图表绘制都在这个框架下完成。microeco这个国产包也是一个很好用的选择,它对中文资料支持更好,教程系统化,适合刚入门的用户。
另一条路线是QIIME 2 View(https://view.qiime2.org)直接在线可视化,零代码操作,适合快速查看单个产物。但如果要做出版级别的图,还是需要用R重绘。Python生态里scikit-bio和q2-diversity的可视化能力相对偏底层,除非有特殊需求,一般不用做主力可视化工具。
3.2 核心图表类型与绘制要点
柱状图是最常见的门水平或属水平丰度展示方式。要点在于分组顺序、颜色方案和“Others”类别的合并阈值。一般把相对丰度低于1%的物种归入Others,避免图面过于杂乱。如果需要展示每个样本的具体组成,堆叠柱状图很好用;如果要展示组间平均组成,则建议用分组均值加误差线的形式。颜色方案推荐使用RColorBrewer的Set3或Paired,如果是色盲友好场景可以使用viridis。
热图一般用来展示物种丰度在样本或分组间的聚类关系。数据预处理上,建议先对丰度做log2(x+1)转换或者Z-score标准化,避免高丰度物种主导颜色映射。行和列都做层次聚类能帮助快速发现样本和物种的共分组模式。要注意展示的物种数量不能太多,一般选取丰度最高的30-50个属,否则热图可读性会严重下降。
PCoA和NMDS图是Beta多样性可视化中最常用的两种排序图。PCoA基于距离矩阵的特征值分解,NMDS基于秩次迭代算法,两者在解释率上有差异:PCoA的轴可以解释方差比例,NMDS则用一个Stress值评价降维的可靠性。绘制时建议同时标注样本点、组别置信椭圆(conf.int=0.95)和统计检验结果。如果要展示的系统发育信号更强,可以用PCoA对UniFrac距离做排序。
还有一个经常被忽略但很实用的图是稀释曲线。它直观反映测序深度是否足够。R里rarecurve函数就可以绘制,QIIME 2中可以通过alpha-rarefaction可视化。如果曲线没有达到平台期,说明测序深度不足,后续所有Alpha多样性比较都可能存在偏差。
3.3 交互式可视化的延伸方向
出版图表通常是静态的,但项目内部探索阶段,交互式可视化能明显提高效率。我常用的方案有几种:一是用plotly包把PCoA图做成可旋转的三维或二维交互图,鼠标悬停时查看样本信息和分组;二是用shiny搭建一个简易的在线报告,把核心图表和分析结果整合成可筛选的网页;三是用gganimate制作排序图随时间或分级变量变化的动态演示。
从研讨会的反馈看,越来越多的研究团队开始把交互式图表集成到补充材料或项目网页中,审稿人和读者对这种形式普遍接受度较高。不过要注意,交互式HTML文件通常比较大,上传期刊系统或者作为补充材料时要确认文件大小限制,必要时对数据进行抽稀。
4. 常见问题与排查技巧实录
4.1 三代测序数据量虚高但有效数据少
这是我做三代扩增子项目时踩过最大的坑。PacBio平台下机的原始subreads数据量看着非常大,但经过CCS聚合成HiFi reads后,有效数据可能只剩下20%-30%。这是因为CCS聚合会过滤掉聚合次数不够、准确率不达标的reads。因此三代项目的测序量设计必须预留足够的冗余,一般建议按最终需要的HiFi reads数量的3-5倍进行上机。Nanopore数据则要留意测序过程中reads长度的真实分布,有些短片段reads虽然数量很多但对全长扩增子分析没有价值,必须在长度过滤阶段清零。
4.2 ITS分析中的unclassified比例异常偏高
ITS分析中unclassified序列比例高,首先检查Primer和数据库是否匹配。UNITE数据库的注释格式和Greengenes/SILVA不一样,如果下载动态聚类版本的UNITE,注释结果里会有大量未命名的真菌类群,这其实是正常现象,毕竟真菌中大量物种尚未被培养和描述。其次检查引物去除是否干净,残留的引物序列会导致DADA2把该区域识别为高错误率区域,进而产生偏差。还有一个容易被忽略的点:宿主污染。植物样本的ITS扩增中,叶绿体18S和线粒体区域的非特异性扩增可能导致无参考注释的reads占比居高不下,建议扩增前使用特异性引物阻断或扩增后进行分类学过滤。
4.3 多样性分析结果与预期不符
组间Beta多样性差异不显著,最可能是组内个体差异过大,生物学重复数量不足。QIIME 2的beta-group-significance可以逐组检验,但即使不显著,也不能草率下结论——先检查测序深度、质控后的reads数是否在样本间严重不均衡,这会影响距离计算。其次检查是否有多度异常偏高的优势物种,当某个物种相对丰度超过80%时,Bray-Curtis距离会被这个物种主导,UniFrac距离则会因为物种间的系统发育位置而产生扭曲。这时建议做一次去优势种(de novo)的敏感性分析,或者同时查看Jaccard与Unweighted UniFrac结果。
4.4 关于计算资源与流程复现的避坑建议
扩增子分析虽然不像宏基因组那么吃资源,但在DADA2降噪环节,三代数据的计算压力还是比较大的。建议用带有16核以上CPU、64GB内存的工作站或者云服务器跑三代数据,二代数据16GB内存也能勉强运行。所有工具链务必使用conda或docker管理环境,QIIME 2不同插件版本之间兼容性问题不少,建议锁定一个稳定的QIIME 2发行版本,不要在生产环境里频繁升级。分析记录建议用Nextflow或Snakemake这样的流程管理工具,便于追溯每一步参数,保证研究成果的可复现性,这也是当前很多期刊的隐性要求。
4.5 一个典型的完整分析命令序列参考
这里给出一个我最常用的QIIME 2分析模板,针对二代PE300的16S V3-V4扩增子数据:
# 导入数据(使用Fastq manifest格式) qiime tools import \ --type 'SampleData[PairedEndSequencesWithQuality]' \ --input-path manifest.tsv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33V2 # 质控与可视化 qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv # 去除引物 qiime cutadapt trim-paired \ --i-demultiplexed-sequences demux.qza \ --p-front-f CCTACGGGNGGCWGCAG \ --p-front-r GACTACHVGGGTATCTAATCC \ --p-error-rate 0.1 \ --o-trimmed-sequences trimmed.qza \ --verbose # 降噪(DADA2) qiime dada2 denoise-paired \ --i-demultiplexed-sequences trimmed.qza \ --p-trim-left-f 0 --p-trim-left-r 0 \ --p-trunc-len-f 280 --p-trunc-len-r 220 \ --p-max-ee-f 2 --p-max-ee-r 2 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats stats.qza # 物种注释(以SILVA 138为例) qiime feature-classifier classify-sklearn \ --i-classifier silva-138-99-nb-classifier.qza \ --i-reads rep-seqs.qza \ --p-confidence 0.7 \ --o-classification taxonomy.qza # 多样性分析 qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 30000 \ --m-metadata-file metadata.tsv \ --output-dir core-metrics-results这套命令的关键参数值得复盘一下:--p-trunc-len的数值是依据demux.qzv里的质量曲线确定的,如果300bp读段在290bp之后质量还是很高,可以适当延长;如果F/R两端质量不对称,就要分开设置,不能图省事用同一个值。--p-sampling-depth(抽样深度)不能设置得比最低样本的feature总数还高,否则样本会被抽空,一般选择样本reads数中位数的90%。
5. 从数据分析到生物学故事:研讨会上大家最关心的几个延伸问题
5.1 三代数据时代还需要做OTU聚类吗
在完整16S和ITS序列基础上,很多三代数据分析流程默认直接生成ASV表,再把ASV按物种分类名称合并到属或种水平进行下游分析。那传统的97% OTU聚类还有必要吗?我个人的看法是:如果目标是跨研究比较,或者要和过去发表的大量基于V3-V4的历史数据对接,OTU聚类仍有价值。但如果项目只关注当前数据集的组成和差异,三代全长ASV的分析粒度已经足够,做OTU聚类反而会在“同一物种内不同菌株”这个层面抹掉真实差异。
5.2 多组学整合分析中的扩增子定位
扩增子数据提供的只是微生物群落的“物种名录”和“相对丰度画像”,它不直接告诉我们这些微生物在功能上做了什么。现在越来越多的项目在做扩增子+宏基因组的联合分析:扩增子负责回答“谁在”,宏基因组负责回答“能干什么”。代谢物检测、宿主转录组等信息与菌群组成整合后,相关网络分析可以找到核心物种和关键代谢通路。视觉呈现上,这类整合分析常用网络图、和弦图和具有通路层级的圈图,R包如ggraph、circlize是主力工具。研讨会现场讨论最热烈的话题之一,就是如何避免多组学关联分析中过度解释的陷阱——相关不等于因果,扩增子得出的“标志物种”必须用功能层面的证据去验证。
5.3 机器学习和扩增子数据的组合
基于扩增子数据做随机森林分类或者构建疾病预测模型,这几年论文里很常见,实操上也有一些值得注意的地方。16S数据的稀疏性、高维度和成分性(compositionality)特征,决定了常规机器学习流程里很多默认参数并不适用。建议使用vegan包中专门的变换方法(如CLR变换)先对特征表做处理,再进行特征筛选和模型构建。训练集和测试集的划分必须在样本层面而不是ASV层面进行,否则会严重高估模型性能。交叉验证策略和置换检验是扩增子机器学习模型评估中的标配,汇报模型结果时一定要包含这些指标。
写在最后:关于这场2026年4月研讨会的个人感受
这次研讨会筹备期间,我回头把自己过去五年做过的扩增子项目重新梳理了一遍,最大的感触是:分析工具在不断更新,但数据分析的底层逻辑其实没有变过——问清楚每个样本代表什么,控制住每个环节的误差来源,选择合适的粒度回答自己的生物学问题。二代和三代测序的博弈,其实给了研究者更多选择空间,而不是替你做选择。QIIME 2、phyloseq、R语言生态这些工具本身就代表了一套非常成熟的方法论,真正拉开项目差距的,往往是对数据的理解和细节的处理,而不是工具的先进程度。
最后再分享一个我自己的小习惯:每次分析到一个关键节点,比如完成质控、生成ASV表、拿到PCoA结果时,都养成把图表和主要参数保存到一个时间戳文件夹的习惯。这个习惯在项目复盘和论文返稿时能节省大量时间。扩增子分析最怕的不是踩坑,而是踩过坑之后没有留下脚印。研讨会结束后,如果你有任何具体的分析问题,欢迎对照自己项目的数据慢慢调参、多跑几遍试试看,真实的经验永远是在自己的数据上磨出来的。