引入clusterProfiler與注釋數據
library(clusterProfiler)
library(DOSE)
library(DO.db)
library(org.Mm.eg.db)
GO(gene ontology)分析
GO，Gene Ontology，是基因功能國際標准分類體系。它旨在建立一個適用於各種物種的，對基因和蛋白質功能進行限定和描述的，並能隨着研究不斷深入而更新的語言詞匯標准。GO分為分子功能（Molecular Function）、生物過程（Biological Process）、和細胞組成（Cellular Component）三個部分。
ID轉換
org.Mm.eg.db數據對象里面包含着各大主流數據庫的數據，一般人都比較熟悉的entrez ID 和ensembl 數據庫的ID。
#查看ID類型
keytypes(org.Mm.eg.db)
gene <- row.names(diff_gene_deseq2)
#使用select函數進行表格拼接，這是來自數據庫的函數
#keys是原始的ID，columns是轉換之后的ID，keytype是要指定的原始ID類型
tansid <- select(org.Mm.eg.db,keys = gene,columns = c("SYMBOL","ENTREZID","GENENAME"),keytype = "ENSEMBL")
可視化
#GO分析及畫圖
ego <- enrichGO(gene = tansid$ENSEMBL,
keyType = "ENSEMBL",
OrgDb = "org.Mm.eg.db",
ont = "BP"
)
# 氣泡圖
dotplot(ego, font.size=10)
#網絡圖,但是找不到函數
enrichMap(ego, vertex.label.cex=1.2, layout=igraph::layout.kamada.kawai)
#GO圖
plotGOgraph(ego)
橫軸為GeneRatio, 代表該GO term下的差異基因個數占差異基因總數的比例，縱軸為富集到的GO Terms的描述信息
有向無環。矩形代表富集到的top10個GO terms, 顏色從黃色過濾到紅色，對應p值從大到小。橢圓形代表未富集到前十的GOterms。
##GSEA
# 按照log2FC從大到小排序基因列表
genelist <- diff_gene_deseq2$log2FoldChange
names(genelist) <- rownames(diff_gene_deseq2)
genelist <- sort(genelist, decreasing = TRUE)
#gsea
gsemf <- gseGO(genelist,
OrgDb = org.Mm.eg.db,
keyType = "ENSEMBL",
ont="MF"
)
head(gsemf)
# 畫出GSEA圖
gseaplot(gsemf, geneSetID="GO:0000977")
GSEA的輸入是一個基因表達量矩陣，其中的樣本分成了A和B兩組，首先對所有基因進行排序，在之前的文章中也有提到排序的標准，這里簡單理解就是foldchange, 用來表示基因在兩組間表達量的變化趨勢。排序之后的基因列表其頂部可以看做是上調的差異基因，其底部是下調的差異基因。
GSEA分析的是一個基因集下的所有基因是否在這個排序列表的頂部或者底部富集，如果在頂部富集，我們可以說，從總體上看，該基因集是上調趨勢，反之，如果在底部富集，則是下調趨勢。
- RANK IN GENE LIST代表該基因在排序號的列表中的位置
- RANK METRIC SCORE對應排序，比如foldchange值
- RUNNIG ES代表累計的Enrichment score, CORE ENRICHMENT代表是否屬於核心基因，即對該基因集的Enerchment score做出了主要貢獻的基因。出現基因集中的基因就增加ES值，不出現就減少（上調）。
https://blog.csdn.net/weixin_43569478/article/details/83745105
https://cloud.tencent.com/developer/article/1426130
KEGG
# id轉換
x=bitr(rownames(diff_gene_deseq2),fromType = "ENSEMBL",toType = "ENTREZID", OrgDb = "org.Mm.eg.db")
# 獲取keggID
kegg <- x[,2]
# KEGG分析，小鼠mmu，http://www.genome.jp/kegg/catalog/org_list.html
ekk <- enrichKEGG(kegg, keyType = "kegg",organism = "mmu", pvalueCutoff = 0.05, pAdjustMethod = "BH", qvalueCutoff = 0.1)
head(summary(ekk))
# 氣泡圖
barplot(ekk)
cnetplot(ekk, showCategory = 5)
browseKEGG(ekk, "mmu04742")
# 將GO/KEGG結果轉換成CSV格式輸出
write.csv(as.data.frame(ekk),"KEGG-enrich.csv",row.names =F)
write.csv(as.data.frame(ego),"GO-enrich.csv",row.names =F)
- 橫軸為該pathway的差異基因個數，縱軸為富集到的pathway的描述信息
- 灰色的點代表基因，黃色的點代表富集到的pathways, pathways節點的大小對應富集到的基因個數。
- 富集到的差異基因會用紅色方框表示,但是沒有一個基因比對上去。
https://www.jianshu.com/p/380bc3508d72
