參考這篇
ID轉(zhuǎn)換不用怕，clusterProfiler幫你忙
這篇詳細(xì)講了無參考基因組該怎么辦豹休，值得一看

簡(jiǎn)單介紹一下幾種常用的ID：
Ensemble id：一般以ENSG開頭炊昆，后邊跟11位數(shù)字。如TP53基因：ENSG00000141510
Entrez id：通常為純數(shù)字威根。如TP53基因：7157
Symbol id：為我們常在文獻(xiàn)中報(bào)道的基因名稱凤巨。如TP53基因的symbol id為TP53
Refseq id：NCBI提供的參考序列數(shù)據(jù)庫：可以是NG、NM洛搀、NP開頭敢茁，代表基因，轉(zhuǎn)錄本和蛋白質(zhì)留美。如TP53基因的某個(gè)轉(zhuǎn)錄本信息可為NM_000546

1 轉(zhuǎn)換ID

1.1 載入包library(clusterProfiler)

library(clusterProfiler)
library(DOSE)
library(org.Hs.eg.db)

org.Hs.eg.db是人類基因組注釋彰檬，如果需要其他物種的注釋信息，可以去bioconductor這里獲得

1.2. keytypes()查看注釋包中支持的ID轉(zhuǎn)換類型谎砾，基本包括了所以常用的ID類型

keytypes()

1.3. bitr()函數(shù)轉(zhuǎn)換ID

函數(shù)全部?jī)?nèi)容如下

bitr(geneID, fromType, toType, OrgDb, drop = TRUE

geneID：輸入待轉(zhuǎn)換的geneID
fromType：輸入的ID類型
toType：希望輸出的ID類型
OrgDb：注釋對(duì)象的信息
drop：drop = FALSE 保留空值

以TP53為例僧叉，希望輸出'ENTREZID','ENSEMBL','REFSEQ'

my_id <- c("TP53")

output <- bitr(my_id,
               fromType = 'SYMBOL',
               toType = c('ENTREZID','ENSEMBL','REFSEQ'),
               OrgDb = 'org.Hs.eg.db')

轉(zhuǎn)換后的結(jié)果

1.4. bitr_kegg()函數(shù)進(jìn)行基因ID與蛋白質(zhì)ID的轉(zhuǎn)換

函數(shù)全部?jī)?nèi)容如下，與bitr()類似

bitr_kegg(geneID, fromType, toType, organism, drop = TRUE)

??????在參數(shù)這里與bitr()有些區(qū)別棺榔！
geneID：輸入待轉(zhuǎn)換的geneID
fromType：輸入的ID類型只能是kegg（與entrez id相同）, ncbi-geneid, ncbi-proteinid or uniprot中的一個(gè)
toType：希望輸出的ID類型瓶堕，也只能是kegg（與entrez id相同）, ncbi-geneid, ncbi-proteinid or uniprot中的一個(gè)
organism：hsa，代表人類症歇，其他的物種名稱可以參考kegg的網(wǎng)站
drop：去除空值與否

TP53的entrez id為7157

kegg <- bitr_kegg(7157,
                   fromType = 'kegg',
                   toType = c('ncbi-geneid', 'ncbi-proteinid', 'uniprot'),
                   organism = 'hsa')

報(bào)錯(cuò)了郎笆，所以一次應(yīng)該只能進(jìn)行一種轉(zhuǎn)換谭梗！

改成這樣↓

kegg <- bitr_kegg(7157,
                  fromType = 'kegg',
                  toType = 'ncbi-proteinid',
                  organism = 'hsa')

轉(zhuǎn)換為uniprot

kegg <- bitr_kegg(7157,
                  fromType = 'kegg',
                  toType = 'uniprot',
                  organism = 'hsa')

出現(xiàn)了三個(gè)uniprot

P04637

K7PPA8顯示的是轉(zhuǎn)錄本水平的表達(dá)，級(jí)別是Unreviewed凄吏，就是其來源為UniProtKB/TrEMBL

K7PPA8

Q53GA5顯示的是轉(zhuǎn)錄本水平的表達(dá)远舅，級(jí)別是Unreviewed，就是其來源為UniProtKB/TrEMBL

