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