在人医学遗传学研究中，SNP 与疾病相关性一直是广泛的研究课题。很多人在刚开始接触课题，在没有前期研究基础指示的目的基因时，都会选择从公共数据库中寻找与疾病相关的基因或者SNP进行研究。本文小编带你学习如何从数据库中挖掘基因，并聚焦疾病相关的重要通路。

技术路线

1.疾病/复杂性状相关数据库

百度搜索可以获得很多人疾病相关数据库的使用说明，在此不再赘述。我们常用的人类疾病数据库是 DisGeNET（https://www.disgenet.org/search ），常用 GWAS 数据库是 EMBL-EBI 的 GWAS catalog（https://www.ebi.ac.uk/gwas/ ）。

2.疾病相关最全基因 list

通过数据库下载疾病相关基因列表；相关 SNP 也可以在 VEP 在线注释工具（http://asia.ensembl.org/Multi/Tools/VEP ）注释其所在的基因。将两部分基因合并，获得基本相关较为全面的基因列表。

3.富集分析----聚焦疾病/复杂性状相关通路

富集分析使用 R 语言的 clusterProfiler 程序包。即使不会 R 语言，不懂编程，一样可以完成分析。

安装 R
百度搜索 R，找到合适的下载源；也可直接点击链接https://cran.dcc.uchile.cl/，选择合适的版本下载。

#设置工作目录
运行 R 后，在《文件》菜单下选择《改变工作目录》；
将基因名替换好以后，直接复制下列代码到 R，回车运行即可，除此之外可以不做任何改动。如果已经安装 clusterProfiler 程序包，请从基因编号转换开始。

#安装 clusterProfiler 程序包，此种安装方法适合 R3.5.2 及以下版本，R3.6.0 以上版本请参考文后补充说明。
source("https://bioconductor.org/biocLite.R")
biocLite("clusterProfiler")

#安装 KEGG.db
biocLite("kegg.db")

#安装人 org.db 数据库
biocLite(“org.Hs.eg.db”)

#基因编号转换

#将基因名称列表复制给（或者任何你喜欢的文件名），如需分析自己特定的基因集，可替换括号内容，每个基因名称，用””,隔开。
library(clusterProfiler)
yh <- c("GPX3",  "GLRX",   "LBP",   "CRYAB", "DEFB1", "HCLS1",   "SOD2",   "HSPA2",
       "ORM1",  "IGFBP1", "PTHLH", "GPC3",  "IGFBP3","TOB1",    "MITF",   "NDRG1",       "NR1H4", "FGFR3", "PVR", "IL6", "PTPRM", "ERBB2",   "NID2",   "LAMB1",       "COMP",  "PLS3",   "MCAM",  "SPP1",  "LAMC1", "COL4A2", "COL4A1", "MYOC",       "ANXA4", "TFPI2",  "CST6",  "SLPI",  "TIMP2", "CPM",     "GGT1",   "NNMT",       "MAL",   "EEF1A2", "HGD",   "TCN2",  "CDA",   "PCCA",    "CRYM",   "PDXK",       "STC1",  "WARS",  "HMOX1", "FXYD2", "RBP4",   "SLC6A12", "KDELR3", "ITM2B")

#利用 cluterProfiler 内置的 bitr 函数进行基因编号转换，并将转换后的信息存储在 gene 文件中
gene <- bitr(yh, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb="org.Hs.eg.db")
head(gene)

##提取 gene 数据中的 ENTREZID 列，并赋值给 DE_list
DE_list <- gene$ENTREZID

#去除重复值
DE_list[duplicated(DE_list)]
integer(0)

#调用 org.Hs.eg.db，并查看文件的各列名称信息
library(org.Hs.eg.db)
columns(org.Hs.eg.db)

