天堂精品a,午夜在线一区二区

Introduction

功能富集分析是一種用于分析基因集合或基因組數(shù)據(jù)中功能模式富集程度的計(jì)算方法。它可以幫助揭示在特定生物學(xué)背景下，哪些功能模塊、代謝通路、基因家族等在統(tǒng)計(jì)上富集或顯著過(guò)表示。

功能富集分析通常包括以下步驟：

數(shù)據(jù)預(yù)處理：根據(jù)研究問(wèn)題和數(shù)據(jù)類(lèi)型，選擇適當(dāng)?shù)幕蚣匣蚧蚪M數(shù)據(jù)，例如基因表達(dá)數(shù)據(jù)、基因注釋數(shù)據(jù)或基因列表。
注釋和功能分類(lèi)：將基因集合與已知的功能注釋數(shù)據(jù)庫(kù)進(jìn)行比較，例如基因本體論（Gene Ontology）、KEGG（Kyoto Encyclopedia of Genes and Genomes）通路數(shù)據(jù)庫(kù)等。這一步將基因與特定的功能或生物過(guò)程相關(guān)聯(lián)。
統(tǒng)計(jì)分析：使用合適的統(tǒng)計(jì)方法，如超幾何分布、Fisher’s exact test、GSEA（基因集富集分析）等，評(píng)估每個(gè)功能的富集程度。這些方法會(huì)計(jì)算一個(gè)得分或P值，用于判斷功能是否在給定基因集合中富集。
結(jié)果解釋和可視化：根據(jù)統(tǒng)計(jì)分析的結(jié)果，識(shí)別在給定條件下顯著富集的功能模塊，并將結(jié)果進(jìn)行解釋和可視化。這可以幫助研究人員理解基因集合或基因組數(shù)據(jù)中的生物學(xué)特征和功能。

功能富集分析可應(yīng)用于多個(gè)研究領(lǐng)域，如基因表達(dá)分析、蛋白質(zhì)組學(xué)、微生物組學(xué)等。它可以幫助研究人員理解基因集合的生物學(xué)意義，從而揭示生物過(guò)程、代謝通路、細(xì)胞組分等在特定條件下的調(diào)控機(jī)制，并為進(jìn)一步的實(shí)驗(yàn)設(shè)計(jì)和研究提供有價(jià)值的指導(dǎo)。

$P=1-\sum_{i=0}^{m-1}\frac{C_M^iC_{N-M}^{n-i}}{C_N^n}$

P：某pathway的富集顯著性；N：注釋上KEGG的所有基因的數(shù)量；n：所有顯著差異的基因數(shù)量；M：所有基因中注釋到某pathway的基因數(shù)量；m：所有差異基因中注釋到某pathway的基因數(shù)量

R函數(shù)phyper:

1-phyper(k-1,m, N-m, n,)
phyper(k-1,M, N-M, n, lower.tail=F)

Reporter score是一種改良的微生物富集分析的新方法，此方法最初是為了揭示代謝網(wǎng)絡(luò)中的轉(zhuǎn)錄調(diào)控模式而開(kāi)發(fā)的，目前已被引入微生物研究中進(jìn)行功能富集分析。

Method

Reporter score算法最初由Patil和Nielsen于2005年開(kāi)發(fā)，用于識(shí)別代謝調(diào)節(jié)熱點(diǎn)的代謝物 (1)。

應(yīng)用于宏基因組分析，則是基于基因的KO（KEGG orthology，同源基因）注釋?zhuān)@得KO的差異信息，再”上升”至KEGG pathway的功能層面。主要步驟如下：

