最近，我開發(fā)的R包ReporterScore（GRSA方法）發(fā)表在Briefings in Bioinformatics上了，這是一種靈活的，可用于復(fù)雜多組學(xué)數(shù)據(jù)的功能富集新方法。

R包主頁(yè)：https://github.com/Asa12138/ReporterScore

GRSA具有比GSEA等常用富集方法更好的sensitivity以獲取更多的testable hypotheses。
GRSA可以直接適用于multi-group and longitudinal 實(shí)驗(yàn)設(shè)計(jì)，而不需要重復(fù)做兩兩比較分析，可以直接對(duì)復(fù)雜實(shí)驗(yàn)設(shè)計(jì)進(jìn)行富集分析。
我們將該GRSA應(yīng)用于轉(zhuǎn)錄組、代謝組、微生物組等數(shù)據(jù)并有了新的發(fā)現(xiàn)。
我們預(yù)計(jì)該富集方法可以替代常用富集分析方法用在各類組學(xué)數(shù)據(jù)上。
這個(gè)R包內(nèi)置了各類主流可視化方法以及主要功能學(xué)數(shù)據(jù)庫(kù)，屬于一站式解決方案。歡迎大家下載使用，提出寶貴意見。

下面對(duì)這篇文章進(jìn)行簡(jiǎn)要地介紹：

Paper Info

• 標(biāo)題：Generalized reporter score-based enrichment analysis for omics data
• 譯名：基于廣義報(bào)告評(píng)分的富集分析可用于各類組學(xué)數(shù)據(jù)
• 期刊：Briefings in Bioinformatics
  IF：9.5
• 發(fā)表時(shí)間：2024年3月27日
• 鏈接：https://doi.org/10.1093/bib/bbae116

引言

功能富集分析是一種常用的生物信息學(xué)方法，用于理解大型組學(xué)數(shù)據(jù)集（如轉(zhuǎn)錄組、宏基因組和代謝組數(shù)據(jù)）的生物學(xué)意義。通過識(shí)別富集的功能類別（如基因本體術(shù)語(yǔ)或生物通路），我們可以洞察底層的生物過程和功能，并提出下游實(shí)驗(yàn)研究的假設(shè)。功能富集分析方法可以根據(jù)統(tǒng)計(jì)方法大致分為三類：（i）過表示分析（ORA）、（ii）功能類別評(píng)分（FCS）和（iii）通路拓?fù)洌≒T）。Goeman和Bühlmann將富集分析方法根據(jù)底層零假設(shè)歸類為“競(jìng)爭(zhēng)性”或“自包含”。在“競(jìng)爭(zhēng)性”方法中，將基因集與不在該集合中的所有基因背景進(jìn)行比較，以評(píng)估統(tǒng)計(jì)差異水平是否超過背景水平；而“自包含”方法則分析每個(gè)基因集的孤立情況。

基于報(bào)告分?jǐn)?shù)的分析（RSA）算法最初由Patil和Nielsen于2005年開發(fā)，用于識(shí)別與代謝網(wǎng)絡(luò)的調(diào)控?zé)狳c(diǎn)相關(guān)的代謝物。RSA近年來(lái)因其在微生物組研究中功能富集分析方面的擴(kuò)展應(yīng)用而備受青睞。RSA是一種基于報(bào)告分?jǐn)?shù)的競(jìng)爭(zhēng)性FCS方法，基于選擇的統(tǒng)計(jì)分析的P值解析，沒有先驗(yàn)截止值（無(wú)閾值）。其基本原理是，P值可以被視為標(biāo)準(zhǔn)化的統(tǒng)計(jì)量，反映了不同基因或特征之間的差異，而不受平均表達(dá)值的影響。與背景P值分布相比，具有顯著較低P值的通路是富集的。

然而，由于缺乏特定工具和對(duì)算法的系統(tǒng)理解，RSA經(jīng)常被誤用。此外，經(jīng)典RSA中每條通路的報(bào)告分?jǐn)?shù)的符號(hào)（正負(fù)號(hào)）并不表示通路表達(dá)趨勢(shì)的增加或減少；而是報(bào)告分?jǐn)?shù)（包括負(fù)值）低于指定閾值表示相應(yīng)的通路未顯著富集。這經(jīng)常導(dǎo)致結(jié)果的錯(cuò)誤解釋。

