★ FEATURED ARTICLE 简介一份面向微生物组研究者的QIIME 2流程操作指南适合有一定测序数据基础、希望系统掌握16S rRNA基因扩增子分析方法的科研人员与初学者。文档基于《Bio-protocol》发表的标准实验流程整理内容涵盖QIIME 2在Linux服务器、Windows 10子系统、VirtualBox虚拟机及Mac系统四种环境下的安装部署方案与选择建议并完整介绍了数据导入、特征表生成、α多样性分析、β多样性分析、物种组成分析、差异物种分析及数据可视化等标准分析步骤同时对各环节的关键参数和常见报错给出了详细解读与优化思路。资源为单个docx文档压缩包约578KB便于查阅与本地保存目前已有2796人学习下载。文档还提供了配套分析代码、测序数据与预期结果的获取线索以及计算硬件配置建议可帮助缺乏生物信息学背景的研究者快速上手并完整复现微生物组分析流程。1. 先搞懂一件事16S扩增子测序数据到底在分析什么第一次拿到16S rRNA基因扩增子测序数据的人大多是懵的测序公司给了一堆几十G的压缩文件文件名像乱码解压出来全是ATCG碱基你真正想知道的“样品里有哪些微生物、各组之间差在哪”在这些序列里完全看不出来。QIIME 2 就是当前微生物组扩增子分析里用得最广的开源流程它用 Python 3 重写把导入数据、降噪、多样性分析、物种注释这些步骤串成一条标准流水线你照着顺序执行就能把原始序列变成特征表、多样性指数和能放进论文的图。适合谁懂实验设计、但没系统学过生物信息的人尤其是刚接触微生物组的新手研究生——这篇文章就是为你准备的从安装到跑通第一个例子一步步来。2. 为什么不直接学老版 QIIME以及怎么从零装好它2.1 老版 QIIME 很好但已经追不上现在的数据量老版 QIIME 在扩增子分析里统治了将近十年引用次数早就过了两万很多人至今还在用。但它的底层设计是单机图形界面加脚本堆叠跑现在的数据量越来越吃力交互式图表少想换一个算法就得等开发者适配。QIIME 2 把这一切推倒重来它用 Python 3 编写每个分析步骤都是一个独立插件比如 DADA2 负责降噪、feature-classifier 负责物种注释、diversity 负责多样性分析插件可以单独更新和扩展。这意味着新算法出来以后你不用等整个软件重装装一个新插件就能用。更关键的是可重复性。微生物组分析步骤多、参数杂老版流程跑完以后别人很难知道你每一步到底用了什么参数。QIIME 2 的每个结果文件都自带 provenance 记录也就是“这份结果是怎么一步步生成的”完整家谱对发表论文和团队协作非常友好。我第一次学的时候就是被这一点打动的——做实验翻车还能重来分析翻车连参数都记不住才叫真的头疼。2.2 安装前先想明白你准备在哪台机器上跑它QIIME 2 不支持 Windows 直接安装这是新手劝退的第一道坎。实际可选的方案有四种按推荐程度排Linux 服务器强烈推荐适合大项目Windows 10 自带的 Linux 子系统也就是 WSL 里的 Ubuntu 20.04推荐适合小数据Windows 里用 VirtualBox 开一个 Ubuntu 虚拟机能用但效率低Mac 直接装兼容性问题比较多不建议折腾如果你是课题组里已经有计算服务器的直接走方案一用 XShell 这类终端工具远程连上去操作分析速度最快也不用折腾子系统。如果只是个人电脑想先跑通流程方案二最省心在 Microsoft Store 里搜索 Ubuntu 20.04 LTS安装后打开就是一个 Linux 命令行数据放在 C 盘就对应/mnt/c/路径操作起来和真 Linux 几乎一样。注意虚拟机方案虽然看起来“最像 Linux”但 IO 效率很低跑 DADA2 这种密集计算任务会明显变慢新手不建议选。2.3 用 Miniconda 安装 QIIME 2一条命令创建一个干净环境装 QIIME 2 之前先装 Miniconda它是一个软件管理器作用是帮你创建互相隔离的 Python 环境。为什么要隔离因为 QIIME 2 依赖特定版本的 Python 库装到系统全局很容易和别的软件冲突用 Conda 单独建一个环境装上 QIIME 2不影响系统里其他任何东西。安装 Miniconda 的三步命令# 1. 从 Miniconda 官网下载 Linux 版安装脚本 wget -c Miniconda3-latest-Linux-x86_64.sh 的下载链接 # 2. 运行安装脚本全程按回车和 yes 即可 bash Miniconda3-latest-Linux-x86_64.sh # 3. 让 Conda 命令初始化到当前终端 ~/miniconda3/bin/conda init完成后关掉终端重开一次conda命令就能用了。接着下载 QIIME 2 这个版本的软件列表创建环境并激活# 下载 QIIME 2 2020.6 的软件依赖清单 wget -c qiime2-2020.6-py36-linux-conda.yml 的下载链接 # 根据清单创建名为 qiime2-2020.6 的新环境 conda env create -n qiime2-2020.6 --file qiime2-2020.6-py36-linux-conda.yml # 每次分析前先激活这个环境 conda activate qiime2-2020.6这里-n后面是环境名你可以改成自己喜欢的名字但建议保留版本号方便以后区分。激活成功以后命令行前面会出现(qiime2-2020.6)前缀说明你已经进入 QIIME 2 的专属环境了。2.4 下载并导入数据库这一步决定物种注释准不准QIIME 2 本身不带物种数据库你需要单独下载参考数据库并训练一个“分类器”。类比一下分类器就像一本词典后面做物种注释时拿着序列去词典里查最像哪个物种。16S 扩增子分析最常用的参考数据库之一是 GreenGenes 13.8按 99% 相似度聚类的那套序列最常用。下载解压后先导入参考序列再导入物种分类信息# 下载数据库压缩包并解压 wget -c gg_13_8_otus.tar.gz 的下载链接 tar -zxvf gg_13_8_otus.tar.gz # 导入参考序列fasta 格式为 QIIME 2 的数据文件 qza qiime tools import \ --type FeatureData[Sequence] \ --input-path gg_13_8_otus/rep_set/99_otus.fasta \ --output-path 99_otus.qza # 导入每个序列对应的物种分类信息 qiime tools import \ --type FeatureData[Taxonomy] \ --input-format HeaderlessTSVTaxonomyFormat \ --input-path gg_13_8_otus/taxonomy/99_otu_taxonomy.txt \ --output-path ref-taxonomy.qza然后用这两份导入好的数据训练分类器# 训练全长分类器耗时约半小时 time qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads 99_otus.qza \ --i-reference-taxonomy ref-taxonomy.qza \ --o-classifier classifier_gg_13_8_99.qza如果你扩增的是 16S 的特定区域比如 V5-V7 区还可以从参考序列里先把对应引物区间提取出来再基于这个区间训练一个区域特异分类器精度会更高。提取序列用extract-reads把前后引物序列作为参数传给插件qiime feature-classifier extract-reads \ --i-sequences 99_otus.qza \ --p-f-primer AACMGGATTAGATACCCKG \ --p-r-primer ACGTCATCCCCACCTTCC \ --o-reads ref-seqs.qza # 基于提取出的目标区段训练耗时约 8 分钟 qiime feature-classifier fit-classifier-naive-bayes \ --i-reference-reads ref-seqs.qza \ --i-reference-taxonomy ref-taxonomy.qza \ --o-classifier classifier_gg_13_8_99_V5-V7.qza看到qza文件生成环境就彻底准备好了。3. 跑通第一个分析从原始测序数据到特征表3.1 分析前的三样东西项目目录、元数据、序列文件正式开始分析前先建一个干净的目录组织所有文件。我习惯在 C 盘建一个qiime2文件夹然后进入 Linux 子系统的/mnt/c/qiime2路径操作所有中间结果都放在这里不乱跑mkdir -p qiime2 cd qiime2 conda activate qiime2-2020.6接下来需要准备两份关键文件。第一份是元数据metadata.txt它是每个样本的“身份证”至少有四列样本名、分组信息以及你关心的其他变量比如采样部位、处理条件。第二份是测序数据本身。通常测序公司会返回每个样本一对双端测序文件比如sampleA_1.fq.gz和sampleA_2.fq.gz。如果公司只给了混池的原始数据需要先找公司要拆分后的单样本文件QIIME 2 默认从单样本文件开始分析。为了让你能复现我用的是一份公开的植物根系微生物组数据的子集3 个组、每组 6 个生物学重复、共 18 个样本每个样本随机抽取了 5 万对序列体积小、跑得快适合用来练习完整流程。3.2 写一个 manifest 文件告诉 QIIME 2 每个文件在哪QIIME 2 不像老版那样靠目录结构猜文件位置它要求显示提供一份“清单”说明每个样本的序列文件具体路径。这份清单可以用一行 awk 命令从元数据里自动生成awk NR1{print sample-id\tforward-absolute-filepath\treverse-absolute-filepath} \ NR1{print $1\t$PWD/seq/$1_1.fq.gz\t$PWD/seq/$1_2.fq.gz} \ metadata.txt manifest这段 awk 的作用是第一行打印表头之后的每一行取第 1 列为样本名按$PWD/seq/样本名_1.fq.gz拼出正反向文件的绝对路径。这里的$PWD是当前目录会自动展开成完整路径保证了后面导入时能找到文件。生成后可以cat manifest看一眼确认路径没问题再继续。3.3 数据导入把 fastq 变成 QIIME 2 的 qza所有原始数据都要被转成 QIIME 2 的内部格式.qza。它看起来像一个文件实际上相当于一个压缩包里面有数据本身和所有处理记录。执行导入命令qiime tools import \ --type SampleData[PairedEndSequencesWithQuality] \ --input-path manifest \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33V2参数含义从左到右--type声明数据类型这里指定的是“带质量值的双端测序数据”--input-path指向刚才生成的 manifest 清单--input-format告诉 QIIME 2 manifest 文件的格式规范Phred33 是当前 Illumina 测序最通用的质量编码方式。如果压缩后的 fastq 每个样本不到 100M导入只要几十秒如果是未压缩的 1G 文件可能需要几分钟耐心等。提示导入成功后可以用qiime tools peek demux.qza查看这个文件的内部信息和类型确认没导错。3.4 DADA2 降噪把百万条序列整理成一张特征表导入之后就要去噪。原始序列里有很多测序错误、引物残留和嵌合体DADA2 是 QIIME 2 里默认的降噪方法它不按 97% 相似度聚类 OTU而是直接推断每个真实的生物序列变异也就是 ASV。这一步输出的特征表是后续所有多样性分析的基础表里每一行是一个 ASV每一列是一个样本单元格里的数字是该样本中这个 ASV 的序列数。qiime dada2 denoise-paired \ --i-demultiplexed-seqs demux.qza \ --p-n-threads 8 \ --p-trim-left-f 29 --p-trim-left-r 18 \ --p-trunc-len-f 0 --p-trunc-len-r 0 \ --o-table dada2-table.qza \ --o-representative-sequences dada2-rep-seqs.qza \ --o-denoising-stats denoising-stats.qza--p-trim-left-f和--p-trim-left-r是从左端切除的长度作用是去掉引物和低质量碱基这里的 29 和 18 对应本次数据使用的引物长度换成你自己的实验时这里要改成你的引物长度。--p-trunc-len是从右端截断的位置填 0 表示不截断。这两个参数是 DADA2 分析里最影响结果的“玄学参数”跑之前务必看测序报告里每个碱基位置的质量分数。降噪完成后把 DADA2 的结果复制成主流程的标准命名方便后面命令统一引用cp dada2-table.qza table.qza cp dada2-rep-seqs.qza rep-seqs.qza3.5 查看特征表统计抽平阈值到底怎么定特征表生成后先看数据概况再决定多样性分析时的抽平深度sampling depth。这一步很多人直接抄别人的数字最后结果出问题都不知道原因。执行统计命令qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv \ --m-sample-metadata-file metadata.txt生成的.qzv文件是可视化报表把它上传到 QIIME 2 官网提供的网页预览工具中打开最需要关注的是左边的分位数表和样本量分布图。我这份练习数据的分位数长这样指标数值最小值Minimum frequency27060第一分位数1st quartile28581中位数Median frequency30867第三分位数3rd quartile32860最大值Maximum frequency34663平均值Mean frequency30779抽平阈值的经典选法是“最小值优先”选 27060这样 18 个样本全部保留。如果最小值和第一分位数差距特别大说明个别样本测序深度明显不够这时要在两者之间选一个折中值宁可少保留那个极低深度的样本也要保证多数样本能用上足够的序列。低于阈值的样本会在多样性分析时自动丢弃所以别为了多留一个样本把阈值压得太低。提示现在 Illumina 测序数据量普遍很大抽平阈值小于 5000 基本没有意义推荐从 1 万起步越大越好。4. 多样性分析和物种注释从数字里读出生物学结论4.1 先建一棵进化树多样性分析的前置动作Beta 多样性里的 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、保守区过滤mask和建树FastTree三个阶段。中途出现告警不用慌只要最终生成rooted-tree.qza就能继续。4.2 Alpha 和 Beta 多样性一个看组内一个看组间多样性分析最核心的一步是core-metrics-phylogenetic它一次会计算多种 index输出到一个目录里。执行时把前面定的抽平阈值填进去qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 27060 \ --m-metadata-file metadata.txt \ --output-dir core-metrics-results运行完core-metrics-results目录里会出现四类结果Alpha 多样性指数向量faith_pd、shannon、observed_features、evennessBeta 多样性距离矩阵unweighted_unifrac、weighted_unifrac、bray_curtis、jaccardPCoA 可视化文件以及对应的组间显著性检验结果。Alpha 多样性回答的问题是“每个样本内部微生物种类有多丰富”Beta 多样性回答“样本之间群落结构差异有多大”。看组间差异是否显著用这两条命令# 检验 Alpha 多样性指数在各组间是否有显著差异 qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/observed_features_vector.qza \ --m-metadata-file metadata.txt \ --o-visualization core-metrics-results/observed_features-group-significance.qzv # 检验 Beta 多样性用加权 UniFrac 距离在各组间是否有差异 qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/weighted_unifrac_distance_matrix.qza \ --m-metadata-file metadata.txt \ --m-metadata-column Group \ --o-visualization core-metrics-results/weighted_unifrac-Group-significance.qzv \ --p-pairwiseAlpha 结果里会给出基于 Kruskal-Wallis 检验的 p-value 和 q-valueq-value 是校正过多重比较后的版本看结果时优先看 q-value。Beta 结果里是 PERMANOVA 置换检验的结果同样关注 q-value 是否小于 0.05。我这份练习数据里三个组两两之间 q-value 都小于 0.05说明组间群落结构差异显著可以继续做物种组成找差异来源。4.3 物种注释把 ASV 编号变成看得懂的生物学名字特征表里的 ASV 是一串 32 位的哈希字符串光看它什么结论都得不出要用分类器给每个 ASV 注释上从界到属的物种信息qiime feature-classifier classify-sklearn \ --i-classifier classifier_gg_13_8_99_V5-V7.qza \ --i-reads rep-seqs.qza \ --o-classification taxonomy.qza # 把注释结果转成可视化的交互表格 qiime metadata tabulate \ --m-input-file taxonomy.qza \ --o-visualization taxonomy.qzv打开taxonomy.qzv你会看到一个表格每一行是一个 ASV列分别是特征 ID、注释的物种分类串和置信度。分类串的形式类似k__Bacteria; p__Proteobacteria; c__Betaproteobacteria; ...k是界、p是门、c是纲、o是目、f是科、g是属。如果只注释到属就停了说明数据库里没有更细的参考序列这种情况很常见不一定是你分析错了。4.4 堆叠柱状图展示每个样本的物种组成物种组成最直观的展示方式是堆叠柱状图每个柱子是一个样本每段颜色代表一个分类群qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file metadata.txt \ --o-visualization taxa-bar-plots.qzv在生成的交互式图表里你可以切换分类级别到门Level 2或者属Level 6调整配色按分组排序。如果想自己用 R 画更个性化的图可以在界面里把对应的 csv 数据导出来再丢进 ggplot2。4.5 差异物种分析找到组间显著变化的类群堆叠柱状图能看出“长得不一样”但到底哪些物种在组间显著变化要用 ANCOM 这类差异分析方法。ANCOM 要求数据先加一个伪计数规避零值问题# 为特征表添加伪计数 qiime composition add-pseudocount \ --i-table table.qza \ --o-composition-table comp-table.qza # 在 ASV 水平做差异分析 qiime composition ancom \ --i-table comp-table.qza \ --m-metadata-file metadata.txt \ --m-metadata-column Group \ --o-visualization ancom-Group.qzvANCOM 刚跑完时你可能会有点失望火山图上可能只有零星几个点显著。这很正常ASV 粒度太细检验功效有限。更常用的做法是先按属水平合并特征表再做 ANCOM# 按属水平合并level 6 对应属分类层级 qiime taxa collapse \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --p-level 6 \ --o-collapsed-table table-l6.qza # 对合并后的表加伪计数并运行 ANCOM qiime composition add-pseudocount \ --i-table table-l6.qza \ --o-composition-table comp-table-l6.qza qiime composition ancom \ --i-table comp-table-l6.qza \ --m-metadata-file metadata.txt \ --m-metadata-column Group \ --o-visualization ancom-Group-l6.qzv属水平的结果里差异物种的名字能直接对应到生物学讨论比如哪个属在实验组显著升高再去查文献找机制逻辑就顺了。5. 避坑手册安装和使用中的 5 个高频问题5.1 wget 下载进度条半天不动现象下载软件列表或数据库时进度条一直停在 0%最后超时报错。原因大部分时候是网络连接不稳定也有少数情况是官网服务器响应慢。解决先用wget -c断点续传中断了重新执行同一条命令就能继续如果多次重试还是不行换一个时间再试或者用浏览器直接下载这个文件再手动上传到服务器曲线救国。5.2 conda activate 后提示 CommandNotFoundError现象输入conda activate qiime2-2020.6报错说找不到 activate 命令。原因Miniconda 装完后没有执行conda initShell 还不知道去哪找 Conda 的脚本。解决执行一次~/miniconda3/bin/conda init然后关掉终端重新打开再执行conda activate qiime2-2020.6。如果还不行手动执行source ~/miniconda3/etc/profile.d/conda.sh也能临时解决。5.3 在 Windows 子系统里访问不到 C 盘文件现象在 Ubuntu 子系统里cd /mnt/c/qiime2提示目录不存在。原因路径写错了或者数据不是放在 C 盘根目录。解决先ls /mnt/c/看看 C 盘有没有正确挂载然后确认你的数据在 C 盘哪个目录目录名注意大小写和空格比如C:\qiime2 data对应的是/mnt/c/qiime2 data空格要保留。5.4 DADA2 报错说序列长度太短或引物切完没剩多少现象降噪过程报错提示 reads 长度不足或者处理之后大部分 reads 被丢弃。原因--p-trim-left-f 29 --p-trim-left-r 18切得太多把整条序列切没了或者你的测序读长本来就很短切完引物后只剩几十个碱基。解决先看测序报告的碱基质量图确定引物实际长度不要把质量很差的尾部强行保留。对读长特别短的数据考虑改用单端模式或者降低--p-trunc-len截断长度。5.5 物种注释结果全是 k__Bacteria;p__注释不到属现象taxonomy.qzv里很多 ASV 只能注释到门或纲再往下全是空的。原因参考数据库的覆盖范围有限你的样本里含有数据库里没有的近缘种或者训练分类器时用的引物区域和实际扩增区域不一致。解决优先确认扩增区域和训练分类器时提取的区段是否一致不一致就重新训练区域特异分类器数据库方面可以换成 SILVA 或 Greengenes 其他版本对比看看。注意有些 ASV 注释不到属是客观事实不代表你跑错了。5.6 抽平阈值选了最小值后样本还是被丢现象core-metrics-phylogenetic跑完结果里样本数少了几个。原因某些样本的实际序列数低于你设定的--p-sampling-depth。解决回头看table.qzv的分位数表把阈值改到“除去极端低深度样本后能保留尽量多样本”的数值。注意低于阈值的样本会直接从多样性分析里剔除写论文时要如实说明样本剔除标准。6. 下一步从“跑通流程”到“跑出好结果”第一个流程跑通以后你已经不是那个对着 fastq 文件发懵的新手了。接下来我建议按这个顺序往前走。第一为自己的实验重新训练一个区域特异分类器。数据库部署那一步里extract-reads的引物一定要改成你自己实验的真实引物序列这是提升物种注释精度的性价比最高的一步比换数据库更有效。第二学会从 qzv 交互图表里提取信息。每次分析完不要只看一眼 p-value把table.qzv的分位数表、alpha-rarefaction.qzv的稀释曲线都点开看看这些图能帮你判断数据质量也能帮你更合理地设置抽平阈值。第三把特征表和物种注释结果导出来学一点 R 语言用 ggplot2 重画堆叠柱状图和 PCoA 图。QIIME 2 自带的图表适合快速查看但发表级的图通常要自己用 R 精细调排版。第四有余力可以了解功能预测分析比如用 PICRUSt2 根据 16S 序列预测样本的潜在功能基因虽然不是实测但能为课题提供假设。我自己刚开始学的时候跑完 DADA2 就直接把表拿去画图结果发现抽平阈值根本没调组间比较全是假的差异白白浪费了两周。从那以后我每次拿到测序数据都强制自己先看特征表统计确定阈值再往下走顺序不乱结果才不乱。希望这套流程和踩坑经验能帮到你剩下的就是下载数据、打开终端、一步步执行跑起来你就知道它没那么可怕。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
