簡介圍繞16S擴(kuò)增子功能預(yù)測與元素循環(huán)解讀的PDF講義聚焦FAPROTAX將OTU分類表與人工整理的功能注釋數(shù)據(jù)庫關(guān)聯(lián)用于海洋、湖泊等環(huán)境樣本中碳、氫、氮、磷、硫循環(huán)的功能預(yù)測。內(nèi)容適合具備一定微生物組分析基礎(chǔ)的學(xué)生與研究者用于理解FAPROTAX的數(shù)據(jù)庫來源、預(yù)測邊界、結(jié)果矩陣形式以及功能豐度熱圖、環(huán)境因子相關(guān)性泡泡圖等呈現(xiàn)方式。包內(nèi)共1個(gè)PDF文件約2.47MB篇幅集中便于作為專題課配套資料查閱。資源還結(jié)合哺乳動(dòng)物與人類皮膚微生物、水稻氮利用等案例討論功能預(yù)測在差異挖掘與基因/代謝通路關(guān)聯(lián)中的用法并提及ImageGP一鍵分析與Linux本地分析兩條路徑的適用場景及依賴、輸入表等常見問題。已有1148人學(xué)習(xí)適合需要快速建立FAPROTAX知識(shí)框架、補(bǔ)充元素循環(huán)分析思路的讀者。1. 16S 數(shù)據(jù)做完物種注釋之后為什么我先補(bǔ)一層元素循環(huán)功能預(yù)測擴(kuò)增子跑完OTU 表和物種注釋擺在面前緊接著一定會(huì)有人問一句這些菌到底在干什么碳氮磷硫循環(huán)里誰在干活。PICRUSt2 給的是一堆 KO 編號(hào)Tax4Fun 給的是通路百分比而 FAPROTAX 直接把 OTU 翻譯成「硝化」「反硝化」「甲烷氧化」「硫酸鹽還原」這種能直接寫進(jìn)結(jié)論的功能條目。它源自 2016 年那篇把海洋微生物組功能與分類解耦的工作數(shù)據(jù)庫靠人工從可培養(yǎng)菌文獻(xiàn)整理全部代碼就一個(gè)collapse_table.py。覆蓋度確實(shí)不高但每條注釋都能追到文獻(xiàn)做環(huán)境樣本、農(nóng)田土壤、湖泊沉積物的元素循環(huán)題目時(shí)比一堆 KO 好解釋得多。2. 拆開 FAPROTAX 的黑匣子一張二值表加一個(gè) python 腳本動(dòng)手之前先把這套東西的內(nèi)部結(jié)構(gòu)看清楚后面報(bào)錯(cuò)時(shí)才知道該往哪兒找。FAPROTAX 沒有模型、沒有隱層、沒有任何機(jī)器學(xué)習(xí)成分它就是一次查表加一次求和把所有不確定性都?jí)涸凇笖?shù)據(jù)庫里到底有沒有你這個(gè)屬」這一件事上。2.1 數(shù)據(jù)庫長什么樣菌與功能的二值映射FAPROTAX.txt是純文本、制表符分隔。第一列是分類單元名絕大多數(shù)落在屬一級(jí)少量到種或科后面每一列是一個(gè)功能條目交叉位置用 1/0 標(biāo)記該分類單元是否被文獻(xiàn)支持具備這項(xiàng)功能。整張表本質(zhì)是一個(gè)二部圖一個(gè)屬可以同時(shí)掛上硝化、亞硝酸鹽氧化、好氧氨氧化好幾個(gè)標(biāo)簽一個(gè)功能條目下面也拴著幾十上百個(gè)屬。腳本做的工作就是把這個(gè)二部圖和你的 OTU 表做乘法——OTU 按注釋歸到屬屬映射到功能同一功能下所有 OTU 的豐度加起來。這里有個(gè)必須提前知道的細(xì)節(jié)功能條目的粒度差別很大。nitrification是粗線條的匯總條目下面還掛著aerobic_ammonia_oxidation、aerobic_nitrite_oxidation這些細(xì)分條目跑完的結(jié)果里兩類會(huì)同時(shí)出現(xiàn)。如果你直接拿去畫熱圖同一個(gè) OTU 會(huì)被數(shù)兩遍圖上呈現(xiàn)的就是自相關(guān)審稿人一眼能看出來。我一般會(huì)在出圖前把匯總級(jí)別的條目單獨(dú)拎出來只留細(xì)分條目做展示。還有一類條目要特別小心像fermentation、chemoheterotrophy、aerobic_chemoheterotrophy這種范圍極寬的幾乎每個(gè)樣本都有、豐度還都不低做組間差異分析永遠(yuǎn)不顯著。它們不是錯(cuò)只是信息量太低屬于「知道就行、別往結(jié)論里寫」的那一檔。2.2 和 KO 通路、表型預(yù)測工具的分工同一個(gè) OTU 表可以喂給好幾個(gè)預(yù)測工具選哪個(gè)取決于你要回答什么問題。下面這張表是我自己在項(xiàng)目里做選型時(shí)用的對(duì)照參數(shù)和輸出形態(tài)都按實(shí)際跑出來的結(jié)果填工具輸出內(nèi)容參考依據(jù)適合回答的問題主要短板FAPROTAX 1.2.4具體元素循環(huán)功能條目人工整理的可培養(yǎng)菌文獻(xiàn)碳氮磷硫循環(huán)里誰在動(dòng)覆蓋度低注釋率常在 20%~60%PICRUSt2KO / EC / 代謝通路基因組數(shù)據(jù)庫外推通路層面的整體代謝潛力中間步驟多結(jié)果難落到具體過程Tax4Fun2KO 通路百分比16S 參考序列庫快速通路比較準(zhǔn)確率高度依賴參考庫覆蓋表型預(yù)測類工具革蘭氏、氧氣利用、致病性已有表型數(shù)據(jù)庫群落表型的粗略畫像表型條目有限粒度粗選型邏輯其實(shí)很樸素題目落在元素生物地球化學(xué)循環(huán)上樣本來自海洋、湖泊、沉積物、濕地、農(nóng)田這類環(huán)境FAPROTAX 是對(duì)口的題目落在宿主相關(guān)樣本上也能用但覆蓋度會(huì)明顯掉一截因?yàn)樗牡鬃邮强膳囵B(yǎng)的環(huán)境菌。有研究把它用在哺乳動(dòng)物和人類皮膚菌群的比較上也有人把它和秈稻根系氮利用的表型關(guān)聯(lián)起來這兩類工作能發(fā)出來靠的不是覆蓋度而是「命中的那部分功能確實(shí)對(duì)得上已知機(jī)制」。2.3 輸入輸出的形狀四個(gè)文件之間的關(guān)系輸入側(cè)只有兩樣?xùn)|西一份不帶注釋的 OTU 豐度表轉(zhuǎn)成 biom 格式和一份把 OTUID 對(duì)上 taxonomy 字符串的元數(shù)據(jù)文件。腳本默認(rèn)取注釋字符串的最后一級(jí)去和數(shù)據(jù)庫匹配通常就是屬名所以如果你的注釋只到科匹配率會(huì)掉得很難看。輸出側(cè)是兩個(gè)文件。faprotax.txt是功能 × 樣本的豐度矩陣行是功能條目列是樣本數(shù)值是落在該功能下的所有 OTU 豐度之和。faprotax_report.txt是明細(xì)報(bào)告每個(gè)功能條目下列出命中了哪些 OTU末尾還有一行統(tǒng)計(jì)格式類似Assigned 979 records to groups, 1924 records were leftovers。別小看這個(gè) report后面寫討論、做圖注、解釋「為什么某個(gè)功能突然冒起來」全靠它把功能追回到具體序列。中間還有一個(gè)容易被忽略的產(chǎn)物biom 里掛上去的 observation metadata。它既不是輸入也不是輸出但整條鏈路成不成就看它第 5 章里一半的坑都出在這一步。3. 本地環(huán)境搭起來python3、numpy、biom 與 1.2.4 版腳本自檢這套腳本沒有網(wǎng)絡(luò)版依賴純本地跑好處是數(shù)據(jù)不出內(nèi)網(wǎng)、參數(shù)完全可控壞處是環(huán)境得自己收拾。1.2.4 版是 2020 年 11 月底更新的數(shù)據(jù)庫做了擴(kuò)充腳本本身也從 python2 遷到了 python3這個(gè)遷移是很多人翻車的源頭。3.1 下載與解包先認(rèn)清目錄里有什么壓縮包解壓后是一個(gè)以版本號(hào)命名的目錄里面至少要有三樣?xùn)|西collapse_table.py主腳本、FAPROTAX.txt數(shù)據(jù)庫、以及一份說明文件。三者的相對(duì)位置不要?jiǎng)幽_本按相對(duì)路徑去找數(shù)據(jù)庫是常見寫法挪了位置就得在命令行里顯式指定。# 假設(shè)壓縮包已經(jīng)下載到當(dāng)前目錄先解壓 unzip FAPROTAX_1.2.4.zip -d . # 看一眼目錄結(jié)構(gòu)確認(rèn)關(guān)鍵文件都在 ls -lh FAPROTAX_1.2.4/ # 預(yù)期看到 collapse_table.py、FAPROTAX.txt、README 之類操作本身沒什么難度容易出問題的是「版本混用」網(wǎng)上 2019 年以前流傳的教程命令里帶的是 python2 的寫法直接抄到 python3 環(huán)境里會(huì)報(bào)語法錯(cuò)?？吹綀?bào)錯(cuò)信息里出現(xiàn)print語句相關(guān)的SyntaxError基本可以斷定你抄了老教程或者手里拿的是老版本腳本。判斷依據(jù)很簡單python --version輸出 3.x腳本也必須是 1.2.4 這一版兩者對(duì)齊才談得上后面的事。3.2 conda 環(huán)境與依賴安裝為什么不建議裝在 base 里依賴只有兩個(gè)大頭numpy和biom。聽起來簡單但這兩者在不同渠道的包名和版本約束不一樣裝在 base 環(huán)境里很容易和你已有的其他生信工具打架。我一般會(huì)單開一個(gè)環(huán)境用完就扔出問題直接刪掉重建比修依賴快得多。# 新建一個(gè)干凈的 python3 環(huán)境名字隨意 conda create -n faprotax python3.9 -y conda activate faprotax # 裝依賴biom 在 conda 里的包名是 biom-format conda install -c conda-forge numpy biom-format -y # 確認(rèn)解釋器路徑落在新建的環(huán)境下避免裝完還是用系統(tǒng) python which pythonconda create里指定 python 版本這一步別省直接繼承系統(tǒng) python 有時(shí)候會(huì)帶進(jìn)一堆歷史包袱。biom-format這個(gè)包名值得記一下很多教程里寫的是conda install biom能不能裝上取決于頻道里有沒有同名包用-c conda-forge加全名是最穩(wěn)的寫法。如果你機(jī)器上已經(jīng)有一個(gè)裝好 QIIME2 的環(huán)境里面自帶 numpy 和 biom-format也可以直接激活它來跑省一次安裝。代價(jià)是 QIIME2 的依賴樹比較深偶爾會(huì)有版本鎖死的情況真出問題時(shí)排查成本比自己建環(huán)境高。選哪條路取決于你手頭有沒有現(xiàn)成的環(huán)境沒有就從零建。3.3 自檢一行命令彈出幫助才算裝好裝完不要急著上真實(shí)數(shù)據(jù)先用不帶參數(shù)的方式跑一次腳本。這類命令行工具通常在不給參數(shù)時(shí)打印幫助信息然后退出能正常彈出幫助就說明解釋器、依賴、腳本三樣都通了。# 不帶任何參數(shù)運(yùn)行預(yù)期打印 usage 幫助信息 python FAPROTAX_1.2.4/collapse_table.py如果這一步報(bào)ModuleNotFoundError: No module named numpy說明當(dāng)前激活的環(huán)境里沒有 numpy回到上一步執(zhí)行conda install -c conda-forge numpy -y。報(bào)No module named biom同理裝biom-format。報(bào)cant open file則是路徑寫錯(cuò)了用ls確認(rèn)腳本的實(shí)際位置。還有一種情況是命令能跑但立刻退出、什么也不打印這多半是腳本讀參數(shù)時(shí)提前 return 了把-h顯式加上再試一次。自檢這一步花三十秒能省掉后面半小時(shí)的無效排查屬于典型的后悔藥。4. 從 otutab.txt 到 faprotax.txt三步轉(zhuǎn)換加一次預(yù)測整條鏈路只有三段文本表轉(zhuǎn) biom、把物種注釋掛上去、跑 collapse。每一段都有它的格式約束任何一環(huán)不合規(guī)最壞的結(jié)果不是報(bào)錯(cuò)而是靜默地給你一張全零矩陣。4.1 txt 轉(zhuǎn) biom json兩個(gè)格式約束先滿足原始 OTU 表通常是制表符分隔的純文本第一行是#OTU ID加一串樣本名下面是數(shù)字。轉(zhuǎn)換前先做兩件事把最后一列如果有的話——通常是 taxonomy 字符串——單獨(dú)拆出去存成另一個(gè)文件確認(rèn)所有單元格是純數(shù)字不帶千分位逗號(hào)、不帶空格、不帶NA。# 文本 OTU 表轉(zhuǎn) biom 的 json 格式 biom convert -i otutab_rare.txt -o otutab_rare.biom \ --table-typeOTU table --to-json--table-type這個(gè)參數(shù)必須寫寫錯(cuò)了后面的工具鏈會(huì)按錯(cuò)誤的類型解析--to-json生成的是 json 版本而不是 HDF5好處是純文本、能直接head看內(nèi)容、出問題時(shí)肉眼可查對(duì)調(diào)試非常友好。輸入表里的樣本順序和分組表順序不要求一致biom 靠列名匹配但列名里的空格、中文、特殊符號(hào)建議提前清掉這些在后續(xù)各種工具里都是雷。轉(zhuǎn)換完成后建議立刻驗(yàn)證一遍用biom summarize-table -i otutab_rare.biom看每個(gè)樣本的序列數(shù)、總 OTU 數(shù)是否和原表對(duì)得上。這一步不做等到預(yù)測結(jié)果全零時(shí)再回頭找就得從頭查起了。4.2 add-metadata 把 taxonomy 掛上去表頭和分隔符是關(guān)鍵這一步是整個(gè)流程里最容易翻車的地方。需要準(zhǔn)備一個(gè)兩列的元數(shù)據(jù)文件第一列是 OTUID必須和 biom 里的 observation ID 逐字符一致第二列是 taxonomy 字符串各級(jí)之間用分號(hào)分隔。# 把物種注釋作為 observation metadata 掛到 biom 上 biom add-metadata -i otutab_rare.biom \ --observation-metadata-fp taxonomy2.txt \ -o otutab_rare_tax.biom \ --sc-separated taxonomy \ --observation-header OTUID,taxonomy--sc-separated taxonomy告訴程序這一列內(nèi)部是按分號(hào)切分的多級(jí)結(jié)構(gòu)切完會(huì)在 biom 里存成列表--observation-header OTUID,taxonomy是元數(shù)據(jù)文件的表頭聲明這兩個(gè)詞必須一字不差。寫成OTU_ID、#OTUID、OTUID,taxonomy_name都不會(huì)報(bào)錯(cuò)但掛上去的 taxonomy 是空的直到你跑完預(yù)測看到一張全零的矩陣才會(huì)發(fā)現(xiàn)。我見過不止一次有人在這上面耗掉一整天。掛完之后同樣要驗(yàn)biom table -i otutab_rare_tax.biom --observations隨便看幾條確認(rèn)每個(gè) OTU 后面都跟著一段注釋字符串。另外注意注釋字符串里最后一級(jí)是不是屬名中間有沒有混進(jìn)g__、unclassified、norank這類前綴或占位詞這些在匹配階段都是純噪聲。4.3 collapse_table.py 的參數(shù)逐個(gè)說清參數(shù)不多但每一個(gè)都直接決定結(jié)果對(duì)不對(duì)值得逐條過一遍參數(shù)含義常見取值與注意事項(xiàng)-i輸入 biom必須是掛了 taxonomy 的那個(gè)版本-g數(shù)據(jù)庫文件指向FAPROTAX.txt路徑寫錯(cuò)會(huì)直接報(bào)錯(cuò)或零注釋-o功能豐度矩陣行是功能、列是樣本-r明細(xì)報(bào)告每個(gè)功能下命中的 OTU 列表與統(tǒng)計(jì)行--collapse_by_metadata用哪一列做歸并寫taxonomy必須和上一步的列名一致-v打印過程信息建議開著能看到進(jìn)度和警告--force覆蓋已有輸出重跑時(shí)必加否則遇到同名文件直接退出對(duì)應(yīng)的完整命令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跑完之后終端會(huì)打印一行統(tǒng)計(jì)形如Assigned 979 records to groups, 1924 records were leftovers。前半句是成功映射到功能的記錄數(shù)后半句是沒匹配上的。這個(gè)比例就是你的注釋率它決定了后面所有分析的可信度看到數(shù)字的第一反應(yīng)應(yīng)該是判斷它合不合理而不是直接進(jìn)入畫圖環(huán)節(jié)。環(huán)境樣本的注釋率能到 40%~60% 算正常宿主樣本掉到 20% 以下也不奇怪如果環(huán)境樣本只有個(gè)位數(shù)百分比那一定是前面某一步出了問題。4.4 讀 report 判斷質(zhì)量先看統(tǒng)計(jì)行再看命中列表faprotax_report.txt的價(jià)值被嚴(yán)重低估。它按功能條目分塊每個(gè)塊下列出命中的 OTU 及其歸屬末尾帶統(tǒng)計(jì)。拿它做三件事判斷注釋率是否合理、看某個(gè)關(guān)鍵功能是被哪些 OTU 撐起來的、以及在討論部分解釋「這個(gè)功能的升高來自哪幾個(gè)屬」。如果要把 report 整理成可用的矩陣常見做法是先把每個(gè)功能塊的首行和它下面的 ASV 行篩出來去掉分隔線再按功能分組展開。下面的腳本做的是同一件事用 python 解析比記一堆 grep 參數(shù)更省心import collections # faprotax_report.clean 是篩掉分隔線后的報(bào)告功能名行 縮進(jìn)的 OTU 行交替出現(xiàn) func2otu collections.OrderedDict() cur None with open(faprotax_report.clean) as fh: for line in fh: line line.rstrip(\n) if not line.strip(): continue if line.startswith(ASV_) or line.startswith(OTU_): if cur: # 當(dāng)前功能下追加一條命中記錄 func2otu[cur].append(line.split(\t)[0]) else: cur line.strip() # 遇到功能名行切換當(dāng)前功能 func2otu.setdefault(cur, []) # 輸出「功能 - 命中OTU數(shù)」的速查表用來判斷哪些功能是靠單個(gè)OTU撐起來的 for f, otus in func2otu.items(): print(f, len(otus))這段腳本的輸出可以直接當(dāng)質(zhì)控表用。命中 OTU 數(shù)只有一兩個(gè)的功能條目在組間差異里極不穩(wěn)定換個(gè)抽平深度可能就消失了寫結(jié)論時(shí)要留余地。反過來命中幾十個(gè) OTU 的功能更穩(wěn)健但往往也就是那些寬泛條目差異不顯著。兩頭都要具體分析不能一刀切。5. 五條踩坑記錄注釋率、報(bào)錯(cuò)與結(jié)果對(duì)不上號(hào)的排查路徑下面這些是我和周圍人真實(shí)踩過、并且在別人數(shù)據(jù)上也復(fù)現(xiàn)過的問題。每條按「現(xiàn)象→原因→解決」寫遇到對(duì)得上的直接照方抓藥。5.1 現(xiàn)象Assigned 只有一百多條leftovers 兩三千原因通常有三個(gè)。最常見的是 taxonomy 只注到科腳本取最后一級(jí)拿到的是科名數(shù)據(jù)庫里根本沒有這一級(jí)其次是屬名帶了一堆前綴和占位符比如g__、unclassified、norank、uncultured這些字符串拿去匹配當(dāng)然一條都中不了第三是生境本身特殊數(shù)據(jù)庫的底子是可培養(yǎng)的環(huán)境菌你的樣本來自極端環(huán)境或宿主腸道庫里確實(shí)沒有對(duì)應(yīng)條目。解決順序建議從便宜到貴先寫個(gè)腳本把所有注釋字符串的最后一級(jí)抽出來統(tǒng)計(jì)有多少是unclassified之類的占位詞再把g__前綴統(tǒng)一去掉、屬名首字母大寫規(guī)范化如果注釋只到科換一套能到屬的參考庫重新注釋經(jīng)驗(yàn)上 132 版的 SILVA 和這套數(shù)據(jù)庫配合得比較好都做完注釋率還是低那就是生境問題接受它并在方法里如實(shí)寫明覆蓋度。5.2 現(xiàn)象跑腳本報(bào) ModuleNotFoundError明明剛裝過原因幾乎都是解釋器和環(huán)境不匹配。裝依賴時(shí)裝在 A 環(huán)境跑腳本時(shí)激活的是 B 環(huán)境或者干脆用了系統(tǒng)自帶的 python。conda 環(huán)境切換不徹底是很常見的操作失誤。解決運(yùn)行which python和which conda確認(rèn)前者路徑里帶envs/你的環(huán)境名。如果不對(duì)conda deactivate之后重新conda activate再不行就直接用環(huán)境里的絕對(duì)路徑調(diào)用解釋器類似~/miniconda3/envs/faprotax/bin/python collapse_table.py。這條路徑寫死在命令行里雖然難看但絕對(duì)不會(huì)跑錯(cuò)環(huán)境。5.3 現(xiàn)象預(yù)測跑完了faprotax.txt 全是 0原因八成在上一步的 add-metadata。要么--observation-header寫錯(cuò)導(dǎo)致 taxonomy 列沒掛上要么元數(shù)據(jù)文件里的 OTUID 和 biom 里的不一致——多了個(gè)#、少了ASV_前綴、大小寫不同、行尾帶\r任何一種都能讓匹配全軍覆沒而且整個(gè)過程不報(bào)錯(cuò)。解決先用biom table -i otutab_rare_tax.biom --observations | head看前幾條確認(rèn) taxonomy 字段真的有內(nèi)容。如果為空回去查表頭如果有內(nèi)容但預(yù)測結(jié)果還是全零就對(duì)比兩邊的 ID我一般會(huì)cut -f1 taxonomy2.txt | sort a.txt、把 biom 的 observation ID 也導(dǎo)出來排序diff一下差異一目了然。Windows 上生成的文件記得dos2unix。5.4 現(xiàn)象換一套物種注釋重跑兩次結(jié)果相關(guān)性只有 0.6原因不是腳本不穩(wěn)定而是輸入變了。不同參考庫對(duì)同一個(gè) OTU 的屬級(jí)歸屬本來就可能不同能命中數(shù)據(jù)庫的 OTU 集合隨之改變功能豐度自然整體漂移。如果兩次流程里還有一次做了抽平、另一次沒做差異會(huì)更大。解決同一批數(shù)據(jù)從頭到尾鎖死一套注釋流程中途不換庫確實(shí)需要比較時(shí)先把兩邊都轉(zhuǎn)成相對(duì)豐度再算相關(guān)更嚴(yán)謹(jǐn)?shù)淖龇ㄊ前褍纱?report 里命中的 OTU 列表取交集只在交集上比較功能豐度把「注釋差異」這個(gè)變量控制住。5.5 現(xiàn)象某個(gè)功能在所有樣本里豐度幾乎一模一樣原因有兩種。一種是這個(gè)功能命中的 OTU 極少且豐度很低抽平之后被隨機(jī)性主導(dǎo)看起來就像常數(shù)另一種是命中的是那幾個(gè)無處不在的寬泛條目本來就沒什么組間變化。解決先用 4.4 那段腳本看該功能命中了幾個(gè) OTU只有一兩個(gè)的直接從圖里刪掉別硬畫如果是寬泛條目保留在表里但在正文里不展開。低豐度 OTU 在抽平前先做一次過濾比如保留在至少 10% 樣本中出現(xiàn)的 OTU能明顯減少這類噪聲。6. 拿到功能矩陣之后歸一化、追到 OTU、用實(shí)測指標(biāo)兜底faprotax.txt里是絕對(duì)豐度加和直接跨樣本比 reads 數(shù)沒有意義測序深度和抽平方式都會(huì)影響結(jié)果。第一步永遠(yuǎn)是轉(zhuǎn)相對(duì)豐度把每個(gè)樣本的功能豐度除以該樣本的總和再乘 100import pandas as pd df pd.read_csv(faprotax.txt, sep\t, index_col0) # 行是功能列是樣本 df df.div(df.sum(axis0), axis1) * 100 # 按樣本列歸一化成百分比 df df.loc[df.sum(axis1) 0.01] # 剔掉幾乎不出現(xiàn)的功能 df.to_csv(faprotax_rel.txt, sep\t)歸一化之后才談得上熱圖、氣泡圖、組間差異。做環(huán)境因子相關(guān)的泡泡圖時(shí)功能這一側(cè)建議只保留與碳氮磷硫直接相關(guān)的細(xì)分條目把寬泛條目擋在外面否則圖上會(huì)擠滿一堆不顯著的點(diǎn)。第二步是把功能追回到 OTU。這一步很多人跳過結(jié)果討論部分只能寫「硝化功能在組間顯著升高」說不清是誰升高的。用第 4.4 節(jié)的腳本把 report 拆成「功能—OTU」對(duì)應(yīng)表再和 OTU 表的物種注釋一join就能落到具體的屬。我做氮循環(huán)相關(guān)題目時(shí)這一步幾乎必做因?yàn)樗軒湍闩袛嘤^察到的升高是真信號(hào)還是某個(gè)低豐度 OTU 的隨機(jī)波動(dòng)。第三步也是最能提升說服力的一步用實(shí)測指標(biāo)做錨點(diǎn)。功能預(yù)測終究是預(yù)測如果條件允許同一批樣本補(bǔ)測幾個(gè)關(guān)鍵基因的 qPCR硝化補(bǔ) amoA、反硝化補(bǔ) nirK/nirS、甲烷氧化補(bǔ) pmoA把預(yù)測豐度和實(shí)測拷貝數(shù)做個(gè)相關(guān)。方向一致后面的機(jī)制討論才站得住方向相反那就是個(gè)明確的信號(hào)說明這套數(shù)據(jù)庫在你的生境里不適用趁早換工具比硬寫強(qiáng)。沒有 qPCR 條件的退一步用環(huán)境因子的實(shí)測值做相關(guān)也行銨態(tài)氮、硝態(tài)氮、有機(jī)碳這些理化指標(biāo)很多項(xiàng)目本來就測了。血淚教訓(xùn)是我曾經(jīng)拿著一份注釋率只有 12% 的結(jié)果寫了一整節(jié)機(jī)制討論投稿后被審稿人一句「你這條通路只由一個(gè)未培養(yǎng)屬的單個(gè) OTU 支撐」直接打回。從那以后我每跑完一次collapse_table.py都強(qiáng)制先把 Assigned 比例、每個(gè)功能的命中 OTU 數(shù)、以及關(guān)鍵功能的屬級(jí)歸屬這三樣抄進(jìn)分析日志達(dá)不到預(yù)期就回頭查輸入絕不帶著可疑的結(jié)果往下走。希望幫到你。本文還有配套的精品資源點(diǎn)擊獲取