受經(jīng)典RSA啟發(fā)，我們開發(fā)了改進(jìn)的廣義報(bào)告分?jǐn)?shù)分析（GRSA）方法，并在R包ReporterScore中實(shí)現(xiàn)了該方法，還包括全面的可視化方法和通路數(shù)據(jù)庫(kù)。GRSA是一種無(wú)閾值方法，適用于所有類型的生物醫(yī)學(xué)特征，如基因、化合物和微生物物種。GRSA可在mixed模式（經(jīng)典RSA）和directed模式（增強(qiáng)RSA）下工作。directed模式使用報(bào)告分?jǐn)?shù)的符號(hào)區(qū)分上調(diào)或下調(diào)的通路，更直觀。重要的是，GRSA支持多組和縱向?qū)嶒?yàn)設(shè)計(jì)，因?yàn)樗硕嘟M兼容的統(tǒng)計(jì)方法。ReporterScore包還支持自定義的分層和關(guān)系數(shù)據(jù)庫(kù)，為高級(jí)用戶提供額外的靈活性。在本研究中，我們描述了GRSA的全面實(shí)用性。我們?cè)诙鄠€(gè)數(shù)據(jù)集上將GRSA與其他流行的富集方法進(jìn)行了基準(zhǔn)測(cè)試，并展示了GRSA在各種組學(xué)數(shù)據(jù)集上的應(yīng)用。

方法

圖1:ReporterScore包中GRSA的整體工作流程。

GRSA算法簡(jiǎn)要流程如下（細(xì)節(jié)請(qǐng)參見原文）：

(1) 計(jì)算P值

使用統(tǒng)計(jì)方法（Wilcox, t-test, ANOVA, pearson,...）獲得特征（feature，即,可以是基因，代謝物，KO等等）的P值。這里簡(jiǎn)單地用KO表示不同的特征。

(2) 將P值轉(zhuǎn)換為Z分?jǐn)?shù)

對(duì)于經(jīng)典的mxied RSA，使用逆正態(tài)累積分布函數(shù) () 將KO的P值轉(zhuǎn)換為Z分?jǐn)?shù)。因此，假設(shè)在隨機(jī)數(shù)據(jù)假設(shè)范圍內(nèi)，P值均勻分布，得到的Z分?jǐn)?shù)將遵循標(biāo)準(zhǔn)正態(tài)分布。

對(duì)于新的directed RSA，首先將P值除以2，將P值范圍從轉(zhuǎn)換為。然后使用逆正態(tài)累積分布函數(shù)將KO的P值轉(zhuǎn)換為Z分?jǐn)?shù)。當(dāng)P值為0.5時(shí)，轉(zhuǎn)換后的Z分?jǐn)?shù)等于0。

確定KO是上調(diào)還是下調(diào)，并計(jì)算。

對(duì)于差異豐度分析（兩組設(shè)計(jì)）：

對(duì)于相關(guān)性分析（兩組、多組和縱向設(shè)計(jì)）：

給每個(gè)Z分?jǐn)?shù)分配正負(fù)號(hào)，如果小于0，則Z分?jǐn)?shù)取負(fù)值；否則取正值。

(3) 評(píng)分通路

使用KO的Z分?jǐn)?shù)對(duì)通路進(jìn)行評(píng)分。首先選擇通路數(shù)據(jù)庫(kù)作為參考。對(duì)于選定數(shù)據(jù)庫(kù)中的每個(gè)通路，計(jì)算通路的Z分?jǐn)?shù)。接下來(lái)，對(duì)Z分?jǐn)?shù)進(jìn)行校正，以評(píng)估富集的顯著性。最后，計(jì)算通路的報(bào)告分?jǐn)?shù)和相關(guān)的P值。