使用Wilcoxon秩和檢驗(yàn)獲得兩分組間每個(gè)KO差異顯著性的P值（即 $P_{koi}$ ，i代表某個(gè)KO）；
采用逆正態(tài)分布，將每個(gè)KO的P值轉(zhuǎn)化為Z值（ $Z_{koi}$ ）,公式： $Z_{koi}=\theta ^{-1}(1-P_{koi})$ ；
將KO”上升”為pathway： $Z_{koi}$ ，計(jì)算通路的Z值, $Z_{pathway}=\frac{1}{\sqrt{k}}\sum Z_{koi}$ ，其中k表示對(duì)應(yīng)通路共注釋到k個(gè)KO；
評(píng)估顯著程度：置換（permutation）1000次，獲得 $Z_{pathway}$ 的隨機(jī)分布，公式： $Z_{adjustedpathway}=(Z_{pathway}-\mu _k)/\sigma _k$ ， $μ_k$ 為隨機(jī)分布的均值， $σ_k$ 為隨機(jī)分布的標(biāo)準(zhǔn)差。

最終獲得的 $Z_{adjustedpathway}$ ，即為每條代謝通路富集的Reporter score值，Reporter score是非方向性的，Reporter score越大代表富集越顯著，但不能指示通路的上下調(diào)信息。

Misuse

最近有一篇文章就討論了reporter-score的正負(fù)號(hào)誤用問(wèn)題 (2)：

主要結(jié)論是 reporter score算法（上述）是一種忽略通路中KOs上/下調(diào)節(jié)信息的富集方法，直接將reporter score的符號(hào)視為通路的調(diào)節(jié)方向是不正確的。

但是我們應(yīng)該可以將其改為能夠考慮通路內(nèi)KO上下調(diào)的方式，我稱(chēng)為directed 模式, 參考自https://github.com/wangpeng407/ReporterScore。

具體步驟如下:

使用Wilcoxon秩和檢驗(yàn)或者t.test獲得兩分組間每個(gè)KO差異顯著性的P值（即 $P_{koi}$ ，i代表某個(gè)KO），再將P值除以2，即將(0,1]的范圍變?yōu)?0,0.5]， $P_{koi}=P_{koi}/2$ ；
采用逆正態(tài)分布，將每個(gè)KO的P值轉(zhuǎn)化為Z值（ $Z_{koi}$ ）,公式： $Z_{koi}=\theta ^{-1}(1-P_{koi})$ ，由于上述P值小于0.5，則Z值將全部大于0；
考慮每個(gè)KO是上調(diào)還是下調(diào)，計(jì)算 $\Delta KO_i$ ，

$\Delta KO_i=\overline {KO_{i_{g1}}}-\overline {KO_{i_{g2}}}$

其中， $\overline {KO_{i_{g1}}}$ 是組1的 $KO_i$ 的平均豐度, $\overline {KO_{i_{g2}}}$ 是組2的 $KO_i$ 的平均豐度，然后：

$Z_{koi} = \begin{cases} -Z_{koi}, & (\Delta KO_i<0) \\ Z_{koi}, & (\Delta KO_i \ge 0) \end{cases}$

這樣的話 $Z_{koi}$ 大于0為上調(diào)， $Z_{koi}$ 小于0為下調(diào)。

將KO”上升”為pathway： $Z_{koi}$ ，計(jì)算通路的Z值, $Z_{pathway}=\frac{1}{\sqrt{k}}\sum Z_{koi}$ ，其中k表示對(duì)應(yīng)通路共注釋到k個(gè)KO；
評(píng)估顯著程度：置換（permutation）1000次，獲得 $Z_{pathway}$ 的隨機(jī)分布，公式： $Z_{adjustedpathway}=(Z_{pathway}-\mu _k)/\sigma _k$ ， $μ_k$ 為隨機(jī)分布的均值， $σ_k$ 為隨機(jī)分布的標(biāo)準(zhǔn)差。

最終獲得的 $Z_{adjustedpathway}$ ，即為每條代謝通路富集的Reporter score值，在這種模式下，Reporter score是方向性的，更大的正值代表顯著上調(diào)富集，更小的負(fù)值代表顯著下調(diào)富集。

但是這種方法的缺點(diǎn)是當(dāng)一條通路顯著上調(diào)KO和顯著下調(diào)KO差不多時(shí)，最終的Reporter score絕對(duì)值可能會(huì)趨近0，成為沒(méi)有被顯著富集的通路。

