生信實操|如何利用R語言進行GSEA分析

專注生物分析最前沿

定期解讀生信文章

提供生信分析思路和套路

看圖說話欄目曾介紹過GSEA的原理（看圖說話|GSEA分析--教你解鎖高級的富集分析），今天我們來看一下如何利用R語言進行GSEA分析。

如果你有RNA-seq的數據，就可以這樣做，先把數據整成這樣，一共兩列，一列是SYMBOL，一列是foldChange。

然後輸入代碼：

> setwd("E:\\ 20201214GSEA分析怎麼做")
> #if (!requireNamespace("BiocManager", quietly = TRUE))
> #  install.packages("BiocManager")
> #BiocManager::install("clusterProfiler")
> library(clusterProfiler)#主要用到clusterProfiler包
> df = read.table("easy_input.txt",header = T,as.is = T)
> head(df)#查看前面幾行
    SYMBOL  foldChange
1      A2M -0.50361737
2     NAT1 -3.25012047
3     NAT2 -0.28492113
4 SERPINA3 -0.63137249
5    AADAC -0.18896102
6     AAMP -0.03605823
> dim(df)#數據總共幾行幾列
[1] 12444     2
>

可以看出，輸入的數據總共12444個基因，當然，這並不僅僅是是差異基因。
在分析之前，還是要把SYMBOL轉換成ENTREZ ID，之前講過，ENTREZ ID在染色體上是有編號的，一個編號對應一個基因，錯不了，但是如果是gene name，很容易出錯，另外也不好計算。
這裡要用到clusterProfiler包裡面的bitr函數，這個函數很方便的就可以轉換ENTREZID和SYMBOL。

