富集分析的原理與實現

TOP生物資訊發表於2021-10-29

一般做完差異分析都會做這一步，目的是找到差異基因富集到的通路，進而與生物學意義聯絡起來。具體的統計方法很簡單，這篇筆記裡面的程式碼可以從零搭建一個富集分析工具。

後臺回覆20211007獲取本文的測試資料和程式碼，以及（單細胞）轉錄組分析中可能用到的GO KEGG富集分析程式碼（這部分本文不演示）。

關於Gene Ontology (GO), KEGG這些背景就不講了，網上很多資料。

將富集分析中的問題抽象出來，其實就是下圖的“摸球”問題。

藍色方框中的球是所有的基因【共N個】，在探究某個特定通路P時，通路里面涉及到的基因用紅色表示【共M個】。綠色圓圈是一次摸球事件，用來表示做了一次差異分析得到的基因【共n個】，這些基因中，有屬於通路P的(紅色球)【共k個】，有不屬於的(黑色球)。

用摸球問題中的語言再描述一次：袋子中共有黑球和紅球N個，其中紅球M個。某次抽樣中，一共摸球n個，其中紅球k個，問在這次摸球中，紅球的佔比是否顯著高於袋子中紅球的佔比？
（以前學摸球問題/超幾何分佈的時候，可能只求概率，沒有進一步到這個統計檢驗）

回答這個問題，需要求出問題中這個事件的概率以及更極端事件的概率之和，也就是p值，小於0.05或者0.01就能認為是顯著了。

1. 一個通路，計算p值

接下來以一個GO term (GO_extracellular_matrix_organization)為例，計算p值。
用的通路基因集是小鼠 GO BP的基因集，差異基因集是單細胞轉錄組分析中一個cluster的高表達基因

library(tidyverse)
library(clusterProfiler)

gmt.df=read.gmt("Mm.c5.bp.v7.1.SYMBOL.gmt")
deg=read.table("test_deg.txt",header = T,sep = "\t",stringsAsFactors = F)
deg=deg[deg$gene %in% gmt.df$gene,]

這種情況下，前面說的幾個引數的值如下：

  • 全部球的個數/全部基因數：
    N=length(unique(gmt.df$gene))
  • 全部紅球的個數/通路基因集的基因數：
    one.set=gmt.df[ gmt.df$term %in% c("GO_extracellular_matrix_organization") ,] M=length(one.set$gene)
  • 摸球數/差異基因數：
    n=length(deg$gene)
  • 摸球中紅球的個數/差異基因中屬於這個通路的基因數：
    k=sum(deg$gene %in% one.set$gene)

N M n k的值分別為: 23210, 271, 47, 6

求p值之前，先看一下滿足N M n這三個引數的超幾何分佈，在不同的k值之下的概率：

df1=data.frame(x=1:47,y=dhyper(x=1:47, M, N-M, n))
df1%>%ggplot(aes(x,y))+geom_point()+
  geom_line()+
  geom_vline(xintercept=k,color="red",linetype=5)+  #這裡k等於6，其作為閾值
  theme_bw()+
  theme(panel.grid = element_blank())

dhyper()用來求概率，四個引數分別是：摸球中紅球個數的向量，袋中紅球數，袋中黑球數，摸球數

我們要計算的就是紅線以及右邊那些紅球數對應的概率之和，如下：

> phyper(k-1,M, N-M, n, lower.tail=FALSE)
[1] 1.722659e-05

當lower.tail=FALSE，計算的是P[X > x]，即大於第一個引數的概率之和。上面的程式碼第一個引數寫的是k-1，因為我們需要求k以及k右邊的概率之和。

以上是對一個通路求p值

2. 多個通路，依次計算p值

如果是多個通路，需要迴圈操作，依次對每個通路進行富集分析。
下面的演示用到的差異基因集和GO BP基因集同上

分析哪些pathway？要滿足兩個條件：

  • 通路里面基因的數量滿足一定要求
  • 至少和deg有基因交集

下面的程式碼就是對通路做過濾的

bp.stat=as.data.frame(table(gmt.df$term))
colnames(bp.stat)[1]="pathway"
bp.stat=bp.stat%>%filter(Freq >= 2 & Freq <= 2000)

tmp.df=gmt.df
tmp.df$TF=tmp.df$gene %in% deg$gene
tmp.stat=as.data.frame(tmp.df %>% dplyr::group_by(term) %>% dplyr::summarize(counts=sum(TF)))
tmp.stat=tmp.stat%>%filter(counts > 0)
keep.pw=sort(intersect(bp.stat$pathway,tmp.stat$term))

下面就是迴圈求p值了

N=length(unique(gmt.df$gene))
n=length(deg$gene)
term=c()
pvalue=c()

for (i in keep.pw) {
  one.set=gmt.df[ gmt.df$term %in% i ,]
  M=length(one.set$gene)
  k=sum(deg$gene %in% one.set$gene)
  one.pvalue=phyper(k-1,M, N-M, n, lower.tail=FALSE)
  term=append(term,i)
  pvalue=append(pvalue,one.pvalue)
}

my_go_res=data.frame(term=term,pvalue=pvalue)
my_go_res=my_go_res%>%arrange(pvalue)

到這會兒還沒有結束，還差一個FDR

3. 計算矯正p值

FDR, False Discovery Rates。為什麼要控制FDR，降低假陽性。
這裡用到的是The Benjamini-Hochberg method

The Benjamini-Hochberg method

假設我們對10個通路做了富集分析，我們會先得到10個p值：

  1. 將這10個p值從小到大排序
  2. 從1到10給這些p值排序
  3. 最大的FDR adjusted p value（第10位）等於原來最大的那個p值
  4. 第9位的FDR adjusted p value等於這兩個值中的較小值：
    ①前一位矯正的p值；
    ②當前未矯正的p值 * （p值總個數/當前位數）
  5. 重複第4步，直到第1位

程式碼如下：

fdr=c()
for (i in dim(my_go_res)[1]:1) {
  if (i==dim(my_go_res)[1]) {
    tmpfdr=my_go_res$pvalue[i]
  }else{
    tmpfdr=min(tmpfdr,my_go_res$pvalue[i] * (dim(my_go_res)[1] / i))
  }
  fdr=append(fdr,tmpfdr)
}
my_go_res$p.adj=rev(fdr)

到這兒富集分析的完整流程才算結束

4. 輪子有現成的

當然，這個演算法已經非常常見了，clusterProfiler的enricher()就能夠自定義基因集做富集分析。使用如下：

deg_gmt=clusterProfiler::enricher(deg$gene,TERM2GENE = gmt.df,minGSSize = 2,maxGSSize = 2000)
go_res=deg_gmt@result

和上面的結果是一模一樣的。

今天的內容就到這裡，後臺回覆20211007獲取本文的測試資料和程式碼，以及（單細胞）轉錄組分析中可能用到的GO KEGG富集分析程式碼（這部分本文不演示）。

ref

因水平有限，有錯誤的地方，歡迎批評指正！

相關文章