公司简介
首页 > 关于翼和 > 公司动态 > (会议)【技术贴】如何从数据库挖掘基因并筛选TagSNP(医学篇)(1)

【技术贴】如何从数据库挖掘基因并筛选TagSNP(医学篇)(1)

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

技术路线

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

百度搜索可以获得很多人疾病相关数据库的使用说明,在此不再赘述。我们常用的人类疾病数据库是DisGeNEThttps://www.disgenet.org/search ),常用GWAS数据库是EMBL-EBIGWAS cataloghttps://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”)

 

#基因编号转换

#将基因名称列表复制给sxl(或者任何你喜欢的文件名),如需分析自己特定的基因集,可替换括号内容,每个基因名称,用””,隔开。

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 则展示SYMBOLFALSE 则展示原来的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 则展示SYMBOLFALSE 则展示原来的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 则展示SYMBOLFALSE 则展示原来的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)

#基因和富集排名第1pathway对应关系

cnetplot(ekp, showCategory = 1)

#输出pathway富集结果,可以用excel打开查看

write.table(ekp, file = "ekp.txt",

              sep = "\t", quote = F, row.names = T)

#查看通路

browseKEGG(ekp,'hsa04512')

 

补充说明

#R3.6.0以上版本安装方法不同于R3.5.2及以下版本,biocManager安装方法如下:

If(!requireNamespace(“BiocManager”,quietly=TRUE))

Install.packages(“BiocManager”)

BiocManager::install(“clusterProfiler”,version = “3.8”)

 

过程文件展示

MF_barplot

 

MF_dotplot


KEGG_barplot

 

KEGG_dotplot

 

 

KEGG enrichment pathway browse