這部分開始進(jìn)行基本的富集分析,，兩類

• A：差異基因富集分析（不需要表達(dá)值,，只需要gene name）
• B: 基因集(gene set)富集分析（不管有無差異,，需要全部genes表達(dá)值）

############################################################

A：差異基因富集分析（不需要表達(dá)值，只需要gene name）

############################################################

-----------先說富集什么-----------

• 最常用的基因注釋工具是GO和KEGG注釋,，這基本上是差異基因分析一定做的兩件事,。GO可以在GO:BP(生物過程），GO:MF(分子功能）,，GO:CC(細(xì)胞組分）三個(gè)方面分別進(jìn)行注釋,，用的比較多的是GO:BP，但其他兩方面也很重要,。
• 外還有一個(gè)軟件不得不提,，那就是IPA（Ingenuity pathway analysis），這是一個(gè)收費(fèi)軟件,，有基本版和高級版,。我個(gè)人覺得它的upstream regulator analysis還是很不錯(cuò)的。分子激活功能等也可以用用,。另外一個(gè)就是它內(nèi)置的熱圖功能,。高級版我沒用過，但是知道可以導(dǎo)出一些數(shù)據(jù)等,。

-----------什么是富集（原理）-----------

富集的統(tǒng)計(jì)學(xué)基礎(chǔ)是超幾何分布,，簡單來說根據(jù)下面的Fisher精確檢驗(yàn)（Fisher exact test）公式，對每個(gè)GO或KEGG term計(jì)算一個(gè)p值 p=(M/K)[(N-M)/(n-k)]/(N/n),，其中 N：所有g(shù)ene總數(shù) n：N中差異表達(dá)gene的總數(shù) M：N中屬于某個(gè)GO term的gene個(gè)數(shù) k: n中屬于某個(gè)GO term的gene個(gè)數(shù) p：表示差異表達(dá)gene富集到這個(gè)GO term上的可信程度

• 當(dāng)p<0.05或0.01,，則認(rèn)為差異表達(dá)gene顯著到這個(gè)GO term上（自己定義p值）
• 意義：提供的信息更集中，更有意義

---------------拿什么來富集---------------

得到的差異表達(dá)基因列表就可以,，也就是說不需要其他的值

---------------用什么工具富集--------------

只能說實(shí)在是太多太多了,。。,。,。但是用的時(shí)候要小心，因?yàn)槟愣嘤脦讉€(gè)工具，即使設(shè)定同樣的p值也會發(fā)現(xiàn)結(jié)果有出入,，有時(shí)還差異挺大,。

1 按使用方式來說（簡單度）有3種

• （1）在線版：最主流的就是DAVID,各種級別雜志總見其身影，使用非常簡單,，不再贅述,。另外還有Gather，GOrilla,revigo,還有很多很多我就不在貼了,。網(wǎng)頁版有網(wǎng)頁版的好處,，可以先大概看下自己篩選的genes。另外很多工具有很好的可視化功能,，自己一一去探索吧,。
• （2）客戶端版：IPA（IPA不是用的GO和KEGG數(shù)據(jù)庫）和FUNRICH，后者更新速度很慢,，但里面有好玩又實(shí)用的功能,，并且可以加載自己的數(shù)據(jù)。
• （3）R包：介紹一個(gè)就行了,，那就是Y叔的clusterProfiler,，我論文中的富集功能很多都是用這個(gè)包做的（還有的用了IPA）。 ##########################################################

B: 基因集(gene set)富集分析（不管有無差異,，需要全部genes表達(dá)值）

##########################################################

• 好處：可以發(fā)現(xiàn)被差異基因舍棄的genes可能參與了某重要生理過程或信號通路(稍后我會把以前手寫的GSEA原理和結(jié)果意義解讀發(fā)上來)
• 工具：GSEA
• 使用方法：R（還是clusterProfiler）或客戶端

-------------------具體操作---------------------

#enrichment analysis using clusterprofiler package created by yuguangchuang
library(clusterProfiler)
library(DOSE)
library(org.Mm.eg.db)
#get the ENTREZID for the next analysis
setwd("F:/rna_seq/data/matrix")
sig.gene<-read.csv(file="DEG_treat_vs_control.csv")
head(sig.gene)
gene<-sig.gene[,1]
head(gene)
gene.df<-bitr(gene, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = org.Mm.eg.db)

head(gene.df)

GO富集分析：

#Go classification
#Go enrichment
ego_cc<-enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Mm.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "CC",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
ego_bp<-enrichGO(gene       = gene.df$ENSEMBL,
                 OrgDb      = org.Mm.eg.db,
                 keyType    = 'ENSEMBL',
                 ont        = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)
barplot(ego_bp,showCategory = 18,title="The GO_BP enrichment analysis of all DEGs ")+ 
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(ego_bp) str_wrap(ego_bp,width = 25))

gobp.jpeg

#KEGG enrichment
install.packages("stringr")
library(stringr)
kk<-enrichKEGG(gene      =gene.df$ENTREZID,
               organism = 'mmu',
               pvalueCutoff = 0.05)
kk[1:30]
barplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))
dotplot(kk,showCategory = 25, title="The KEGG enrichment analysis of all DEGs")+
  scale_size(range=c(2, 12))+
  scale_x_discrete(labels=function(kk) str_wrap(kk,width = 25))

kegg.jpeg

keggdot.jpeg

# Gene Set Enrichment Analysis（GSEA）
# 獲取按照log2FC大小來排序的基因列表
genelist <- sig.gene$log2FoldChange
names(genelist) <- sig.gene[,1]
genelist <- sort(genelist, decreasing = TRUE)
# GSEA分析
gsemf <- gseGO(genelist,
               OrgDb = org.Mm.eg.db,
               keyType = "ENSEMBL",
               ont="BP"
)
# 查看信息
head(gsemf)
# 畫出GSEA圖
gseaplot(gsemf, geneSetID="GO:0001819")

gsea.jpeg

后記：做完這部分富集分析,，接著按我的流程進(jìn)入下一部分分析RNA-seq(10):KEGG通路可視化，因?yàn)橹苯佑玫竭@部分?jǐn)?shù)據(jù),，

參考Y叔的包說明,，里面寫的特別詳細(xì) 還有lxmic的簡書