#GO_MF 富集,基于基因数目，如果使用的是个人电脑，配置不高，为防止程序卡死，建议 MF\CC\BP 单个来运行，生成的图片也逐个生成保存后再运行下一个。
MF <- enrichGO(gene          = DE_list, #差异基因 vector
                keyType       ="ENTREZID",
OrgDb         = org.Hs.eg.db, #对应的OrgDb
                ont           = "MF", #GO 分类名称，CC BP MF
                pAdjustMethod = "BH", #Pvalue 矫正方法
                pvalueCutoff  = 0.05, #Pvalue 阈值
                qvalueCutoff  = 0.05, #qvalue 阈值
                readable      = TRUE) #TRUE 则展示SYMBOL，FALSE 则展示原来的ID

#将 MF 对象转换为 dataframe，新版本可以用 as.data.frame(MF)
MF_results<-summary(MF)
#生成 barplot PDF 格式，x 轴为 GeneRatio，展示前 20 富集的 GO，数字可以调整
pdf(file = "MF_barplot.pdf")
barplot(MF, showCategory=20, x = "GeneRatio")
dev.off()

#生成 MF 气泡图
dotplot(MF)

#GO_CC 富集,基于基因数目
CC <- enrichGO(gene          = DE_list, #差异基因 vector
                keyType       ="ENTREZID",
OrgDb         = org.Hs.eg.db, #对应的OrgDb
                ont           = "CC", #GO 分类名称，CC BP MF
                pAdjustMethod = "BH", #Pvalue 矫正方法
                pvalueCutoff  = 0.05, #Pvalue 阈值
                qvalueCutoff  = 0.05, #qvalue 阈值
                readable      = TRUE) #TRUE 则展示SYMBOL，FALSE 则展示原来的ID

#将 CC 对象转换为 dataframe，新版本可以用 as.data.frame(CC)
CC_results<-summary(CC)

#生成 barplot PDF 格式，x 轴为 GeneRatio，展示前 20 富集的 GO，数字可调整
pdf(file = "CC_barplot.pdf")
barplot(CC, showCategory=20, x = "GeneRatio")
dev.off()

#生成 CC 气泡图
dotplot(CC)

#GO_BP 富集,基于基因数目
BP <- enrichGO(gene          = DE_list, #差异基因 vector
                keyType       ="ENTREZID",
OrgDb         = org.Hs.eg.db, #对应的OrgDb
                ont           = "BP", #GO 分类名称，CC BP MF
                pAdjustMethod = "BH", #Pvalue 矫正方法
                pvalueCutoff  = 0.05, #Pvalue 阈值
                qvalueCutoff  = 0.05, #qvalue 阈值
                readable      = TRUE) #TRUE 则展示SYMBOL，FALSE 则展示原来的ID

#将 BP 对象转换为 dataframe，新版本可以用 as.data.frame(BP)
BP_results<-summary(BP)

#生成 barplot PDF 格式，x 轴为 GeneRatio，展示前 20 富集的 GO，数字可调整
pdf(file = "BP_barplot.pdf")
barplot(BP, showCategory=20, x = "GeneRatio")
dev.off()

#生成 BP 气泡图
dotplot(BP)

#KEGG pathway 富集
ekp <- enrichKEGG(gene         = DE_list,
                 keyType = "kegg",
                 organism     = 'hsa',
                 pvalueCutoff = 0.05)
ekp_results <- summary(ekp)

#生成 KEGG 富集分析的 barplot 图，数字可调整
barplot(ekp, showCategory=20, x = "GeneRatio")

#生成气泡图
dotplot(ekp)

#基因和富集排名第 1 的pathway对应关系
cnetplot(ekp, showCategory = 1)

#输出 pathway 富集结果，可以用 excel 打开查看
write.table(ekp, file = "ekp.txt",
              sep = "\t", quote = F, row.names = T)

#查看通路
browseKEGG(ekp,'hsa04512')

过程文件展示

MF_barplot

MF_dotplot

KEGG_barplot

KEGG_dotplot

KEGG enrichment pathway browse


补充说明

#R3.6.0 以上版本安装方法不同于 R3.5.2 及以下版本，biocManager 安装方法如下：
If(!requireNamespace(「BiocManager」,quietly=TRUE))
Install.packages(「BiocManager」)
BiocManager::install(「clusterProfiler」,version = 「3.8」)