|
作者：MC學(xué)公衛(wèi) 來(lái)源：百味科研芝士 GeneOntology富集分析是高通量數(shù)據(jù)分析的標(biāo)配，不管是轉(zhuǎn)錄組、甲基化、ChIP-seq還是重測(cè)序，都會(huì)用到對(duì)一個(gè)或多個(gè)集合的基因進(jìn)行功能富集分析。分析結(jié)果可以指示這個(gè)集合的基因具有什么樣的功能偏好性，進(jìn)而據(jù)此判斷相應(yīng)的生物學(xué)意義。 GO富集分析是先篩選差異基因，再判斷差異基因在哪些注釋的通路存在富集。
1. 安裝r包c(diǎn)lusterProfiler 在本地：
> source('http:///biocLite.R')
> biocLite(\\'clusterProfiler\\')
2.基因列表Entrez ID與FC值的準(zhǔn)備 利用一些差異表達(dá)的基因作為輸入的基因，之前Seurat包的 FindAllMarkers有一批這些基因，但是沒(méi)有EntrezID，需要進(jìn)行一些處理。總的來(lái)說(shuō)，就是需要整理出需要作富集分析的基因列表 先查看一下：
spleen.markers <- FindAllMarkers(object = spleen, only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25)
print(x= head(x=spleen.markers,n = 10))
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
MS4A1 2.410289e-189 1.592266 0.974 0.258 3.773307e-185 0 MS4A1
CD79A 1.245084e-177 1.485304 0.971 0.270 1.949180e-173 0 CD79A
HLA-DRA 1.792014e-151 1.370242 1.000 0.556 2.805397e-147 0 HLA-DRA
CD79B 2.280181e-149 1.265725 0.938 0.315 3.569624e-145 0 CD79B
CD74 3.217361e-144 1.256090 0.997 0.871 5.036778e-140 0 CD74
HLA-DPB1 4.238869e-137 1.154889 0.997 0.616 6.635949e-133 0 HLA-DPB1
HLA-DRB5 6.911780e-136 1.123850 0.990 0.437 1.082039e-131 0 HLA-DRB5
HLA-DQB1 4.004476e-134 1.121013 0.940 0.366 6.269007e-130 0 HLA-DQB1
HLA-DRB1 2.026996e-133 1.115709 0.992 0.452 3.173262e-129 0 HLA-DRB1
HLA-DPA1 3.492170e-130 1.099676 0.995 0.570 5.466992e-126 0 HLA-DPA1
spleen.markers一共有3,108個(gè)基因marker。 獲取基因列表 genelist <- spleen.markers$gene
head(genelist)
[1] 'MS4A1' 'CD79A' 'HLA-DRA' 'CD79B' 'CD74' 'HLA-DPB1'
加載包，下載人類(lèi)基因組注釋的數(shù)據(jù)org.Hs.eg.db library(clusterProfiler)
BiocInstaller::biocLite('org.Hs.eg.db')
library(org.Hs.eg.db)
安裝與細(xì)節(jié)可在此網(wǎng)站查看Genome wide annotation for Human（含引用文獻(xiàn)）： 轉(zhuǎn)化為EntrezID s.EntrezID <- bitr(genelist, fromType='SYMBOL', toType='ENTREZID', OrgDb='org.Hs.eg.db')
head(s.EntrezID)
SYMBOL ENTREZID
1 MS4A1 931
2 CD79A 973
3 HLA-DRA 3122
4 CD79B 974
5 CD74 972
6 HLA-DPB1 3115
查看行數(shù)： [1] 2043 2
原本是3,108個(gè)，現(xiàn)在剩下2043個(gè)基因，原因是去除了重復(fù)的個(gè)數(shù)。  接下來(lái)進(jìn)行富集分析： s.ego <- enrichGO(gene = s.EntrezID.df$ENTREZID, OrgDb = \\'org.Hs.eg.db\\', ont = \\'ALL\\', pAdjustMethod = \\'BH\\',pvalueCutoff = 0.05, qvalueCutoff = 0.2, keyType = \\'ENTREZID\\')
head(s.ego)
篩選的p值為0.05，q值為0.2 ONTOLOGY ID Description GeneRatio
GO:0042119 BP GO:0042119 neutrophil activation 211/1931
GO:0002283 BP GO:0002283 neutrophil activation involved in immune response 208/1931
GO:0043312 BP GO:0043312 neutrophil degranulation 207/1931
GO:0002446 BP GO:0002446 neutrophil mediated immunity 208/1931
GO:0006613 BP GO:0006613 cotranslational protein targeting to membrane 77/1931
GO:0045047 BP GO:0045047 protein targeting to ER 78/1931
BgRatio pvalue p.adjust qvalue
GO:0042119 496/17381 2.402001e-74 1.352326e-70 9.264138e-71
GO:0002283 486/17381 7.667076e-74 2.158282e-70 1.478535e-70
GO:0043312 485/17381 3.218981e-73 6.040954e-70 4.138367e-70
GO:0002446 499/17381 2.372136e-71 3.338782e-68 2.287239e-68
GO:0006613 96/17381 5.653664e-56 6.366026e-53 4.361058e-53
GO:0045047 100/17381 5.659096e-55 5.310118e-52 3.637706e-52
geneID
GO:0042119 HVCN1/EEF2/SELL/IMPDH2/PTPN6/CTSH/CCL5/B2M/HSPA8/PTPRC/RAB27A/STOM/SURF4/HSP90AA1/BIN2/HSPA1A/HSPA6/HSPA1B/ATP6V0A1/CTSZ/RNASET2/HSP90AB1/SERPINA1/FCN1/CST3/CD68/CFP/FGL2/S100A12/FPR1/LYZ/CD14/CD36/MNDA/CFD/FCGR2A/MGST1/TLR2/CD93/C5AR1/TYROBP/FCER1G/LGALS3/
………………
………………
Count
GO:0042119 211
GO:0002283 208
GO:0043312 207
GO:0002446 208
GO:0006613 77
GO:0045047 78
剩下1931個(gè)基因： dim(s.ego)
[1] 1884 10
> dim(s.ego[s.ego$ONTOLOGY==\\'BP\\',])
[1] 1507 10
> dim(s.ego[s.ego$ONTOLOGY==\\'CC\\',])
[1] 215 10
> dim(s.ego[s.ego$ONTOLOGY==\\'MF\\',])
[1] 162 10
看來(lái)大部分的基因富集在BP中。 這里科普一下BP是指Biological Process，CC是指Cellular Component，MF是指Molecular Function?？梢愿鶕?jù)自己需要各取所需。 可以作柱狀圖
s.barplot <- barplot(s.ego, showCategory = 20, drop = T)
s.barplot
s.barplot <- barplot(s.ego, showCategory = 25, drop = T)
s.barplot
可以作點(diǎn)圖，作用與柱狀圖類(lèi)似
dotplot(s.ego)
也可以作富集網(wǎng)絡(luò)
enrichMap(s.ego)
這里的GO較多，若要清晰展示排名前10的GO可加上參數(shù)n=10，即enrichMap(s.ego, n = 10)再進(jìn)行繪圖。
|
