别人拿到测序公司交付的fastq数据包解压后几百个.gz文件第一反应通常是找教程从导入数据开始跑。但QIIME2和传统分析工具最大的不同是它不直接吃fastq你得先理解它的artifact机制否则会在导入这一步反复报错连数据都进不了门。这篇博文我按一条龙流程完整走一遍——导入数据、质控去噪、多样性分析、分类学注释、差异丰度与功能预测——每一步我都会给出命令、讲解参数取舍顺便把容易坑人的细节一次性交代清楚。博文正文1. 开工之前先理解QIIME2的数据组织和运行环境1.1 artifact不是格式洁癖是为了数据追溯先不谈命令我把话放前面QIIME2发音是“chime two”之所以要引入.qza这种封装格式不是故意增加学习门槛而是要解决生物信息项目里最让人头疼的“溯源问题”。举个实际例子我把一个table.qza发给合作者里面嵌入了完整的产生历史包含上游去了多少噪音、用了哪些参数、输入文件是什么谁都可以沿着provenance数据谱系往前追溯。这一点在写论文Methods部分时特别爽——你不需要单独维护一张Excel记录每一步参数qiime tools界面里直接能查看完整操作日志。理解了这层逻辑开篇那些看似繁琐的概念就顺了.qza是“分析工件”的容器.qzv是可视化的容器两者类似一个带说明书的数据包。每次操作都会生成新的对象而你之前的原始数据始终不会被覆盖中间每一步都有迹可循。所以这里的建议是给项目建一个清晰的目录结构建议按input、scripts、output分成三层避免后续分析时被一堆.qza淹没。1.2 环境准备conda、版本选择与资源预判QIIME2官方推荐用Miniconda安装避免污染系统Python环境这一点我强烈建议照做。下面是安装的大致流程wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh conda update conda conda create -n qiime2-2023.5 python3.8 conda activate qiime2-2023.5这里的版本号重要。每个版本对应的Python版本和插件版本捆绑较死我不建议你随便用pip去补装插件尽量整体安装官方release。以我常用的2023.5版本为例整个环境装完约占用几个GB空间跑dada2去噪时内存消耗看样本量常规几十个样本8-16GB内存也能跑但如果数据量上了百级PE reads建议准备32GB以上并配合多核CPU用。conda install -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/qiime2/ qiime2如果在国内环境把官方conda源替换成清华镜像会省很多时间但要注意镜像列表里必须包含qiime2这个channel。另外一定记得确认环境中没有同时存在另一个Python版的q2cli否则导入命令会提示找不到可执行文件。1.3 拿到测序数据后的第一件事摸清R1/R2与目录结构Illumina双端测序下机后每个样品会有R1和R2两个fastq.gz文件R1是正向读长R2是反向读长二者中间是插入片段。很多第一次做16S的人容易犯一个错误直接把两个文件单独导入当作单端数据处理这样会把双端信息割裂后续dada2拼接就没有意义。正确思路是一对R1/R2在导入后代表同一个样本的两个方向样本数量应该是文件对数而不是文件个数。先解压并看一眼目录ls rawdata/ | head正常情况你会看到Sample001_R1.fastq.gz、Sample001_R2.fastq.gz这种成对文件第三方公司还常带barcode和adapter信息。这里务必先把样本名统一只保留字母、数字和下划线。如果样本名下划线太多后续分组信息会写得很痛苦。我最开始跑项目时公司文件叫A-1_S1_L001_R1_001.fastq.gz光是把名字改成A1.R1.fastq.gz就折腾了半天。2. 数据导入manifest文件、常见报错与验证2.1 两种导入姿势Casava目录格式与manifest清单QIIME2支持多种导入格式最常见的是CasavaOneEightSingleLanePerSampleDirFmt和manifest。Casava格式适合测序公司交付的原始目录还保留着标准命名的场景——目录层级固定为SampleID_SampleNumber_L001_R1_001.fastq.gz此时可以直接qiime tools import \ --type CasavaOneEightSingleLanePerSampleDirFmt \ --input-path rawdata \ --output-path demux.qza但我遇到更多的情况是公司交付的文件名已经改成了自己的命名习惯或样本来自多个批次。这时老老实实用manifest清单。先在文本编辑器里生成一个manifest.tsv三列sample-id、forward-absolute-filepath、reverse-absolute-filepath。sample-id forward-absolute-filepath reverse-absolute-filepath A1 /absolute/path/rawdata/A1_R1.fastq.gz /absolute/path/rawdata/A1_R2.fastq.gz A2 /absolute/path/rawdata/A2_R1.fastq.gz /absolute/path/rawdata/A2_R2.fastq.gz B1 /absolute/path/rawdata/B1_R1.fastq.gz /absolute/path/rawdata/B1_R2.fastq.gz然后导入qiime tools import \ --type SampleData[PairedEndSequencesWithQuality] \ --input-path manifest.tsv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33V2这里有个细节需要说明Phred33和Phred64两种格式容易混淆。现在主流Illumina下机数据基本都是Phred33默认选PairedEndFastqManifestPhred33V2不会错。如果分析老数据不确定时用seqkit head检查质量值字符范围ASCII码43以下对应Phred3343以上则是Phred64再对应修改格式。2.2 manifest文件的隐形坑绝对路径、制表符和CRLF很多人在导入时报File not found或解析错误多半不是数据问题而是路径问题。manifest里必须写绝对路径而且要保证路径里的目录和文件真的有读权限。服务器环境下如果是从Windows下编辑的manifest记事本默认换行符是CRLF而Linux下解析时会残留\r字符导致样本ID后面多个不可见符号。解决办法是sed -i s/\r$// manifest.tsv还有一种坑是字段分隔符到底用tab还是逗号。PairedEndFastqManifestPhred33V2严格识别TSV制表符分隔而旧的PairedEndFastqManifestPhred33虽也是TSV但如果从Excel复制粘贴会把tab换成空格导致解析字段数不对。我的习惯是用VSCode或Notepad先开启“显示所有字符”来确认分隔符再上传服务器。2.3 导入后第一件事查看每个样本的reads数导入完成不代表数据好用你还需要马上验证qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv然后把demux.qzv拖到QIIME2 View网页里查看。这里主要看两个信息一是每个样本的序列条数判断有没有样本明显偏少二是质量分数图。如果看到某个样本reads数只有几千而其他样本十几万尽早记录样本名后续多样性分析时要么剔除要么做好心理准备。这一步不需要复杂统计纯粹是“提前认清数据”。质量分数图同样值得看坐标横轴是读长位置纵轴是质量分数会分别显示forward和reverse。我通常在导入后只判断一件事总体质量是不是断崖式下降如果是下一步DADA2的截断参数就要重点照顾这段位置。3. 质量控制与DADA2去噪参数别乱抄要看图3.1 质量图到底怎么读剪哪里留哪里从demux.qzv里可以看到两条曲线ForwardR1和ReverseR2。R1的质量通常前几十个碱基不太好测序仪启动不稳定中间稳定末端开始下降R2整体质量低于R1且末端下降更早更快。你的任务就是确定左端要剪掉多少trim-left右端从哪个位置开始切trunc-len。不要照抄网上的--p-trunc-len-f 240 --p-trunc-len-r 200。不同引物、不同测序平台和不同插入片段长度会导致最优参数差很多。我自己的判断逻辑是观察质量曲线找到质量分数中位数开始低于30的那个位置。给“重叠区”留足余量。双端reads后续dada2要合并如果R1/R2截得太短两个方向没有足够重叠会直接导致拼接失败。一般要求截断后至少留有20-50bp重叠长度低于100bp会使dada2默认拒绝合并。不要把曲线尾部硬保质量差的碱基会显著增加错误ASV。3.2 DADA2参数选择trim-left与trunc-len的执行逻辑执行去噪的命令如下qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux.qza \ --p-trim-left-f 0 \ --p-trim-left-r 0 \ --p-trunc-len-f 240 \ --p-trunc-len-r 200 \ --p-n-threads 0 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats denoising-stats.qzatrim-left-f和trim-left-r是把每个read前几个碱基剪掉目的是去除引物序列。如果测序完成后公司已经去引物那这里就用0。如果没去引物按你所用引物的长度来设置比如515F/806R这对V4区引物通常正向剪掉19bp、反向剪掉20bp但具体以引物序列的实际长度为准。注意这里有个原则剪错左端比截错右端更难受因为左端是结构性序列剪多了会整体偏移后续比对注释全乱。trunc-len-f和trunc-len-r是从read的3‘端截断到指定长度。这个参数越接近质量图显示的“下降点”越好。比如forward在240bp后质量明显下滑就设240reverse在200bp后下滑就设200。设得太大会保留大量错误碱基设得太小又会损失重叠区信息。我见过一个项目reverse只有180bp有效长度有人硬按教程设了250结果合并成功率只有20%重跑后设置成180成功率回到85%。3.3 去噪后必须看的两个数量指标DADA2跑完会产出三个文件table.qza特征表、rep-seqs.qza代表序列、denoising-stats.qza统计信息。先查看denoising-statsqiime metadata tabulate \ --m-input-file denoising-stats.qza \ --o-visualization denoising-stats.qzv这个表会列出每个样本的input输入reads数、filtered过滤后、denoised去噪后、merged合并后、non-chimeric去嵌合体后。我一般会算一个比例non-chimeric除以input正常情况在50%-80%之间。如果低于40%先考虑截断参数是否太保守再考虑数据质量是不是本身有问题。紧接着查看特征表qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv这里最关键是看每个样本的“Frequency per sample”分布。如果样本之间测序深度差一个数量级比如最少的才3000最多的30万那下一步多样性分析的采样深度选择会变得很难受。我会在这个环节记录两个数最小样本reads数和所有样本频率分布的中位数后续core-metrics选择采样深度就用得上。4. Alpha与Beta多样性分析样品组的差异从哪里看出来4.1 先建系统发育树再做多样性不是多此一举有些做过旧流程的人会疑惑为什么QIIME2的多样性分析之前要先跑q2-phylogeny因为UniFrac距离需要系统发育树提供物种间的亲缘信息单纯Bray-Curtis不需要树但要做加权/非加权UniFrac就必须有树。一条龙流程里这一步几乎是必跑的qiime phylogeny align-to-tree-mafft-fasttree \ --i-sequences rep-seqs.qza \ --o-alignment aligned-rep-seqs.qza \ --o-masked-alignment masked-aligned-rep-seqs.qza \ --o-unrooted-tree unrooted-tree.qza \ --o-rooted-tree rooted-tree.qza这里我不建议换参数。MAFFT用来多序列比对FastTree用来构建进化树对于16S这种保守区域组合这组工具配合得最稳。构建完成后rooted-tree.qza在UniFrac计算中会被用到而unrooted-tree.qza在需根分析中同样有用。全部保留即可。4.2 Alpha多样性香农指数、Chao1与抽平深度的逻辑执行核心多样性pipelineqiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 10000 \ --m-metadata-file sample-metadata.tsv \ --output-dir core-metrics-results这里最容易犯错的是sampling-depth——抽平深度。这个值的意思是为了让每个样本的测序深度可比把每个样本随机抽到同样数量的reads。如果设成10000那所有低于10000的样本都会被丢弃。选太高样本量骤减选太低数据浪费严重。我的基本原则是先看table.qzv里的频率分布取一个能保留80%以上样本的深度值。如果样本量少或者实验设计不允许丢样本宁可选低一点的深度也要保留更多生物学重复。Alpha多样性结果重点关注这几个指标observed_features观察到的ASV数、shannon香农指数、faith_pd基于系统发育的多样性、evenness均匀度。我一般先看shannon因为它兼顾物种丰富度和均匀度。如果想看每组之间的差异是否显著可以继续用Kruskal-Wallis检验qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/shannon_vector.qza \ --m-metadata-file sample-metadata.tsv \ --o-visualization alpha-shannon-group-significance.qzv4.3 Beta多样性PCoA图怎么讲出生物学故事Beta多样性的核心是样本间差异距离。QIIME2默认输出四个距离矩阵unweighted_unifrac只看物种有无差异、weighted_unifrac同时考虑物种相对丰度、bray_curtis、jaccard。实操里我更关注weighted和unweighted两个UniFrac前者反映群落组成的量变后者反映群落组成的质变。如果两组分开但两个指标结论不一致往往说明差异来自丰度低的稀有物种而不是优势物种变化。直接用主坐标分析PCoA画图太常见但你得会看qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \ --m-metadata-file sample-metadata.tsv \ --m-metadata-column Group \ --o-visualization beta-weight-unifrac-group-significance.qzv在QIIME2 View里看PCoA图时我常做三件事第一按分组着色判断组间是否形成明显聚类第二旋转不同主坐标轴不能只盯着PC1/PC2有些数据在PC3/PC4上分离得更明显第三配合PERMANOVA检验也就是beta-group-significance输出的p-valuep0.05才谈得上组间显著差异。如果PCoA图上看着有分离趋势但p值不显著可能是样本量太少或组内异质性太大不能简单下结论说“没有差异”。5. 分类学注释数据库选择与分类器训练的取舍5.1 三个主流数据库怎么选分类学注释就是把ASV代表序列与参考数据库比对猜出它们属于哪个门纲目科属种。目前最常用的三个16S数据库数据库特点适用场景Greengenes 13_8历史最悠久2013年后停滞更新很多早期文献用它需要与前人结果直接对比的纵向研究SILVA 138覆盖古菌、细菌、真核更新及时分类层级完整通用推荐是目前大多数新研究的选择GTDB以基因组数据为框架重排分类更贴近基因组分类学注重准确分类、愿意接受新分类体系的团队我个人的建议如果课题组之前大量使用了Greengenes为了批阅对比一致性继续用Greengenes也说得通但新项目建议直接用SILVA或GTDB。尤其临床样本或环境样本中有不少未培养类群SILVA在这块的注释能力更完备。5.2 预训练分类器 vs 自己训练QIIME2官方提供了基于Greengenes的预训练分类器515F/806R区域但针对SILVA或GTDB预训练分类器有时不是专门针对你引物区域训练的这时我更建议自己训练。过程分两步第一步按引物区域从参考序列里提取片段qiime feature-classifier extract-reads \ --i-sequences silva-138-99-seqs.qza \ --p-f-primer GTGCCAGCMGCCGCGGTAA \ --p-r-primer GGACTACHVGGGTWTCTAAT \ --p-trunc-len 300 \ --o-reads ref-seqs.qza第二步训练朴素贝叶斯分类器qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads ref-seqs.qza \ --i-reference-taxonomy ref-taxonomy.qza \ --o-classifier silva-138-515-806-classifier.qza自己训练的耗时通常几十分钟到一两个小时完全可以接受。训练好的分类器要妥善保存同一个引物区域的多个项目都能复用。如果你用的是公司交付的通用V3-V4区引物比如341F/806R提取片段时一定换成对应的正向和反向引物序列不要拿V4区的515F/806R硬套。5.3 注释结果的可视化与分层柱状图分类注释的核心命令qiime feature-classifier classify-sklearn \ --i-classifier silva-138-515-806-classifier.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv注释结果里你会看到每个ASV对应的一系列分类层级以及置信度。实操里我有一个检查习惯看属水平上有没有“未分类”的比例过高的样本如果某个样本Unclassified占比超过30%先不要急着怀疑注释回去看看该样本reads数是不是太少或者质控后留下的序列是不是质量偏低。从这里还能生成分类学堆叠柱状图qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file sample-metadata.tsv \ --o-visualization taxa-bar-plots.qzv在View里可以切换不同分类层级门/纲/目/科/属按分组着色查看相对丰度结构。我一般会先用“门水平”看宏观趋势再用“属水平”看具体差异类群。6. 下游差异分析与R生态衔接一条龙的最后一段路程6.1 ANCOM与ANCOM-BC的选型逻辑差异丰度分析用来回答“哪些物种在组间显著差异”。QIIME2内置了ANCOM命令是qiime composition ancom \ --i-table table.qza \ --m-metadata-file sample-metadata.tsv \ --m-metadata-column Group \ --o-visualization ancom.qzv注意ANCOM有严格的内部假设它假设绝大多数分类单元在组间没有差异通常少于25%的OTU发生差异。如果实际数据里差异物种占比很高ANCOM的结果会偏保守很多真正有差异的物种被检测不出来。因此如果样本量允许我更推荐用较新版本的QIIME2中提供的q2-composition扩展或直接在R里跑ANCOM-BC。ANCOM-BC加入了偏差校正可以输出完整的p值和置信区间处理稀疏数据的能力也更强。实操时我还有一个提醒差异分析之前一定要把特征表里丰度太低、出现于太少样本的ASV过滤掉。我个人常用过滤规则至少在某一个分组内超过20%的样本中出现且总丰度大于等于5。过滤命令qiime feature-table filter-features \ --i-table table.qza \ --p-min-frequency 5 \ --p-min-samples 2 \ --o-filtered-table filtered-table.qza不事先过滤少数ASV的极端丰度变化会干扰整体分析的稳健性。6.2 功能预测PICRUSt2注意它不在内置插件里16S本身只提供物种组成信息功能预测则是利用参考基因组推断群落可能具备的代谢通路。PICRUSt2是目前的主流工具。需要注意PICRUSt2已经不再作为QIIME2内置插件发布需单独安装并作为独立命令行使用conda install -c bioconda -c conda-forge picrust2使用流程大致是把QIIME2的table和rep-seqs导出为PICRUSt2能接受的BIOM格式和fasta格式qiime tools export --input-path table.qza --output-path picrust-input qiime tools export --input-path rep-seqs.qza --output-path picrust-input biom convert -i picrust-input/feature-table.biom \ -o picrust-input/feature-table.tsv --to-tsv然后运行PICRUSt2的一个核心命令picrust2_pipeline.py \ -s picrust-input/dna-sequences.fna \ -i picrust-input/feature-table.biom \ -o picrust-output \ -p 4生成的预测结果包含KEGG Orthologs、MetaCyc pathway等。做这一步前务必想清楚16S的功能预测是基于系统发育推断的“潜在功能”不是实测的宏基因组结果。你可以在讨论里描述为“群落代谢潜能的初步推测”不能直接当宏转录组或宏基因组结论使用。6.3 把数据导出给R生态phyloseq与ggplot2的衔接不少下游可视化最终要离开QIIME2界面去R里完成比如个性化的PCoA图、热图或LEfSe分析。最常用的做法是导出BIOM和分类学表qiime tools export --input-path filtered-table.qza --output-path export/table qiime tools export --input-path taxonomy.qza --output-path export/taxonomy qiime tools export --input-path rooted-tree.qza --output-path export/tree导出后的feature-table.biom用R的phyloseq读取library(phyloseq) ps - import_biom(export/table/feature-table.biom) tax - read.delim(export/taxonomy/taxonomy.tsv, row.names 1)然后再做alpha多样性箱线图、NMDS/PCA散点图都相对自由。这里我想额外提一句在QIIME2 View里直接导出的PCoA图适合快速看趋势但论文投稿建议还是用R重新绘制ggplot2对点的大小、颜色、标签、置信椭圆的控制要灵活得多。6.4 LEfSe之类工具的定位很多教程会把LEfSe线性判别分析效应量也塞进QIIME2流程但严格说LEfSe是独立工具通常需要从Galaxy或本地运行。它适合筛选组间生物学差异的标记类群输出LDA得分与柱状图。实际项目中我会先用ANCOM或ANCOM-BC做初步筛选再用LEfSe丰富结果的可视化。如果你有分组且样本量足够LEfSe对故事线的构建很有帮助但它同样对样本量敏感样本太小时结果容易被个别极端值驱动。说到最后我还是想强调“一条龙”并不等于“一路默认参数跑到底”。从导入数据开始每一步的取舍实际上都在影响最后的生物学结论。最典型的就是DADA2截断长度和抽平深度选择这两个参数在不同项目里必须分别判断。作为从业者我建议你每个环节都留下日志、保存好qzv可视化文件并在Methods里准确记录每个参数。这样就算半年后审稿人要求补分析你也能迅速复现整个流程而不是对着文件夹里几百个.qza发呆。如果在实战中遇到导入报错、注释结果异常或者多样性图与预期不符欢迎在评论区带上你的具体报错信息或qzv截图来讨论我看到后会结合自己的项目经验帮你一起排查。16S分析这条路上每一个参数背后都是可积累的细节愿你少走弯路。
QIIME2实战：16S扩增子测序从数据导入到差异分析一条龙
别人拿到测序公司交付的fastq数据包解压后几百个.gz文件第一反应通常是找教程从导入数据开始跑。但QIIME2和传统分析工具最大的不同是它不直接吃fastq你得先理解它的artifact机制否则会在导入这一步反复报错连数据都进不了门。这篇博文我按一条龙流程完整走一遍——导入数据、质控去噪、多样性分析、分类学注释、差异丰度与功能预测——每一步我都会给出命令、讲解参数取舍顺便把容易坑人的细节一次性交代清楚。博文正文1. 开工之前先理解QIIME2的数据组织和运行环境1.1 artifact不是格式洁癖是为了数据追溯先不谈命令我把话放前面QIIME2发音是“chime two”之所以要引入.qza这种封装格式不是故意增加学习门槛而是要解决生物信息项目里最让人头疼的“溯源问题”。举个实际例子我把一个table.qza发给合作者里面嵌入了完整的产生历史包含上游去了多少噪音、用了哪些参数、输入文件是什么谁都可以沿着provenance数据谱系往前追溯。这一点在写论文Methods部分时特别爽——你不需要单独维护一张Excel记录每一步参数qiime tools界面里直接能查看完整操作日志。理解了这层逻辑开篇那些看似繁琐的概念就顺了.qza是“分析工件”的容器.qzv是可视化的容器两者类似一个带说明书的数据包。每次操作都会生成新的对象而你之前的原始数据始终不会被覆盖中间每一步都有迹可循。所以这里的建议是给项目建一个清晰的目录结构建议按input、scripts、output分成三层避免后续分析时被一堆.qza淹没。1.2 环境准备conda、版本选择与资源预判QIIME2官方推荐用Miniconda安装避免污染系统Python环境这一点我强烈建议照做。下面是安装的大致流程wget https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh conda update conda conda create -n qiime2-2023.5 python3.8 conda activate qiime2-2023.5这里的版本号重要。每个版本对应的Python版本和插件版本捆绑较死我不建议你随便用pip去补装插件尽量整体安装官方release。以我常用的2023.5版本为例整个环境装完约占用几个GB空间跑dada2去噪时内存消耗看样本量常规几十个样本8-16GB内存也能跑但如果数据量上了百级PE reads建议准备32GB以上并配合多核CPU用。conda install -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/qiime2/ qiime2如果在国内环境把官方conda源替换成清华镜像会省很多时间但要注意镜像列表里必须包含qiime2这个channel。另外一定记得确认环境中没有同时存在另一个Python版的q2cli否则导入命令会提示找不到可执行文件。1.3 拿到测序数据后的第一件事摸清R1/R2与目录结构Illumina双端测序下机后每个样品会有R1和R2两个fastq.gz文件R1是正向读长R2是反向读长二者中间是插入片段。很多第一次做16S的人容易犯一个错误直接把两个文件单独导入当作单端数据处理这样会把双端信息割裂后续dada2拼接就没有意义。正确思路是一对R1/R2在导入后代表同一个样本的两个方向样本数量应该是文件对数而不是文件个数。先解压并看一眼目录ls rawdata/ | head正常情况你会看到Sample001_R1.fastq.gz、Sample001_R2.fastq.gz这种成对文件第三方公司还常带barcode和adapter信息。这里务必先把样本名统一只保留字母、数字和下划线。如果样本名下划线太多后续分组信息会写得很痛苦。我最开始跑项目时公司文件叫A-1_S1_L001_R1_001.fastq.gz光是把名字改成A1.R1.fastq.gz就折腾了半天。2. 数据导入manifest文件、常见报错与验证2.1 两种导入姿势Casava目录格式与manifest清单QIIME2支持多种导入格式最常见的是CasavaOneEightSingleLanePerSampleDirFmt和manifest。Casava格式适合测序公司交付的原始目录还保留着标准命名的场景——目录层级固定为SampleID_SampleNumber_L001_R1_001.fastq.gz此时可以直接qiime tools import \ --type CasavaOneEightSingleLanePerSampleDirFmt \ --input-path rawdata \ --output-path demux.qza但我遇到更多的情况是公司交付的文件名已经改成了自己的命名习惯或样本来自多个批次。这时老老实实用manifest清单。先在文本编辑器里生成一个manifest.tsv三列sample-id、forward-absolute-filepath、reverse-absolute-filepath。sample-id forward-absolute-filepath reverse-absolute-filepath A1 /absolute/path/rawdata/A1_R1.fastq.gz /absolute/path/rawdata/A1_R2.fastq.gz A2 /absolute/path/rawdata/A2_R1.fastq.gz /absolute/path/rawdata/A2_R2.fastq.gz B1 /absolute/path/rawdata/B1_R1.fastq.gz /absolute/path/rawdata/B1_R2.fastq.gz然后导入qiime tools import \ --type SampleData[PairedEndSequencesWithQuality] \ --input-path manifest.tsv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33V2这里有个细节需要说明Phred33和Phred64两种格式容易混淆。现在主流Illumina下机数据基本都是Phred33默认选PairedEndFastqManifestPhred33V2不会错。如果分析老数据不确定时用seqkit head检查质量值字符范围ASCII码43以下对应Phred3343以上则是Phred64再对应修改格式。2.2 manifest文件的隐形坑绝对路径、制表符和CRLF很多人在导入时报File not found或解析错误多半不是数据问题而是路径问题。manifest里必须写绝对路径而且要保证路径里的目录和文件真的有读权限。服务器环境下如果是从Windows下编辑的manifest记事本默认换行符是CRLF而Linux下解析时会残留\r字符导致样本ID后面多个不可见符号。解决办法是sed -i s/\r$// manifest.tsv还有一种坑是字段分隔符到底用tab还是逗号。PairedEndFastqManifestPhred33V2严格识别TSV制表符分隔而旧的PairedEndFastqManifestPhred33虽也是TSV但如果从Excel复制粘贴会把tab换成空格导致解析字段数不对。我的习惯是用VSCode或Notepad先开启“显示所有字符”来确认分隔符再上传服务器。2.3 导入后第一件事查看每个样本的reads数导入完成不代表数据好用你还需要马上验证qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzv然后把demux.qzv拖到QIIME2 View网页里查看。这里主要看两个信息一是每个样本的序列条数判断有没有样本明显偏少二是质量分数图。如果看到某个样本reads数只有几千而其他样本十几万尽早记录样本名后续多样性分析时要么剔除要么做好心理准备。这一步不需要复杂统计纯粹是“提前认清数据”。质量分数图同样值得看坐标横轴是读长位置纵轴是质量分数会分别显示forward和reverse。我通常在导入后只判断一件事总体质量是不是断崖式下降如果是下一步DADA2的截断参数就要重点照顾这段位置。3. 质量控制与DADA2去噪参数别乱抄要看图3.1 质量图到底怎么读剪哪里留哪里从demux.qzv里可以看到两条曲线ForwardR1和ReverseR2。R1的质量通常前几十个碱基不太好测序仪启动不稳定中间稳定末端开始下降R2整体质量低于R1且末端下降更早更快。你的任务就是确定左端要剪掉多少trim-left右端从哪个位置开始切trunc-len。不要照抄网上的--p-trunc-len-f 240 --p-trunc-len-r 200。不同引物、不同测序平台和不同插入片段长度会导致最优参数差很多。我自己的判断逻辑是观察质量曲线找到质量分数中位数开始低于30的那个位置。给“重叠区”留足余量。双端reads后续dada2要合并如果R1/R2截得太短两个方向没有足够重叠会直接导致拼接失败。一般要求截断后至少留有20-50bp重叠长度低于100bp会使dada2默认拒绝合并。不要把曲线尾部硬保质量差的碱基会显著增加错误ASV。3.2 DADA2参数选择trim-left与trunc-len的执行逻辑执行去噪的命令如下qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux.qza \ --p-trim-left-f 0 \ --p-trim-left-r 0 \ --p-trunc-len-f 240 \ --p-trunc-len-r 200 \ --p-n-threads 0 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats denoising-stats.qzatrim-left-f和trim-left-r是把每个read前几个碱基剪掉目的是去除引物序列。如果测序完成后公司已经去引物那这里就用0。如果没去引物按你所用引物的长度来设置比如515F/806R这对V4区引物通常正向剪掉19bp、反向剪掉20bp但具体以引物序列的实际长度为准。注意这里有个原则剪错左端比截错右端更难受因为左端是结构性序列剪多了会整体偏移后续比对注释全乱。trunc-len-f和trunc-len-r是从read的3‘端截断到指定长度。这个参数越接近质量图显示的“下降点”越好。比如forward在240bp后质量明显下滑就设240reverse在200bp后下滑就设200。设得太大会保留大量错误碱基设得太小又会损失重叠区信息。我见过一个项目reverse只有180bp有效长度有人硬按教程设了250结果合并成功率只有20%重跑后设置成180成功率回到85%。3.3 去噪后必须看的两个数量指标DADA2跑完会产出三个文件table.qza特征表、rep-seqs.qza代表序列、denoising-stats.qza统计信息。先查看denoising-statsqiime metadata tabulate \ --m-input-file denoising-stats.qza \ --o-visualization denoising-stats.qzv这个表会列出每个样本的input输入reads数、filtered过滤后、denoised去噪后、merged合并后、non-chimeric去嵌合体后。我一般会算一个比例non-chimeric除以input正常情况在50%-80%之间。如果低于40%先考虑截断参数是否太保守再考虑数据质量是不是本身有问题。紧接着查看特征表qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv这里最关键是看每个样本的“Frequency per sample”分布。如果样本之间测序深度差一个数量级比如最少的才3000最多的30万那下一步多样性分析的采样深度选择会变得很难受。我会在这个环节记录两个数最小样本reads数和所有样本频率分布的中位数后续core-metrics选择采样深度就用得上。4. Alpha与Beta多样性分析样品组的差异从哪里看出来4.1 先建系统发育树再做多样性不是多此一举有些做过旧流程的人会疑惑为什么QIIME2的多样性分析之前要先跑q2-phylogeny因为UniFrac距离需要系统发育树提供物种间的亲缘信息单纯Bray-Curtis不需要树但要做加权/非加权UniFrac就必须有树。一条龙流程里这一步几乎是必跑的qiime phylogeny align-to-tree-mafft-fasttree \ --i-sequences rep-seqs.qza \ --o-alignment aligned-rep-seqs.qza \ --o-masked-alignment masked-aligned-rep-seqs.qza \ --o-unrooted-tree unrooted-tree.qza \ --o-rooted-tree rooted-tree.qza这里我不建议换参数。MAFFT用来多序列比对FastTree用来构建进化树对于16S这种保守区域组合这组工具配合得最稳。构建完成后rooted-tree.qza在UniFrac计算中会被用到而unrooted-tree.qza在需根分析中同样有用。全部保留即可。4.2 Alpha多样性香农指数、Chao1与抽平深度的逻辑执行核心多样性pipelineqiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 10000 \ --m-metadata-file sample-metadata.tsv \ --output-dir core-metrics-results这里最容易犯错的是sampling-depth——抽平深度。这个值的意思是为了让每个样本的测序深度可比把每个样本随机抽到同样数量的reads。如果设成10000那所有低于10000的样本都会被丢弃。选太高样本量骤减选太低数据浪费严重。我的基本原则是先看table.qzv里的频率分布取一个能保留80%以上样本的深度值。如果样本量少或者实验设计不允许丢样本宁可选低一点的深度也要保留更多生物学重复。Alpha多样性结果重点关注这几个指标observed_features观察到的ASV数、shannon香农指数、faith_pd基于系统发育的多样性、evenness均匀度。我一般先看shannon因为它兼顾物种丰富度和均匀度。如果想看每组之间的差异是否显著可以继续用Kruskal-Wallis检验qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/shannon_vector.qza \ --m-metadata-file sample-metadata.tsv \ --o-visualization alpha-shannon-group-significance.qzv4.3 Beta多样性PCoA图怎么讲出生物学故事Beta多样性的核心是样本间差异距离。QIIME2默认输出四个距离矩阵unweighted_unifrac只看物种有无差异、weighted_unifrac同时考虑物种相对丰度、bray_curtis、jaccard。实操里我更关注weighted和unweighted两个UniFrac前者反映群落组成的量变后者反映群落组成的质变。如果两组分开但两个指标结论不一致往往说明差异来自丰度低的稀有物种而不是优势物种变化。直接用主坐标分析PCoA画图太常见但你得会看qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \ --m-metadata-file sample-metadata.tsv \ --m-metadata-column Group \ --o-visualization beta-weight-unifrac-group-significance.qzv在QIIME2 View里看PCoA图时我常做三件事第一按分组着色判断组间是否形成明显聚类第二旋转不同主坐标轴不能只盯着PC1/PC2有些数据在PC3/PC4上分离得更明显第三配合PERMANOVA检验也就是beta-group-significance输出的p-valuep0.05才谈得上组间显著差异。如果PCoA图上看着有分离趋势但p值不显著可能是样本量太少或组内异质性太大不能简单下结论说“没有差异”。5. 分类学注释数据库选择与分类器训练的取舍5.1 三个主流数据库怎么选分类学注释就是把ASV代表序列与参考数据库比对猜出它们属于哪个门纲目科属种。目前最常用的三个16S数据库数据库特点适用场景Greengenes 13_8历史最悠久2013年后停滞更新很多早期文献用它需要与前人结果直接对比的纵向研究SILVA 138覆盖古菌、细菌、真核更新及时分类层级完整通用推荐是目前大多数新研究的选择GTDB以基因组数据为框架重排分类更贴近基因组分类学注重准确分类、愿意接受新分类体系的团队我个人的建议如果课题组之前大量使用了Greengenes为了批阅对比一致性继续用Greengenes也说得通但新项目建议直接用SILVA或GTDB。尤其临床样本或环境样本中有不少未培养类群SILVA在这块的注释能力更完备。5.2 预训练分类器 vs 自己训练QIIME2官方提供了基于Greengenes的预训练分类器515F/806R区域但针对SILVA或GTDB预训练分类器有时不是专门针对你引物区域训练的这时我更建议自己训练。过程分两步第一步按引物区域从参考序列里提取片段qiime feature-classifier extract-reads \ --i-sequences silva-138-99-seqs.qza \ --p-f-primer GTGCCAGCMGCCGCGGTAA \ --p-r-primer GGACTACHVGGGTWTCTAAT \ --p-trunc-len 300 \ --o-reads ref-seqs.qza第二步训练朴素贝叶斯分类器qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads ref-seqs.qza \ --i-reference-taxonomy ref-taxonomy.qza \ --o-classifier silva-138-515-806-classifier.qza自己训练的耗时通常几十分钟到一两个小时完全可以接受。训练好的分类器要妥善保存同一个引物区域的多个项目都能复用。如果你用的是公司交付的通用V3-V4区引物比如341F/806R提取片段时一定换成对应的正向和反向引物序列不要拿V4区的515F/806R硬套。5.3 注释结果的可视化与分层柱状图分类注释的核心命令qiime feature-classifier classify-sklearn \ --i-classifier silva-138-515-806-classifier.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv注释结果里你会看到每个ASV对应的一系列分类层级以及置信度。实操里我有一个检查习惯看属水平上有没有“未分类”的比例过高的样本如果某个样本Unclassified占比超过30%先不要急着怀疑注释回去看看该样本reads数是不是太少或者质控后留下的序列是不是质量偏低。从这里还能生成分类学堆叠柱状图qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file sample-metadata.tsv \ --o-visualization taxa-bar-plots.qzv在View里可以切换不同分类层级门/纲/目/科/属按分组着色查看相对丰度结构。我一般会先用“门水平”看宏观趋势再用“属水平”看具体差异类群。6. 下游差异分析与R生态衔接一条龙的最后一段路程6.1 ANCOM与ANCOM-BC的选型逻辑差异丰度分析用来回答“哪些物种在组间显著差异”。QIIME2内置了ANCOM命令是qiime composition ancom \ --i-table table.qza \ --m-metadata-file sample-metadata.tsv \ --m-metadata-column Group \ --o-visualization ancom.qzv注意ANCOM有严格的内部假设它假设绝大多数分类单元在组间没有差异通常少于25%的OTU发生差异。如果实际数据里差异物种占比很高ANCOM的结果会偏保守很多真正有差异的物种被检测不出来。因此如果样本量允许我更推荐用较新版本的QIIME2中提供的q2-composition扩展或直接在R里跑ANCOM-BC。ANCOM-BC加入了偏差校正可以输出完整的p值和置信区间处理稀疏数据的能力也更强。实操时我还有一个提醒差异分析之前一定要把特征表里丰度太低、出现于太少样本的ASV过滤掉。我个人常用过滤规则至少在某一个分组内超过20%的样本中出现且总丰度大于等于5。过滤命令qiime feature-table filter-features \ --i-table table.qza \ --p-min-frequency 5 \ --p-min-samples 2 \ --o-filtered-table filtered-table.qza不事先过滤少数ASV的极端丰度变化会干扰整体分析的稳健性。6.2 功能预测PICRUSt2注意它不在内置插件里16S本身只提供物种组成信息功能预测则是利用参考基因组推断群落可能具备的代谢通路。PICRUSt2是目前的主流工具。需要注意PICRUSt2已经不再作为QIIME2内置插件发布需单独安装并作为独立命令行使用conda install -c bioconda -c conda-forge picrust2使用流程大致是把QIIME2的table和rep-seqs导出为PICRUSt2能接受的BIOM格式和fasta格式qiime tools export --input-path table.qza --output-path picrust-input qiime tools export --input-path rep-seqs.qza --output-path picrust-input biom convert -i picrust-input/feature-table.biom \ -o picrust-input/feature-table.tsv --to-tsv然后运行PICRUSt2的一个核心命令picrust2_pipeline.py \ -s picrust-input/dna-sequences.fna \ -i picrust-input/feature-table.biom \ -o picrust-output \ -p 4生成的预测结果包含KEGG Orthologs、MetaCyc pathway等。做这一步前务必想清楚16S的功能预测是基于系统发育推断的“潜在功能”不是实测的宏基因组结果。你可以在讨论里描述为“群落代谢潜能的初步推测”不能直接当宏转录组或宏基因组结论使用。6.3 把数据导出给R生态phyloseq与ggplot2的衔接不少下游可视化最终要离开QIIME2界面去R里完成比如个性化的PCoA图、热图或LEfSe分析。最常用的做法是导出BIOM和分类学表qiime tools export --input-path filtered-table.qza --output-path export/table qiime tools export --input-path taxonomy.qza --output-path export/taxonomy qiime tools export --input-path rooted-tree.qza --output-path export/tree导出后的feature-table.biom用R的phyloseq读取library(phyloseq) ps - import_biom(export/table/feature-table.biom) tax - read.delim(export/taxonomy/taxonomy.tsv, row.names 1)然后再做alpha多样性箱线图、NMDS/PCA散点图都相对自由。这里我想额外提一句在QIIME2 View里直接导出的PCoA图适合快速看趋势但论文投稿建议还是用R重新绘制ggplot2对点的大小、颜色、标签、置信椭圆的控制要灵活得多。6.4 LEfSe之类工具的定位很多教程会把LEfSe线性判别分析效应量也塞进QIIME2流程但严格说LEfSe是独立工具通常需要从Galaxy或本地运行。它适合筛选组间生物学差异的标记类群输出LDA得分与柱状图。实际项目中我会先用ANCOM或ANCOM-BC做初步筛选再用LEfSe丰富结果的可视化。如果你有分组且样本量足够LEfSe对故事线的构建很有帮助但它同样对样本量敏感样本太小时结果容易被个别极端值驱动。说到最后我还是想强调“一条龙”并不等于“一路默认参数跑到底”。从导入数据开始每一步的取舍实际上都在影响最后的生物学结论。最典型的就是DADA2截断长度和抽平深度选择这两个参数在不同项目里必须分别判断。作为从业者我建议你每个环节都留下日志、保存好qzv可视化文件并在Methods里准确记录每个参数。这样就算半年后审稿人要求补分析你也能迅速复现整个流程而不是对着文件夹里几百个.qza发呆。如果在实战中遇到导入报错、注释结果异常或者多样性图与预期不符欢迎在评论区带上你的具体报错信息或qzv截图来讨论我看到后会结合自己的项目经验帮你一起排查。16S分析这条路上每一个参数背后都是可积累的细节愿你少走弯路。