Rpackage

因?yàn)槲铱茨壳皼](méi)有現(xiàn)成的工具完成Reporter Score分析（除了一些云平臺(tái)，但可能不太方便），所以我參考https://github.com/wangpeng407/ReporterScore 寫(xiě)了一個(gè)R包幫助分析（雖然也不是特別復(fù)雜）

地址：https://github.com/Asa12138/ReporterScore

安裝方法：

install.packages("devtools")
devtools::install_github('Asa12138/ReporterScore',dependencies=T)

使用方法：

library(ReporterScore)
library(dplyr)
library(ggplot2)
#準(zhǔn)備KO豐度表和實(shí)驗(yàn)metadata
data(KO_test)

head(KO_abundance)

##                CG1         CG2         CG3         CG4         CG5         CG6
## K03169 0.002653545 0.005096380 0.002033923 0.000722349 0.003468322 0.001483028
## K07133 0.000308237 0.000280458 0.000596527 0.000859854 0.000308719 0.000878098
## K03088 0.002147068 0.002030742 0.003797459 0.004161979 0.002076596 0.003091182
##                CG7         CG8         CG9        CG10        CG11        CG12
## K03169 0.002261685 0.004114644 0.002494258 0.002793671 0.004053729 0.002437170
## K07133 0.000525566 0.000356138 0.000445409 0.000268306 0.000293546 0.000465780
## K03088 0.003098506 0.002558730 0.002896506 0.002618472 0.002367986 0.002082786
##               CG13        CG14        CG15         EG1         EG2         EG3
## K03169 0.002187500 0.001988374 0.002304885 0.003317368 0.001150671 0.002610814
## K07133 0.000507992 0.000409447 0.000327910 0.002916018 0.004820742 0.001973789
## K03088 0.002680792 0.003066870 0.002975895 0.002257401 0.002889640 0.001997586
##                EG4         EG5         EG6         EG7         EG8         EG9
## K03169 0.000900673 0.001545374 0.001640295 0.001445024 0.001096728 0.001026556
## K07133 0.003359927 0.001913932 0.001384079 0.001321643 0.002376473 0.004391014
## K03088 0.002613441 0.002803388 0.002251835 0.002981244 0.002944061 0.003113215
##               EG10        EG11        EG12        EG13        EG14        EG15
## K03169 0.001513195 0.001812732 0.003256782 0.006723067 0.001769819 0.001233307
## K07133 0.002479040 0.003484868 0.000790457 0.000127818 0.000634529 0.004746572
## K03088 0.003177522 0.002790092 0.001607913 0.002574928 0.001662157 0.002614489
##  [ reached 'max' / getOption("max.print") -- omitted 3 rows ]

head(Group_tab)

##     Group
## CG1    CG
## CG2    CG
## CG3    CG
## CG4    CG
## CG5    CG
## CG6    CG

分組檢驗(yàn)獲得P值,threads多線程可加速

ko_pvalue=ko_test(KO_abundance,"Group",Group_tab,threads = 1,verbose = F)

## Compared groups: CG and EG
## Total KO number: 4535
## Time use: 1.348

head(ko_pvalue)

##    KO_id       avg_CG        sd_CG       avg_EG        sd_EG      diff_mean
## 1 K03169 0.0026728975 0.0011094132 0.0020694937 0.0014922004  0.00060340387
## 2 K07133 0.0004554658 0.0001951678 0.0024480601 0.0014916566 -0.00199259427
## 3 K03088 0.0027767713 0.0006253559 0.0025519275 0.0004966801  0.00022484380
## 4 K03530 0.0005779169 0.0008952163 0.0005197504 0.0001435245  0.00005816647
## 5 K06147 0.0020807307 0.0007731661 0.0014321838 0.0004273716  0.00064854693
## 6 K05349 0.0021064422 0.0005243558 0.0017419317 0.0005382770  0.00036451047
##         p.value
## 1 0.03671754164
## 2 0.00002654761
## 3 0.48636476395
## 4 0.01125600770
## 5 0.00987482265
## 6 0.12614740102

