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