我們使用了50個(gè)基準(zhǔn)數(shù)據(jù)集（41個(gè)人類基因表達(dá)數(shù)據(jù)集和9個(gè)小鼠基因表達(dá)數(shù)據(jù)集）比較了GRSA和六種常用的富集分析方法（Fisher’s exact test，GSEA，SAFE，GSA，PADOG，GSVA）。我們使用了不同的統(tǒng)計(jì)方法對(duì)特征進(jìn)行P值計(jì)算，并使用Benjamini和Hochberg方法進(jìn)行了P值調(diào)整。我們?cè)u(píng)估了每種方法在確定目標(biāo)通路方面的性能。

GRSA可應(yīng)用于多組和縱向組學(xué)數(shù)據(jù)

GRSA的一個(gè)重要特征是新開發(fā)的directed模式。directed模式與mixed模式（經(jīng)典RSA）的主要區(qū)別在于，在directed模式中，報(bào)告分?jǐn)?shù)的正負(fù)符號(hào)表示通路的增長(zhǎng)或減少趨勢(shì)。相比之下，在mixed模式中，報(bào)告分?jǐn)?shù)的符號(hào)不表示通路的趨勢(shì)。我們?cè)诠驳膃x_KO_profile數(shù)據(jù)集上以兩種模式（directed模式和mixed模式）進(jìn)行了GRSA。在directed模式中富集的通路中，大多數(shù)KO具有相同的趨勢(shì)。具有一致增加（減少）KO的通路將獲得比背景更大（更?。┑木酆蟌分?jǐn)?shù)。相反，如果通路內(nèi)的KO具有相反的趨勢(shì)，那么帶符號(hào)的Z分?jǐn)?shù)將互相抵消，導(dǎo)致結(jié)果不顯著。與之相比，在mixed模式中，無(wú)法確定富集通路的增長(zhǎng)和減少趨勢(shì)。因此，directed模式有助于發(fā)現(xiàn)具有一致變化KO的通路。一些先前的研究旨在使用directed模式的結(jié)果，但錯(cuò)誤地使用了經(jīng)典RSA（mixed模式）。

GRSA的另一個(gè)主要優(yōu)勢(shì)是支持多組和縱向組學(xué)數(shù)據(jù)。ReporterScore包使用差異豐度分析或相關(guān)性分析計(jì)算組間每個(gè)特征的P值。Kruskal–Wallis檢驗(yàn)或ANOVA評(píng)估特征豐度是否在多個(gè)組之間顯著變化。默認(rèn)的相關(guān)性分析將組分配視為有序的（例如，組‘G1’、‘G2’和‘G3’將轉(zhuǎn)換為1、2和3），因此相關(guān)性分析可以評(píng)估特征豐度是否呈線性增加或減少。此外，ReporterScore包還支持任何指定的模式。例如，如果期望出現(xiàn)指數(shù)增長(zhǎng)趨勢(shì)，可以將組‘G1’、‘G2’和‘G3’設(shè)置為1、10和100。為了探索數(shù)據(jù)中的潛在模式，可以使用聚類方法，例如C均值聚類。

作為一般規(guī)則，用戶必須確保所選的統(tǒng)計(jì)方法適用于數(shù)據(jù)集和實(shí)驗(yàn)設(shè)計(jì)。我們?cè)诙鄠€(gè)基準(zhǔn)數(shù)據(jù)集上使用不同的統(tǒng)計(jì)方法應(yīng)用了GRSA。對(duì)于經(jīng)典的兩組設(shè)計(jì)，參數(shù)方法的Jaccard相似度超過了0.84，非參數(shù)方法超過了0.78，但參數(shù)方法與非參數(shù)方法之間的Jaccard相似度低于0.63。主要差異主要源于參數(shù)方法與非參數(shù)方法之間的差異。對(duì)于多組數(shù)據(jù)，如果目標(biāo)是在組間顯著改變的通路中富集，用戶可以選擇差異豐度分析。如果目標(biāo)是富集顯示一致增加或減少模式的通路，則相關(guān)性分析是首選。最后，GRSA還支持其他統(tǒng)計(jì)測(cè)試，例如‘DESeq2’、‘Edger’、‘Limma’、‘ALDEX’和‘ANCOM’，以計(jì)算報(bào)告分?jǐn)?shù)。

