1. 非模式生物

非模式生物需要提前從AnnotationHub抓取目標(biāo)物種的OrgDb數(shù)據(jù)庫，這個庫里包括Uniport，ENSEMBL，ENTREZID，REFSEQ，GO...。擁有它，你就可以實現(xiàn)GENENAME，ENSEMBL，ENTREZID等多種格式的自由轉(zhuǎn)換。不過這個功能一些在線網(wǎng)站也能做到，如果只有少量數(shù)據(jù)轉(zhuǎn)換，DAVID/gProfiler還是挺香的。下面以油菜(Brassica_napus)為例介紹OrgDb：

#R script
setwd("C:\\Users\\Desktop\\Transcriptome\\class\\class3")

####非模式生物通過AnnotationHub在線檢索并抓取OrgDb
BiocManager::install("AnnotationHub")
library(AnnotationHub)
hub<-AnnotationHub::AnnotationHub()
query(hub,"Brassica") #抓取包含Brassica字符串的物種OrgDb包，并獲取其編號
bna <- hub[["AH75769"]] 
#玉米query(hub, "zea");maize <- hub[['AH55736']]

若你想深入了解OrgDb數(shù)據(jù)庫，記住5個函數(shù)function():

• 1 columns() 顯示括號內(nèi)對象擁有的數(shù)據(jù)名稱, 每一類數(shù)據(jù)為key
• 2 keys() 返回當(dāng)前對象的keys
• 2 keytypes() 字義理解keys的種類，用在select和mapIds中作參數(shù)轉(zhuǎn)化
• 3 select() 輸入一組key，對應(yīng)其column，可將其轉(zhuǎn)化為目標(biāo)keytype。輸出的數(shù)據(jù)框同時包含轉(zhuǎn)化前后數(shù)據(jù)。
• 4 mapId() 輸入輸入一組key，對應(yīng)其keytype，可將其轉(zhuǎn)化為目標(biāo)column中的key。輸出的數(shù)據(jù)只轉(zhuǎn)化后數(shù)據(jù)(即column)。
  這樣的文字介紹還是抽象，不好理解。上機操作：

columns(bna)
# [1] "ACCNUM"  "ALIAS"  "CHR"  "ENTREZID" "EVIDENCE"   
# [6] "EVIDENCEALL" "GENENAME"   "GID"  "GO"  "GOALL"      
# [11] "ONTOLOGY"    "ONTOLOGYALL" "PMID"  "REFSEQ" "SYMBOL"     
# [16] "UNIGENE" 
keytypes(bna)
# [1] "ACCNUM"      "ALIAS"       "ENTREZID"    "EVIDENCE"    "EVIDENCEALL"
# [6] "GENENAME"    "GID"         "GO"          "GOALL"       "ONTOLOGY"   
# [11] "ONTOLOGYALL" "PMID"        "REFSEQ"      "SYMBOL"      "UNIGENE" 
head(keys(bna,keytype="SYMBOL"))
#[1] "BNAA01G00220D" "BNAA01G00960D" "BNAA01G01590D" "BNAA01G01860D"
#[5] "BNAA01G01890D" "BNAA01G02080D"
genename<-select(bna,keys=c("BNAA01G00220D","BNAA01G02080D"),column=c("GENENAME"),keytype = "SYMBOL")
head(genename)
#         SYMBOL                      GENENAME
# 1 BNAA01G00220D uncharacterized BNAA01G00220D
# 2 BNAA01G02080D uncharacterized BNAA01G02080D
entrezid<-mapIds(bna,keys=c("BNAA01G00220D","BNAA01G02080D"), column = c("ENTREZID"), keytype = "SYMBOL")
head(entrezid)
# BNAA01G00220D BNAA01G02080D 
# "106357927"   "106370211" 

模式動物基因功能KEGG/GO注釋(ClusterProfiler)

模式動物有現(xiàn)成的OrgDb包，安裝加載目標(biāo)物種的R package就可以了。

####模式動物(人為例)使用clusterProfiler進行GO/KEGG分析
#Bioconductor 安裝 clusterProfiler
BiocManager::install("clusterProfiler")
BiocManager::install("org.Hs.eg.db")
install.packages("org.Hs.eg.db")
#加載clusterProfiler
library(clusterProfiler)
library(ggplot2)
library(org.Hs.eg.db)

hsa<-org.Hs.eg.db
keytypes(hsa)
head(keys(hsa, keytype="SYMBOL"))