Q53GA5

一般ID轉(zhuǎn)換僅僅為開始的準(zhǔn)備工作痕钢，將自己的數(shù)劇轉(zhuǎn)換好之后可以進(jìn)行后續(xù)的分析图柏。
另外，利用clusterProfiler包可以進(jìn)行許多豐富的下游分析任连，比如GO分析蚤吹、KEGG分析等等

2 富集分析

參考這篇
最常用的基因注釋工具是GO和KEGG注釋，這基本上是和差異基因分析一樣随抠，必做的事裁着，GO可在BP(生物過程），MF(分子功能）拱她，CC(細(xì)胞組分）三方面進(jìn)行注釋

GO中的注釋信息

KEGG數(shù)據(jù)庫全部的物種及其簡(jiǎn)寫（三個(gè)字母）

有很多在線（DAVID二驰、Gather、GOrilla,revigo）或者客戶端版的工具（IPA（IPA不是用的GO和KEGG數(shù)據(jù)庫）和FUNRICH）可以用來分析差異基因秉沼，我主要實(shí)踐一下clusterProfiler

??clusterProfiler提供有兩種富集方式??：

① ORA（Over-Representation Analysis）差異基因富集分析（不需要表達(dá)值诸蚕，只需要gene name）

② FCS，它的代表就是GSEA氧猬，即基因集(gene set)富集分析（不管有無差異背犯，需要全部genes表達(dá)值）

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

1) 先讀進(jìn)來差異基因

sig.gene <- read.delim("你的路徑/final_result.txt（之前計(jì)算出的差異基因文件）",  sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)
gene<-sig.gene[,1]
head(gene)

#轉(zhuǎn)換ID
gene.df<-bitr(gene, fromType = "ENSEMBL", 
              toType = c("SYMBOL","ENTREZID"),
              OrgDb = org.Hs.eg.db,
              drop = FALSE)

head(gene.df)

2) GO富集

ego_cc<-enrichGO(gene = gene.df$ENSEMBL,
                 OrgDb  = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont  = "CC",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)

ego_bp<-enrichGO(gene = gene.df$ENSEMBL,
                 OrgDb = org.Hs.eg.db,
                 keyType = 'ENSEMBL',
                 ont  = "BP",
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.05)

參數(shù)意義：
keyType 指定基因ID的類型盅抚，默認(rèn)為ENTREZID
OrgDb 指定該物種對(duì)應(yīng)的org包的名字
ont 代表GO的3大類別漠魏，BP, CC, MF，也可是全部ALL

pAdjustMethod 指定多重假設(shè)檢驗(yàn)矯正的方法妄均，有“ holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”中的一種
cufoff 指定對(duì)應(yīng)的閾值
readable=TRUE 代表將基因ID轉(zhuǎn)換為gene symbol柱锹。

barplot() #富集柱形圖
dotplot() #富集氣泡圖
cnetplot() #網(wǎng)絡(luò)圖展示富集功能和基因的包含關(guān)系
emapplot() #網(wǎng)絡(luò)圖展示各富集功能之間共有基因關(guān)系
heatplot() #熱圖展示富集功能和基因的包含關(guān)系

3) barplot()作圖顯示，依賴于ggplot2包

library(ggplot2)
barplot(ego_bp,showCategory = 18,title="The GO_BP enrichment analysis of all DEGs ")  

ggplot2參數(shù)意義具體可以看這篇

barplot()

dotplot()

emapplot()

4) 或者是KEGG分析

輸入的gene要為vector格式的entrezid丰包，可以用is.vector()判斷是否是vector格式

library(stringr)

kk<-enrichKEGG(gene = gene.df$ENTREZID,
               organism = "hsa",
               keyType = "kegg",          
               pAdjustMethod = "BH",
               pvalueCutoff = 0.01)

KEGG

library(DOSE)

data(geneList, package = "DOSE")
gene <- names(geneList)[abs(geneList) > 2]

kk<-enrichKEGG(gene = gene,
               organism = "hsa",
               keyType = "kegg",
               pAdjustMethod = "BH",
               pvalueCutoff = 0.01)

就可以了瞧毙，繼續(xù)往下進(jìn)行