GRSA具有較高的靈敏度

圖2:GRSA與其他常用富集分析方法的比較。(A - B)箱型圖顯示了在24個(gè)基因表達(dá)數(shù)據(jù)集上通過四種方法得到的目標(biāo)通路的排名(A)和調(diào)整后的p值(B)。數(shù)字表示每種方法的中位數(shù)。(C - D)箱形圖顯示了四種方法在9個(gè)野生型/敲除基因表達(dá)數(shù)據(jù)集上的敏感性(C)和特異性(D)。數(shù)字表示每種方法的中位數(shù)。

接下來(lái)，我們?cè)u(píng)估了GRSA的性能，并將其與其他常用的富集分析方法在幾個(gè)基準(zhǔn)數(shù)據(jù)集上進(jìn)行了比較。在某些情況下，基于通路拓?fù)浣Y(jié)構(gòu)的方法可能比非基于拓?fù)浣Y(jié)構(gòu)的方法更好地識(shí)別具有生物意義的通路。然而，基于拓?fù)浣Y(jié)構(gòu)的方法需要通路的全面結(jié)構(gòu)，限制了它們?cè)谄渌侨祟惿镏械膽?yīng)用。因此，我們將重點(diǎn)放在與非基于拓?fù)浣Y(jié)構(gòu)的富集分析方法的比較上。Nguyen等人提出了幾種比較富集方法的方法，我們采用了他們的方法，并使用相同的通路數(shù)據(jù)庫(kù)（KEGG v109.0）評(píng)估了GRSA與其他流行的富集分析方法的性能。

首先，我們比較了不同方法在識(shí)別與已知人類疾病相關(guān)的24個(gè)基因表達(dá)數(shù)據(jù)集中的目標(biāo)疾病通路的能力。由于每個(gè)數(shù)據(jù)集都與特定的疾病相關(guān)聯(lián)的KEGG通路（即目標(biāo)通路），因此一個(gè)優(yōu)化的富集分析方法應(yīng)該將目標(biāo)通路排名在所有342個(gè)通路的前列，并且以較小的調(diào)整后P值富集目標(biāo)通路。結(jié)果表明，在給定一個(gè)較小的排名給目標(biāo)通路方面，PADOG、GSA、GRSA、GSEA和SAFE的表現(xiàn)相似，因?yàn)樗鼈兊闹形粩?shù)排名都落在前20%內(nèi)。此外，GRSA實(shí)現(xiàn)了目標(biāo)通路的最低中位數(shù)調(diào)整后P值。我們還使用了由‘GSEABenchmarkeR’包提供的GEO列表中的數(shù)據(jù)集進(jìn)行進(jìn)一步的基準(zhǔn)測(cè)試，并發(fā)現(xiàn)了類似的結(jié)果?？傮w而言，GRSA在無(wú)閾值的FCS方法中表現(xiàn)良好，優(yōu)于傳統(tǒng)的ORA方法。

接下來(lái)，我們?cè)u(píng)估了不同方法在檢測(cè)基因敲除實(shí)驗(yàn)中受擾動(dòng)通路的能力。在基因敲除實(shí)驗(yàn)中，敲除基因是已確認(rèn)的擾動(dòng)源。在這些假設(shè)下，我們可以計(jì)算方法的敏感性和特異性。GRSA在考慮的方法中顯示了最高的中位數(shù)敏感性，盡管其特異性略低于其他方法。我們優(yōu)先考慮方法的敏感性，因?yàn)閷?duì)于包含敲除基因的通路，刪除它應(yīng)該會(huì)對(duì)通路產(chǎn)生相當(dāng)大的影響；然而，對(duì)于不包含敲除基因的通路，鑒于通路和基因數(shù)據(jù)庫(kù)的可能不完整性，僅將這些富集的通路歸因?yàn)榧訇?yáng)性可能并不總是恰當(dāng)?shù)摹?/p>