#讀入基因數(shù)據(jù)
gene<-read.csv("DEGs_identified.csv")
head(gene)
# Ensembl        padj log2FoldChange
# 1 ENSG00000001084 1.06000e-32       1.371851
# 2 ENSG00000002587 5.12000e-57       2.516021
# 3 ENSG00000003096 9.67315e-04      -1.088586
# 4 ENSG00000003137 6.01000e-18       1.170100
# 5 ENSG00000003400 8.16000e-07      -1.092483
# 6 ENSG00000004487 1.04000e-24       1.061878
dim(gene)
# [1] 934   3

#篩選|FC|>=2且FDR<=0.05的為顯著差異基因
gene_ensembl<-na.omit(gene[(abs(gene$log2FoldChange)>= 1)&(gene$padj<= 0.05),1])
head(gene_ensembl);length(gene_ensembl)
#[1] ENSG00000001084 ENSG00000002587 ENSG00000003096 ENSG00000003137 ENSG00000003400 ENSG00000004487
# 934 Levels: ENSG00000001084 ENSG00000002587 ENSG00000003096 ENSG00000003137 ENSG00000003400 ... ENSG00000276600
# [1] 934
#通過ENSEMBL，獲取ENTREZID，GENENAME
genelist<- bitr(gene_ensembl,fromType="ENSEMBL",toType=c("ENTREZID","GENENAME"),OrgDb=hsa)
#'select()' returned 1:many mapping between keys and columns
#Warning message:
#In bitr(gene_ensembl, fromType = "ENSEMBL", toType = c("ENTREZID",  :
#0.86% of input gene IDs are fail to map...

值得一提的是，除了bitr()還有個bitr_kegg()函數(shù)，也是轉(zhuǎn)換用的，不過看名字就知道，這個函數(shù)服務(wù)于kegg。bitr_kegg中的轉(zhuǎn)化參數(shù)有“Path”, “Module”, “ncbi-proteinid”, “ncbi-geneid”, “uniprot”, “kegg”。提供其中一項可以互相轉(zhuǎn)化，舉個栗子??:

kegg<-bitr_kegg(genelist$ENTREZID,fromType="ncbi-geneid",toType="kegg",organism='hsa')
head(kegg)
#  ncbi-geneid      kegg
# 1   100131827 100131827
# 2   100132074 100132074

Path<-bitr_kegg(genelist$ENTREZID,fromType="ncbi-geneid",toType="Path",organism='hsa')
head(Path)
#

人類KEGG庫中對應(yīng)編碼hsa, 如果你不知道目標(biāo)物種對應(yīng)的縮寫organism=?，請點擊kegg_organism_list。而'kegg'和‘Path‘都是取自KEGG PATHWAY數(shù)據(jù)庫，kegg表示KEGG PATHWAY數(shù)據(jù)庫中基因命名方式，大多數(shù)采用ENTREZID，也用采用特殊方式，比如擬南芥采用'TAIR'。如果你不確定目標(biāo)物種的kegg命名，可以使用ncbi-geneid(ENTREZID)轉(zhuǎn)化為kegg。而Path則為目標(biāo)基因?qū)?yīng)的kegg pathway的編號。

在使用bitr_kegg時常會報錯，比如

如果你確定物種縮寫沒有錯，keyType使用正確，那就不是你的問題，可能是網(wǎng)絡(luò)的問題，多跑幾遍???♂?。

差異基因功能GO分析

GO/kegg分析的在線工具很多，比如gProfiler。但是如果需要批量處理幾組差異基因分析，在線工具還是挺費事的。你沒時間研究R的話，用網(wǎng)頁點點點也是錯不了。

GO(GENE ONTOLOGY)對功能的注釋分為三大類:分別是BP，MF，和CC。通過對ont的選擇(ont='BP','MF','CC')就能分門別類地分析，但如果你想一次解決，直接令ont="ALL"。

下面進入ClusterProfiler的高光時刻，走起：

go.all<- enrichGO(gene=genelist$ENTREZID, 
                  OrgDb = org.Hs.eg.db, 
                  ont='ALL', #ont='BP'
                  pAdjustMethod = 'BH',
                  pvalueCutoff = 0.05, 
                  qvalueCutoff = 0.2,
                  keyType = 'ENTREZID')
