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