最后，我們?cè)u(píng)估了不同方法富集具有生物意義的通路的能力。我們比較了GRSA、競(jìng)爭(zhēng)工具和兩者都富集到的通路的比例，以所有顯著通路的數(shù)量作為分母。在這些數(shù)據(jù)集中，GRSA一致識(shí)別出比ORA方法更多的通路，并且與GSEA的重疊程度很大。例如，在腎細(xì)胞癌數(shù)據(jù)集中，僅有GRSA富集到與細(xì)胞因子-細(xì)胞因子受體相互作用、IL-17信號(hào)和PI3K-Akt信號(hào)相關(guān)的通路。因此，GRSA具有識(shí)別出與研究的疾病相關(guān)的更多生物學(xué)上相關(guān)通路的潛力。

組學(xué)數(shù)據(jù)實(shí)例

接下來(lái)，我們展示了GRSA在不同類型組學(xué)數(shù)據(jù)中的通用應(yīng)用。

案例研究1:皮膚微生物群的功能分析和與年齡相關(guān)的動(dòng)態(tài)

圖3:GRSA在IHSMGC數(shù)據(jù)集皮膚微生物組中的應(yīng)用。(A) KO-Module網(wǎng)絡(luò)富集于m型(綠色)和c型(藍(lán)色)。只顯示了與維生素生物合成相關(guān)的模塊。大點(diǎn)代表模塊;小點(diǎn)代表ko。小點(diǎn)的顏色代表原型。陰影表示參與同一種維生素生物合成的模塊。深淺顏色表示富含m型剪切型(綠色)或富含c型剪切型(藍(lán)色)的模塊。(B - C)模塊“M00866”(B)和“M00061”(C)不同年齡的箱型圖。線條的顏色代表了KO在模塊中相對(duì)豐度的趨勢(shì)?！癕00866”的陽(yáng)性報(bào)告得分最大(增加)，而“M00061”的陰性報(bào)告得分絕對(duì)值最大(減少)。(D)柱狀圖顯示隨著年齡的增長(zhǎng)，模塊顯著豐富;報(bào)告者得分閾值為2.5，對(duì)應(yīng)的置信度約為0.995，這些模塊根據(jù)KEGG進(jìn)行分組。顏色表示隨著年齡增長(zhǎng)而上調(diào)(紅色)或下調(diào)(綠色)的模塊。

案例研究2:心肌細(xì)胞分化過程中的功能轉(zhuǎn)錄動(dòng)力學(xué)

圖4:GRSA在心肌細(xì)胞分化過程轉(zhuǎn)錄組數(shù)據(jù)集中的應(yīng)用。(A)四個(gè)分化階段基因豐度譜的c均值聚類結(jié)果。每條線的alpha(透明度)與其成員得分的值相關(guān)，y軸表示標(biāo)準(zhǔn)化豐度。(B)箱形圖顯示了“GO:0003015”(心臟過程)在四個(gè)時(shí)間點(diǎn)上的基因豐度;線的顏色表示每個(gè)基因與聚類6在GO項(xiàng)內(nèi)的相關(guān)顯著性。'GO:0003015 '是簇6的代表性術(shù)語(yǔ)。(C)柱狀圖顯示了每個(gè)與分化階段相對(duì)應(yīng)的聚類模式的GO項(xiàng)顯著豐富。條形圖的顏色表示聚類信息，并顯示每個(gè)聚類中報(bào)告者得分較高的代表性GO術(shù)語(yǔ)。左側(cè)的文本標(biāo)簽根據(jù)具有最高表達(dá)的階段著色。一般來(lái)說(shuō)，簇2對(duì)應(yīng)第0天，簇4和5對(duì)應(yīng)第2天，簇1對(duì)應(yīng)第4天，簇3和6對(duì)應(yīng)CM。

案例研究3:母體系統(tǒng)性代謝組學(xué)變化與胎齡相關(guān)