> df.id<-bitr(df$SYMBOL, #轉換的列是df數據框中的SYMBOL列
+             fromType = "SYMBOL",#需要轉換ID類型
+             toType = "ENTREZID",#轉換成的ID類型
+             OrgDb = "org.Hs.eg.db")#對應的物種，小鼠的是org.Mm.eg.db
'select()' returned 1:many mapping between keys and columns
Warning message:
In bitr(df$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db") :
  1.28% of input gene IDs are fail to map
> head(df.id)
    SYMBOL ENTREZID
1      A2M        2
2     NAT1        9
3     NAT2       10
4  SERPINA3       12
5    AADAC       13
6     AAMP       14
> dim(df.id)
[1] 12287     2
>

代碼說有1.28%沒有map，所以，剛才是12444個基因，現在就變成了了12287個基因。
接下來，將最開始的文件和轉成ENTREZID的文件合併，用到merge函數。

> #讓基因名、ENTREZID、foldchange對應起來
> easy.df<-merge(df,df.id,by="SYMBOL",all=F)
> head(easy.df)
  SYMBOL  foldChange ENTREZID
1   A1CF -0.13606162    29974
2    A2M -0.50361737        2
3 A4GALT  0.22796773    53947
4  A4GNT -0.08822687    51146
5   AAAS -0.12161476     8086
6   AACS  0.23240224    65985

為了做GSEA分析，先要對所有的基因進行排序，就是大家常看到的那個蝴蝶圖，這裡用foldchange作為量化的標準，將基因從大到小的排序。

> #按照foldchange排序
> sortdf<-easy.df[order(easy.df$foldChange, decreasing = T),]#降序排序
> head(sortdf)
     SYMBOL foldChange ENTREZID
6398   MMP1   4.572613     4312
1724  CDC45   4.514594     8318
7042    NMU   4.418218    10874
1731  CDCA8   4.144075    55143
6190  MCM10   3.876258    55388
1706  CDC20   3.677857      991
> gene.expr = sortdf$foldChange#把foldchange按照從大到小提取出來
> head(gene.expr)
[1] 4.572613 4.514594 4.418218 4.144075 3.876258 3.677857
> names(gene.expr) <- sortdf$ENTREZID#給上面提取的foldchange對應上ENTREZID
> head(gene.expr)
    4312     8318    10874    55143    55388      991 
4.572613 4.514594 4.418218 4.144075 3.876258 3.677857

有了這個ENTREZID和foldchange信息的gene.expr之後，就可以使用clusterProfiler進行GSEA分析了，進行GSEA分析。
這裡就很簡單，GSEA就一句代碼：

kk <- gseKEGG(gene.expr, organism = "hsa")

是不是出奇的簡單。

> kk <- gseKEGG(gene.expr, organism = "hsa")
preparing geneSet collections...
GSEA analysis...
leading edge analysis...
done
Warning messages:
1: In preparePathwaysAndStats(pathways, stats, minSize, maxSize, gseaParam,  :
  There are ties in the preranked stats (0.02% of the list).
The order of those tied genes will be arbitrary, which may produce unexpected results.
2: In serialize(data, node$con) :
  'package:stats' may not be available when loading
3: In serialize(data, node$con) :
  'package:stats' may not be available when loading
4: In serialize(data, node$con) :
  'package:stats' may not be available when loading
5: In serialize(data, node$con) :
  'package:stats' may not be available when loading
6: In fgseaMultilevel(...) :
  For some pathways, in reality P-values are less than 1e-10. You can set the `eps` argument to zero for better estimation.
> head(kk)
               ID                  Description setSize enrichmentScore      NES
hsa04110 hsa04110                   Cell cycle     114       0.6698529 2.740968
hsa03050 hsa03050                   Proteasome      43       0.7099227 2.425653
hsa05169 hsa05169 Epstein-Barr virus infection     193       0.4337991 1.932303
hsa04657 hsa04657      IL-17 signaling pathway      85       0.5624088 2.172820
hsa03030 hsa03030              DNA replication      33       0.7233029 2.309149
hsa01230 hsa01230  Biosynthesis of amino acids      62       0.5891945 2.178783
               pvalue     p.adjust      qvalues rank                   leading_edge
hsa04110 1.000000e-10 3.240000e-08 2.431579e-08 1210 tags=40%, list=10%, signal=37%
hsa03050 1.358588e-08 2.200913e-06 1.651757e-06 2469 tags=65%, list=20%, signal=52%
hsa05169 8.736724e-08 9.435662e-06 7.081345e-06 2770 tags=39%, list=23%, signal=31%
hsa04657 1.413224e-07 1.144711e-05 8.590914e-06 2830 tags=49%, list=23%, signal=38%
hsa03030 2.526900e-07 1.637431e-05 1.228871e-05 1867 tags=64%, list=15%, signal=54%
hsa01230 1.176431e-06 6.352728e-05 4.767642e-05 2868 tags=55%, list=23%, signal=42%
                                                                                                                                                                                                                                                                                                                                                                               core_enrichment
hsa04110                                                                                                                                                         8318/991/9133/890/983/4085/7272/1111/891/4174/9232/4171/993/990/5347/701/9700/898/23594/4998/9134/4175/4173/10926/6502/994/699/4609/5111/1869/1029/8317/4176/2810/3066/1871/1031/9088/995/1019/4172/5885/11200/7027/1875/7534
hsa03050                                                                                                                                                                                                                                       5688/5709/5698/5693/3458/5713/11047/5721/5691/5685/5690/5684/5686/5695/10213/23198/7979/5699/5714/5702/5708/5692/5704/5683/5694/5718/51371/5682
hsa05169 3627/890/6890/9636/898/9134/6502/6772/3126/3112/4609/917/5709/1869/3654/919/915/4067/4938/864/4940/5713/5336/11047/3066/54205/1871/578/1019/637/916/3383/4939/10213/23586/4793/5603/7979/7128/6891/930/5714/3452/6850/5702/4794/7124/3569/7097/5708/2208/8772/3119/5704/7186/5971/3135/1380/958/5610/4792/10018/8819/3134/10379/9641/1147/5718/6300/3109/811/5606/2923/3108/5707/1432
hsa04657                                                                                                                                                                 4312/6280/6279/6278/3627/2921/6364/8061/4318/3576/3934/6347/727897/1051/6354/3458/6361/6374/2919/9618/5603/7128/1994/7124/3569/8772/5743/7186/3596/6356/5594/4792/9641/1147/2932/6300/5597/27190/1432/7184/64806/3326
hsa03030                                                                                                                                                                                                                                                                           4174/4171/4175/4173/2237/5984/5111/10535/1763/5427/23649/4176/5982/5557/5558/4172/5424/5983/5425/54107/6119
hsa01230                                                                                                                                                                                                             29968/26227/875/445/5214/440/65263/6472/7086/5723/3418/7167/586/2597/5230/2023/5223/5831/6888/50/5315/3419/10993/2805/22934/3421/5634/3417/5232/221823/2027/5211/384/5832

返回的kk文件很複雜,總共48行，11列。

> dim(kk)
[1] 48 11

1，代表的KEGG中的信號通路，比如hsa04110就是KEGG中的Cell cycle的信號通路
2，對信號通路的描述，也就是信號通路的名字
3，該信號通路的基因個數
4，很熟悉的富集分數，就是ES
5，這是標準化以後的ES，全稱是normalized enrichment score
6，P值
7，矯正以後的p值
8，Q值，也有的寫FDR q-val（false discovery rate），表示錯誤發現率
9，這裡的rank就是ES在頂點時候，那個最幸運基因的位置，或者說是排名
10，這裡的leading_edge有點複雜，對富集貢獻最大的基因成員,即領頭亞集，用於定義Leading-edge subset的參數有：Tags，List，Signal，對於一個基因集而言，定義其中對Enrichment score貢獻最大的基因為核心基因，也稱之為leading edge subset，也就是ES頂點的之前的基因。tags表示核心基因佔該基因集基因總數的比例，而list表示核心基因佔所有基因總數的比例，signal利用這兩個指標計算得到，公式如下：

N代表所有基因的數目，Nh代表該基因集下的基因總數。對於一個基因集而言，當核心基因的數目和該基因集下的基因總數相同，signal取值最大，當該基因集的基因數目和所有基因數目接近時，signal的取值接近於0。當然，我們希望的是signal越大越好。
11，這個core_enrichment就是主要富集的基因，下面的數字其實是ENTREZID，就是一個一個基因，換句話說，就是做GSEA分析的基因列表中，Hit到該目的通路的基因列表。

接下裡，需要按照ES的大小，給通路排序，然後保存起來，後面用來可視化。

#按照enrichment score從高到低排序，便於查看富集通路
sortkk<-kk[order(kk$enrichmentScore, decreasing = T),]
head(sortkk)
dim(sortkk)
write.table(sortkk,"gsea_output2.txt",sep = "\t",quote = F,col.names = T,row.names = F)

通過查看gsea_output2.txt,選擇你想畫的pathway，記下通路的term ID
開始畫圖,不需要另外準備輸入文件，clusterProfiler直接可以用gseaplot出圖：

gseaplot(kk, "hsa04510")#很簡單的代碼是吧

可以改變參數，更改線條的顏色。

gseaplot(kk, "hsa04510", color.line='steelblue')

在這裡，我們看到ES是負值，ES是負值，就是表示，功能集中在尾部，就是表示實驗組低表達，對照組高表達。
這裡用到clusterProfiler畫圖，如果覺得不夠美觀，還可以用enrichplot的包來畫圖

#另一種出圖方式是模擬Broad institute的GSEA軟體的圖，可以用gseaplot2出圖：
library(enrichplot)
gseaplot2(kk, "hsa04510")
gseaplot2(kk, "hsa04510", color = "firebrick", rel_heights=c(1, .2, .6))#改變更多參數，為了美觀

是不是感覺有點Nature的味道了。
它還有更外一個功能，就是支持同時展示多個pathways的結果。

gseaplot2(kk, row.names(sortkk)[1:4])

上圖用的是ES排名前4個畫圖，也可以用你自己感興趣的通路畫圖，自己在剛才保存的txt文件裡挑選就行。

paths <- c("hsa04510", "hsa04512", "hsa04974", "hsa05410")
gseaplot2(kk, paths)

這裡的GSEA分析其實由三個圖構成，GSEA分析的runningES折線圖，還有下面基因的位置圖，還有所謂的蝴蝶圖。如果不想同時展示，還可以通過subplots改變。

gseaplot2(kk, paths, subplots=1)#只要第一個圖
gseaplot2(kk, paths, subplots=1:2)#只要第一和第二個圖
gseaplot2(kk, paths, subplots=c(1,3))#只要第一和第三個圖

如果想把p值標上去，也是可以的。

gseaplot2(kk, paths, color = colorspace::rainbow_hcl(4), pvalue_table = TRUE)

到目前為止，基本上用KEGG庫做注釋的GSEA分析就做的差不多了，還想了解更多的，可能需要對clusterProfiler ，enrichplot兩個包好好研讀一下，說不定有驚喜。
最後甩出一個代碼：

gseaplot2(kk,#數據
          row.names(sortkk)[39],#畫那一列的信號通路
          title = "Prion disease",#標題
          base_size = 15,#字體大小
          color = "green",#線條的顏色
          pvalue_table = TRUE,#加不加p值
          ES_geom="line")#是用線，還是用d點