這部分開始進(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ì)胞組分）三個方面分別進(jìn)行注釋设拟，用的比較多的是GO:BP慨仿，但其他兩方面也很重要。
- 外還有一個軟件不得不提纳胧，那就是IPA（Ingenuity pathway analysis）镶骗，這是一個收費軟件，有基本版和高級版躲雅。我個人覺得它的upstream regulator analysis還是很不錯的。分子激活功能等也可以用用骡和。另外一個就是它內(nèi)置的熱圖功能相赁。高級版我沒用過，但是知道可以導(dǎo)出一些數(shù)據(jù)等慰于。
-----------什么是富集（原理）-----------
富集的統(tǒng)計學(xué)基礎(chǔ)是超幾何分布钮科，簡單來說根據(jù)下面的Fisher精確檢驗（Fisher exact test）公式，對每個GO或KEGG term計算一個p值
p=(M/K)[(N-M)/(n-k)]/(N/n)婆赠，其中
N：所有g(shù)ene總數(shù)
n：N中差異表達(dá)gene的總數(shù)
M：N中屬于某個GO term的gene個數(shù)
k: n中屬于某個GO term的gene個數(shù)
p：表示差異表達(dá)gene富集到這個GO term上的可信程度
- 當(dāng)p<0.05或0.01绵脯，則認(rèn)為差異表達(dá)gene顯著到這個GO term上（自己定義p值）
- 意義：提供的信息更集中，更有意義
---------------拿什么來富集---------------
得到的差異表達(dá)基因列表就可以，也就是說不需要其他的值
---------------用什么工具富集--------------
只能說實在是太多太多了蛆挫。赃承。。悴侵。但是用的時候要小心瞧剖，因為你多用幾個工具，即使設(shè)定同樣的p值也會發(fā)現(xiàn)結(jié)果有出入可免，有時還差異挺大抓于。
1 按使用方式來說（簡單度）有3種
- （1）在線版：最主流的就是DAVID,各種級別雜志總見其身影，使用非常簡單浇借，不再贅述捉撮。另外還有Gather，GOrilla,revigo,還有很多很多我就不在貼了妇垢。網(wǎng)頁版有網(wǎng)頁版的好處巾遭，可以先大概看下自己篩選的genes。另外很多工具有很好的可視化功能修己，自己一一去探索吧恢总。
- （2）客戶端版：IPA（IPA不是用的GO和KEGG數(shù)據(jù)庫）和FUNRICH，后者更新速度很慢睬愤，但里面有好玩又實用的功能片仿，并且可以加載自己的數(shù)據(jù)。
- （3）R包：介紹一個就行了尤辱，那就是Y叔的clusterProfiler砂豌，我論文中的富集功能很多都是用這個包做的（還有的用了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