go.all
# over-representation test
#...@organism    Homo sapiens 
#...@ontology    GOALL 
#...@keytype     ENTREZID 
#...@gene        chr [1:929] "2729" "9957" "90293" "56603" "843" "23028" "115703" "26224" "114881" "51330" "80853" "55200" ...
#...pvalues adjusted by 'BH' with cutoff <0.05 
#...627 enriched terms found
#'data.frame':   627 obs. of  10 variables:
# $ ONTOLOGY   : Factor w/ 2 levels "BP","MF": 1 1 1 1 1 1 1 1 1 1 ...
# $ ID         : chr  "GO:0001763" "GO:0061138" "GO:0045766" "GO:0048754" ...
# $ Description: chr  "morphogenesis of a branching structure" "morphogenesis of a branching epithelium" "positive regulation of angiogenesis" "branching morphogenesis of an epithelial tube" ...
# $ GeneRatio  : chr  "33/832" "30/832" "31/832" "26/832" ...
# $ BgRatio    : chr  "196/18670" "182/18670" "204/18670" "150/18670" ...
# $ pvalue     : num  4.02e-11 5.27e-10 2.13e-09 2.49e-09 2.89e-09 ...
# $ p.adjust   : num  2.10e-07 1.38e-06 3.03e-06 3.03e-06 3.03e-06 ...
# $ qvalue     : num  1.62e-07 1.06e-06 2.33e-06 2.33e-06 2.33e-06 ...
# $ geneID     : chr  "54903/6591/639/59277/2294/2138/374/7422/7475/1746/5228/53335/157506/6662/650/54567/10253/2247/1969/9353/170690/"| __truncated__ "54903/6591/59277/2294/2138/374/7422/7475/5228/157506/6662/650/54567/10253/2247/1969/9353/170690/133/10481/57216"| __truncated__ "9734/5743/23411/3162/3696/5054/7422/3552/2152/5228/3553/79924/694/2113/9314/9982/7057/1545/2247/133/7424/51738/"| __truncated__ "54903/2294/2138/374/7422/7475/5228/157506/6662/650/54567/10253/2247/1969/9353/170690/57216/64783/4824/2637/2303"| __truncated__ ...
# $ Count      : int  33 30 31 26 33 44 48 47 44 21 ...

#隨后對富集結(jié)果進行總覽，查看BP，CC，MF的個數(shù)
dim(go.all[go.all$ONTOLOGY=='BP',]);dim(go.all[go.all$ONTOLOGY=='CC',]);dim(go.all[go.all$ONTOLOGY=='MF',])
#保存結(jié)果
write.csv(go.all@result,'DEG_go.all.result.csv',row.names=F)

用enrichplot畫個圖看看：

dotplot(go.all,showCategory = 10, size = NULL,font.size = 10, title = "GO enrichment", split = "ONTOLOGY") + facet_grid(ONTOLOGY ~ ., scales = "free")
barplot(go.all,showCategory = 10, size = NULL,font.size = 10, title = "GO enrichment", split = "ONTOLOGY") + facet_grid(ONTOLOGY ~ ., scales = "free")

如果要將ont='BP'篩選出來，并且對P值設(shè)置更為嚴格閾值，有兩種方式可以直接篩選:1.直接篩選；2.利用clusterProfiler.dplyr。

##1. 直接篩選，MF沒有富集到terms，不作演示
enrich.BP<-go.all[go.all@result$pvalue<0.01&go.all$ONTOLOGY=='BP',]
enrich.CC<-go.all[go.all@result$pvalue<0.01&go.all$ONTOLOGY=='CC',]
dim(enrich.BP);dim(enrich.CC)

還可以安裝個為cluster profiler服務(wù)的數(shù)據(jù)處理包clusterProfiler.dplyr。然后一條命令就解決了。

install.packages("devtools")
devtools::install_github("YuLab-SMU/clusterProfiler.dplyr")
BiocManager::install('clusterProfiler.dplyr')
library(clusterProfiler.dplyr)
enrich.BP<-filter(go.all,ONTOLOGY=='BP',p.adjust <.01, qvalue < 0.1)
enrich.BP@ontology<-"BP"
dim(enrich.BP)
BP.bar<-barplot(enrich.BP,showCategory=20);BP.dot<-dotplot(enrich.BP,showCategory=20)
BP.bar
BP.dot

通常得到的GOterms中有一些相似度較高的term，我們可以去除這些冗余terms，讓結(jié)果更加簡潔。