圖5:GRSA在丹麥妊娠隊(duì)列代謝數(shù)據(jù)中的應(yīng)用及IHSMGC數(shù)據(jù)集的物種富集分析。(A)珠包圖顯示了代謝組學(xué)研究中GRSA發(fā)現(xiàn)的顯著富集通路的層次關(guān)系。圓圈的大小表示報(bào)告者得分的絕對(duì)值，圓圈的顏色表示報(bào)告者得分的符號(hào)。陽(yáng)性報(bào)告因子評(píng)分表示該通路增加(橙色)，陰性報(bào)告因子評(píng)分表示該通路減少(紫色)。(B)熱圖顯示了“類固醇激素生物合成”途徑中代謝物的豐度。列是按胎齡增加排序的樣本。(C)柱狀圖顯示C-cutotype和M-cutotype的屬明顯富集。(D)網(wǎng)絡(luò)圖顯示了g_Moraxella和g_Cutibacterium中富含M-cutotype(綠色)或C-cutotype(藍(lán)色)的物種。

案例研究4:GRSA在功能富集分析之外的應(yīng)用

GRSA算法表明，任何組織在層次關(guān)系中的特征都可以用作富集數(shù)據(jù)庫(kù)。例如，我們可以利用微生物的系統(tǒng)發(fā)育關(guān)系，如屬-種關(guān)系，進(jìn)行物種分類富集分析。我們使用了IHSMGC數(shù)據(jù)集的物種豐度表，并查找了在兩種cutotype中富集的屬。我們發(fā)現(xiàn)在M-cutotype中富集的有Psychrobacter、Paracoccus、Chryseobacterium、Elizabethkingia、Deinococcus和Microbacterium，而在C-cutotype中富集的有Acidipropionibacterium、Staphylococcus、Corynebacterium和Cutibacterium（圖5C），其中一些與前期研究中通過共存網(wǎng)絡(luò)發(fā)現(xiàn)的差異物種模塊高度一致。然而，我們額外發(fā)現(xiàn)了一些屬，如Brevundimonas和Rhodobacter，在M-cutotype中富集，而Pahexavirus（丙酸桿菌和皮脂桿菌的噬菌體）在C-cutotype中富集（圖5C），可能是由于GRSA的更高靈敏度。

在先前的研究中，使用了兩種物種Moraxella osloensis和Cutibacterium acnes來(lái)定義cutotype。有趣的是，雖然Cutibacterium屬在cutotype之間是一個(gè)很好的生物標(biāo)志物，但Moraxella屬不是，因?yàn)榘奈锓N沒有共享相同的趨勢(shì)（圖5D）。因此，除了功能富集分析外，GRSA還可以擴(kuò)展到任何層次關(guān)系數(shù)據(jù)結(jié)構(gòu)。

討論與結(jié)論

我們開發(fā)了ReporterScore軟件包，并展示了GRSA富集分析的廣泛應(yīng)用。我們改進(jìn)了經(jīng)典的RSA方法，以便更容易解釋Reporter Score的加減符號(hào)。更重要的是，我們將GRSA的范圍從兩組設(shè)計(jì)擴(kuò)展到多組和縱向設(shè)計(jì)。我們用metagenomic、transcriptomic和metabolomic數(shù)據(jù)進(jìn)行了案例研究驗(yàn)證了先前的關(guān)鍵發(fā)現(xiàn)，并獲得了新的生物學(xué)見解。總之，GRSA和ReporterScore軟件包可以極大地促進(jìn)多種類型的組學(xué)數(shù)據(jù)的功能富集分析，具有更高的靈敏度、與多組和縱向設(shè)計(jì)的兼容性，以及對(duì)定制數(shù)據(jù)庫(kù)的靈活性。

蔣超實(shí)驗(yàn)室博士生彭晨為論文的第一作者，蔣超研究員為論文的通訊作者，參與該研究的還有蔣超實(shí)驗(yàn)室博士后陳瓊，華大研究員譚上進(jìn)和斯坦福大學(xué)博士后申小濤。該研究得到了國(guó)家自然科學(xué)基金等項(xiàng)目資助，以及生命科學(xué)研究院NECHO高性能計(jì)算集群的支持。