簡介一份聚焦微生物群落功能預(yù)測的專題 PDF圍繞 FAPROTAX 在 16S 擴(kuò)增子分析中的應(yīng)用展開面向微生物生態(tài)、環(huán)境科學(xué)方向的研究生與科研人員幫助理解環(huán)境樣本中碳、氫、氮、磷、硫等元素循環(huán)的功能注釋思路也適合作為組會匯報與課題設(shè)計時的概念參考。內(nèi)容從數(shù)據(jù)庫的文獻(xiàn)構(gòu)建依據(jù)講起說明 OTU 分類表與功能注釋的關(guān)聯(lián)方式并延伸到功能—樣本矩陣、分組功能豐度熱圖、環(huán)境因子相關(guān)性泡泡圖等結(jié)果呈現(xiàn)同時對比 ImageGP 一鍵分析與 Linux 本地運(yùn)行兩種實(shí)現(xiàn)路徑也提示了預(yù)測準(zhǔn)確度與覆蓋度的適用邊界。壓縮包僅含 1 個 PDF 文件大小約 2.47MB輕量便于查閱。已有 1148 人學(xué)習(xí)讀者可據(jù)此建立 FAPROTAX 預(yù)測流程的整體框架并通過哺乳動物皮膚微生物、水稻氮利用等案例把握功能預(yù)測在差異挖掘與機(jī)制關(guān)聯(lián)中的用法。1. 拿到 16S 結(jié)果之后功能預(yù)測這一步怎么落地擴(kuò)增子測序做到最后ASV 表和物種注釋都齊了合作方或者審稿人往往不會滿足于“這批樣本里擬桿菌門多、變形菌門少”他們要的是——這些菌在干哪些活氮循環(huán)、硫循環(huán)、碳降解到底有沒有差異。FAPROTAXFunctional Annotation of Prokaryotic Taxa就是為這個問題準(zhǔn)備的它不訓(xùn)練模型也不依賴參考基因組而是用一份人工從文獻(xiàn)里整理出來的“分類名—功能”對照表把 16S 的屬種注釋直接映射成元素循環(huán)相關(guān)的功能組豐度。它最適合環(huán)境、農(nóng)業(yè)、水體沉積物這類關(guān)注生物地球化學(xué)循環(huán)的樣本對腸道、皮膚這類宿主相關(guān)樣本也能用但解釋要謹(jǐn)慎。想預(yù)測整條 KO 通路的人應(yīng)該去看 PICRUSt2想把“這些菌可能參與哪幾步氮轉(zhuǎn)化”講清楚的人繼續(xù)往下看。2. FAPROTAX 在算什么物種名到元素循環(huán)功能的映射表2.1 數(shù)據(jù)庫不是模型推導(dǎo)而是文獻(xiàn)里“抄”出來的對照表理解 FAPROTAX 的關(guān)鍵是先接受它“簡單得有點(diǎn)粗暴”的設(shè)計。數(shù)據(jù)庫文件 FAPROTAX.txt 本質(zhì)上是一張兩列的映射表左邊是分類單元名屬名或種名右邊是用制表符分隔的若干功能標(biāo)簽比如nitrogen_fixation、nitrification、sulfate_respiration、methanotrophy、cellulolysis。這些標(biāo)簽不是算法推斷出來的而是整理者逐篇翻閱已發(fā)表的純培養(yǎng)、富集培養(yǎng)和生理實(shí)驗(yàn)文獻(xiàn)確認(rèn)“這個屬確實(shí)具備這個代謝能力”之后才寫進(jìn)去的。這樣做的好處很直接每一條注釋背后都有一篇可追溯的文獻(xiàn)預(yù)測出來的功能不是統(tǒng)計相關(guān)性而是生理學(xué)證據(jù)。代價同樣明顯只有被研究過、被培養(yǎng)過、被驗(yàn)證過的類群才會進(jìn)表冷門門類和大量未培養(yǎng)微生物基本是空白。所以它在環(huán)境樣本上的表現(xiàn)取決于樣本里有多少類群落在“已培養(yǎng)、有文獻(xiàn)”的范圍里。這也解釋了它為什么和 2016 年那篇海洋微生物組研究綁在一起那項(xiàng)工作想說明海洋里功能組成和物種組成并不是一一對應(yīng)的同一個功能可以被很多不同類群分擔(dān)而 FAPROTAX 提供了把物種翻譯成功能的工具讓這種“解耦”能在 16S 數(shù)據(jù)上被觀察。經(jīng)過幾年迭代數(shù)據(jù)庫在 2020 年 11 月隨 1.2.4 版做了一次較大更新注釋比例比早期版本高不少腳本運(yùn)行環(huán)境也從 Python 2 遷到了 Python 3。提示FAPROTAX 的功能標(biāo)簽是“能力”而不是“活性”。一個屬被標(biāo)了denitrification只說明它有這個基因通路不代表樣本里真的在反硝化。寫文章時把措辭放在“潛在功能”“功能潛力”上比直接寫“反硝化強(qiáng)度”穩(wěn)得多。2.2 和 PICRUSt2、Tax4Fun 的邊界覆蓋度換準(zhǔn)確度同樣是 16S 功能預(yù)測幾條技術(shù)路線的取舍完全不同選錯工具比參數(shù)調(diào)錯更致命。PICRUSt2 走的是“參考基因組樹 基因家族推斷”的路線把 ASV 序列放進(jìn)參考樹推測基因拷貝數(shù)再歸到 KO 通路。它的覆蓋度很高幾乎每個 ASV 都能落進(jìn)去輸出的是完整的 KEGG 通路譜適合做全局的功能差異分析。但它推斷的是“基因組潛力”鏈條長、假設(shè)多落在具體某一步氮轉(zhuǎn)化上時中間誤差會被放大。Tax4Fun 系含 Tax4Fun2基于 SILVA 注釋結(jié)果做最近鄰匹配輸出的是百分比形式的功能組成實(shí)現(xiàn)簡單但對參考庫依賴強(qiáng)遇到數(shù)據(jù)庫里沒有的類群就會失真。FAPROTAX 的定位正好卡在另一端覆蓋度最低但每一條映射都有生理學(xué)依據(jù)且直接對應(yīng)到“固氮、硝化、反硝化、硫酸鹽還原、甲烷氧化、纖維素降解”這種可以直接和實(shí)測理化指標(biāo)對話的功能組。樣本是水稻根際、湖泊水體、海洋沉積物、活性污泥研究問題是“氮素轉(zhuǎn)化菌群隨施肥怎么變”FAPROTAX 的說服力比通路預(yù)測強(qiáng)樣本是人體腸道研究問題是“整體代謝通路差異”就該換 PICRUSt2。我的習(xí)慣是兩條線并跑FAPROTAX 出元素循環(huán)的功能組用來和理化數(shù)據(jù)做關(guān)聯(lián)、畫機(jī)制圖PICRUSt2 出 KO 通路用來做無假設(shè)的差異篩選。兩者結(jié)論互相印證時故事才立得住只有 FAPROTAX 結(jié)果時討論部分一定要寫清覆蓋度這個前提。2.3 輸入的硬性要求帶物種注釋的 OTU/ASV 表FAPROTAX 不吃普通的文本豐度表它要的是 biom 格式而且物種注釋必須寫進(jìn) biom 的 observation metadata 里。這一步是絕大多數(shù)人第一次翻車的地方。具體來說輸入需要滿足三條第一豐度表是“行 OTU/ASV列 樣本”的矩陣第一列是 OTU ID表頭是樣本名值是非負(fù)整數(shù)或抽平后的數(shù)。第二每個 OTU 必須有分類字符串且字符串的最后一級落在屬或種。FAPROTAX 是靠屬名/種名去查表的如果你的注釋只到門p__Proteobacteria或者最后一級是g__uncultured它就只能算作 leftovers。第三分類字符串的層級順序要和數(shù)據(jù)庫預(yù)期一致通常是界門綱目科屬種依次排列用分號分隔。前綴寫不寫k__、p__無所謂腳本會做清理但順序不能亂也不能中間插別的字段。把這三條對齊后面的命令基本就是流水線。對不齊就會看到Assigned 12 records to groups這種讓人心里一涼的結(jié)果——我見過最慘的一次兩千多條 ASV 只匹配上十幾條原因就是上游注釋表的屬名那一級全帶著g__前綴加空格。3. 本地跑通 1.2.4從 txt 豐度表到功能矩陣3.1 建一個干凈的 Python3 環(huán)境別和 QIIME2 混用FAPROTAX 1.2.4 的腳本本身只依賴 numpy但整條流程要用 biom 命令行工具做格式轉(zhuǎn)換而 biom 的版本和 QIIME2 環(huán)境里的依賴經(jīng)常打架。血淚經(jīng)驗(yàn)是單獨(dú)建一個環(huán)境別圖省事直接在 QIIME2 的 conda 環(huán)境里跑。# 建一個只裝必要依賴的干凈環(huán)境 conda create -n faprotax python3.8 -y conda activate faprotax # numpy 是腳本硬依賴biom-format 提供 biom 命令行 conda install -c conda-forge numpy biom-format -y # 驗(yàn)證 biom 命令可用正常會打印版本號 biom --version環(huán)境建好后去項(xiàng)目主頁下載FAPROTAX_1.2.4.zip并解壓目錄里至少要有collapse_table.py和FAPROTAX.txt兩個文件。先跑一次幫助能彈出參數(shù)列表就說明解釋器路徑和依賴都沒問題。# 只看幫助不做任何計算 python FAPROTAX_1.2.4/collapse_table.py -h如果這一步報ModuleNotFoundError: No module named numpy說明 conda 環(huán)境沒激活對或者你在別的 shell 里跑報biom相關(guān)錯誤則是后面章節(jié)才會遇到的問題此處不涉及。先把-h跑通能省掉后面一半的排查時間。3.2 txt 轉(zhuǎn) biom 再掛物種注釋兩條命令假設(shè)手上是抽平后的otutab_rare.txt和一份單獨(dú)的注釋文件taxonomy2.txt。注釋文件的格式必須是兩列、制表符分隔、帶表頭表頭寫作OTUID和taxonomy順序不能反OTUID taxonomy ASV_1 k__Bacteria;p__Proteobacteria;c__Gammaproteobacteria;o__Nitrosomonadales;f__Nitrosomonadaceae;g__Nitrosomonas ASV_2 k__Bacteria;p__Proteobacteria;c__Alphaproteobacteria;o__Rhizobiales;f__Rhizobiaceae;g__Rhizobium轉(zhuǎn)換分兩步先把文本表轉(zhuǎn)成 biom JSON再把注釋作為 observation metadata 掛上去。# 第一步txt 轉(zhuǎn) biomJSON 格式注意表類型寫 OTU table biom convert -i otutab_rare.txt -o otutab_rare.biom \ --table-typeOTU table --to-json # 第二步把物種注釋寫入 observation metadata # --sc-separated 表示 taxonomy 這一列按分號拆成多值 # --observation-header 指定注釋文件的表頭字段順序 biom add-metadata -i otutab_rare.biom \ --observation-metadata-fp taxonomy2.txt \ -o otutab_rare_tax.biom \ --sc-separated taxonomy \ --observation-header OTUID,taxonomy # 第三步抽查前幾條 observation確認(rèn) taxonomy 真的進(jìn)去了 biom summarize-table -i otutab_rare_tax.biom | head -20 biom convert -i otutab_rare_tax.biom -o - --to-tsv | head -5這三條命令里--sc-separated是最容易被忽略的參數(shù)。它告訴 biom這一列的值不是普通字符串而是用分號分隔的多值字段后續(xù) FAPROTAX 才能按層級取出屬名。如果漏寫taxonomy 會被當(dāng)成一個整串存進(jìn)去腳本取最后一段時就會拿到整個字符串匹配全軍覆沒。--observation-header的順序也要和文件實(shí)際列順序嚴(yán)格對應(yīng)。寫成taxonomy,OTUID而文件是OTUID,taxonomybiom 不會報錯但會把 ASV ID 塞進(jìn) taxonomy 字段結(jié)果是注釋全錯、匹配為零屬于最隱蔽的一類錯誤。3.3 collapse_table.py 的參數(shù)逐條拆輸入準(zhǔn)備好之后一條命令出結(jié)果python FAPROTAX_1.2.4/collapse_table.py \ -i otutab_rare_tax.biom \ -g FAPROTAX_1.2.4/FAPROTAX.txt \ --collapse_by_metadata taxonomy \ -v --force \ -o faprotax.txt \ -r faprotax_report.txt逐個說清楚-i是上一步產(chǎn)出的、帶注釋的 biom 表。注意不能用原始 txt也不能用沒掛 metadata 的 biom。-g指向數(shù)據(jù)庫文件FAPROTAX.txt路徑寫錯會直接報文件不存在。建議用絕對路徑或確認(rèn)當(dāng)前工作目錄。--collapse_by_metadata taxonomy告訴腳本按 observation metadata 里名為taxonomy的字段來做歸并。這個字段名必須和biom add-metadata時寫的一致大小寫敏感。-v輸出詳細(xì)過程信息會打印每個功能組匹配到多少條記錄、哪些記錄被歸為 leftovers。排查階段強(qiáng)烈建議開著。--force強(qiáng)制覆蓋已有輸出文件批處理時必備否則第二次運(yùn)行會因?yàn)槲募嬖诙袛唷?o是功能豐度矩陣行是功能組列是樣本值是該功能組在所有匹配 OTU 上的豐度加和。-r是明細(xì)報告記錄每個功能組由哪些 OTU/ASV 貢獻(xiàn)這是后面做證據(jù)鏈的關(guān)鍵文件。運(yùn)行完屏幕上會出現(xiàn)類似Assigned 979 records to groups, 1924 records were leftovers的統(tǒng)計。這個比例要盯住匹配上的記錄數(shù)除以總記錄數(shù)就是這次預(yù)測的有效覆蓋度。經(jīng)驗(yàn)上環(huán)境樣本常見三到六成落到一兩成就得回頭查注釋格式別急著往下分析。腳本還支持用--column之類的參數(shù)指定從分類字符串里取第幾段作為匹配鍵不同版本參數(shù)名略有差別跑一次-h確認(rèn)比照著舊教程硬套可靠。3.4 兩個輸出文件分別看什么faprotax.txt是主結(jié)果第一列是功能組名后面每列一個樣本。它可以直接讀進(jìn) R 或 Python 做統(tǒng)計但要注意兩點(diǎn)一是行數(shù)遠(yuǎn)少于 KO 表通常幾十到一百多行這是覆蓋度決定的二是同一個 OTU 可能被同時計入多個功能組所以列求和會大于總 reads 數(shù)這是“一對多”映射的正?，F(xiàn)象不是數(shù)據(jù)錯了。faprotax_report.txt是黑匣子鑰匙。它按功能組組織每個組下列出歸屬到該組的 OTU/ASV 及其分類信息。想看“固氮功能到底是誰貢獻(xiàn)的”直接在這里翻審稿人質(zhì)疑“你這個固氮功能是不是被某個占主導(dǎo)的屬帶偏了”也在這里查。很多人只拿faprotax.txt畫圖把 report 丟在一邊等被問細(xì)節(jié)時才發(fā)現(xiàn)答不上來。4. report 拆成三類表讓結(jié)果能進(jìn)熱圖和分組比較4.1 先把 report 清洗成“功能組 記錄”兩行一組report 文件的結(jié)構(gòu)是每個功能組先占一行組名下面跟著若干行以 ASV/OTU ID 開頭的記錄行組與組之間用分隔線隔開。這種嵌套結(jié)構(gòu)不適合直接做表需要先壓平成“組名—記錄”交替的形式。# 抓取 ASV 行以及它的上一行即所屬功能組名去掉純分隔線 grep ASV_ -B 1 faprotax_report.txt | grep -v -P ^--$ faprotax_report.clean-B 1的作用是連帶輸出匹配行的前一行這樣每條記錄都跟著它的功能組-v -P ^--$把分隔線剔掉避免后面解析時出現(xiàn)空組。如果你的 OTU ID 不是ASV_開頭把 grep 的模式換成實(shí)際前綴比如OTU_。這一步看著土但對格式規(guī)整的 report 極其有效比寫正則解析整個文件穩(wěn)。4.2 用腳本拆出 mat / otu_func / func_otu清洗完的文件是交替的兩行結(jié)構(gòu)用一段短腳本拆成三張表OTU×功能的有無矩陣、OTU 到功能的列表、功能到 OTU 的列表。三種表對應(yīng)三種下游用途缺一不可。# report2tables.py # 讀入清洗后的 report輸出三類整理表 import sys clean, prefix sys.argv[1], sys.argv[2] pairs [] # (功能組, ASV) cur_group None with open(clean) as fh: for line in fh: line line.rstrip(\n) if not line.strip(): continue if not line.startswith(prefix): # 不以 ASV_ 開頭的行是功能組名 cur_group line.strip() else: asv line.split(\t)[0].strip() # 記錄行首個字段是 OTU ID if cur_group: pairs.append((cur_group, asv)) # 表一OTU - 功能一行一個 OTU功能用分號連接 otu2func {} for g, a in pairs: otu2func.setdefault(a, set()).add(g) with open(faprotax_report.otu_func, w) as out: out.write(OTUID\tfunctions\n) for a, gs in otu2func.items(): out.write(f{a}\t{;.join(sorted(gs))}\n) # 表二功能 - OTU用于反查某個功能由誰貢獻(xiàn) func2otu {} for g, a in pairs: func2otu.setdefault(g, set()).add(a) with open(faprotax_report.func_otu, w) as out: out.write(function\tOTUIDs\n) for g, as_ in func2otu.items(): out.write(f{g}\t{;.join(sorted(as_))}\n) # 表三OTU x 功能 0/1 矩陣直接拿去畫熱圖 funcs sorted({g for g, _ in pairs}) otus sorted(otu2func) with open(faprotax_report.mat, w) as out: out.write(OTUID\t \t.join(funcs) \n) for a in otus: row [1 if g in otu2func[a] else 0 for g in funcs] out.write(a \t \t.join(row) \n) print(f{len(otus)} OTUs, {len(funcs)} functions, {len(pairs)} pairs)跑法是python report2tables.py faprotax_report.clean ASV_第二個參數(shù)就是 grep 用的前綴兩邊保持一致。最后打印的三個數(shù)字可以當(dāng)作快速自檢OTU 數(shù)、功能組數(shù)、配對關(guān)系數(shù)。配對關(guān)系數(shù)應(yīng)該明顯大于功能組數(shù)如果兩者接近說明大部分功能組只對應(yīng)一條記錄樣本的功能譜很可能被少數(shù)類群撐起來解讀時要收著點(diǎn)。.mat是 OTU 層級的有無矩陣適合展示“某個功能由哪些 ASV 分擔(dān)”.otu_func適合在補(bǔ)充材料里做逐條溯源.func_otu適合在正文表格里寫“固氮功能由 12 個 ASV 貢獻(xiàn)歸屬 3 個屬”。4.3 下游三種圖的數(shù)據(jù)準(zhǔn)備主結(jié)果faprotax.txt進(jìn)分析之前要做標(biāo)準(zhǔn)化。因?yàn)楣δ芙M豐度受測序深度影響樣本間直接比大小沒有意義。常見做法是轉(zhuǎn)成相對豐度每列除以該列總和或者用 TSStotal sum scaling也有人用 centered log-ratio 處理組成型數(shù)據(jù)。抽平后的數(shù)據(jù)我一般先用相對豐度熱圖里再做行方向的 z-score視覺上更好讀。熱圖的行是功能組列是樣本按實(shí)驗(yàn)分組排序配上分組色條就能看到“某個處理下硝化相關(guān)功能整體抬升”這種模式。分組比較用 STAMP 或者 R 里的 Wilcoxon BH 校正注意功能組數(shù)量少多重檢驗(yàn)校正后很容易全不顯著這時候更該看效應(yīng)量和具體功能組的生物學(xué)含義而不是死盯 q 值。和環(huán)境因子的關(guān)聯(lián)分析是 FAPROTAX 的高光場景功能組豐度矩陣與 pH、總氮、有機(jī)質(zhì)這些理化指標(biāo)做 Spearman 相關(guān)畫成氣泡圖橫縱軸是功能組和環(huán)境因子氣泡大小表示相關(guān)系數(shù)絕對值顏色表示方向。因?yàn)楣δ芙M數(shù)量有限一張圖能塞下全部比 KO 通路那種幾百行的矩陣清爽得多。做相關(guān)時樣本量小的項(xiàng)目要謹(jǐn)慎n 小于 10 時相關(guān)系數(shù)的可信區(qū)間非常寬圖可以畫結(jié)論別寫死。5. 常見問題排查五個把我卡住的現(xiàn)場5.1 匹配數(shù)為零或個位數(shù)現(xiàn)象命令跑完屏幕打印Assigned 5 records to groups兩千多條 ASV 幾乎全是 leftovers。原因八成出在分類字符串上。最常見的是最后一級不是屬/種而是g__uncultured、f__Xanthobacteraceae這類其次是從上游帶下來的g__前綴后面跟了空格或不可見字符還有一種是分類順序被打亂屬名跑到了門的前面。解決先抽幾條看實(shí)際字符串a(chǎn)wk -F\t NR3{print $2} taxonomy2.txt直接肉眼確認(rèn)再統(tǒng)計最后一級的分布如果大量卡在科級說明上游注釋器只給到科FAPROTAX 無能為力要么換注釋工具重跑要么在文章里如實(shí)寫“僅部分類群可歸到屬級”?？崭窈吞厥庾址胹ed s/ //g清一遍再重新掛注釋。5.2 taxonomy 掛上去了但腳本報找不到字段現(xiàn)象collapse_table.py提示 metadata 里沒有taxonomy字段或者歸并結(jié)果為空。原因biom add-metadata的--observation-header順序和注釋文件列順序不一致導(dǎo)致字段名和內(nèi)容錯位或者--sc-separated拼寫成了別的字段名。解決用biom convert -i xxx.biom -o - --to-tsv | head -3把 biom 打回文本直接看 observation metadata 那幾列到底叫什么、值長什么樣。這一步比猜快得多。5.3 biom 命令報版本或依賴錯誤現(xiàn)象biom convert拋ImportError或者與 QIIME2 的命令互相覆蓋。原因QIIME2 環(huán)境里的 biom 是綁定版本和獨(dú)立安裝的 biom-format 沖突或者在一個 Python 2 的舊環(huán)境里跑 1.2.4 的腳本。解決按 3.1 的方式單獨(dú)建環(huán)境全程只在這個環(huán)境里操作。轉(zhuǎn)換和預(yù)測分兩個環(huán)境跑也可以中間用文件傳遞不要試圖在一個環(huán)境里同時滿足所有工具的依賴。5.4 把功能組豐度當(dāng)絕對量解讀現(xiàn)象某樣本固氮功能豐度是另一組的 3 倍直接寫“固氮能力增強(qiáng) 3 倍”。原因功能組豐度是匹配到的 OTU 豐度加和既受該類群相對豐度影響也受測序深度和抽平方式影響還受覆蓋度影響——不同樣本匹配上的類群集合可能根本不一樣。解決統(tǒng)一抽平深度后再比較報告相對豐度而非絕對 reads同時打開 report 看兩組匹配到的 OTU 是否同源如果一組匹配到 20 個 ASV、另一組只有 3 個這個倍數(shù)沒有意義。表述改成“固氮相關(guān)功能組相對豐度顯著升高”并附上貢獻(xiàn) OTU 列表。5.5 覆蓋度低就宣布“樣本沒有該功能”現(xiàn)象結(jié)果里沒有methanotrophy直接得出“該環(huán)境不存在甲烷氧化菌”的結(jié)論。原因FAPROTAX 是陽性證據(jù)庫沒匹配上只代表“沒有已知的、被文獻(xiàn)收錄的類群”不代表功能不存在。未培養(yǎng)微生物、新屬、數(shù)據(jù)庫未收錄的類群全都落在 leftovers 里。解決先報覆蓋度再談陰性結(jié)果。表述改成“在 FAPROTAX 可注釋的范圍內(nèi)未檢出”并在討論里說明未培養(yǎng)類群的可能貢獻(xiàn)。有條件的話用 qPCR 或宏基因組對關(guān)鍵功能基因做一次驗(yàn)證把話說圓。6. 讓結(jié)果經(jīng)得起追問抽平、注釋庫與證據(jù)鏈功能預(yù)測這類工作最容易被問的不是“你怎么跑的”而是“你憑什么信”。三個動作能把結(jié)果的可信度抬一個臺階。第一個是抽平的一致性。所有樣本必須抽到同一深度再算功能否則高深度樣本匹配到的 OTU 更多功能組豐度天然偏高組間差異會摻進(jìn)測序深度的假信號。抽平深度取樣本 reads 數(shù)的最小值或去掉極端低值后的分位數(shù)抽三次看功能組排序是否穩(wěn)定穩(wěn)了再定稿。第二個是注釋庫的選擇。FAPROTAX 不限注釋庫但經(jīng)驗(yàn)上不同庫的屬名寫法、層級完整性差別很大同一批數(shù)據(jù)換庫重跑匹配率能差出一大截。我一般用 SILVA 132 打底它在屬級的信息完整度和 FAPROTAX 的命名習(xí)慣比較對得上也見過用 RDP 注釋后匹配率明顯更高的項(xiàng)目這跟樣本類群有關(guān)沒有普適答案。穩(wěn)妥做法是拿一小批樣本分別用兩個庫跑一遍比一下 Assigned 比例再決定用哪個別直接上全量。第三個是證據(jù)鏈。每一張功能熱圖、每一個顯著差異的功能組都要能順著 report 反查到具體 ASV 和分類。寫補(bǔ)充材料時把func_otu表附上審稿人問“這個功能由誰貢獻(xiàn)”一查就有。下面這張表是我常用的自檢清單檢查項(xiàng)看什么合格線經(jīng)驗(yàn)值匹配覆蓋度Assigned / 總記錄數(shù)環(huán)境樣本 30%~60%低于 15% 需復(fù)查注釋抽平一致性三次抽平后功能組排序相關(guān)性Spearman 0.9功能組冗余每個功能組對應(yīng)的 ASV 數(shù)多數(shù)功能組應(yīng)由多個 ASV 貢獻(xiàn)單一 ASV 主導(dǎo)需說明差異結(jié)果穩(wěn)健性換分組檢驗(yàn)方法Wilcoxon / t 檢驗(yàn)主要結(jié)論方向一致外部驗(yàn)證與理化指標(biāo)或 qPCR 的相關(guān)方向關(guān)鍵功能組方向不矛盾還有一個具體技巧值得單獨(dú)說做組間比較前先把faprotax.txt里在所有樣本中相對豐度都低于 0.01% 的功能組過濾掉。這些低豐度功能組本身由一兩條 ASV 撐著稍一抽平就消失留著只會制造假陽性讓多重檢驗(yàn)校正后的結(jié)果更難看。過濾閾值別拍腦袋按你數(shù)據(jù)的總豐度分布來定寫進(jìn)方法里就行。最后提一句和實(shí)測數(shù)據(jù)的交叉驗(yàn)證。如果項(xiàng)目里同時測了總氮、銨態(tài)氮、pH把硝化、反硝化、固氮這幾個功能組單獨(dú)拎出來做散點(diǎn)或相關(guān)分析方向?qū)Φ蒙瞎适戮烷]環(huán)了方向?qū)Σ簧舷葎e急著刪數(shù)據(jù)回頭查抽平和注釋庫很多時候問題就出在這兩處。從那以后我每次跑 FAPROTAX都強(qiáng)制先看 Assigned 比例、再抽平三次自檢、最后反查一遍 report 里的貢獻(xiàn) OTU這三步不走完后面的圖一張都不出。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取