#bar圖
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))
#點(diǎn)圖
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))

barplot()

dotplot()

cnetplot()

heatplot()

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

這篇為網(wǎng)頁版操作
好處：可以發(fā)現(xiàn)被差異基因舍棄的genes可能參與了某重要生理過程或信號(hào)通路
要點(diǎn)：用所有差異基因來跑GSEA，而不是按篩選條件篩選出來的宙彪，否則GSEA就失去了意義

1）讀入全部差異基因數(shù)據(jù)

library(dplyr)
all_deseq <- read.delim("你的路徑/DESeq2-res.txt",  sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

部分內(nèi)容

2） 提取出gene_id 和 log2FoldChange列

ensemblid_log2fdc <-dplyr::select(all_deseq, gene_id, log2FoldChange)

3) 可以看到這時(shí)gene_id還不是整數(shù)形式矩动，所以對(duì)這列gsub()取整，命名為symbol列添加到ensemblid_log2fdc中

ensemblid_log2fdc$symbol <- gsub("\\.\\d*", "", all_deseq[,2])

4）轉(zhuǎn)換id：將ENSEMBLID轉(zhuǎn)為ENTREZID

entrezid <-bitr(ensemblid_log2fdc$symbol, 
              fromType = "ENSEMBL", 
              toType = c("ENTREZID"),
              OrgDb = org.Hs.eg.db,
              drop = TRUE)

5）因?yàn)?）中merge需要相同列名释漆，所以將symbol列列名改為ENSEMBL

names(ensemblid_log2fdc) <- c("gene_id", "log2FoldChange", "ENSEMBL")

6）ensemblid_log2fdc 和 entrezid拼接到一起悲没，依靠的是它們具有相同列名的列——ENSEMBL

final <- merge(ensemblid_log2fdc, entrezid)

7) 提取entrezid 和log2fdc列

library(dplyr)    
entrezid_log2fdc <- dplyr::select(final, ENTREZID, log2FoldChange)

8) 用arrange()按log2fdc排序，desc()表示降序

final.sort <- arrange(entrezid_log2fdc, desc(log2FoldChange))

9) 構(gòu)建一下要GSEA分析的數(shù)據(jù)

genelist <- final.sort$log2FoldChange
names(genelist) <- final.sort[,1]

10) 進(jìn)行GSEA分析男图，需要entrezid示姿，因此我前面才會(huì)進(jìn)行ID轉(zhuǎn)化

gsemf <- gseGO(genelist,
           OrgDb = org.Hs.eg.db,
           keyType = "ENTREZID",
           ont="BP",
           pvalueCutoff = 1)
head(gsemf)

11）畫出GSEA圖

gseaplot(gsemf, 1)

GSEA圖

3 可視化

1.GO富集分析結(jié)果可視化

barplot

barplot(ego, showCategory = 10) #默認(rèn)展示顯著富集的top10個(gè)，即p.adjust最小的10個(gè)

dotplot

dotplot(ego, showCategory = 10)

DAG有向無環(huán)圖

plotGOgraph(ego) #矩形代表富集到的top10個(gè)GO terms, 顏色從黃色過濾到紅色逊笆，對(duì)應(yīng)p值從大到小栈戳。

igraph布局的DAG

goplot(ego)

GO terms關(guān)系網(wǎng)絡(luò)圖（通過差異基因關(guān)聯(lián)）

emapplot(ego, showCategory = 30)

GO term與差異基因關(guān)系網(wǎng)絡(luò)圖

cnetplot(ego, showCategory = 5)

2.Pathway富集分析結(jié)果可視化

barplot

barplot(kk, showCategory = 10)

dotplot

dotplot(kk, showCategory = 10)

pathway關(guān)系網(wǎng)絡(luò)圖（通過差異基因關(guān)聯(lián)）

emapplot(kk, showCategory = 30)

pathway與差異基因關(guān)系網(wǎng)絡(luò)圖

cnetplot(kk, showCategory = 5)

pathway映射

browseKEGG(kk, "hsa04934") #在pathway通路圖上標(biāo)記富集到的基因，會(huì)鏈接到KEGG官網(wǎng)