1. 從一堆陌生基因序列說起COG注釋到底在解決什么問題做過微生物基因組或者宏基因組項(xiàng)目的人大概率都經(jīng)歷過這樣一個(gè)場(chǎng)景測(cè)序公司交付了一堆結(jié)果里面有個(gè)文件叫“基因預(yù)測(cè)結(jié)果”打開一看幾萬個(gè)基因ID排得整整齊齊每個(gè)ID后面跟著一段A、T、C、G組成的序列。你盯著這些序列心里只有一個(gè)念頭——這些基因到底是干什么的COG注釋分析就是回答這個(gè)問題的核心手段之一。COG全稱是Clusters of Orthologous Groups中文一般叫“直系同源基因簇”。它的核心邏輯并不復(fù)雜把來自不同物種的、推測(cè)具有共同祖先的基因歸為一類每一類就是一個(gè)COG每個(gè)COG代表一個(gè)蛋白家族對(duì)應(yīng)某種特定的生物學(xué)功能。你把自己的基因序列拿去和這個(gè)數(shù)據(jù)庫比對(duì)就能知道你的基因大概屬于哪個(gè)功能家族進(jìn)而推斷它的功能。這件事的價(jià)值在于它把“未知序列”變成了“已知功能類別”。比如你做了一個(gè)土壤微生物組的宏基因組項(xiàng)目注釋完之后發(fā)現(xiàn)大量基因集中在“氨基酸轉(zhuǎn)運(yùn)與代謝”“能量產(chǎn)生與轉(zhuǎn)換”“轉(zhuǎn)錄調(diào)控”這幾個(gè)COG類別里那你就能初步判斷這個(gè)土壤微生物群落的功能偏好。再比如你做的是一個(gè)極端環(huán)境樣本注釋結(jié)果里“防御機(jī)制”“翻譯后修飾”類COG顯著富集那可能暗示這個(gè)環(huán)境對(duì)微生物施加了某種選擇壓力。適合看這篇內(nèi)容的人我大致分三類第一類是剛接觸生物信息學(xué)的研究生手里有數(shù)據(jù)但不知道怎么下手做注釋和畫圖第二類是做微生物組或者基因組項(xiàng)目的從業(yè)者之前可能用過在線工具一鍵出結(jié)果但想搞清楚背后的邏輯和參數(shù)細(xì)節(jié)第三類是需要把COG注釋結(jié)果整理成論文圖表或者項(xiàng)目報(bào)告的人想知道怎么把注釋結(jié)果做得既專業(yè)又好看。我自己的經(jīng)驗(yàn)是COG注釋這件事工具用對(duì)只是第一步真正拉開差距的是對(duì)注釋結(jié)果的理解和可視化呈現(xiàn)。同樣一份數(shù)據(jù)有人做出來就是一張干巴巴的餅圖有人能做出層次分明、信息密度高的組合圖審稿人和項(xiàng)目評(píng)審的觀感完全不一樣。接下來我會(huì)從方案設(shè)計(jì)、數(shù)據(jù)庫選擇、實(shí)操流程、可視化技巧到問題排查把整個(gè)鏈路拆開講清楚。2. 方案選型與整體設(shè)計(jì)為什么COG仍然是功能注釋的常青樹2.1 COG、KEGG、GO三者的定位差異剛?cè)腴T的人最容易困惑的一點(diǎn)是功能注釋的數(shù)據(jù)庫那么多COG、KEGG、GO、Pfam、Swiss-Prot到底該用哪個(gè)我的建議是不要把它們當(dāng)成互斥選項(xiàng)而是理解各自擅長(zhǎng)的層面。COG的強(qiáng)項(xiàng)在于功能分類的粗粒度歸納。它把基因按功能大類歸攏比如“碳水化合物轉(zhuǎn)運(yùn)與代謝”“細(xì)胞壁/膜/包被生物合成”“復(fù)制、重組與修復(fù)”等一共二十多個(gè)大類。這種粗粒度看起來不夠精細(xì)但恰恰適合做宏觀功能概覽。你拿到一個(gè)陌生基因組先跑一遍COG能快速知道這個(gè)物種或者群落在功能層面的整體傾向。KEGG的強(qiáng)項(xiàng)在于代謝通路和信號(hào)通路的映射。它能把基因定位到具體的通路圖上比如糖酵解、TCA循環(huán)、雙組分系統(tǒng)等。如果你關(guān)心的是“這個(gè)群落能不能降解某種污染物”“有沒有完整的固氮通路”KEGG更合適。GO的強(qiáng)項(xiàng)在于功能描述的標(biāo)準(zhǔn)化。它從分子功能、生物學(xué)過程、細(xì)胞組分三個(gè)維度給基因打標(biāo)簽適合做富集分析和跨物種比較。實(shí)際操作中我通常的做法是COG做宏觀概覽和分類統(tǒng)計(jì)KEGG做通路級(jí)深入分析GO做富集驗(yàn)證。三者互相補(bǔ)充而不是只用一個(gè)。COG之所以在這些年里一直是功能注釋的常青樹核心原因是它的數(shù)據(jù)庫結(jié)構(gòu)清晰、注釋速度快、結(jié)果解讀門檻低特別適合作為功能分析的第一站。2.2 數(shù)據(jù)庫版本選擇NOG、COG、eggNOG的區(qū)別這里有一個(gè)很多人踩過的坑網(wǎng)上搜“COG數(shù)據(jù)庫下載”出來的結(jié)果五花八門有人下的是NCBI的老版COG有人下的是eggNOG的NOG還有人下的是KOG。這幾個(gè)東西名字像但適用范圍差別很大。NCBI的原始COG數(shù)據(jù)庫主要基于細(xì)菌、古菌和真核生物的完整基因組構(gòu)建更新頻率較低但勝在經(jīng)典、穩(wěn)定很多老項(xiàng)目的注釋結(jié)果都是基于它。eggNOG則是升級(jí)版覆蓋的物種范圍更廣除了COG之外還有KOG真核生物直系同源組和NOG更寬泛的同源組。eggNOG的更新頻率高注釋信息也更豐富目前主流做法是用eggNOG的COG子集來做注釋。我的建議是如果你做的是細(xì)菌或古菌項(xiàng)目用eggNOG的COG子集就夠了覆蓋度和更新性都更好。如果你做的是真核生物比如真菌或者小型真核生物可以考慮KOG。如果你需要和之前發(fā)表的老項(xiàng)目做對(duì)比那就用NCBI原始COG保持一致。數(shù)據(jù)庫版本這件事沒有絕對(duì)的對(duì)錯(cuò)關(guān)鍵是要和你項(xiàng)目的分析目標(biāo)以及對(duì)比對(duì)象匹配。2.3 注釋工具選型本地比對(duì)還是在線服務(wù)工具層面常見的選擇有幾種。第一種是NCBI的CD-Search在線服務(wù)適合少量序列快速驗(yàn)證但批量處理幾萬個(gè)基因就不現(xiàn)實(shí)了。第二種是本地跑BLAST或者DIAMOND比對(duì)這是最主流的方式靈活性和可控性最高。第三種是用eggNOG-mapper這類封裝好的工具一條命令搞定比對(duì)和注釋適合不想折騰參數(shù)的人。我自己最常用的是DIAMOND加本地COG數(shù)據(jù)庫的組合。原因很簡(jiǎn)單DIAMOND的比對(duì)速度比傳統(tǒng)BLAST快幾個(gè)數(shù)量級(jí)而且對(duì)短序列的敏感度也夠用。具體參數(shù)上我一般用--evalue 1e-5作為閾值--max-target-seqs 1只保留最佳比對(duì)結(jié)果--outfmt 6輸出制表符分隔格式方便后續(xù)處理。如果你的基因數(shù)量在幾萬條以內(nèi)用BLAST也不是不行但時(shí)間成本會(huì)高很多。在線服務(wù)的好處是省事但有兩個(gè)風(fēng)險(xiǎn)一是數(shù)據(jù)隱私二是服務(wù)穩(wěn)定性。如果你的數(shù)據(jù)涉及未發(fā)表成果我強(qiáng)烈建議本地跑。本地跑還有一個(gè)好處是你可以隨時(shí)調(diào)整參數(shù)重新注釋不用等在線隊(duì)列。3. 核心細(xì)節(jié)解析從序列到功能類別的完整鏈路3.1 輸入數(shù)據(jù)的準(zhǔn)備與格式要求COG注釋的輸入通常是蛋白序列而不是核酸序列。這一點(diǎn)很多人一開始會(huì)搞混。原因在于COG數(shù)據(jù)庫本身是蛋白層面的同源組比對(duì)也是在蛋白層面進(jìn)行的。所以你需要先把基因預(yù)測(cè)得到的核酸序列翻譯成蛋白序列。翻譯這一步常用的工具是Prodigal或者GeneMark。Prodigal的優(yōu)勢(shì)是速度快、對(duì)原核生物基因預(yù)測(cè)準(zhǔn)確率高而且可以直接輸出蛋白序列。我一般用Prodigal的-p meta模式處理宏基因組數(shù)據(jù)用默認(rèn)模式處理單基因組數(shù)據(jù)。翻譯完成后檢查一下序列文件是否符合FASTA格式每條序列以開頭后面跟序列ID和描述信息然后是序列本身。序列ID最好保持和原始基因ID一致方便后續(xù)回溯。還有一個(gè)細(xì)節(jié)容易被忽略如果蛋白序列里存在終止密碼子或者內(nèi)部終止子比對(duì)結(jié)果會(huì)受影響。Prodigal輸出的蛋白序列一般已經(jīng)處理過這個(gè)問題但如果你是自己用其他工具翻譯的建議用seqkit或者biopython檢查一下把含有內(nèi)部終止子的序列過濾掉。3.2 比對(duì)參數(shù)的選擇邏輯比對(duì)參數(shù)直接決定了注釋結(jié)果的靈敏度和特異性。我拿DIAMOND舉例幾個(gè)關(guān)鍵參數(shù)的選擇邏輯如下。--evalue控制的是比對(duì)結(jié)果的統(tǒng)計(jì)顯著性。值越小結(jié)果越嚴(yán)格假陽性越少但可能漏掉一些真實(shí)同源但序列差異較大的基因。我通常用1e-5作為起點(diǎn)如果注釋率偏低可以放寬到1e-3試試。但要注意放寬閾值會(huì)引入更多低置信度的注釋后續(xù)分析時(shí)要留意。--max-target-seqs控制每條序列保留多少個(gè)比對(duì)結(jié)果。設(shè)為1表示只保留最佳比對(duì)適合做COG分類統(tǒng)計(jì)。如果你想看一個(gè)基因可能屬于多個(gè)COG的情況可以設(shè)為5或者10但后續(xù)統(tǒng)計(jì)時(shí)要決定怎么處理多映射。--query-cover和--subject-cover控制比對(duì)覆蓋度。我一般要求覆蓋度不低于50%否則即使E值顯著也可能只是局部同源不能代表整個(gè)基因的功能。--id控制序列一致性。對(duì)于跨物種的COG注釋一致性閾值不宜設(shè)得太高30%左右是比較常用的起點(diǎn)。設(shè)太高會(huì)導(dǎo)致注釋率大幅下降設(shè)太低會(huì)引入噪聲。這些參數(shù)沒有一套放之四海而皆準(zhǔn)的數(shù)值核心原則是先跑一版默認(rèn)參數(shù)看注釋率和結(jié)果分布再根據(jù)項(xiàng)目需求微調(diào)。我習(xí)慣在項(xiàng)目記錄里把每版參數(shù)和對(duì)應(yīng)的注釋率都記下來方便回溯和對(duì)比。3.3 COG功能大類的映射與統(tǒng)計(jì)比對(duì)完成后你得到的是每條基因?qū)?yīng)的COG ID。但COG ID本身只是一串編號(hào)比如COG0001、COG0002直接看沒有意義。你需要把它映射到功能大類上。COG數(shù)據(jù)庫提供了一個(gè)功能分類表把每個(gè)COG ID歸入一個(gè)大類用單個(gè)字母表示。比如J代表翻譯、核糖體結(jié)構(gòu)與生物合成K代表轉(zhuǎn)錄L代表復(fù)制、重組與修復(fù)D代表細(xì)胞周期控制、細(xì)胞分裂、染色體分割等等。一共二十多個(gè)大類每個(gè)大類下面還有更細(xì)的功能描述。統(tǒng)計(jì)這一步我通常做兩個(gè)層面的匯總。第一個(gè)層面是大類層面的計(jì)數(shù)每個(gè)功能大類里有多少條基因占總基因數(shù)的百分比是多少。這個(gè)結(jié)果適合做餅圖或者柱狀圖給人一個(gè)宏觀印象。第二個(gè)層面是具體COG層面的計(jì)數(shù)每個(gè)COG ID對(duì)應(yīng)多少條基因按數(shù)量排序取前20或者前30做展示。這個(gè)結(jié)果適合做條形圖能看出哪些具體功能家族在樣本中富集。這里有一個(gè)實(shí)操心得大類層面的統(tǒng)計(jì)建議同時(shí)輸出絕對(duì)數(shù)量和百分比。因?yàn)椴煌瑯颖镜幕蚩倲?shù)可能差異很大只看百分比會(huì)丟失規(guī)模信息只看絕對(duì)數(shù)量又不好跨樣本比較。兩個(gè)都給讀者自己判斷。4. 實(shí)操過程從原始序列到可視化圖表的完整復(fù)現(xiàn)4.1 環(huán)境準(zhǔn)備與數(shù)據(jù)庫下載我假設(shè)你用的是Linux環(huán)境這是生物信息分析的主流平臺(tái)。先建一個(gè)工作目錄把原始數(shù)據(jù)、數(shù)據(jù)庫、中間文件、結(jié)果文件分開放避免文件混亂。數(shù)據(jù)庫下載這一步eggNOG的官網(wǎng)提供了預(yù)構(gòu)建的DIAMOND數(shù)據(jù)庫文件直接下載解壓就能用省去了自己建庫的時(shí)間。下載完成后用diamond makedb命令把蛋白序列文件轉(zhuǎn)成DIAMOND格式的數(shù)據(jù)庫。這一步只需要做一次后續(xù)所有項(xiàng)目都可以復(fù)用。工具安裝方面DIAMOND可以用conda直接裝Prodigal也是。如果你不想折騰環(huán)境用conda創(chuàng)建一個(gè)獨(dú)立環(huán)境是最省事的做法。我一般會(huì)固定工具版本比如DIAMOND 2.1.x和Prodigal 2.6.x避免不同版本之間參數(shù)行為差異導(dǎo)致結(jié)果不一致。4.2 基因預(yù)測(cè)與蛋白序列提取假設(shè)你拿到的是組裝好的基因組或者宏基因組contig文件。第一步是用Prodigal做基因預(yù)測(cè)prodigal -i assembly.fasta -a proteins.faa -d genes.fna -o genes.gbk -p meta-a輸出蛋白序列-d輸出核酸序列-o輸出完整的基因預(yù)測(cè)報(bào)告。-p meta表示宏基因組模式如果是單基因組就把這個(gè)參數(shù)去掉。跑完之后檢查一下proteins.faa文件里有多少條序列。如果序列數(shù)量和你預(yù)期的基因數(shù)量差距很大可能是組裝質(zhì)量或者預(yù)測(cè)參數(shù)的問題。我遇到過幾次因?yàn)閏ontig太短導(dǎo)致Prodigal預(yù)測(cè)不出基因的情況后來把最短contig長(zhǎng)度閾值調(diào)到500bp以上就正常了。4.3 DIAMOND比對(duì)與結(jié)果過濾比對(duì)命令如下diamond blastp -d eggnog_cog.dmnd -q proteins.faa -o blast_results.tsv --evalue 1e-5 --max-target-seqs 1 --outfmt 6 --query-cover 50 --subject-cover 50 --threads 8--threads根據(jù)你的機(jī)器配置調(diào)整一般設(shè)成CPU核心數(shù)的80%左右比較穩(wěn)妥留一些資源給系統(tǒng)。比對(duì)完成后blast_results.tsv里每行是一條基因的比對(duì)結(jié)果包含基因ID、COG ID、一致性、覆蓋度、E值等信息。接下來需要把COG ID映射到功能大類。eggNOG提供了一個(gè)cog_category的映射文件格式是兩列COG ID和功能大類字母。用join或者awk做映射就行。我一般會(huì)寫一個(gè)簡(jiǎn)單的Python腳本來處理這一步因?yàn)橐瑫r(shí)做幾件事過濾低質(zhì)量比對(duì)、映射功能大類、統(tǒng)計(jì)每個(gè)大類的基因數(shù)量、輸出多個(gè)格式的結(jié)果文件。腳本邏輯不復(fù)雜但手寫一遍比每次用命令行拼湊更可靠。4.4 可視化圖表的制作要點(diǎn)COG注釋結(jié)果的可視化常見的圖表類型有幾種。餅圖適合展示大類層面的占比但缺點(diǎn)是當(dāng)類別超過8個(gè)時(shí)小扇區(qū)會(huì)擠在一起看不清。柱狀圖適合展示具體COG的豐度排序橫軸是COG ID或者功能描述縱軸是基因數(shù)量。堆疊柱狀圖適合做多樣本比較每個(gè)樣本一根柱子不同顏色代表不同功能大類。我個(gè)人的偏好是大類層面用橫向柱狀圖而不是餅圖因?yàn)闄M向柱狀圖的標(biāo)簽更容易閱讀排序也更直觀。具體COG層面用條形圖取Top 20或者Top 30其余歸為“其他”。如果是多樣本比較用堆疊柱狀圖但顏色不要超過8種否則辨識(shí)度會(huì)下降。配色方面我建議用色盲友好的調(diào)色板比如ColorBrewer的Set2或者Paired。避免用紅綠對(duì)比因?yàn)橛幸徊糠秩舜嬖诩t綠色覺障礙。圖表標(biāo)題和坐標(biāo)軸標(biāo)簽要寫清楚單位要標(biāo)明。如果圖是給論文用的字體大小和分辨率要符合期刊要求一般300dpi起步。還有一個(gè)細(xì)節(jié)COG功能大類的名稱通常比較長(zhǎng)比如“翻譯后修飾、蛋白質(zhì)周轉(zhuǎn)、伴侶蛋白”直接放在坐標(biāo)軸上會(huì)占很多空間。我一般會(huì)縮寫或者用字母代號(hào)然后在圖注里給出完整名稱。這樣圖面干凈信息也不丟失。5. 常見問題與排查技巧實(shí)錄5.1 注釋率偏低怎么辦注釋率偏低是新手最常遇到的問題。所謂注釋率就是成功比對(duì)到COG數(shù)據(jù)庫的基因數(shù)占總基因數(shù)的比例。一般來說細(xì)菌基因組的注釋率在70%到85%之間算正常宏基因組的注釋率可能低一些50%到70%也常見。如果你的注釋率明顯低于這個(gè)范圍可以從幾個(gè)方向排查。第一檢查輸入序列是不是蛋白序列。如果誤把核酸序列當(dāng)?shù)鞍仔蛄腥ケ葘?duì)結(jié)果會(huì)慘不忍睹。第二檢查數(shù)據(jù)庫是否完整下載和解壓。有時(shí)候下載中斷導(dǎo)致數(shù)據(jù)庫文件不完整比對(duì)結(jié)果會(huì)異常。第三嘗試放寬E值閾值和覆蓋度閾值看注釋率是否明顯提升。如果放寬后提升很大說明你的序列和數(shù)據(jù)庫的差異較大可能需要考慮用更寬泛的NOG數(shù)據(jù)庫。第四檢查基因預(yù)測(cè)是否合理。如果預(yù)測(cè)出的蛋白序列普遍偏短可能是基因預(yù)測(cè)參數(shù)不合適。5.2 多映射與結(jié)果沖突的處理有些基因會(huì)比對(duì)到多個(gè)COG上而且這些COG可能屬于不同的功能大類。這種情況在宏基因組數(shù)據(jù)里尤其常見因?yàn)楹昊蚪M里混雜了來自不同物種的序列同源關(guān)系更復(fù)雜。處理多映射常見策略有三種。第一種是只保留最佳比對(duì)也就是E值最小、一致性最高的那個(gè)。這是最簡(jiǎn)單也最常用的做法適合做宏觀統(tǒng)計(jì)。第二種是保留所有比對(duì)但在統(tǒng)計(jì)時(shí)按權(quán)重分配比如一個(gè)基因比對(duì)到三個(gè)COG每個(gè)COG計(jì)0.33。這種做法更精細(xì)但解釋起來復(fù)雜。第三種是只保留一致性超過某個(gè)閾值的比對(duì)低于閾值的丟棄。我一般用第一種策略做常規(guī)分析用第二種策略做深入分析。關(guān)鍵是要在方法部分寫清楚你用了哪種策略因?yàn)椴煌呗詴?huì)導(dǎo)致結(jié)果差異。5.3 圖表信息密度與可讀性的平衡做可視化的時(shí)候很容易陷入一個(gè)誤區(qū)想把所有信息都塞進(jìn)一張圖里。結(jié)果就是圖面擁擠、標(biāo)簽重疊、顏色混亂讀者根本看不懂。我的經(jīng)驗(yàn)是一張圖只講一件事。大類占比就只講大類占比不要同時(shí)疊加具體COG的細(xì)節(jié)。具體COG的豐度排序就只講排序不要同時(shí)展示多個(gè)樣本的對(duì)比。如果確實(shí)需要展示多個(gè)層面的信息那就拆成多張圖或者用分面圖facet的方式組織。另外圖表的注釋文字要克制。不要在圖上寫大段解釋把解釋放在圖注或者正文里。圖本身要干凈讓讀者一眼能看出主要趨勢(shì)。5.4 常見問題速查表問題現(xiàn)象可能原因排查方向解決建議注釋率低于50%輸入序列格式錯(cuò)誤檢查是否為蛋白序列重新翻譯核酸序列注釋率低于50%數(shù)據(jù)庫不完整檢查數(shù)據(jù)庫文件大小重新下載解壓注釋率低于50%閾值過嚴(yán)放寬E值和覆蓋度逐步調(diào)整參數(shù)比對(duì)結(jié)果為空數(shù)據(jù)庫路徑錯(cuò)誤檢查-d參數(shù)路徑確認(rèn)數(shù)據(jù)庫文件存在多映射嚴(yán)重宏基因組復(fù)雜度高查看比對(duì)結(jié)果分布只保留最佳比對(duì)圖表標(biāo)簽重疊類別過多檢查類別數(shù)量合并小類別為“其他”圖表顏色難辨配色不友好檢查色盲友好性換用Set2或Paired配色6. 結(jié)果解讀與后續(xù)分析方向6.1 從功能大類分布看樣本特征COG注釋結(jié)果出來之后怎么解讀是一門功夫。大類層面的分布能給你很多線索。比如“氨基酸轉(zhuǎn)運(yùn)與代謝”類占比高說明樣本中蛋白質(zhì)合成和降解活動(dòng)活躍?！澳芰慨a(chǎn)生與轉(zhuǎn)換”類占比高說明樣本的代謝活性強(qiáng)?！胺烙鶛C(jī)制”類占比高可能暗示環(huán)境中存在選擇壓力?！耙苿?dòng)基因組”類占比高比如轉(zhuǎn)座子、質(zhì)粒相關(guān)基因可能說明樣本中存在水平基因轉(zhuǎn)移。但要注意這些解讀都是概率性的不是絕對(duì)的。一個(gè)功能大類占比高可能是因?yàn)闃颖局写_實(shí)有大量相關(guān)基因也可能是因?yàn)閿?shù)據(jù)庫對(duì)這個(gè)大類的注釋覆蓋度更高。解讀時(shí)要結(jié)合樣本背景和其他分析結(jié)果不要單憑COG分布下結(jié)論。6.2 與KEGG、GO結(jié)果的交叉驗(yàn)證COG注釋的結(jié)果最好和KEGG、GO的結(jié)果交叉驗(yàn)證。如果COG顯示“碳水化合物代謝”類富集KEGG也顯示糖酵解和TCA循環(huán)通路完整那這個(gè)結(jié)論就比較可靠。如果兩者矛盾就需要深入排查原因可能是注釋閾值不同也可能是數(shù)據(jù)庫覆蓋度差異。交叉驗(yàn)證還有一個(gè)好處是能發(fā)現(xiàn)新的線索。比如COG注釋顯示某個(gè)功能大類富集但KEGG通路分析沒有顯著結(jié)果那可能說明這個(gè)大類里的基因還沒有被映射到已知通路上值得進(jìn)一步挖掘。6.3 多組比較與差異功能分析如果你有多個(gè)樣本或者多個(gè)處理組COG注釋結(jié)果可以做差異功能分析。基本思路是先統(tǒng)計(jì)每個(gè)樣本在每個(gè)功能大類上的基因數(shù)量或者百分比然后做組間比較找出顯著差異的功能大類。統(tǒng)計(jì)方法上如果樣本量小可以用簡(jiǎn)單的倍數(shù)變化加卡方檢驗(yàn)。如果樣本量大可以考慮用DESeq2或者edgeR這類專門做差異分析的工具把功能大類當(dāng)成“基因”來處理。不過要注意COG大類層面的計(jì)數(shù)是匯總數(shù)據(jù)直接套用基因?qū)用娴牟町惙治龉ぞ呖赡懿煌耆线m結(jié)果解釋要謹(jǐn)慎。我自己的做法是先做描述性統(tǒng)計(jì)看組間分布差異再用統(tǒng)計(jì)檢驗(yàn)確認(rèn)顯著性最后結(jié)合生物學(xué)知識(shí)判斷哪些差異是真正有意義的。統(tǒng)計(jì)顯著不等于生物學(xué)顯著這一點(diǎn)在功能分析里尤其重要。7. 我踩過的坑和幾條實(shí)用建議第一個(gè)坑是數(shù)據(jù)庫版本混用。有一次我做一個(gè)對(duì)比項(xiàng)目?jī)蓚€(gè)樣本分別用了不同版本的COG數(shù)據(jù)庫注釋結(jié)果功能大類分布差異很大后來發(fā)現(xiàn)是數(shù)據(jù)庫更新導(dǎo)致某些COG的分類變了。從那以后我所有對(duì)比項(xiàng)目都固定用同一個(gè)版本的數(shù)據(jù)庫并且在方法里寫清楚版本號(hào)。第二個(gè)坑是忽略序列ID的對(duì)應(yīng)關(guān)系。Prodigal預(yù)測(cè)基因時(shí)會(huì)自動(dòng)生成ID如果你后續(xù)用其他工具處理過序列ID可能會(huì)變。一旦ID對(duì)應(yīng)不上注釋結(jié)果就沒法回溯到原始基因。我的做法是從基因預(yù)測(cè)開始所有中間文件的ID都保持一致不做重命名。第三個(gè)坑是圖表配色。早期我做堆疊柱狀圖用了默認(rèn)的彩虹配色結(jié)果打印出來是黑白的完全分不清。后來改用灰度加紋理的方案或者用色盲友好的配色問題就解決了。如果你的圖要投稿提前確認(rèn)期刊對(duì)彩色的要求。第四個(gè)坑是注釋結(jié)果的過度解讀。COG注釋給的是功能類別的歸屬不是功能的直接證據(jù)。一個(gè)基因被注釋到“轉(zhuǎn)錄調(diào)控”大類不代表它一定是一個(gè)轉(zhuǎn)錄因子只是說它和已知的轉(zhuǎn)錄調(diào)控相關(guān)基因有同源性。結(jié)論要留有余地不要說得太絕對(duì)。最后分享一個(gè)小技巧如果你要做大量樣本的COG注釋建議把比對(duì)和統(tǒng)計(jì)步驟腳本化用Snakemake或者Nextflow做流程管理。這樣不僅省時(shí)間還能保證每次運(yùn)行的參數(shù)一致減少人為錯(cuò)誤。我早期手動(dòng)跑流程的時(shí)候經(jīng)常因?yàn)閰?shù)記錯(cuò)或者文件路徑寫錯(cuò)導(dǎo)致結(jié)果異常后來改成流程化管理之后這類問題基本消失了。