##去除冗余terms
go.filter<-simplify(enrich.BP,cutoff = 0.7,
                    by = "p.adjust",
                    select_fun = min)
dim(go.filter)
#將結(jié)果數(shù)據(jù)框中的ENTREZID轉(zhuǎn)化為SYMBOL
y <-setReadable(go.all,OrgDb = org.Hs.eg.db,keyType="ENTREZID")
head(y)

接下來是KEGG，KEGG需要注意的是organism，keyType，use_internal_data設(shè)置。

enrich.kegg <- enrichKEGG(gene =genelist$ENTREZID,
                          organism ="hsa",
                          keyType = "kegg",
                          pvalueCutoff = 1,
                          pAdjustMethod = "BH",
                          minGSSize = 10,
                          maxGSSize = 500,
                          qvalueCutoff = 1,
                          use_internal_data =FALSE)
dim(enrich.kegg)

minGSSize：minimal size of genes annotated for testing，用于測試基因集的最小數(shù)目。
maxGSSize：maximal size of genes annotated for testing，用于測試的基因注釋最大數(shù)目。
這兩個參數(shù)常常被忽略。
use_internal_data : 是否使用內(nèi)部數(shù)據(jù)，當(dāng)use_internal_data =FALSE時，爬取在線KEGG PATHWAY數(shù)據(jù)庫，特點就是實時更新，而內(nèi)部數(shù)據(jù)就有滯后性了。

需要注意pvalue或p.adjust設(shè)置

sig.kegg<-filter(enrich.kegg,pvalue<.05,qvalue<0.2)
dim(sig.kegg)
kegg.bar<-barplot(sig.kegg,showCategory=20,color = "pvalue")
kegg.dot<-dotplot(sig.kegg,showCategory=20,color = "pvalue")
db<-plot_grid(kegg.bar,kegg.dot,ncol=2)
db
ggsave("kegg_bar_dot1.png",db,width=14,height=6)

kegg_bar_dot1.png

問題一：橫坐標(biāo)，氣泡溢出

p1<-kegg.dot + xlim(NA,0.085) #數(shù)值>橫坐標(biāo)顯示值
p1
db1<-plot_grid(kegg.bar,p1,ncol=2)
db1
ggsave("kegg_bar_dot2.png",db1,width=14,height=6)

kegg_bar_dot2.png

問題二：注釋文字太長，截斷成兩行

library(stringr)
p2<- kegg.dot + scale_y_discrete(labels=function(y) stringr::str_wrap(y,width=36))
p2

#####問題三：改變顏色
p30 <-p2 + scale_color_continuous(low='purple',high='green')+ xlim(NA,0.084)
p30
p31 <-p2 + scale_color_gradientn(colours = rainbow(5))+ xlim(NA,0.084)
p31
p32 <-p2+scale_color_gradient(low="yellow", high="green")+ xlim(NA,0.084)
p32
db<-plot_grid(p30,p31,p32,ncol=3)
db
ggsave("kegg_bar_dot3.png",db,width=16,height=6)

kegg_bar_dot3.png

問題四：改變形狀

p3 <-p31 + aes(shape=GeneRatio > 0.05) 
p3

問題五：氣泡變大

p6<- p31 + scale_size(range=c(2,12))
p6

問題六：橫坐標(biāo)不是GeneRatio,而是richFactor或者FoldEnrichment

y1 <- mutate(sig.kegg, richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))
#y2 <- mutate(sig.kegg, FoldEnrichment = parse_ratio(GeneRatio) / parse_ratio(BgRatio))
y1 

#若橫軸是FoldEnrichment將y1替換成y2,
library(enrichplot)
library(forcats)
rf<-ggplot(y1, showCategory = 20, 
       aes(richFactor, fct_reorder(Description,richFactor))) + 
    geom_point(aes(color=pvalue, size = Count)) +
    scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
    scale_size_continuous(range=c(2,10)) +
    #theme_minimal() + 
    xlab("Rich Factor") +
    ylab(NULL) + 
    ggtitle("Top 20 KEGG Enrichment")+
    theme_bw()
ggsave("kegg_bar_dot4.png",rf,width=7,height=5)

kegg_bar_dot4.png

如果你還想了解更多，點擊下方鏈接，直接進入clusterProfiler官方解說
clusterProfiler 官方document
一本關(guān)于clusterProfiler的小冊子 by Y叔