1. 非模式生物
非模式生物需要提前從AnnotationHub抓取目標物種的OrgDb數(shù)據(jù)庫到逊，這個庫里包括Uniport滤钱，ENSEMBL，ENTREZID件缸，REFSEQ他炊，GO...争剿。擁有它痊末，你就可以實現(xiàn)GENENAME凿叠，ENSEMBL，ENTREZID等多種格式的自由轉換盒件。不過這個功能一些在線網站也能做到，如果只有少量數(shù)據(jù)轉換锰茉，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() 顯示括號內對象擁有的數(shù)據(jù)名稱, 每一類數(shù)據(jù)為key
- 2 keys() 返回當前對象的keys
- 2 keytypes() 字義理解keys的種類全谤，用在select和mapIds中作參數(shù)轉化
- 3 select() 輸入一組key，對應其column补憾，可將其轉化為目標keytype。輸出的數(shù)據(jù)框同時包含轉化前后數(shù)據(jù)腾务。
- 4 mapId() 輸入輸入一組key削饵，對應其keytype，可將其轉化為目標column中的key启昧。輸出的數(shù)據(jù)只轉化后數(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包田炭，安裝加載目標物種的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ù)惭蹂，也是轉換用的，不過看名字就知道盾碗，這個函數(shù)服務于kegg廷雅。bitr_kegg中的轉化參數(shù)有“Path”, “Module”, “ncbi-proteinid”, “ncbi-geneid”, “uniprot”, “kegg”京髓。提供其中一項可以互相轉化商架，舉個栗子??:
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庫中對應編碼hsa, 如果你不知道目標物種對應的縮寫organism=?，請點擊kegg_organism_list诚些。而'kegg'和‘Path‘都是取自KEGG PATHWAY數(shù)據(jù)庫皇型，kegg表示KEGG PATHWAY數(shù)據(jù)庫中基因命名方式砸烦，大多數(shù)采用ENTREZID幢痘，也用采用特殊方式，比如擬南芥采用'TAIR'颜说。如果你不確定目標物種的kegg命名门粪，可以使用ncbi-geneid(ENTREZID)轉化為kegg。而Path則為目標基因對應的kegg pathway的編號玄妈。
在使用bitr_kegg時常會報錯拟蜻，比如
如果你確定物種縮寫沒有錯，keyType使用正確诡必，那就不是你的問題搔扁，可能是網絡的問題，多跑幾遍???♂?碳抄。
差異基因功能GO分析
GO/kegg分析的在線工具很多场绿，比如gProfiler。但是如果需要批量處理幾組差異基因分析璧尸，在線工具還是挺費事的。你沒時間研究R的話爷光，用網頁點點點也是錯不了蛀序。
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 ...
#隨后對富集結果進行總覽崇堵，查看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',])
#保存結果
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值設置更為嚴格閾值幔摸，有兩種方式可以直接篩選: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服務的數(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咧七，讓結果更加簡潔。
##去除冗余terms
go.filter<-simplify(enrich.BP,cutoff = 0.7,
by = "p.adjust",
select_fun = min)
dim(go.filter)
#將結果數(shù)據(jù)框中的ENTREZID轉化為SYMBOL
y <-setReadable(go.all,OrgDb = org.Hs.eg.db,keyType="ENTREZID")
head(y)
接下來是KEGG，KEGG需要注意的是organism废酷，keyType澈蟆，use_internal_data設置。
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 : 是否使用內部數(shù)據(jù)，當use_internal_data =FALSE時彩掐，爬取在線KEGG PATHWAY數(shù)據(jù)庫灰追，特點就是實時更新狗超，而內部數(shù)據(jù)就有滯后性了谐檀。
需要注意pvalue或p.adjust設置
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
問題一：橫坐標桐猬，氣泡溢出
p1<-kegg.dot + xlim(NA,0.085) #數(shù)值>橫坐標顯示值
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
問題六：橫坐標不是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
一本關于clusterProfiler的小冊子 by Y叔
