簡介面向微生物組研究者和生物信息學(xué)入門人員的一份技術(shù)文檔聚焦于16S rRNA基因擴(kuò)增子測序數(shù)據(jù)的QIIME 2分析流程。文檔基于Bio-protocol論文格式整理系統(tǒng)介紹軟件在Linux服務(wù)器、Windows 10子系統(tǒng)、VirtualBox虛擬機(jī)及Mac等多種環(huán)境下的安裝部署與數(shù)據(jù)庫配置并詳細(xì)拆解從測序數(shù)據(jù)導(dǎo)入、特征表生成、α/β多樣性分析到物種組成與差異物種分析、可視化的完整標(biāo)準(zhǔn)操作文中穿插參數(shù)優(yōu)化思路與常見報(bào)錯(cuò)的排錯(cuò)方法并涉及GreenGenes分類數(shù)據(jù)庫部署、Miniconda環(huán)境管理等實(shí)用細(xì)節(jié)幫助缺乏生信背景的科研人員縮短上手周期。資源為單個(gè)docx文檔大小578KB包含圖文說明、操作要點(diǎn)以及配套視頻指引、分析代碼與參考數(shù)據(jù)來源可結(jié)合原文鏈接對照學(xué)習(xí)。當(dāng)前已有2793人學(xué)習(xí)瀏覽適合需要系統(tǒng)性掌握QIIME 2主流分析流程的微生物生態(tài)學(xué)與醫(yī)學(xué)研究者參考。1. 用 QIIME 2 分析 16S rRNA 擴(kuò)增子數(shù)據(jù)為什么它值得成為你的第一套流程拿到公司返回的 16S rRNA 基因擴(kuò)增子測序數(shù)據(jù)之后多數(shù)人不是死在算法選型上而是死在第一步一堆 fastq.gz 文件名長得不一樣雙端怎么拼低質(zhì)量區(qū)切多少引物要不要去掉切完之后兩個(gè) read 還剩多少重疊區(qū)域。如果每一步都靠手搓臨時(shí)腳本換一個(gè)項(xiàng)目就要重來一遍而且后續(xù)審稿人問你要“可重復(fù)的分析流程”時(shí)根本拿不出東西。QIIME 2 之所以在微生物組分析里被當(dāng)默認(rèn)選項(xiàng)就是因?yàn)樗言紲y序數(shù)據(jù)到特征表、分類學(xué)注釋、alpha/beta 多樣性這套鏈條統(tǒng)一成一套帶類型檢查的命令行框架背后跑的是 DADA2、sklearn、MAFFT/FastTree 這類成熟工具但暴露給用戶的是 q2cli 里統(tǒng)一風(fēng)格的參數(shù)。它能幫你把“我有哪些菌、組間差多少、誰導(dǎo)致差異”這三件事一次做完還能留出完整的日志和中間產(chǎn)物。這篇文章沒有太多教科書概念按我實(shí)際分析一份 16S 擴(kuò)增子數(shù)據(jù)時(shí)踩過的路徑來寫先講數(shù)據(jù)導(dǎo)入和環(huán)境搭建再把 DADA2 去噪、分類器和多樣性分析的關(guān)鍵參數(shù)說透最后列出最容易翻車的五個(gè)排查現(xiàn)場。新手可以跟著命令一步一步跑熟手可以直接跳到參數(shù)邊界和避坑章節(jié)對照自己的報(bào)錯(cuò)找答案。2. 數(shù)據(jù)輸入與最小環(huán)境導(dǎo)入前把測序數(shù)據(jù)裝成 .qza 工件2.1 先認(rèn)清三種常見格式別在導(dǎo)入階段就選錯(cuò)輸入類型QIIME 2 不像傳統(tǒng)腳本那樣直接讀 fastq.gz它要求所有輸入先變成自己的二進(jìn)制工件格式 .qza。這意味著第一步不是跑算法而是把測序公司交付的原始文件導(dǎo)入成 QIIME 2 認(rèn)識的類型。常見的到貨格式有三種我按實(shí)際使用頻率從高到低排格式典型特征QIIME 2 導(dǎo)入類型適合場景雙端 manifest一個(gè) csv 標(biāo)明每個(gè)樣本的正反向文件路徑SampleData[PairedEndSequencesWithQuality]幾乎所有新項(xiàng)目文件命名不規(guī)范時(shí)最穩(wěn)妥Casava OneStep文件名按樣本ID_S1_L001_R1_001.fastq.gz規(guī)律排列SampleData[PairedEndSequencesWithQuality] Casava 格式公司嚴(yán)格按 Illumina 下機(jī)命名時(shí)單端 manifest只有單端 fastqSampleData[SequencesWithQuality]已經(jīng)合并或?qū)ｉT做單端 16S 的舊項(xiàng)目對于 16S 擴(kuò)增子我默認(rèn)推薦雙端 manifest。原因很簡單雙端 manifest 不依賴文件名規(guī)律只要 CSV 里寫清楚文件路徑就能導(dǎo)入哪怕是公司把 R1 和 R2 文件放在不同子目錄也只要把路徑寫對就行。一個(gè)標(biāo)準(zhǔn)的雙端 manifest 長這樣sample-id,forward-absolute-filepath,reverse-absolute-filepath Ctrl-1,/home/user/project/data/Ctrl-1_S1_L001_R1_001.fastq.gz,/home/user/project/data/Ctrl-1_S1_L001_R1_002.fastq.gz Treat-1,/home/user/project/data/Treat-1_S1_L001_R1_001.fastq.gz,/home/user/project/data/Treat-1_S1_L001_R1_002.fastq.gz Ctrl-2,/home/user/project/data/Ctrl-2_S1_L001_R1_001.fastq.gz,/home/user/project/data/Ctrl-2_S1_L001_R1_002.fastq.gz有幾個(gè)細(xì)節(jié)直接影響導(dǎo)入成敗。第一路徑必須是絕對路徑不能寫~/project/data/這種帶波浪號的寫法也不能寫相對路徑第二CSV 編碼必須是 UTF-8 無 BOMWindows 下用記事本編輯出來的 CSV 經(jīng)常帶 BOM導(dǎo)入時(shí)會(huì)報(bào)“第一列 sample-id 不存在”處理方式是用dos2unix轉(zhuǎn)一次或者用 VS Code 重新保存第三列名一個(gè)字母都不能錯(cuò)forward-absolute-filepath和reverse-absolute-filepath是固定的別寫成forward_path或R1。最后樣本 ID 里盡量只保留字母、數(shù)字、下劃線和連字符不要有空格不要用中文否則后面畫 PCoA 圖時(shí)會(huì)出現(xiàn)無數(shù)個(gè)你意想不到的坑。另外要確認(rèn)一個(gè)技術(shù)細(xì)節(jié)Phred 質(zhì)量編碼。絕大多數(shù)公司現(xiàn)在交付的都是 Phred33但如果你拿到的測序數(shù)據(jù)來自比較舊的平臺或者公司沒有做轉(zhuǎn)換也有可能是 Phred64。判斷方法很簡單用zcat解壓一個(gè) fastq.gz看質(zhì)量行字符的 ASCII 碼范圍。質(zhì)量行里如果出現(xiàn)大量B、C、D這類大于I的字母大概率是 Phred64如果幾乎都是#、5、9、?這類字符基本就是 Phred33。判斷錯(cuò)了導(dǎo)入階段不一定報(bào)錯(cuò)但后續(xù) DADA2 會(huì)把高質(zhì)量和低質(zhì)量完全搞反得到的結(jié)果整批報(bào)廢。我一般會(huì)在導(dǎo)入前用seqkit stats或者直接解壓看一個(gè)文件確認(rèn)無誤再寫 manifest。2.2 用 conda 搭建最小環(huán)境隔離是避免插件沖突的唯一出路QIIME 2 的安裝沒有太多捷徑官方文檔會(huì)根據(jù)發(fā)布季提供對應(yīng)的 conda 環(huán)境文件常見做法是下載那個(gè).yml文件然后執(zhí)行conda env create -n qiime2 --file qiime2-2024.xx.yml。為什么不推薦自己conda install qiime2因?yàn)?QIIME 2 的插件依賴互相之間版本要求很苛刻比如 q2-dada2、q2-feature-classifier、q2-phylogeny 各自鎖定了不同的 numpy / pandas / scikit-learn 版本手動(dòng)安裝很容易出現(xiàn)“主程序能跑但插件加載失敗”的尷尬局面。官方打包的 yml 文件已經(jīng)把這一套依賴鎖好是最省事的路。如果你實(shí)在找不到合適的 yml也可以通過 mamba 從 conda-forge 和 bioconda 頻道創(chuàng)建基礎(chǔ)環(huán)境。命令類似mamba create -n qiime2 -c conda-forge -c bioconda qiime2 q2cli q2-dada2 q2-feature-classifier q2-phylogeny但頻道優(yōu)先級必須對齊否則裝出來的插件版本可能對不上。這里我想強(qiáng)調(diào)一個(gè)經(jīng)驗(yàn)QIIME 2 不要裝在 conda 的 base 環(huán)境里也不要和 R、Python 數(shù)據(jù)分析環(huán)境混裝。我之前在一個(gè)項(xiàng)目里圖省事把 QIIME 2 和 R 的很多包裝在同一個(gè)環(huán)境結(jié)果某次pip install 某個(gè)包把環(huán)境里的 scikit-learn 升級了之后classify-sklearn直接報(bào)ModuleNotFoundError最后只能重建環(huán)境重跑。這個(gè)教訓(xùn)后來成了我的默認(rèn)習(xí)慣QIIME 2 一個(gè)環(huán)境R 統(tǒng)計(jì)一個(gè)環(huán)境互不干擾。創(chuàng)建好環(huán)境后先做一次完整性檢查conda activate qiime2 qiime --help qiime infoqiime info會(huì)列出當(dāng)前環(huán)境識別到的所有插件。如果這里缺了 diversity、dada2 或 feature-classifier說明安裝源或依賴有問題不要急著往后跑先把環(huán)境補(bǔ)好。安裝過程如果特別慢常見原因是 conda 默認(rèn)源在國外可以考慮給 conda 配一個(gè)國內(nèi)鏡像源但注意只改 conda 的 channel 配置不要去動(dòng) QIIME 2 內(nèi)部的依賴。2.3 導(dǎo)入為工件用 qiime tools import 生成 demux.qza準(zhǔn)備就緒后把 manifest 里指向的 fastq.gz 導(dǎo)入成 QIIME 2 的雙端測序工件這一步用的是qiime tools import。我的標(biāo)準(zhǔn)命令是這樣的qiime tools import \ --type SampleData[PairedEndSequencesWithQuality] \ --input-path manifest.csv \ --output-path demux.qza \ --input-format PairedEndFastqManifestPhred33V2--type聲明這是一份雙端帶質(zhì)量的測序數(shù)據(jù)。如果項(xiàng)目是單端類型要改成SampleData[SequencesWithQuality]導(dǎo)入格式也要換成SingleEndFastqManifestPhred33V2。--input-format里的Phred33對應(yīng)質(zhì)量編碼如果你的數(shù)據(jù)是 Phred64就要寫成PairedEndFastqManifestPhred64V2。--output-path生成的 demux.qza 是后面所有質(zhì)控和去噪步驟的輸入。導(dǎo)入完成后需要立刻做一次“看圖”操作驗(yàn)證導(dǎo)入是否成功同時(shí)為下一步的截?cái)鄥?shù)做準(zhǔn)備qiime demux summarize \ --i-data demux.qza \ --o-visualization demux.qzvdemux.qzv 是一個(gè)可視化文件用 QIIME 2 View 打開后能看到每個(gè)樣本的序列數(shù)和前幾條序列的質(zhì)量分布。這一步我不只是看一下序列數(shù)更關(guān)鍵的是確認(rèn)樣本數(shù)量和 metadata 對得上以及雙端 read 的長度分布是否正常。如果導(dǎo)入時(shí)路徑寫錯(cuò)或文件損壞這里會(huì)直接暴露。還有個(gè)小技巧查看 demux.qzv 時(shí)留意樣本名是否帶有多余空格有些公司交付的 sample-id 里會(huì)偷偷帶一個(gè)換行符或空格后期做樣本分組統(tǒng)計(jì)時(shí)會(huì)變成“兩個(gè)不同樣本”這種問題越早發(fā)現(xiàn)越容易處理。3. 去噪與質(zhì)控DADA2 是這套流程里最值得調(diào)參數(shù)的一步3.1 先解讀質(zhì)量圖再定截?cái)辔恢脙蓚€(gè)人的差距從這一步拉開DADA2 的核心思想不是把序列按相似度聚成 OTU而是通過誤差模型把測序錯(cuò)誤從真實(shí)序列中分離出來得到單核苷酸分辨率的 ASV。這個(gè)目標(biāo)決定了它對輸入質(zhì)量比傳統(tǒng) OTU 聚類更敏感。因此拿到的 fastq 不能一股腦全跑進(jìn)去必須先做質(zhì)量截?cái)唷＝財(cái)鄥?shù)怎么定基本都來自 demux.qzv 里的質(zhì)量圖。打開 demux.qzvInteractive Quality Plot 會(huì)畫出 R1 和 R2 每個(gè)堿基位置的 Boxplot。橫軸是 read 位置縱軸是 Phred 質(zhì)量分?jǐn)?shù)箱體的中位數(shù)一般會(huì)在某個(gè)位置明顯往下掉。比如常見的雙端 250 bp 測序R1 質(zhì)量往往在 230 bp 以后開始下滑R2 則在 200 bp 以后就開始下滑。截?cái)嚅L度要選在質(zhì)量明顯下降之前的位置同時(shí)兼顧雙端合并時(shí)需要的重疊長度。16S 擴(kuò)增子測序通常讀長在 250 bp 左右兩個(gè) read 合并后需要至少 20 bp 以上重疊DADA2 才能可靠地合并出完整的擴(kuò)增子序列。如果截?cái)嗵瘫热?R1 只留 150 bp、R2 只留 150 bp重疊區(qū)域可能不足大量 read 會(huì)在 merge 階段別丟棄最后得到的特征表會(huì)非常稀疏。我一般會(huì)根據(jù)質(zhì)量圖做一張簡單的參數(shù)表再?zèng)Q定read 方向質(zhì)量明顯下降位置推薦截?cái)嚅L度備注R1240 bp 左右230–250若目標(biāo)區(qū)域較長保留更多長度R2210 bp 左右190–220R2 質(zhì)量通常比 R1 差截短一點(diǎn)不心疼注意這里有一個(gè)常見誤用有人把截?cái)嚅L度當(dāng)成“統(tǒng)一切到相同長度”這是不對的。DADA2 的trunc-len是兩個(gè) read 分別設(shè)置的R1 和 R2 可以不同。實(shí)際項(xiàng)目中 R1 截 240、R2 截 200 十分常見。另外如果測序片段本身很短比如做 V4 區(qū)域用 PE150 測序那就要確保兩端截完后重疊區(qū)還有 20 bp 以上而不是盲目按質(zhì)量圖把低質(zhì)量區(qū)全切掉。3.2 引物處理先 cutadapt再接 DADA2很多公司交付的數(shù)據(jù)里其實(shí)已經(jīng)去掉了引物和 barcode但也有相當(dāng)一部分只做了 demultiplex沒有切引物。判斷方法很簡單用seqkit head或直接解壓看幾條 read如果 5 端開頭十幾個(gè)堿基和你的引物序列一致就需要先切。QIIME 2 里習(xí)慣用 q2-cutadapt 插件完成這一步。qiime cutadapt trim-paired \ --i-demultiplexed-sequences demux.qza \ --p-front-f GTGYCAGCMGCCGCGGTAA \ --p-front-r GGACTACNVGGGTWTCTAAT \ --p-cores 4 \ --o-trimmed-sequences trimmed.qza--p-front-f是正向引物序列--p-front-r是反向引物序列。對于 V4 區(qū)常用的 515F/806R 引物序列就是上面這兩個(gè)。--p-cores控制并行核心數(shù)這里設(shè) 4 就夠因?yàn)?cutadapt 本身很快。切完之后原來 read 長度會(huì)減少約 20 到 30 bp后面設(shè)置 DADA2 截?cái)嚅L度時(shí)要把這個(gè)損失算進(jìn)去。切完引物后建議再跑一次qiime demux summarize重新看一眼質(zhì)量分布因?yàn)?5 端被切掉之后質(zhì)量圖的起始位置發(fā)生了變化原來決定的截?cái)嚅L度可能需要微調(diào)。如果公司已經(jīng)明確說“數(shù)據(jù)已去除引物”或者你目測 read 開頭不是引物序列可以直接跳過 cutadapt把 demux.qza 作為 DADA2 的輸入。這里不要硬切切多了會(huì)損失可用的序列信息。3.3 DADA2 去噪命令截?cái)?、pooling、chimera 怎么配合接下來是整套流程中計(jì)算密度最高、也最容易翻車的一步。以修剪完引物的雙端數(shù)據(jù)為例我的去噪命令長這樣qiime dada2 denoise-paired \ --i-demultiplexed-seqs trimmed.qza \ --p-trim-left-f 0 --p-trim-left-r 0 \ --p-trunc-len-f 230 --p-trunc-len-r 210 \ --p-pooling-method pseudo \ --p-chimera-method consensus \ --p-n-threads 8 \ --o-table table.qza \ --o-representative-sequences rep-seqs.qza \ --o-denoising-stats stats.qza先解釋最容易理解的三個(gè)參數(shù)。--p-trim-left-f和--p-trim-left-r是切 read 5 端的長度。如果之前已經(jīng)用 cutadapt 去掉引物這里通常設(shè)為 0如果沒用 cutadapt 且引物長度固定也可以在這里直接設(shè)成引物長度比如 V4 引物 20 bp 就設(shè)--p-trim-left-f 20 --p-trim-left-r 20。--p-trunc-len-f和--p-trunc-len-r是 3 端截?cái)嚅L度按質(zhì)量圖定。--p-n-threads控制線程數(shù)默認(rèn) 0 表示使用所有可用核心但這一步對內(nèi)存的影響很大尤其樣本多時(shí)我建議先保守設(shè)成物理核心數(shù)的一半比如 8 核機(jī)器就設(shè) 4后面第 5 章會(huì)詳細(xì)講 OOM 問題。--p-pooling-method是 DADA2 特有的一個(gè)選擇取值有independent、pseudo和pooled。默認(rèn)是independent意思是每個(gè)樣本獨(dú)立學(xué)習(xí)自己的誤差模型速度快、內(nèi)存占用小但低深度樣本會(huì)很吃虧因?yàn)闃颖緝?nèi)部序列太少誤差模型估計(jì)不準(zhǔn)。pseudo是先對所有樣本的序列做一次預(yù)聚類用全局信息輔助誤差學(xué)習(xí)但又不至于像pooled那樣把所有樣本混在一起是深度不均勻項(xiàng)目的推薦選項(xiàng)。如果樣本數(shù)不多且深度都還可以用默認(rèn)independent就夠。pooled最嚴(yán)格也最慢內(nèi)存開銷最大一般不推薦用于幾百個(gè)樣本的大型項(xiàng)目。--p-chimera-method默認(rèn)是consensus這是 DADA2 的推薦模式先在每個(gè)樣本內(nèi)部識別嵌合體再跨樣本確認(rèn)。如果你發(fā)現(xiàn)嵌合體比例特別高比如去噪后非嵌合體序列只剩 50%可以改成pooled更嚴(yán)格地識別但也會(huì)多去掉一些可能真實(shí)的稀有序列。我這里沒有把--p-max-ee之類的參數(shù)寫進(jìn)命令因?yàn)椴煌?QIIME 2 版本暴露的參數(shù)名有差異真需要深調(diào)的可以去查當(dāng)前版本qiime dada2 denoise-paired --help的輸出。去噪完成后會(huì)生成三個(gè)文件table.qza 是特征表行是樣本列是 ASVrep-seqs.qza 是每個(gè) ASV 的代表序列stats.qza 是每一步保留序列數(shù)的統(tǒng)計(jì)。很多新手只看 table.qza 就繼續(xù)分析忽略 stats.qza這是不對的。3.4 去噪結(jié)果怎么看stats 和 table 摘要?jiǎng)e跳過去噪結(jié)束后第一時(shí)間查看去噪統(tǒng)計(jì)。QIIME 2 里可以用 metadata 方式把 stats.qza 轉(zhuǎn)成可視化qiime metadata tabulate \ --m-input-file stats.qza \ --o-visualization stats.qzv打開 stats.qzv每一行是一個(gè)樣本列包含 Input、Filtered、Denoised、Merged、Non-chimeric 等幾個(gè)階段。我判斷數(shù)據(jù)質(zhì)量的標(biāo)準(zhǔn)是從 Input 到 Non-chimeric 的整體保留率低于 50%就要回頭調(diào)截?cái)鄥?shù)高于 70% 說明數(shù)據(jù)質(zhì)量不錯(cuò)可以繼續(xù)。注意 DADA2 的結(jié)果不是“保留得越多越好”因?yàn)樗言S多低質(zhì)量 read 在學(xué)習(xí)誤差模型時(shí)識別為錯(cuò)誤序列并剔除去除錯(cuò)誤也是對后續(xù)分析的保護(hù)。特征表的摘要同樣不能省qiime feature-table summarize \ --i-table table.qza \ --m-sample-metadata-file sample-metadata.tsv \ --o-visualization table.qzv這里提前用到了 sample-metadata.tsv這是后續(xù)所有分組統(tǒng)計(jì)的元數(shù)據(jù)文件格式要求是制表符分隔第一列必須是sample-id第二列開始是你的分組信息、環(huán)境因子等。table.qzv 里能看到每個(gè)樣本的總序列數(shù)和特征數(shù)用于確定后面多樣性分析時(shí)的采樣深度這個(gè)深度選不對后面所有 alpha/beta 多樣性都會(huì)失真。這里只提一個(gè)初步過濾的操作如果某個(gè)樣本總序列數(shù)只有幾百條而其他樣本有幾萬條這類樣本要么做實(shí)驗(yàn)時(shí) DNA 量太低要么建庫失敗建議在進(jìn)入多樣性分析前剔除不要留著當(dāng)極值影響后續(xù)統(tǒng)計(jì)。4. 分類學(xué)注釋、進(jìn)化樹和多樣性讓表格能說“組間差異”4.1 選擇正確的訓(xùn)練分類器區(qū)域匹配比算法更重要特征表和代表序列出來后下一步是讓每個(gè) ASV 獲得一個(gè)物種名。QIIME 2 的 q2-feature-classifier 采用 Naive Bayes 分類器對代表序列進(jìn)行分類但分類器的效果很大程度上取決于訓(xùn)練數(shù)據(jù)。這里的一個(gè)核心原則訓(xùn)練分類器用的參考序列最好來自與你擴(kuò)增區(qū)域相同的高變區(qū)。比如你做的是 V4 區(qū)就應(yīng)該用基于 V4 區(qū)訓(xùn)練的分類器你用全長 16S 訓(xùn)練的分類器去分類 V4 短序列結(jié)果往往會(huì)在屬、種水平上出現(xiàn)大量“張冠李戴”。實(shí)際使用中官方生態(tài)提供多個(gè)預(yù)訓(xùn)練分類器常見的有 Greengenes 13_8 的 99% OTU 全長分類器、Silva 138 的不同區(qū)域分類器等。我一般會(huì)優(yōu)先選 Silva 數(shù)據(jù)庫因?yàn)樗姆诸悓W(xué)層級更新較新對未培養(yǎng)菌群的注釋也更友好。如果你做的是特定區(qū)域最好找與該區(qū)域匹配的預(yù)訓(xùn)練分類器文件否則就要自己訓(xùn)練后者在后面的高級章節(jié)里提到。分類命令如下qiime feature-classifier classify-sklearn \ --i-classifier silva-138-99-v4-classifier.qza \ --i-reads rep-seqs.qza \ --p-confidence 0.7 \ --p-n-jobs 1 \ --o-classification taxonomy.qza--i-classifier指定預(yù)訓(xùn)練分類器的 .qza 文件--i-reads指向去噪得到的代表序列--o-classification輸出分類結(jié)果。需要解釋的關(guān)鍵參數(shù)是--p-confidence默認(rèn)值是 0.7。這個(gè)值表示分類器對某個(gè)分類學(xué)層級的置信度閾值低于這個(gè)值就會(huì)被標(biāo)記為 Unassigned。如果閾值設(shè)太高比如 0.9很多短序列 ASV 會(huì)被標(biāo)記為 unassigned設(shè)太低比如 0.5又會(huì)出現(xiàn)大量置信度很低的“強(qiáng)行注釋”。對于 16S V4 區(qū)數(shù)據(jù)我建議從默認(rèn)的 0.7 開始如果 Unassigned 比例仍然偏高可以降到 0.6但不要低于 0.5。分類之后用qiime taxa tabulate和qiime taxa barplot檢查結(jié)果qiime taxa tabulate \ --i-taxonomy taxonomy.qza \ --o-visualization taxonomy.qzv qiime taxa barplot \ --i-table table.qza \ --i-taxonomy taxonomy.qza \ --m-metadata-file sample-metadata.tsv \ --o-visualization taxa-barplot.qzvtaxa-barplot.qzv 能按樣本顯示各分類學(xué)水平的相對豐度堆疊圖一眼就能看出分組之間門水平的差異方向。如果看到某個(gè)樣本幾乎所有序列都是同一個(gè)門甚至完全是單一屬先不要急著做生物學(xué)解釋回頭檢查是不是樣本標(biāo)簽搞混了或者實(shí)驗(yàn)過程中發(fā)生污染。4.2 建樹對齊與掩蔽別在系統(tǒng)發(fā)育樹上省步驟如果后續(xù)要計(jì)算 Faith‘s PD 這類系統(tǒng)發(fā)育多樣性或者用 Unifrac 距離衡量組間差異就需要建一棵系統(tǒng)發(fā)育樹。QIIME 2 提供了整合的流程qiime phylogeny align-to-tree-mafft-fasttree \ --i-sequences rep-seqs.qza \ --p-n-threads 4 \ --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這個(gè)流程內(nèi)部先用 MAFFT 做多序列比對然后對高變區(qū)進(jìn)行掩蔽再用 FastTree 建樹。輸出中rooted-tree.qza是我們之后做 Unifrac 分析要用的輸入。要特別提醒的是當(dāng) ASV 數(shù)量很多比如超過 5 萬個(gè)時(shí)MAFFT 和 FastTree 會(huì)非常耗時(shí)幾個(gè)小時(shí)的運(yùn)行時(shí)間很常見。如果機(jī)器配置一般建議先通過低頻特征過濾減少 ASV 數(shù)量或者用--p-n-threads適當(dāng)增加并行度。有一個(gè)容易被忽視的點(diǎn)建樹的代表序列最好和分類學(xué)注釋用同一份 rep-seqs.qza不要因?yàn)檎頂?shù)據(jù)時(shí)誤用了另一批序列。輸出 unrooted-tree.qza 和 rooted-tree.qza 的區(qū)別在于根的位置后續(xù)所有計(jì)算系統(tǒng)發(fā)育多樣性的地方都用 rooted-tree別拿錯(cuò)文件。如果項(xiàng)目只需要 Bray-Curtis、Jaccard 這類非系統(tǒng)發(fā)育距離理論上可以不建樹用core-metrics的非系統(tǒng)發(fā)育版本就夠了。但為了讓結(jié)果能橫向比較也為了后續(xù)能算 Faith‘s PD我還是建議跑一次建樹流程一次投入后面不用返工。4.3 alpha/beta 多樣性采樣深度決定你能比較什么多樣性分析里的第一個(gè)重要決策是采樣深度sampling depth。這個(gè)值表示在計(jì)算 alpha 和 beta 多樣性之前每一樣本都抽到相同的測序深度以消除測序深度不同帶來的偏差。很多人直接復(fù)制教程里的--p-sampling-depth 1000或--p-sampling-depth 5000但正確的做法是先看 table.qzv 里每個(gè)樣本的總頻率分布選一個(gè)能保留大部分樣本、又不低得離譜的值。比如有 50 個(gè)樣本最少的樣本有 3000 條序列大部分樣本在 8000 到 12000那么采樣深度設(shè) 3000 到 5000 是合理的低于 3000 的樣本會(huì)被自動(dòng)剔除不會(huì)參與后續(xù)分析。核心多樣性分析命令如下qiime diversity core-metrics-phylogenetic \ --i-phylogeny rooted-tree.qza \ --i-table table.qza \ --p-sampling-depth 5000 \ --m-metadata-file sample-metadata.tsv \ --output-dir core-metrics這條命令會(huì)根據(jù)是否提供樹一次性生成多位 alpha 多樣性和 beta 多樣性結(jié)果。在 core-metrics 目錄里你會(huì)看到 shannon_vector.qza、observed_features_vector.qza、faith_pd_vector.qza、evenness_vector.qza以及 unweighted_unifrac_distance_matrix.qza、weighted_unifrac_distance_matrix.qza、jaccard_distance_matrix.qza、bray_curtis_distance_matrix.qza還有對應(yīng)的 PCoA 可視化文件。這里有一個(gè)常見誤用把--p-sampling-depth理解為“過濾掉低序列樣本的閾值”。它不是過濾閾值而是稀疏化深度。如果設(shè)得比某些樣本總序列數(shù)高那些樣本會(huì)在 alpha 多樣性中被當(dāng)作 NaN 排除如果設(shè)得太低又會(huì)浪費(fèi)大量測序信息。我通常會(huì)把采樣深度設(shè)在所有樣本總頻率的第 10 百分位附近這樣既保留了絕大多數(shù)樣本又不會(huì)讓稀疏化損失太多信息。得到距離矩陣后組間差異顯著性檢驗(yàn)通常用 beta-group-significanceqiime diversity beta-group-significance \ --i-distance-matrix core-metrics/bray_curtis_distance_matrix.qza \ --m-metadata-file sample-metadata.tsv \ --m-metadata-column Group \ --p-pairwise \ --o-visualization beta-bray-curtis.qzv--m-metadata-column指定按 metadata 里的哪個(gè)分組列做檢驗(yàn)比如Group列寫的是 Ctrl 和 Treat。--p-pairwise會(huì)在多組時(shí)輸出兩兩比較的 PERMANOVA 結(jié)果?？唇Y(jié)果時(shí)不要只盯著 p 值還要看排列圖中各組的分離程度。一個(gè)常見的翻車點(diǎn)metadata 里的組名帶有空格或特殊字符QIIME 2 內(nèi)部會(huì)把它對數(shù)值處理導(dǎo)致分組失效。所以建 metadata 時(shí)盡量用簡單的英文單詞比如Group列下寫Ctrl、Treat不要寫含有空格或中文的內(nèi)容。alpha 多樣性的組間比較也是類似邏輯用qiime diversity alpha-group-significance以 shannon_vector.qza 為輸入按分組列計(jì)算 Kruskal-Wallis 檢驗(yàn)。這些結(jié)果放在一起才構(gòu)成論文里常見的“組間 alpha 多樣性差異 beta 多樣性 PCoA 圖 PERMANOVA 檢驗(yàn)”三件套。5. 避坑排查QIIME 2 實(shí)操中 5 個(gè)高頻故障現(xiàn)場5.1 導(dǎo)入一直失敗不是格式錯(cuò)誤而是路徑和 Phred 編碼在作怪現(xiàn)象是最常見的qiime tools import命令回車之后報(bào)錯(cuò)要么是“No samples in manifest file”要么是“Could not find file /home/xx/manifest.csv”。不少人反復(fù)檢查 manifest甚至重新生成好幾次問題依舊。原因通常有兩個(gè)。一是 manifest 里的路徑寫的是相對路徑或者帶~QIIME 2 不會(huì)解釋~必須寫成以/開頭的絕對路徑二是樣本 ID 列名帶了不可見字符尤其是 Windows 下編輯的 CSV 帶 BOM 頭導(dǎo)致第一列實(shí)際變成了\ufeffsample-idQIIME 2 識別不到sample-id這個(gè)列名。Phred 編碼錯(cuò)誤的表現(xiàn)又不一樣它可能順利導(dǎo)入但后續(xù) DADA2 結(jié)果糟糕得離譜或者直接提示質(zhì)量格式不匹配。解決辦法是三步。第一用cat -A manifest.csv查看文件里是否有^M回車符和 BOM有的話用dos2unix manifest.csv或sed -i s/\r// manifest.csv轉(zhuǎn)換。第二確認(rèn) manifest 中路徑為絕對路徑且文件確實(shí)存在不要用編輯器肉眼判斷直接用ls -l逐個(gè)路徑檢查。第三在導(dǎo)入前先解壓一個(gè) fastq.gz 文件觀察質(zhì)量字符范圍確認(rèn) Phred 編碼后再選擇...Phred33...還是...Phred64...格式。這個(gè)故障是最不值得花時(shí)間定位的一次檢查能省半天。5.2 DADA2 保留率暴跌多半是截?cái)嚅L度和引物殘留的問題現(xiàn)象去噪完成后打開 stats.qzv發(fā)現(xiàn)很多樣本從 Input 到 Non-chimeric 的保留率不到 30%甚至有的樣本 Non-chimeric 序列數(shù)為零?？雌饋硐?DADA2 把好序列都丟了實(shí)際上多半是你的截?cái)鄥?shù)和引物處理打架了。如果數(shù)據(jù)里還殘留引物而你沒有用 cutadapt 就去跑 DADA2DADA2 會(huì)把 5 端的引物序列當(dāng)成真實(shí)生物學(xué)序列的一部分導(dǎo)致合并時(shí)無法正確配對大量 read 在 merge 階段失敗。如果截?cái)嚅L度設(shè)置得太短兩條 read 合并后的重疊區(qū)域不足 20 bp同樣會(huì)導(dǎo)致大量 read 無法合并。還有一種情況是 R2 質(zhì)量整體很差而你把 trunc-len-r 設(shè)得和 R1 一樣長導(dǎo)致很多 R2 read 在截?cái)嗪笾皇Ｒ恍《斡行蛄?。解決思路是先回質(zhì)量圖確認(rèn)引物是否已經(jīng)切干凈然后按 R1/R2 各自質(zhì)量分布調(diào)整截?cái)嚅L度確保合并后至少保留 20 bp 重疊。如果項(xiàng)目測序讀長本身就很短比如 PE150 做 V3-V4 區(qū)域重疊不足不是調(diào)參數(shù)能解決的需要換引物或改測序策略。我一般會(huì)在去噪前先跑一次qiime demux summarize不只看質(zhì)量還看 read 長度分布確保長度一致。如果某些樣本 read 長度特別短可以直接先用qiime feature-table filter-samples把它們剔除避免影響整體誤差模型。5.3 分類結(jié)果出現(xiàn)過高 unassigned 或只在門級別打轉(zhuǎn)現(xiàn)象分類學(xué)注釋完成后unassigned 比例超過 40%或者每個(gè) ASV 都只注釋到門、綱根本沒有屬和種的水平。非生物學(xué)原因占了大半。第一你用的預(yù)訓(xùn)練分類器與你的擴(kuò)增區(qū)域不匹配。比如用全長 16S 分類器去分類 V4 短序列因?yàn)槎绦蛄腥鄙偃L序列里的保守位點(diǎn)Naive Bayes 很難給到高置信度。第二置信度閾值設(shè)太高。0.8 或 0.9 時(shí)短 ASV 很難達(dá)到這個(gè)閾值于是大量被記為 Unassigned。第三代表序列里可能混入了嵌合體或宿主序列這些序列本身在參考數(shù)據(jù)庫里找不到近鄰。處理方式按順序來先看 rep-seqs.qzv 里的序列長度分布確認(rèn)平均長度是不是符合你的目標(biāo)區(qū)域然后換一個(gè)與目標(biāo)區(qū)域匹配的分類器最后把--p-confidence降到 0.7 或 0.6 重新分類。如果換了區(qū)域匹配的分類器之后屬水平注釋率仍然很低再考慮是不是數(shù)據(jù)本身有問題這時(shí)候需要挑幾條 unassigned 的代表序列到 NCBI BLAST 里看最接近的是什么不要盲目調(diào)參掩蓋真實(shí)問題。5.4 內(nèi)存不夠不一定是機(jī)器小可能是你開了過多線程或用了 pooled現(xiàn)象qiime dada2 denoise-paired運(yùn)行到一半終端直接報(bào)std::bad_alloc或Segmentation fault更明顯的是用dmesg能看到Out of memory: Kill process的系統(tǒng)日志。很多人第一反應(yīng)是“機(jī)器內(nèi)存太小”但其實(shí) 64G 內(nèi)存跑 200 個(gè)樣本照樣可能崩。DADA2 的算法在執(zhí)行時(shí)需要在內(nèi)存中維護(hù)全體的序列和誤差狀態(tài)它并不是簡單的逐樣本流式處理。--p-pooling-method pooled會(huì)把所有樣本的序列混在一起學(xué)習(xí)誤差模型內(nèi)存開銷最大--p-n-threads 0會(huì)讓 DADA2 把所有核心全部占滿每個(gè)線程都分配獨(dú)立內(nèi)存進(jìn)一步放大了內(nèi)存壓力。這兩個(gè)參數(shù)組合在一起是 OOM 的經(jīng)典原因。解決方式把--p-pooling-method改回pseudo或independent把--p-n-threads從 0 改成 4 或 8先跑通再說。如果還是內(nèi)存不足就要考慮分批次對樣本進(jìn)行去噪再把生成的 table.qza 用qiime feature-table merge合并。分批次去噪在項(xiàng)目設(shè)計(jì)上可行但要注意批次間 DADA2 參數(shù)必須完全一致否則 ASV 邊界可能不一樣。還有一種做法是先用qiime feature-table filter-samples --p-min-frequency把序列數(shù)極少的樣本過濾掉一部分減少整體數(shù)據(jù)量。總之不到萬不得已不要上 pooled內(nèi)存不足時(shí)先降線程再降 pooling 級別。5.5 環(huán)境被搞亂別急著 conda install重建環(huán)境是最穩(wěn)的后悔藥現(xiàn)象某天打開終端qiime --help還能正常顯示但一運(yùn)行qiime diversity就報(bào) “No plugin named diversity”或者qiime feature-classifier報(bào)ModuleNotFoundError: No module named q2_feature_classifier。這時(shí)候如果你繼續(xù)在當(dāng)前環(huán)境里conda install qiime2或conda update某一兩項(xiàng)往往越搞越亂。原因基本只有一個(gè)環(huán)境里的 Python 根據(jù)某個(gè)安裝操作被升級或降級比如你想裝別的包時(shí)用了pip install --upgrade scikit-learn破壞了 QIIME 2 插件所依賴的版本鎖定關(guān)系。QIIME 2 對依賴鎖定非常嚴(yán)格不像普通 Python 包那樣可以隨意升級。正確解法是承認(rèn)環(huán)境已壞不要再修。先用conda list -n qiime2 qiime2-env.txt把當(dāng)前環(huán)境清單導(dǎo)出留檔然后建立一個(gè)新環(huán)境重新按官方 yml 安裝。有人覺得這樣浪費(fèi)時(shí)間但實(shí)際經(jīng)驗(yàn)告訴我修復(fù)一個(gè)依賴混亂的 conda 環(huán)境所花費(fèi)的時(shí)間幾乎總是比重建一個(gè)干凈環(huán)境更多而且還容易殘留隱患。重建之后用qiime info確認(rèn)所有插件都回來了再把之前跑通的 qza 和可視化文件拿到新環(huán)境里繼續(xù)qza 本身是跨版本的舊環(huán)境產(chǎn)物可以繼續(xù)使用。這也是為什么我建議把項(xiàng)目數(shù)據(jù)和 conda 環(huán)境分開存放環(huán)境壞了數(shù)據(jù)還在。6. 進(jìn)階導(dǎo)出到 R 做差異統(tǒng)計(jì)并用“可復(fù)現(xiàn)盒子”守住結(jié)果6.1 導(dǎo)出為 BIOM/TSV交給 phyloseq 繼續(xù)分析QIIME 2 自帶的統(tǒng)計(jì)可視化適合初步探索但到了論文級分析差異菌屬篩選、排序分析和個(gè)性化繪圖往往要交給 R。常見做法是把特征表和分類學(xué)結(jié)果導(dǎo)出成 BIOM 格式再用 phyloseq 讀取。qiime tools export \ --input-path table.qza \ --output-path export-table biom convert \ -i export-table/feature-table.biom \ -o export-table/feature-table.tsv \ --to-tsv導(dǎo)出后的 TSV 文件第一行是# Constructed from biom file第二行是OTU ID和各樣本 ID這個(gè)格式很多人第一次看會(huì)懵它并不是標(biāo)準(zhǔn)表格。我一般在 R 里直接用 phyloseq 讀 BIOM 而不是讀 TSVlibrary(phyloseq) ps - import_biom(export-table/feature-table.biom) tax - read.csv(export-table/taxonomy.tsv, sep \t) # 從 taxonomy.qza 解壓導(dǎo)出后處理更省事的方式是用 qiime2R 包直接讀 qza 原生文件但不同版本的 qiime2R 對 qza 兼容性不一致如果遇到讀取失敗回到 BIOM 導(dǎo)出這條路是最穩(wěn)妥的。導(dǎo)入后就可以做 LEfSe、ANCOM-BC、隨機(jī)森林篩選等下游分析這些已經(jīng)超出 QIIME 2 本身范圍但導(dǎo)出這一步是必備接口。6.2 把每一步命令存成可重跑的腳本并記錄環(huán)境版本我見過太多項(xiàng)目在分析完成兩個(gè)月后想再調(diào)整一個(gè)參數(shù)時(shí)發(fā)現(xiàn)記錄丟了只能靠回憶猜當(dāng)初跑過什么命令。避免這個(gè)黑匣子的方法很簡單從導(dǎo)入開始把每條命令連同參數(shù)寫進(jìn)一個(gè)pipeline.sh并在文件頭部記錄 QIIME 2 版本和插件列表。{ echo $(date) qiime --version conda list | grep q2- pipeline_log.txt qiime dada2 denoise-paired ... } pipeline_log.txt保存命令還不夠環(huán)境版本也要留痕。QIIME 2 同一流程在不同插件版本下DADA2 的 ASV 結(jié)果可能會(huì)有細(xì)微差異。qiime info輸出的版本信息里包含了每個(gè)插件版本號把這部分重定向到日志文件里比在文檔里手寫版本可靠得多。再有一個(gè)習(xí)慣是給輸出文件命名時(shí)帶上關(guān)鍵參數(shù)比如table_trunc230_210_pseudo.qza雖然文件名長一些但以后看到文件名就能想起當(dāng)初用的參數(shù)不用翻日志。6.3 用 tools peek 和 feature-table summarize 做最終驗(yàn)證跑完整套流程后發(fā)布前我會(huì)做最后一項(xiàng)檢查用qiime tools peek逐個(gè)抽查 qza 文件確認(rèn)文件的 type 和 UUID 正常沒有出現(xiàn)損壞或?qū)脲e(cuò)誤。對關(guān)鍵文件比如最終的 table.qza 和 rooted-tree.qza再跑一次qiime feature-table summarize核對特征表行列數(shù)、樣本 ID 集合是否和 metadata 完全一致。這一步能發(fā)現(xiàn)很多潛在問題比如某樣本文檔在 metadata 里存在但在特征表里為空或者樣本名因大小寫不一致被拆成兩列這些都是常見翻車點(diǎn)。我自己現(xiàn)在每次新項(xiàng)目開始都會(huì)固定建三層目錄data/放原始 fastq 和 manifestscripts/放所有命令日志results/放 qza 和 qzv。環(huán)境用一個(gè)獨(dú)立 conda env 管理版本號寫入 scripts 里的 timestamp 文件。上次因?yàn)樨潏D方便在同一個(gè)環(huán)境里跑了個(gè)舊項(xiàng)目的新插件把整體環(huán)境弄壞結(jié)果重跑所有數(shù)據(jù)花了兩天從那以后我再也不混環(huán)境也再不敢省版本記錄這一步。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取