★ FEATURED ARTICLE 简介围绕微生物组16S rRNA基因扩增子分析展开的QIIME 2实操文档面向微生物组研究入门者、缺乏生物信息学背景的实验人员及需要快速搭建分析流程的研究生。内容从Linux服务器、Windows 10子系统Ubuntu与Mac等安装环境选择讲起覆盖数据导入、特征表生成、α与β多样性、物种组成、差异物种分析和可视化等标准环节并对中间步骤与关键参数进行解读帮助读者理解参数优化和结果解释。压缩包仅含1个docx文件约578KB篇幅集中、便于随查随用。文档还总结软件安装与使用中的常见问题及解决策略并给出配套视频、分析代码、测序数据和预期结果线索便于复现完整流程、对照排查报错并形成可重复的分析记录。已有2799人学习适合作为QIIME 2上手与流程规范的参考材料。1. 拿到一摞双端 fq 之后QIIME 2 在 16S 扩增子流程里到底管哪几段测序公司把一摞双端 fq 甩过来很多人的第一反应是装个 QIIME 2 一路往下跑。真跑起来才发现卡住的地方往往不是 16S rRNA 基因扩增子测序数据分析本身而是环境、路径、格式这三件杂事Windows 下装不上、manifest 少一列、分类器和实验引物对不上。这套 QIIME 2 流程把安装、数据导入、特征表生成、Alpha/Beta 多样性、物种注释到差异分析串成一条可复现的线配套给了元数据模板、批处理脚本和示例数据。它适合两类人一是没有生物信息背景、但要独立出微生物组结果的研究生二是跑过 q2 但每次都在抽平阈值和分类器上凭感觉调的分析人员。下面按我实际复现的顺序拆先环境后数据最后落到参数和翻车点。2. 装不上就别谈分析四条部署路径的成本对比与 Conda 环境配置2.1 为什么 Windows 直装必然失败四条路的成本对比QIIME 2 的 conda 包只发布 linux-64 和 osx-64 两个平台没有 win-64。这不是文档写得含糊是根本没编 Windows 版所以任何在 Windows 命令提示符里敲conda install qiime2的尝试都会在解析依赖时失败。绕开的方式只有一条给 Windows 借一个 Linux 环境。原文列了四条路我把它们的实际成本摊开部署方式适合数据量装起来的耗时主要问题Linux 服务器大几十 GB 以上中需要会 ssh 和传文件账号、磁盘配额、排队Win10 子系统 Ubuntu 20.04小几 GB 以内低应用商店点一下路径在 /mnt/c 下 IO 偏慢内存受宿主限制VirtualBox 虚拟机小高装系统加配共享目录效率最低DADA2 多线程跑不满Mac 直装小中依赖兼容问题多出问题不好搜选择逻辑很简单数据超过 10 GB、或者要跑多个项目直接上 Linux 服务器只是复现一次教学数据Windows 子系统足够虚拟机只在既没有服务器又不想动系统设置时用跑 DADA2 的时候你会后悔。Mac 不是不能跑是踩坑之后可参考的经验少新手不建议。提示子系统里访问 Windows 盘符的路径是/mnt/c/...QIIME 2 对路径里的中文和空格很敏感工作目录一律用纯英文。2.2 从 Miniconda 到 qiime2 环境的完整命令链不建虚拟环境直接往 base 里装是后面所有版本冲突的源头。QIIME 2 依赖的 Python 版本、numpy、biom 都是锁死的跟其他生信软件共用一个环境早晚会互相打架。正确做法是按官方发布的 yml 清单单独建一个环境环境名里带上版本号以后同时留两三个版本互不干扰。# 1. 安装 Miniconda3已经装过 conda 的跳过这一步 wget -c https://repo.anaconda.com/miniconda/Miniconda3-latest-Linux-x86_64.sh bash Miniconda3-latest-Linux-x86_64.sh # 把 conda 写进 shell 配置执行完要重开一个终端才生效 ~/miniconda3/bin/conda init # 2. 下载 QIIME 2 2020.6 的环境清单几十个包的版本都锁在里面 wget -c https://data.qiime2.org/distro/core/qiime2-2020.6-py36-linux-conda.yml # 3. 按清单建环境名字带版本号方便以后并存 conda env create -n qiime2-2020.6 --file qiime2-2020.6-py36-linux-conda.yml # 4. 每次开工先激活提示符前面出现环境名才算成功 conda activate qiime2-2020.6 qiime --help这几行里最容易被忽略的是conda init之后的重新登录。它只改了~/.bashrc当前 shell 并不会立刻认识conda activate所以紧接着敲激活命令会报 command not found很多人以为装失败了其实是没重开终端。另一个细节是环境名。写成qiime2-2020.6而不是笼统的qiime2好处是半年后你想复现同一批结果conda env list一眼能看出当时用的哪版不会出现“我记得装过但不知道装的是哪年”的黑匣子。网络慢是这条链上最常见的阻塞点。yml 里几十个包全部从默认通道拉国内直连经常卡在 solving environment 或者下载到一半超时。我的习惯是先给 conda 配一个就近的通道镜像再跑conda env create能少等半小时如果公司网络限制严格就让能上网的机器把包缓存目录打包拷过来用--offline建环境。这一步没有玄学纯粹是网络问题别怀疑 yml 损坏。2.3 数据库落地GreenGenes 导入与分类器训练环境好了下一个必做项是物种分类数据库。QIIME 2 本身不带数据库注释那一步靠的是一个训练好的 naive-bayes 分类器文件。训练分两步先把参考序列和参考分类表导成 qza再用fit-classifier-naive-bayes训练。原文用的是 GreenGenes 13_8 的 99% 聚类版本这个库老是老2013 年定版很多新命名的属查不到但格式规整、体积可控仍然是流程默认。# 下载并解压 GreenGenes 13_899% 聚类版本 wget -c ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_8_otus.tar.gz tar -zxvf gg_13_8_otus.tar.gz # 参考序列导入类型是 FeatureData[Sequence] qiime tools import \ --type FeatureData[Sequence] \ --input-path gg_13_8_otus/rep_set/99_otus.fasta \ --output-path 99_otus.qza # 分类表导入格式名必须写 HeaderlessTSVTaxonomyFormat 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这里有个参数值得单独说--input-format HeaderlessTSVTaxonomyFormat。GreenGenes 的 taxonomy 文件是两列纯文本没有表头第一列是 OTU ID第二列是k__Bacteria; p__...这种分号串。不写这个格式名QIIME 2 会按带表头的 TSV 去解析第一行数据直接被当成列名吃掉后面注释全错但不会报错属于典型的静默失败。全长分类器能覆盖所有 16S 区段代价是参考序列全量参与训练模型大、注释慢。如果实验只扩了某个可变区可以再训一个区段特异的分类器用extract-reads按实际引物把参考序列截取到目标区段再用截取后的序列训练。原文示例用的是 V5-V7 区引物正向AACMGGATTAGATACCCKG反向ACGTCATCCCCACCTTCC注意这里用的是 IUPAC 简并碱基K 代表 G/TM 代表 A/C。# 按实验引物截取参考序列的目标区段引物方向必须和测序一致 time qiime feature-classifier extract-reads \ --i-sequences 99_otus.qza \ --p-f-primer AACMGGATTAGATACCCKG \ --p-r-primer ACGTCATCCCCACCTTCC \ --o-reads ref-seqs.qza # 用截取后的序列再训练文件小、注释更聚焦 time 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两次训练的时间差得不多原文全长约半小时区段特异约 8 分钟但产物大小和注释速度差很多。区段特异的分类器还有个隐性好处参考序列和你的 reads 覆盖同一段区域比对上的位点更集中低置信度的注释会少一些。代价是引物必须和实验完全一致写反一个碱基截取出来的序列集可能就是空的后面训练还能跑完只是分类器里没什么东西注释结果一片k__Bacteria。这个坑在第 5 章还会展开。如果项目对注释精度要求高GreenGenes 覆盖不到的新属可以考虑换 SILVA 或 GTDB但换库不是改个路径那么简单序列版本、分类层级尤其 GTDB 的域级命名和训练参数都要跟着改换之前先用一小批数据对比两套注释结果别在正式数据上试。3. 从双端 fq 到特征表manifest 三列、导入格式和 DADA2 的四个参数3.1 元数据表与 manifest两个最容易被写错的纯文本QIIME 2 的整个流程靠两个纯文本串起来元数据表和 manifest。元数据表描述样本第一列必须是#SampleID后面跟分组、批次这些列manifest 描述文件位置固定三列sample-id、forward-absolute-filepath、reverse-absolute-filepath。这两个文件写错后面所有命令都会连锁失败而且报错信息往往指向别处。# 建工作目录并激活环境 mkdir -p qiime2 cd qiime2 conda activate qiime2-2020.6 # 下载示例元数据第一列 #SampleID后面是 Group 等分组列 wget -c http://example.org/path/metadata.txt # 按元数据里的批次号和样本号批量下载双端数据并改名 mkdir -p seq awk {system(wget -c ftp://download.example.org/gsa/$5/$6/$6_f1.fq.gz -O seq/$1_1.fq.gz)} (tail -n2 metadata.txt) awk {system(wget -c ftp://download.example.org/gsa/$5/$6/$6_r2.fq.gz -O seq/$1_2.fq.gz)} (tail -n2 metadata.txt) # 检查文件大小0 字节或明显偏小的要重下 ls -lsh seq # 生成 manifest路径用 $PWD 拼成绝对路径 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 manifestawk在这两段里干了三件事NR1跳过表头$1、$5、$6取元数据表里的样本名、批次号、样本编号system()把拼好的wget命令交给 shell 执行。(tail -n2 metadata.txt)是进程替换等价于把去掉表头的内容喂给 awk避免表头那一行也触发一次下载。manifest 里的路径必须是绝对路径这一点没有商量余地。QIIME 2 导入时会把这个路径解析成内部引用如果你写的是seq/S1_1.fq.gz在qiime2目录里导入能成功换个目录再跑下游命令就找不到文件了。用$PWD拼是最省事的写法不用手工敲几十行。注意样本名里不要出现下划线以外的特殊字符也别用纯数字。QIIME 2 内部把样本名当字符串键用1和01在某些环节会被当成不同样本导入后对不上元数据。3.2 导入格式选型Phred33 与压缩格式的耗时差异导入这一步的核心是选对--input-format。双端带质量值的数据用PairedEndFastqManifestPhred33V2其中 Phred33 指质量值编码方式老 454 数据才用 Phred64末尾的 V2 对应新版 manifest 的列名V1 的列名是forward-filepath那套写错会直接报列不匹配。# 双端带质量值导入格式名末尾的 V2 对应新版 manifest 列名 qiime tools import \ --type SampleData[PairedEndSequencesWithQuality] \ --input-path manifest \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33V2耗时差异值得单独提一句。原文实测同样一批数据未压缩的 fq 导入要 7 分钟换成 fq.gz 只要 34 秒。原因不复杂QIIME 2 内部本来就按压缩流读写 qza未压缩的文本要先读进来再压一遍等于多做一次全量 IO而 gz 文件是边解压边读解压发生在管道里不落盘。所以测序公司给什么格式就用什么格式不要为了“看起来清爽”先 gunzip 一遍那是纯粹给自己加时间。导入完成后demux.qza里装的是所有样本的序列和质量信息。如果原始数据是混池测序没拆样本QIIME 2 处理不了得先让测序方提供拆分后的单样本 fastq或者用早期 QIIME 的拆分脚本按 barcode 拆完再导入。这一步没有捷径barcode 表错一个碱基拆出来的样本就是错的。3.3 DADA2 去噪trim-left 与 trunc-len 怎么定特征表由 DADA2 去噪产生这一步的参数直接决定后面所有分析的质量。核心参数四个两端各切掉多少trim-left、截断到哪trunc-len、线程数、最大期望错误数。# 去噪输出特征表、代表序列、每步去噪统计三件套 time 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 # 把 dada2 结果接进主流程统一命名方便下游命令 cp dada2-table.qza table.qza cp dada2-rep-seqs.qza rep-seqs.qza参数含义怎么定--p-trim-left-f/r从 5 端切掉的碱基数引物加接头的长度本例正向 29、反向 18--p-trunc-len-f/r在 3 端截断到的位置质量骤降点设 0 表示不截断--p-n-threads线程数按机器核数给8 线程约 77 分钟--p-max-ee-f/r单条序列允许的最大期望错误数默认 2质量差的数据可放宽到 3 到 4trim-left的依据是引物长度。正向引物 29 个碱基、反向 18 个测序读长里包含这些人工序列留着会干扰后续的比对和建树所以整段切掉。注意这里是按固定长度切如果你的引物长度和示例不同必须换成自己的数字照抄 29/18 等于把一部分真实序列也切没了。trunc-len设 0 表示不截断前提是 reads 3 端质量足够好。如果质量曲线在某个位置之后明显下滑就得设一个截断值把低质量尾巴切掉此时trunc-len必须大于trim-left否则切完长度为负DADA2 直接报错。判断依据是demux.qzv里的质量分布图看中位数曲线从哪开始掉。线程数的加速不是线性的。原文实测96 线程 34 分钟、24 线程 44 分钟、8 线程 77 分钟、1 线程 462 分钟。从 1 线程到 8 线程省了六倍时间从 24 到 96 只省了十分钟。所以给 8 到 16 线程是性价比最高的区间剩下的核留给系统和其他任务没必要为了快十分钟把机器跑满。提示denoising-stats.qza要单独看一眼里面记录了每个样本输入多少条、过滤后剩多少条、非嵌合占比多少。某个样本过滤后只剩个位数八成是 trim/trunc 参数和它的读长不匹配别急着往下跑。3.4 特征表统计与抽平阈值的第一次决策特征表出来先别急着做多样性第一步是汇总看清每个样本到底剩多少序列。抽平阈值这个数字一旦定错后面的 Alpha/Beta 全是白算而它只能从汇总结果里读出来。# 汇总特征表每样本测序量、特征数、分位数都在这张表里 qiime feature-table summarize \ --i-table table.qza \ --o-visualization table.qzv \ --m-sample-metadata-file metadata.txt # 代表序列一览每条序列可点开跳 BLAST qiime feature-table tabulate-seqs \ --i-data rep-seqs.qza \ --o-visualization rep-seqs.qzvqza 和 qzv 本质上都是 zip 包qza 里是数据qzv 里是一个网页报告加图表资源。不想装浏览器插件的话直接unzip table.qzv也能看到里面的 html。所有 qzv 都丢到 view.qiime2.org 这个在线查看器里看交互式的分位数表、柱状图、稀释曲线都在网页里。抽平阈值从table.qzv的分位数表读。原文示例的最小值 27060、第一分位 28581.75、中位 30867、最大 34663分布很均匀所以直接取最小值。规则可以记成这样最小值离第一分位不远取最小值保留全部样本最小值明显低于第一分位比如差出一倍就在两者之间取一个数宁可丢掉一两个测序量特别低的样本也别让整体深度被拉到很低。阈值太低会丢掉大量序列信息太高会把样本砍到组内重复不够。顺便说个老标准的问题。抽平阈值下限 1000 是 454 测序时代留下的那时候一个样本几千条序列就算不错。Illumina 通量高一个量级现在最小值一般不低于 5000推荐 1 万以上在样本量允许的前提下越大越好。还在用 1000 当阈值等于主动扔掉九成数据。4. 多样性与物种组成建树、抽平、显著性检验的先后顺序不能乱4.1 建树为什么必须排在多样性之前Alpha 多样性里的 faith_pd、Beta 多样性里的 unifrac 系列都依赖一棵系统发育树。树没建好这两个指数算出来没有意义。所以顺序是固定的先建树再算多样性不能颠倒。# MAFFT 比对 fasttree 建树一次输出对齐、屏蔽对齐、无根树、有根树 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-tree unrooted-tree.qza \ --o-rooted-tree rooted-tree.qza这条命令内部做了两件事MAFFT 做多序列比对fasttree 用近似最大似然建树。输出四个文件里masked-aligned-rep-seqs.qza是把高变区屏蔽之后的比对结果建树用的是它因为高变区比对不可靠留着会把树带偏rooted-tree.qza是定了根的有根树下游core-metrics-phylogenetic要的就是这个。fasttree 是近似算法不是精确最大似然换来的是速度。几万条代表序列分钟级能出树换成 RAxML 跑精确解同一个数据量可能要几小时。对 16S 扩增子这个精度要求fasttree 够用别在这上面纠结。4.2 core-metrics 一次产出四类距离矩阵多样性计算这一步QIIME 2 把它打包成一条命令一次产出四类 Beta 距离矩阵和四个 Alpha 指数省得你挨个跑。# 系统发育多样性一次算完四种距离矩阵和四个 alpha 指数都在输出目录 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输出文件内容*_vector.qza各 Alpha 指数faith_pd、shannon、observed_features、evenness*_distance_matrix.qza各 Beta 距离矩阵*_pcoa_results.qzaPCoA 坐标*_emperor.qzv交互式 PCoA 图rarefied_table.qza抽平后的特征表--p-sampling-depth就是 3.4 节定下来的抽平阈值。抽平是随机无放回抽样随机种子默认固定同样的输入和阈值跑两次结果一致要换种子可以加--p-seed。低于这个深度的样本会被整条丢弃不会补零也不会警告所以跑完先看输出目录里剩了几个样本和元数据表对一下数量。四类距离矩阵的用途不一样bray_curtis只看丰度不看亲缘jaccard只看有无不看丰度weighted_unifrac同时考虑丰度和亲缘unweighted_unifrac只考虑有无和亲缘。写文章时选哪个要跟科学问题对上四个全放上去反而显得没想清楚。4.3 Alpha 显著性与稀释曲线Alpha 多样性回答的是“单个样本内部有多丰富”组间比较靠alpha-group-significance默认用 Kruskal-Wallis 检验加两两比较。# 单个 alpha 指数做组间检验 indexobserved_features qiime diversity alpha-group-significance \ --i-alpha-diversity core-metrics-results/${index}_vector.qza \ --m-metadata-file metadata.txt \ --o-visualization core-metrics-results/${index}-group-significance.qzv # 稀释曲线max-depth 取样本最大测序量异常值时取第三分位与最大值之间 qiime diversity alpha-rarefaction \ --i-table table.qza \ --i-phylogeny rooted-tree.qza \ --p-max-depth 34663 \ --m-metadata-file metadata.txt \ --o-visualization alpha-rarefaction.qzvobserved_features就是老版本的observed_otus数的是样本里出现了多少个特征对低丰度特征最敏感。shannon同时考虑丰富度和均匀度faith_pd把亲缘关系算进去evenness单独看均匀度。四个指数经常给出方向不一致的结论这很正常因为它们衡量的东西不同选一个和你的假设最贴切的当主指标其余当辅助。结果表里同时有 p-value 和 q-value。p-value 是单次比较的显著性q-value 是 FDR 多重检验校正之后的。三组两两比较就是三次检验只看 p-value 会把假阳性放进来正文里报 q-value 更稳妥。原文示例里 KO 对 WT 的 p-value 是 0.010139q-value 是 0.030417校正后仍然显著OE 对 WT 的 p 是 0.054241校正后 q 到 0.081362就只能说趋势。稀释曲线的--p-max-depth取样本最大测序量本例 34663。如果最大值是个异常离群点就在第三分位和最大值之间取一个否则横轴被拉得很长曲线全挤在左边看不清。这条图的作用是验证抽平阈值是否合理曲线在阈值附近已经走平说明再加深也不会有新特征出现抽平没丢信息。4.4 Beta 显著性与 PCoA 出图Beta 多样性比较的是样本之间的距离组间差异用 PERMANOVA 检验可视化用 PCoA。# 指定距离矩阵和分组列pairwise 打开两两比较 distanceweighted_unifrac columnGroup qiime diversity beta-group-significance \ --i-distance-matrix core-metrics-results/${distance}_distance_matrix.qza \ --m-metadata-file metadata.txt \ --m-metadata-column ${column} \ --o-visualization core-metrics-results/${distance}-${column}-significance.qzv \ --p-pairwisePERMANOVA 是基于置换的检验默认置换 999 次特征数多、分组多的时候很耗时所以先指定单一分组列跑别一上来就把所有列都传进去。结果表里三个数字要看懂pseudo-F是组间方差与组内方差的比越大说明分组解释得越好permutations是实际做的置换次数p-value是置换分布里大于等于观测值的比例。原文示例三组两两比较的 pseudo-F 在 2.64 到 4.11 之间p 值都小于 0.05。PCoA 图在weighted_unifrac_emperor.qzv里打开后把颜色映射到分组列再调点的透明度、大小、形状把重叠的样本分开。图上括号里是各轴的解释率Axis1 41.26%、Axis2 20.09%、Axis3 10.89%前两轴加起来六成说明主要变异集中在头两轴分组在 Axis1 上分得开就说明组间差异是主要来源。定好样式之后用右上角 Save plot 导出 SVGLabel会得到主图和图例两个矢量文件投稿拼接时比 PNG 好用得多。4.5 物种注释与堆叠柱状图特征表有了、树有了最后一步是把特征对上物种名。# 用训练好的分类器做注释耗时取决于序列数和分类器大小 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 # 堆叠柱状图可切换分类级别、配色、排序方式 qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file metadata.txt \ --o-visualization taxa-bar-plots.qzvtaxonomy.qzv是三列表特征 ID、分类注释、置信度。置信度低于 0.7 的注释要谨慎对待尤其是只到k__Bacteria就止步的那些多半是参考库里没有对应序列不是数据本身有问题。把这类特征的比例统计一下超过两成就该考虑换数据库或者重训分类器。堆叠柱状图里 Level 2 是门水平Level 6 是属水平。默认配色和排序往往不好看可以在网页里把柱宽调到中等、切到 Level 2、配色换成深色系、按分组排序再按某个优势门降序排列这样组间差异一眼能看出来。除了导 SVG也建议把表格数据导成 csv用 R 的 ggplot2 重画配色和图例位置自由度更大审稿人提修改意见时不用从头再来。5. 避坑排查五类高频翻车现场的现象、原因与解决5.1 安装与环境类现象一conda env create卡在 solving environment或者下载到一半超时。原因yml 里几十个包全部从默认通道拉国内直连不稳定conda 的依赖求解器在老版本上本来就慢包一多容易假死。解决先给 conda 配就近的通道镜像再跑如果网络彻底不通在有网的机器上conda env create成功后把pkgs缓存目录整个拷过来加--offline建环境。判断是不是真卡死看缓存目录里文件大小有没有在涨五分钟没变化就是断了。现象二conda activate qiime2-2020.6报 command not found。原因conda init只改了~/.bashrc当前 shell 还没重新加载或者默认 shell 是 zshinit 写进了 bash 的配置文件。解决source ~/.bashrc或者直接退出重登用 zsh 的话跑一次conda init zsh。还有一种情况是用了sudo装的 miniconda普通用户没权限读这种只能重装到用户目录下。5.2 数据与参数类现象三qiime tools import报 manifest 列名不匹配或者导入成功但样本数为零。原因--input-format写的是 V2manifest 用的却是 V1 那套列名或者路径写成了相对路径导入时能解析换目录就失效。解决manifest 表头严格写成sample-id、forward-absolute-filepath、reverse-absolute-filepath三列用制表符分隔路径用$PWD拼绝对路径。导入后立刻qiime tools peek demux.qza看一眼样本数和元数据行数对不上就别往下跑。现象四extract-reads跑完ref-seqs.qza几乎是空的后面训练照样跑完注释结果一片k__Bacteria。原因引物方向写反或者填了反向互补序列。extract-reads默认只按给定方向在参考序列上匹配方向反了匹配不到它不会报错只会安静地返回一个空集。解决确认引物是 5 到 3 方向、和实验用的完全一致拿不准时先用qiime feature-classifier extract-reads跑一小批用qiime tools peek看输出里序列条数个位数就是匹配失败。这个坑最阴的地方在于流程不中断等注释结果出来才发现中间几小时全白跑。现象五抽平之后样本数变少组内重复不够做统计。原因--p-sampling-depth取高了低于阈值的样本被整条丢弃而且不报警。解决跑core-metrics-phylogenetic之前先看table.qzv的分位数表把阈值定在最小值和第一分位之间如果某个组的样本本来就少宁可把阈值往下调一点牺牲深度保住重复数。跑完立刻ls core-metrics-results/*_vector.qza核对样本数和元数据表逐行对。提示这五类问题里现象四和现象五都属于“不报错但结果错”比直接崩溃更难查。养成每步跑完看一眼输出规模的习惯比事后回溯省时间。6. 差异分析与结果复核ANCOM 的 W 值边界和 qzv 的二次利用差异分析这一步原文用的是 ANCOM。它和 DESeq2 的思路不一样ANCOM 假设“只有少数特征在组间有差异”基于特征间的对数比值做检验对组成型数据的处理更稳。代价是它对输入有要求特征表里不能有零否则对数比值为负无穷所以必须先加伪计数。# ANCOM 要求输入带伪计数否则 log 比值为负无穷 qiime composition add-pseudocount \ --i-table table.qza \ --o-composition-table comp-table.qza # 指定分组列做差异检验置换耗时与特征数相关 columnGroup time qiime composition ancom \ --i-table comp-table.qza \ --m-metadata-file metadata.txt \ --m-metadata-column ${column} \ --o-visualization ancom-${column}.q p a hrefhttps://download.csdn.net/download/weixin_35823067/86347840 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p
阅读完成 · 觉得有帮助？
