|
探索生信之美，解構(gòu)每段代碼背后的故事
▌話說GSEA 通俗來講，GSEA首先將基因在兩種樣品中的差異表達程度（如logFC）或者表型相關(guān)度進行排序（并不是計算差異基因），然后判斷來自功能注釋等預(yù)定義的基因集或自定義的基因集是否傾向于落在有序列表的頂部或底部，也就是說預(yù)定義基因集中大部分的基因是否在樣品中高表達或者低表達，從而判斷功能注釋基因集在不同的樣品中有沒有發(fā)生顯著變化。 ▌GSEA結(jié)果解讀
第1部分是Enrichment Score的折線圖 橫軸排序后的基因，縱軸為對應(yīng)的Running ES, 在折線圖中出現(xiàn)的峰值就是這個基因集的富集分數(shù)（Enrichment Score，ES）。ES是從排序后的表達基因集的第一個基因開始，如果排序表達基因集中的基因出現(xiàn)在功能基因數(shù)據(jù)集中則加分，反之則減分。正值說明在頂部富集，峰值左邊的基因為核心基因，負值則相反。 第2部分為基因位置圖 黑線代表排序后表達數(shù)據(jù)集中的基因存處于當前分析的功能注釋基因集的位置，紅藍相間的熱圖是表達豐度排列，紅色越深的表示該位置的基因logFC越大 ，藍色越深表示logFC越小。如果研究的功能注釋基因集的成員顯著聚集在表達數(shù)據(jù)集的頂部或底部，則說明功能基因數(shù)據(jù)集中的基因在數(shù)據(jù)集中高表達或低表達，若隨機分配，則說明表達數(shù)據(jù)集與該通路無關(guān)。 第3部分表示每個基因?qū)?yīng)的信噪比（Signal2noise） 以灰色面積圖顯展示?；疑幱暗拿娣e比，可以從整體上反映組間的Signal2noise的大小。 除了會看圖形，我們也要關(guān)注一些主要的統(tǒng)計學(xué)指標： NES（校正后的富集分數(shù)） FDR（錯誤發(fā)現(xiàn)率）、p值：FDR<0.25,p<0.05且|NES|>1則表示結(jié)果有意義。 ▌GSEA的優(yōu)點 GSEA是個非常強大的富集分析方法，可以針對多種數(shù)據(jù)庫中的數(shù)據(jù)進行GSEA分析，包括常見的GO數(shù)據(jù)庫，KEGG數(shù)據(jù)庫，Reactome數(shù)據(jù)庫以及MSigDb數(shù)據(jù)庫或其他自定義數(shù)據(jù)集。并且與傳統(tǒng)的GO/KEGG只對顯著的差異基因進行功能注釋相比，GSEA 1. 考慮了基因的表達水平，不需要指定明確的差異基因閾值，檢驗的是基因集合而非單個基因的表達變化，算法會根據(jù)實際數(shù)據(jù)的整體趨勢進行分析。 2. GO和KEGG只能定位到功能，而GSEA通過檢驗預(yù)定義的基因集的基因在某種狀態(tài)下表達水平而提示該通路是否被激活還是抑制。 3. 以待測功能基因集為對象來進行檢驗,使得檢驗結(jié)果針對性和靈敏性提高。 ▌GSEA分析實戰(zhàn) 講了這么多的優(yōu)點，那GSEA分析該怎么做呢？有些經(jīng)驗的小伙伴可能會說GSEA軟件，但是GSEA軟件需要輸入特定的文件格式，圖片也只能是png格式，遠遠達不到發(fā)表級圖片的要求。然而，R語言的clusterprofile包幾條代碼輕松解決所有的難題，還有各種各樣可視化方式，你不來看一看嗎？ #清空環(huán)境變量rm(list = ls())options(stringsAsFactors = F)#加載相關(guān)的包library(data.table)library(clusterProfiler)library(dplyr)library(org.Hs.eg.db)#加載所需物種基因組注釋信息，可有20個物種選擇，http:///packages/release/BiocViews.html#___OrgDblibrary(ggplot2)library(enrichpl)準備文件# 讀入差異分析后文件，包含基因列（gene symbol或entrizid）和logFC列data <- read.table('D:/R study/data_input.txt',sep='\t',header=T,check.names=F)#查看文件內(nèi)容head(data)# SYMBOL logFC#1 FAM182B -5.612055#2 MEP1B -5.134941#3 TMEM167B -3.383565#4 GRIK5 -4.903500#5 XKR9 -5.718563#6 NDST2 -6.843956#本次文件基因名為gene symbol，將gene symbol轉(zhuǎn)換為ENTREZID#如果你的文件是entrizid則跳過該條語句gene.df <- bitr(data$SYMBOL, fromType = 'SYMBOL', toType = c('ENTREZID'), OrgDb = org.Hs.eg.db)#查看前6行數(shù)據(jù)head(gene.df) SYMBOL ENTREZID#1 FAM182B 728882#2 MEP1B 4225#3 TMEM167B 56900#4 GRIK5 2901#5 XKR9 389668#6 NDST2 8509data2=merge(data,gene.df,by.y='SYMBOL')#合并數(shù)據(jù)geneList<-data2$logFC #第二列可以是foldchange，也可以是logFCnames(geneList)=data2$ENTREZID #使用轉(zhuǎn)換好的ENTREZIDgeneList=sort(geneList,decreasing = T) #從高到低排序接下來開始正式的GSEA分析啦#使用GO數(shù)據(jù)庫進行GSEA分析ego <- gseGO( geneList = geneList, OrgDb = org.Hs.eg.db, ont = 'CC', nPerm = 1000, minGSSize = 100, maxGSSize = 500, pvalueCutoff = 0.05, verbose = FALSE, seed = FALSE, by = 'fgsea')##geneList是前面整理好的排序后的ENTREZID##OrgDb選擇物種注釋的數(shù)據(jù)庫##ont同enrichGO分析，可選3大分類BP、CC、MF或ALL##nPerm置換檢驗的次數(shù)，默認為1000##minGSSize富集到某個條目的最小包含基因數(shù)，如果基因數(shù)小于該值則這個條目將被過濾掉，默認為10##maxGSSize富集到某個條目的最大包含基因數(shù)，如果基因數(shù)大于該值則這個條目將被過濾掉，默認為500##verbose是否輸出提示信息##seed是否使結(jié)果具有可重復(fù)性##by選擇使用的統(tǒng)計學(xué)方法，默認為fgseahead(ego)[1:6][1:6]## ID Description setSize enrichmentScore NES##GO:0098687 GO:0098687 chromosomal region 320 0.4187365 1.667091##GO:0000793 GO:0000793 condensed chromosome 209 0.4459184 1.713171##GO:0000775 GO:0000775 chromosome, centromeric region 191 0.4502879 1.722934##GO:0000776 GO:0000776 kinetochore 133 0.4537928 1.677579##GO:0000779 GO:0000779 condensed chromosome, centromeric region 118 0.5091699 1.850798##GO:0072562 GO:0072562 blood microparticle 114 0.4655046 1.680016## pvalue##GO:0098687 0.001052632##GO:0000793 0.001113586##GO:0000775 0.001123596##GO:0000776 0.001154734##GO:0000779 0.001187648##GO:0072562 0.001197605##保存結(jié)果write.table(ego,file='goGSEA.txt',sep=' ',quote=F,row.names = F)#使用KEGG數(shù)據(jù)庫進行GSEA分析kk <- gseKEGG(geneList, organism = 'hsa',pvalueCutoff = 0.2, nPerm = 1000, minGSSize = 10, maxGSSize = 500, verbose = TRUE, seed = FALSE, by = 'fgsea')##參數(shù)同上##按照enrichment score從高到低排序，便于查看富集通路kk_order<-kk[order(kk$enrichmentScore,decreasing=T)]head(kk_order)[1:6][1:6]## ID Description setSize enrichmentScore NES pvalue##hsa03030 hsa03030 DNA replication 36 0.5600885 1.664196 0.008415147##hsa04740 hsa04740 Olfactory transduction 287 0.5005390 1.968686 0.001054852##hsa04110 hsa04110 Cell cycle 124 0.4983236 1.807226 0.001165501##hsa04610 hsa04610 Complement and coagulation cascades 84 0.4868227 1.685794 0.004944376##hsa03008 hsa03008 Ribosome biogenesis in eukaryotes 78 0.4607935 1.581720 0.009962640##hsa05206 hsa05206 MicroRNAs in cancer 160 0.3695695 1.371388 0.019252548##保存結(jié)果write.table(kk_order,file='keggGSEA.txt',sep=' ',quote=F,row.names = F)#使用MsigDb數(shù)據(jù)庫進行GSEA分析##下載離線GMT文件，可以根據(jù)你的需要下載，也可以是全部，網(wǎng)址如下##https://www./gsea/downloads.jspall<-read.gmt('msigdb.v7.4.entrez.gmt') #讀gmt文件gsea<-GSEA(geneList,TERM2GENE = all) ##參數(shù)大致同上，不同是多了TERM2GENE和TERM2NAME##TERM2GENE數(shù)據(jù)集中條目名稱TERM與基因的對應(yīng)關(guān)系,一般是兩列的數(shù)據(jù)框形式##TERM2NAME條目名稱與名稱的對應(yīng)關(guān)系，可用于個性化的GSEA分析，如系定義的基因集，可以為空 write.table(gsea,file='GSEA.txt',sep=' ',quote=F,row.names = F)#使用reactome數(shù)據(jù)庫進行GSEA分析library(ReactomePA)pathway <- gsePathway(geneList)write.table(pathway,file='pathwayGSEA.txt',sep=' ',quote=F,row.names = F)▌GSEA分析可視化 除了上方介紹的經(jīng)典的GSEA結(jié)果展示圖以外，GSEA和GO、KEGG一樣，都有很多讓人眼花繚亂的可視化方法，讓我們來看一看把。
#分面點圖激活和抑制dotplot(kk,split='.sign') facet_grid(~.sign)
# 條形圖##加載需要的R包library(forcats)library(ggstance)##arrange去給NES的絕對值從大到小排序y <- arrange(x, abs(NES)) %>% group_by(sign(NES))##繪圖ggplot(y, aes(NES, fct_reorder(Description, NES), fill=qvalues), showCategory=10) geom_barh(stat='identity') scale_fill_continuous(low='red', high='blue') theme_minimal() ylab(NULL)
#網(wǎng)絡(luò)圖library(enrichplot)ego2 <- pairwise_termsim(ego)##直接運行emapplot會報錯，要加上這條語句emapplot(ego2,showCategory=10,color='pvalue',cex_category=1,layout='kk')##cex_category調(diào)節(jié)節(jié)點的大小##showCategory展示顯示的條目個數(shù)
# 經(jīng)典的gsea結(jié)果圖library(enrichplot)gseaplot2(gsea,1,color='red',pvalue_table = T,title='',base_size=10,ES_geom='line') # 按第一個條目做二維碼圖，并顯示p值#color是enrichment score線的顏色#base_sizexy軸標題字的大小#ES_geom是enrichment score線用線或點顯示
#山脊圖ridgeplot(kk,fill = 'pvalue',5) scale_fill_continuous(type = 'viridis')
好了，今天的分享就到這結(jié)束啦，換上自己的數(shù)據(jù)趁熱趕緊試一下吧！發(fā)表文章時千萬別忘記引用R包文獻哦。
歡迎大家關(guān)注解螺旋生信頻道-挑圈聯(lián)靠公號~ |
|
|
