做微生物組數(shù)據(jù)分析這幾年有一個很深的體會真正決定項(xiàng)目成敗的往往不是高端算法或復(fù)雜模型而是上游那套看起來平平無奇的“數(shù)據(jù)預(yù)處理”流程。從OTU/ASV過濾到相對豐度轉(zhuǎn)換每一個步驟都暗藏著陷阱稍不留神后面的統(tǒng)計(jì)檢驗(yàn)、差異分析、機(jī)器學(xué)習(xí)全都會跑偏。這篇內(nèi)容就圍繞我實(shí)際項(xiàng)目中反復(fù)打磨過的一套標(biāo)準(zhǔn)化流程展開把每一步的“為什么”和“怎么做”都講清楚希望能給正在被微生物組數(shù)據(jù)折磨的人一些參考。這套流程適合誰剛接觸微生物組數(shù)據(jù)分析、用QIIME 2或phyloseq做下游分析的人以及已經(jīng)跑通流程但總覺得結(jié)果不穩(wěn)定、想梳理優(yōu)化分析方案的研究者。文章不會堆砌理論主要講實(shí)操但每個關(guān)鍵環(huán)節(jié)都會解釋背后的統(tǒng)計(jì)學(xué)和生物學(xué)邏輯確保你能理解每一步的意義而不是盲目照搬代碼。1. 內(nèi)容整體設(shè)計(jì)與思路拆解1.1 為什么需要一套標(biāo)準(zhǔn)化的預(yù)處理流程微生物組數(shù)據(jù)最典型的特征就是高維、稀疏、成分性。所謂高維通常一個樣本里有成千上萬個OTU或ASV稀疏指的是絕大多數(shù)樣本中只有少部分微生物能被檢測到很多分類單元在多數(shù)樣本中豐度為零成分性則是說測序得到的豐度只是相對比例而不是絕對數(shù)量各分類單元之間存在“此消彼長”的約束關(guān)系。這三個特征疊加導(dǎo)致如果你不做預(yù)處理后續(xù)分析基本沒法看。舉個例子如果不過濾低豐度的OTU/ASV那些來自測序錯誤或污染物的“幽靈分類群”會占據(jù)大量特征維度稀釋真正有生物學(xué)意義的信號。如果不對樣本做標(biāo)準(zhǔn)化就直接比較豐度測序深度高的樣本天然會“顯得”物種更多、豐度更高但這并不是真實(shí)的生物學(xué)差異只是技術(shù)噪聲。標(biāo)準(zhǔn)化的預(yù)處理流程本質(zhì)上是在做三件事降噪、降維、可比化。降噪是把可疑的測序錯誤和污染信號剔除降維是減少下游多重檢驗(yàn)的負(fù)擔(dān)提升統(tǒng)計(jì)功效可比化則是通過相對豐度等方式讓樣本之間具有可比性。這三件事做好了后續(xù)無論做Alpha多樣性、Beta多樣性、差異豐度分析還是功能預(yù)測結(jié)果都會穩(wěn)健很多。1.2 方案選型OTU還是ASV怎么選才不踩坑開始實(shí)際分析前要明確一個基礎(chǔ)問題使用OTU還是ASV。這不是一個可以隨手選的問題它直接影響后續(xù)所有結(jié)果的分辨率和可比性。OTUOperational Taxonomic Unit操作分類單元是傳統(tǒng)聚類策略的產(chǎn)物通常按97%序列相似度聚類相當(dāng)于把相似的序列“合并同類項(xiàng)”。優(yōu)點(diǎn)是兼容性好很多經(jīng)典數(shù)據(jù)庫和參考文獻(xiàn)都基于OTU構(gòu)建計(jì)算資源消耗相對較低。缺點(diǎn)是聚類閾值是人為定義的不同研究之間、不同聚類算法之間差異明顯而且聚類過程容易把近緣物種合并丟失細(xì)微的生物學(xué)差異。ASVAmplicon Sequence Variant擴(kuò)增子序列變體則走的是“去噪”路子例如DADA2、Deblur這類算法直接把序列解析到單堿基分辨率不需要聚類閾值。優(yōu)點(diǎn)是分辨率高、可重復(fù)性好同一批數(shù)據(jù)用不同版本軟件處理結(jié)果依然高度一致跨研究比較也更可靠。缺點(diǎn)是計(jì)算開銷大對測序質(zhì)量要求高而且在一些數(shù)據(jù)庫中匹配率可能不如OTU理想。我的建議是除非你要跟歷史數(shù)據(jù)集做嚴(yán)格對比、對方明確使用OTU否則新建項(xiàng)目優(yōu)先選ASV。理由很直接——ASV在分辨率和可重復(fù)性上的優(yōu)勢是OTU無法提供的。而且各個主流流程和數(shù)據(jù)庫對ASV的兼容性已經(jīng)非常成熟沒必要為了保守而犧牲數(shù)據(jù)質(zhì)量。1.3 整個流程的核心節(jié)點(diǎn)把標(biāo)準(zhǔn)流程拆開核心節(jié)點(diǎn)大致是這樣特征表生成OTU/ASV表包含每個樣本中每個分類單元的序列數(shù)特征表過濾剔除低豐度、低流行度的可疑特征樣本測序深度檢查與最小深度確認(rèn)標(biāo)準(zhǔn)化與相對豐度轉(zhuǎn)換后續(xù)下游分析前的質(zhì)量核查很多人會把精力全部放在第一步和第五步忽略中間環(huán)節(jié)。實(shí)際上第二步和第四步才是決定分析質(zhì)量的分水嶺。后面我詳細(xì)展開這兩個環(huán)節(jié)的細(xì)節(jié)。2. 核心細(xì)節(jié)解析與實(shí)操要點(diǎn)2.1 過濾前先搞清楚你的數(shù)據(jù)長什么樣拿到特征表后不要著急寫過濾代碼先花點(diǎn)時(shí)間把數(shù)據(jù)的基本盤摸清楚。以R語言phyloseq對象為例我一般會先跑一遍摘要library(phyloseq) ps - readRDS(your_phyloseq_object.rds) ps # phyloseq-class experiment-level object # otu_table() OTU Table: [ 3428 taxa and 120 samples ] # sample_data() Sample Data: [ 120 samples by 8 sample variables ] # tax_table() Tax Table: [ 3428 taxa by 7 taxonomic ranks ]這里能看到總共有多少taxa、多少樣本。接下來要檢查測序深度分布sample_sums(ps) %% summary() # Min. 1st Qu. Median Mean 3rd Qu. Max. # 10532 25634 35120 40211 49220 120453測序深度的分布很關(guān)鍵。如果有些樣本只有幾千條序列有些樣本有十幾萬條差異過大的話后續(xù)標(biāo)準(zhǔn)化壓力會很大而且低深度樣本本身的生物信息可信度存疑。此時(shí)就要考慮是否剔除極低深度的樣本。具體閾值沒有絕對標(biāo)準(zhǔn)但有個經(jīng)驗(yàn)值對16S rRNA基因擴(kuò)增子測序樣本序列數(shù)低于5000甚至10000物種組成的穩(wěn)定性會顯著下降。當(dāng)然這個值還要結(jié)合你研究的環(huán)境類型、預(yù)期物種豐富度來定。比如腸道樣本本身復(fù)雜度高低深度損失的信息更多極端環(huán)境樣本復(fù)雜度低低深度也可能夠用。2.2 低豐度過濾的邊界在哪里低豐度過濾是爭議最多的環(huán)節(jié)因?yàn)椤暗汀钡臉?biāo)準(zhǔn)實(shí)在模糊。我的思路是分層處理而不是一刀切。首先看絕對豐度。如果某個OTU/ASV在整個數(shù)據(jù)集中總讀數(shù)只有個位數(shù)比如總共出現(xiàn)2次、3次這種特征極可能來自測序錯誤或交叉污染。DADA2這類去噪流程已經(jīng)去除了一部分錯誤序列但也不能完全排除。把總讀數(shù)低于某個閾值的特征刪掉是常見做法。閾值可以參考總序列數(shù)的萬分之一也可以根據(jù)負(fù)對照的結(jié)果來確定。其次看流行度prevalence即一個特征在多少比例的樣本中出現(xiàn)。如果一個OTU/ASV只在1個樣本中出現(xiàn)哪怕豐度不低也可能是罕見污染物或個體偶然攜帶對群體水平的分析貢獻(xiàn)有限還會增加多重檢驗(yàn)的負(fù)擔(dān)。實(shí)操中我常用的組合是保留至少在5%~10%樣本中出現(xiàn)的特征同時(shí)去掉絕對總豐度過低的特征。舉個具體例子120個樣本的數(shù)據(jù)集我會先篩掉總豐度少于10的特征再篩掉只在5個以下樣本中出現(xiàn)的特征。這兩個條件組合使用能去掉一大半噪聲特征同時(shí)保留核心微生物群落信號。2.3 相對豐度轉(zhuǎn)換為什么不是簡單的“除以總數(shù)”相對豐度轉(zhuǎn)換的常規(guī)操作是每個樣本中每個特征的讀數(shù)除以該樣本的總讀數(shù)得到0-100%或0-1的比例。這一步本身不復(fù)雜但很多人忽略了背后的兩個問題。第一個問題是“閉合效應(yīng)”closure。一旦轉(zhuǎn)換為相對豐度各特征之間就產(chǎn)生了負(fù)相關(guān)約束一個特征占比上升其他特征必然下降這不是生物學(xué)上的此消彼長而是數(shù)學(xué)上的強(qiáng)制約束。所以在后續(xù)分析中不能直接拿相對豐度去做常規(guī)的相關(guān)系數(shù)分析比如Pearson相關(guān)否則容易得到大量虛假相關(guān)。第二個問題是方差穩(wěn)定性。相對豐度數(shù)據(jù)通常高度偏態(tài)優(yōu)勢物種占比高、稀有物種占比極低直接用原始比例做統(tǒng)計(jì)檢驗(yàn)可能違背正態(tài)性假設(shè)。這時(shí)候需要考慮進(jìn)一步變換比如CLR中心對數(shù)比變換或ILR等距對數(shù)比變換這些會在完整標(biāo)準(zhǔn)化方案中用到。所以相對豐度轉(zhuǎn)換要區(qū)分場景。如果只是做可視化展示展示群落組成堆疊柱狀圖、餅圖直接按總數(shù)縮放沒問題。但如果你要用相對豐度做差異檢驗(yàn)、相關(guān)網(wǎng)絡(luò)或機(jī)器學(xué)習(xí)就要考慮額外處理比如使用CLR變換后的值。3. 實(shí)操過程與核心環(huán)節(jié)實(shí)現(xiàn)3.1 基于phyloseq的特征過濾實(shí)操下面是我在實(shí)戰(zhàn)中反復(fù)打磨的一段過濾代碼兼顧了低豐度過濾和低流行度過濾library(phyloseq) library(tidyverse) # Step 1: 查看每個特征的總豐度確定合適的絕對豐度閾值 taxa_sums - taxa_sums(ps) summary(taxa_sums) # 常見做法過濾掉總豐度為0的特征偶爾會出現(xiàn) ps - prune_taxa(taxa_sums(ps) 0, ps) # Step 2: 按最小總豐度過濾 min_total_abundance - 10 ps_filtered - prune_taxa(taxa_sums(ps) min_total_abundance, ps) # Step 3: 按流行度過濾 min_prevalence - 5 # 至少在5個樣本中檢測到 prev_df - as.data.frame(otu_table(ps_filtered)) prev_count - rowSums(prev_df 0) keep_taxa - prev_count min_prevalence ps_filtered - prune_taxa(keep_taxa, ps_filtered) # 查看過濾前后的變化 ps ps_filtered這段代碼的思路很直白先按絕對豐度粗篩再按流行度細(xì)篩。在實(shí)際項(xiàng)目中我通常會把兩個條件做成循環(huán)或函數(shù)針對不同數(shù)據(jù)集快速調(diào)整參數(shù)。有一點(diǎn)要特別提醒如果設(shè)置了最小流行度強(qiáng)烈建議結(jié)合樣本分組來看。比如研究的是疾病vs對照如果某個特征只在疾病組出現(xiàn)、對照完全沒有那它對差異分析是有意義的不能因?yàn)樗傮w流行度不高就刪掉。所以更穩(wěn)健的做法是分組計(jì)算流行度保留“在任一組中達(dá)到流行度閾值”的特征。3.2 使用microbiome包快速實(shí)現(xiàn)核心特征篩選如果你想更省事可以直接用microbiome包的core()函數(shù)。它的理念就是按流行度和豐度雙重閾值篩選核心微生物組代碼干凈多了library(microbiome) # 將phyloseq對象轉(zhuǎn)換為microbiome可處理的格式 ps_comp - microbiome::transform(ps_filtered, compositional) # 篩選在50%樣本中出現(xiàn)且相對豐度不低于0.1%的特征 ps_core - core(ps_comp, detection 0.001, prevalence 0.5) ps_core這里的detection是相對豐度閾值prevalence是樣本比例閾值。需要注意transform(ps, compositional)已經(jīng)做了相對豐度轉(zhuǎn)換所以detection參數(shù)用的是比例而不是序列數(shù)。這個方案適合快速出結(jié)果但如果你的研究需要精細(xì)控制還是手動過濾更靈活。3.3 相對豐度轉(zhuǎn)換的標(biāo)準(zhǔn)代碼過濾完成后就可以做相對豐度轉(zhuǎn)換了。下面這段代碼會把整數(shù)值的OTU表轉(zhuǎn)成相對豐度表并同步更新phyloseq對象# 方法一phyloseq自帶transform_sample_counts ps_rel - transform_sample_counts(ps_filtered, function(x) x / sum(x)) # 驗(yàn)證轉(zhuǎn)換結(jié)果每個樣本之和應(yīng)該為1 otu_rel - otu_table(ps_rel) colSums(otu_rel) %% head() # Sample1 Sample2 Sample3 Sample4 Sample5 Sample6 # 1 1 1 1 1 1 # 方法二如果需要百分?jǐn)?shù)格式0-100 ps_rel_100 - transform_sample_counts(ps_filtered, function(x) x / sum(x) * 100)如果你用的是microbiome包還有一個專門函數(shù)ps_rel - microbiome::transform(ps_filtered, compositional)transform(ps, compositional)和transform_sample_counts的差異不大選擇哪個順手即可。注意在相對豐度轉(zhuǎn)換前一定要先完成過濾否則那些低豐度噪聲也會被“稀釋”到每個樣本里只是比例變得很小并不會消失。3.4 要不要做rarefaction稀疏化這是整個流程里最讓大家糾結(jié)的問題之一。rarefaction稀疏化就是把所有樣本統(tǒng)一抽平到相同測序深度比如全部抽到10000條序列多的扔掉少的不足則剔除。以前這是標(biāo)準(zhǔn)操作但現(xiàn)在學(xué)術(shù)界爭議很大。反對rarefaction的理由很充分隨機(jī)抽平會丟棄有效數(shù)據(jù)導(dǎo)致統(tǒng)計(jì)功效下降同時(shí)抽樣過程引入新的隨機(jī)性使結(jié)果不穩(wěn)定不同隨機(jī)種子可能給出不同結(jié)論。而且現(xiàn)代標(biāo)準(zhǔn)化方法如DESeq2的幾何均值比率法、edgeR的TMM、以及CLR變換都能更合理地處理測序深度差異。但我個人觀點(diǎn)是rarefaction并沒有完全過時(shí)關(guān)鍵看應(yīng)用場景。如果你的后續(xù)分析手段非常依賴整數(shù)值計(jì)數(shù)比如用vegan做某些多樣性指數(shù)、用機(jī)器學(xué)習(xí)方法對原始計(jì)數(shù)敏感那適當(dāng)rarefaction還能讓結(jié)果更直觀、更容易和舊文獻(xiàn)對比。如果是做差異豐度分析我更推薦用專門的工具處理而不是提前rarefaction。如果決定做rarefaction樣本最低深度要有數(shù)不然會丟掉太多樣本。我一般在過濾后先看深度分布用sample_sums()排序確認(rèn)最小深度和十分位數(shù)把明顯偏離大部隊(duì)的低深度樣本先剔除再對剩余樣本做rarefaction。4. 標(biāo)準(zhǔn)化方案進(jìn)階怎么讓樣本真正可比4.1 相對豐度之外的幾種標(biāo)準(zhǔn)化思路標(biāo)準(zhǔn)化的目的是消除測序深度差異帶來的技術(shù)干擾讓樣本之間可比。相對豐度轉(zhuǎn)換是最樸素的方式但并不是唯一也不是永遠(yuǎn)最優(yōu)。常見的標(biāo)準(zhǔn)化方法大致有幾類總豐度縮放TSS也就是相對豐度、累積和縮放CSS、TMM、DESeq2的幾何均值比率、以及CLR變換。每種方法背后都有不同的統(tǒng)計(jì)假設(shè)適用場景也有差異?？傌S度縮放就是各樣本除以自己的總數(shù)簡單直觀適合粗略比較群落組成但受優(yōu)勢物種影響大。CSS對高豐度物種不敏感它在計(jì)算縮放因子時(shí)會忽略極高豐度分類群的強(qiáng)烈影響適合宏基因組或高動態(tài)范圍的擴(kuò)增子數(shù)據(jù)。TMM方法常用于轉(zhuǎn)錄組它基于加權(quán)截?cái)嗑祵Σ町愇锓N比例較穩(wěn)健在擴(kuò)增子數(shù)據(jù)中使用需要謹(jǐn)慎。CLR則把每個特征除以其在樣本中的幾何均值再取對數(shù)能有效打破閉合效應(yīng)適合后續(xù)做相關(guān)網(wǎng)絡(luò)、PCA這類依賴于歐氏距離或相關(guān)結(jié)構(gòu)的分析。選擇標(biāo)準(zhǔn)化方法的核心原則是先明確你要做什么。如果只是看一眼組成比例TSS足夠如果要找差異物種DESeq2或edgeR更靠譜如果要算樣本間距離、做排序或聚類CLR更加合適。沒有“萬能方法”只有“適合場景的方法”。4.2 CLR變換實(shí)操CLR變換在微生物組分析中越來越常用它的數(shù)學(xué)定義是對每個樣本把每個特征的豐度除以該樣本所有特征的幾何均值然后取自然對數(shù)。用R代碼實(shí)現(xiàn)可以使用compositions包、microbiome包自帶函數(shù)library(microbiome) # 先做相對豐度轉(zhuǎn)換再執(zhí)行CLR ps_clr - microbiome::transform(ps_filtered, clr) # 檢驗(yàn)每個樣本的CLR值之和應(yīng)該接近0因?yàn)閷?shù)比變換的中心化性質(zhì) clr_otu - otu_table(ps_clr) rowSums(clr_otu) %% head()注意CLR變換要求所有值都大于0因?yàn)閷?shù)不能處理0。上游過濾能把一部分零值特點(diǎn)去掉但特征表中依然可能存在大量的0。常用的處理方式有兩種給所有0值加一個小的偽計(jì)數(shù)比如1或者在計(jì)算幾何均值時(shí)只使用非零值。這兩種方法各有缺陷前者會引入偏差后者會改變樣本間可比性。實(shí)際操作中多數(shù)人會選擇用“乘法替換”方法在矩陣中處理零值也就是把零值替換為一個小值同時(shí)調(diào)整其他非零值以保持總數(shù)不變。比如zCompositions包的cmultRepl()函數(shù)就是干這個的。下面是我常用的組合流程library(zCompositions) # 將OTU表轉(zhuǎn)為矩陣并執(zhí)行乘法零值替換 otu_mat - as.matrix(otu_table(ps_filtered)) otu_repl - cmultRepl(otu_mat, method CZM, output p-counts) # 將替換后的表更新回phyloseq對象 otu_table(ps_filtered) - otu_table(otu_repl, taxa_are_rows TRUE) # 然后再執(zhí)行CLR ps_clr - microbiome::transform(ps_filtered, clr)這里有個經(jīng)驗(yàn)點(diǎn)“CZM”方法比簡單的加一偽計(jì)數(shù)更穩(wěn)健它會基于多項(xiàng)分布生成替換值并重復(fù)多次取平均減少單次隨機(jī)性。但也要注意它比較耗時(shí)如果數(shù)據(jù)很大可以先過濾掉極端低豐度特征減少計(jì)算量。4.3 標(biāo)準(zhǔn)化之后的質(zhì)量核查清單轉(zhuǎn)換完成不等于流程結(jié)束我強(qiáng)烈建議在進(jìn)入下游分析前做一輪質(zhì)量核查。我的習(xí)慣列表如下樣本總豐度是否都為1或100如果不是說明轉(zhuǎn)換表/存在非整數(shù)計(jì)數(shù)干擾。聚類/排序前是否重新計(jì)算了距離矩陣記得不要再使用原始計(jì)數(shù)距離。是否檢查了特征表中是否存在替換后的負(fù)值或NaN有時(shí)偽計(jì)數(shù)處理不當(dāng)會產(chǎn)生異常值。樣本標(biāo)簽是否在過濾后仍然正確樣本丟失是常見意外?？梢允褂煤唵蔚目梢暬瘉頇z查批量效應(yīng)和離群樣本library(vegan) # 基于CLR值計(jì)算歐氏距離并做PCA排序 dist_clr - vegdist(otu_table(ps_clr), method euclidean) pcoa - cmdscale(dist_clr, k 2) plot(pcoa)如果樣本點(diǎn)在PCA圖上按照批次聚成清晰的小簇說明批次效應(yīng)明顯需要進(jìn)一步分析是技術(shù)批次還是生物學(xué)分組。5. 常見問題與排查技巧實(shí)錄5.1 過濾后樣本數(shù)量變少或特征數(shù)量變成0這種現(xiàn)象通常有兩種原因。第一種是過濾條件過嚴(yán)比如把流行度閾值設(shè)成50%以上而你的樣本本身來自多個差異巨大的環(huán)境組大量特征只在少數(shù)樣本中出沒過濾后特征數(shù)驟降。第二種是樣本深度過低在過濾低深度樣本后剩余樣本數(shù)不足。排查思路很簡單逐步放寬和收緊參數(shù)做敏感性分析。比如把流行度閾值從20%降到10%、5%看特征數(shù)量變化曲線。如果特征數(shù)在某個閾值附近急劇下降說明數(shù)據(jù)本身就稀疏需要設(shè)定更合理閾值如果特征數(shù)始終很低可能需要回到上游檢查聚類/去噪?yún)?shù)。5.2 相對豐度轉(zhuǎn)換后同一特征在不同樣本中差異巨大這有時(shí)是真實(shí)生物學(xué)差異有時(shí)是測序深度和樣本DNA提取效率導(dǎo)致的偽信號。如果同一批處理的技術(shù)重復(fù)樣本之間相對豐度波動非常大那就要警惕是不是標(biāo)準(zhǔn)化方式選了不合適的。一個常見坑是樣本總讀數(shù)差異特別大比如A樣本總讀數(shù)5萬B樣本總讀數(shù)20萬經(jīng)過TSS轉(zhuǎn)換后A樣本中的稀有特征可能被放大反而造成“假高豐度”。此時(shí)建議考慮不依賴總讀數(shù)的標(biāo)準(zhǔn)化方法例如TMM或CSS。還有一個容易忽略的問題如果你在過濾前用了不同批次的數(shù)據(jù)那么批次間測序深度差異會在標(biāo)準(zhǔn)化后依然殘留。建議在準(zhǔn)備數(shù)據(jù)時(shí)先按批次和樣本類型做好分層或使用R包如sva的ComBat_seq做批次校正再去標(biāo)準(zhǔn)化。5.3 用相對豐度直接做相關(guān)分析得出的“共現(xiàn)網(wǎng)絡(luò)”可能是假的我在實(shí)際項(xiàng)目中踩過最大的坑就是直接用相對豐度計(jì)算Spearman相關(guān)然后構(gòu)建物種共現(xiàn)網(wǎng)絡(luò)。結(jié)果出現(xiàn)了大量正相關(guān)初看很興奮結(jié)果發(fā)現(xiàn)很多相關(guān)性是閉合效應(yīng)逼出來的。比如樣本中物種A占比從40%降到20%那么其他所有物種的總占比自動升高了哪怕它們之間毫無生物學(xué)關(guān)聯(lián)也會產(chǎn)生負(fù)相關(guān)或正相關(guān)。要正確處理應(yīng)該用CLR變換后的數(shù)值做相關(guān)分析至少也要做成分?jǐn)?shù)據(jù)感知的相關(guān)方法比如SparCC或Propr。另外提醒一句CLR或ILR變換后相關(guān)性解釋要更謹(jǐn)慎因?yàn)槠涑叨仁菍?shù)空間數(shù)值大小不能直接按原始豐度理解。5.4 匯總成問題速查表常見問題可能原因推薦處理方式特征數(shù)過低過濾閾值過高、數(shù)據(jù)稀疏逐步降低閾值做敏感性分析樣本總豐度和不為1使用了整數(shù)值表未轉(zhuǎn)換重新執(zhí)行TSS并檢查行列方向CLR后出現(xiàn)NaN存在零值使用乘法替換或偽計(jì)數(shù)處理零值樣本聚類按批次分開批次效應(yīng)明顯分層抽樣、批次校正后再分析相對豐度相關(guān)網(wǎng)絡(luò)假陽性多閉合效應(yīng)使用CLR或SparCC等方法稀有物種對結(jié)果影響過大未過濾低豐度特征先過濾再標(biāo)準(zhǔn)化6. 實(shí)操心得與一個容易被忽略的提升方向最后聊一些不常出現(xiàn)在教程里、但直接影響體驗(yàn)和結(jié)果的事情。第一點(diǎn)是記錄所有參數(shù)。過濾閾值、深度處理規(guī)則、標(biāo)準(zhǔn)化方法、隨機(jī)種子全都要記錄下來。很多項(xiàng)目做到后期需要復(fù)現(xiàn)結(jié)果、回應(yīng)對審稿人的質(zhì)疑如果沒有詳細(xì)參數(shù)記錄返工成本極高。建議在分析項(xiàng)目里建一個parameters.yaml或者簡單的備忘文本把每一步的代碼版本、參數(shù)、日期寫清楚。第二點(diǎn)是所有轉(zhuǎn)換步驟都盡量保留原始數(shù)據(jù)不可變。我習(xí)慣把上游生成的原始特征表單獨(dú)存檔任何過濾和轉(zhuǎn)換都只生成新對象絕不原地覆蓋。這樣如果后續(xù)需要調(diào)整閾值或者換一種標(biāo)準(zhǔn)化方法隨時(shí)可以回到干凈起點(diǎn)。第三點(diǎn)值得一提的高級技巧可以把過濾和標(biāo)準(zhǔn)化封裝成一個函數(shù)輸入原始phyloseq對象和參數(shù)列表輸出預(yù)處理后的對象和摘要報(bào)告。這樣在多組數(shù)據(jù)、多輪分析中保持一致減少人為誤差。下面是我的簡化模板preprocess_ps - function(ps, min_total_abundance 10, min_prevalence_samples 5, norm_method compositional) { # 過濾零豐度 ps - prune_taxa(taxa_sums(ps) 0, ps) # 總豐度過濾 ps - prune_taxa(taxa_sums(ps) min_total_abundance, ps) # 流行度過濾 prev_count - rowSums(otu_table(ps) 0) ps - prune_taxa(prev_count min_prevalence_samples, ps) # 標(biāo)準(zhǔn)化 if (norm_method compositional) { ps - transform_sample_counts(ps, function(x) x / sum(x)) } else if (norm_method clr) { # 這里省略零值處理細(xì)節(jié)實(shí)際需要先替換零值 ps - microbiome::transform(ps, clr) } return(ps) }在實(shí)際操作中我吃過不少虧比如發(fā)現(xiàn)共享同一批數(shù)據(jù)的不同文章因預(yù)處理參數(shù)不同導(dǎo)致重疊物種數(shù)量差異巨大也見過因?yàn)闃?biāo)準(zhǔn)化方式?jīng)]寫清楚審稿人要求重跑所有分析的窘境。所以在這篇實(shí)戰(zhàn)記錄里我盡可能把每一步的細(xì)節(jié)和考慮都攤開來寫。預(yù)處理是微生物組分析中最不顯眼但最值得花心思的環(huán)節(jié)把它做扎實(shí)下游分析才會真正站得住腳。