2.將P值矯正并轉(zhuǎn)為Z-Score，這里提供兩種方法（mixed就是經(jīng)典的方法，另一種是directed方法）

ko_stat=pvalue2zs(ko_pvalue,mode="directed")

## ================================================================================ 
## 
##  Chi-squared test for given probabilities
## 
## data:  up_down_ratio
## X-squared = 21.823, df = 1, p-value = 0.000002991

head(ko_stat)

##    KO_id       avg_CG        sd_CG       avg_EG        sd_EG      diff_mean
## 1 K03169 0.0026728975 0.0011094132 0.0020694937 0.0014922004  0.00060340387
## 2 K07133 0.0004554658 0.0001951678 0.0024480601 0.0014916566 -0.00199259427
## 3 K03088 0.0027767713 0.0006253559 0.0025519275 0.0004966801  0.00022484380
## 4 K03530 0.0005779169 0.0008952163 0.0005197504 0.0001435245  0.00005816647
## 5 K06147 0.0020807307 0.0007731661 0.0014321838 0.0004273716  0.00064854693
## 6 K05349 0.0021064422 0.0005243558 0.0017419317 0.0005382770  0.00036451047
##         p.value sign     type      q.value    Z_score
## 1 0.01835877082    1 Enriched 0.0669268695  1.4990767
## 2 0.00001327381   -1 Depleted 0.0004777517 -3.3033110
## 3 0.24318238198    1 Enriched 0.3058325297  0.5076981
## 4 0.00562800385    1 Enriched 0.0304570375  1.8741186
## 5 0.00493741133    1 Enriched 0.0277461715  1.9150013
## 6 0.06307370051    1 Enriched 0.1607865272  0.9912305

3.將KO”上升”為pathway，計(jì)算ReporterScore：

reporter_s=get_reporter_score(ko_stat)

## =================================Checking file================================= 
## ==================================load KOlist================================== 
## ===================KOlist download time: 2023-05-12 00:07:41=================== 
## ============================Calculating each pathway============================ 
## ID number: 479
## Time use: 4.641

head(reporter_s)

##         ID ReporterScore                              Description K_num
## 1 map00010    -1.6255358             Glycolysis / Gluconeogenesis   106
## 2 map00020    -1.8906022                Citrate cycle (TCA cycle)    67
## 3 map00030    -0.8261448                Pentose phosphate pathway    88
## 4 map00040    -0.6685755 Pentose and glucuronate interconversions    89
## 5 map00051    -2.2494558          Fructose and mannose metabolism   112
## 6 map00052    -0.4196662                     Galactose metabolism    78

結(jié)果進(jìn)行繪圖

plot_report(reporter_s,rs_threshold=c(2,-7),y_text_size=10,str_width=40)+
    labs(title = "CG vs EG")

plot_report(reporter_s,rs_threshold=c(2,-7),mode = 2,y_text_size=10,str_width=40)+
    labs(title = "CG vs EG")

5.挑選一條通路進(jìn)行繪制

plot_KOs_in_pathway(map_id = "map00780",ko_stat = ko_stat)

Reference

K. R. Patil, J. Nielsen, Uncovering transcriptional regulation of metabolism by using metabolic network topology. Proceedings of the National Academy of Sciences of the United States of America. 102, 2685–2689 (2005).

L. Liu, R. Zhu, D. Wu, Misuse of reporter score in microbial enrichment analysis. iMeta. n/a, e95.

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av

Reporter Score 微生物功能富集分析

Reporter Score 微生物功能富集分析

Introduction

Method

Misuse

Rpackage

Reference

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九 欧美,1769亚洲,黄色成人av

Reporter Score 微生物功能富集分析

Introduction

Method

Misuse

Rpackage

Reference

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容

色偷偷精品伊人,欧洲久久精品,欧美综合婷婷骚逼,国产AV主播,国产最新探花在线,九色在线视频一区,伊人大交九欧美,1769亚洲,黄色成人av