• 主页
  • 课程

    关于课程

    • 课程归档
    • 成为一名讲师
    • 讲师信息
    同等学历教学

    同等学历教学

    免费
    阅读更多
  • 特色
    • 展示
    • 关于我们
    • 问答
  • 事件
  • 个性化
  • 博客
  • 联系
  • 站点资源
    有任何问题吗?
    (00) 123 456 789
    weinfoadmin@weinformatics.cn
    注册登录
    恒诺新知
    • 主页
    • 课程

      关于课程

      • 课程归档
      • 成为一名讲师
      • 讲师信息
      同等学历教学

      同等学历教学

      免费
      阅读更多
    • 特色
      • 展示
      • 关于我们
      • 问答
    • 事件
    • 个性化
    • 博客
    • 联系
    • 站点资源

      未分类

      • 首页
      • 博客
      • 未分类
      • R|clusterProfiler-富集分析

      R|clusterProfiler-富集分析

      • 发布者 一览
      • 分类 未分类
      • 日期 2020年3月8日
      • 评论 0评论

      简单总结clusterProfiler包进行GO、KEGG的富集分析方法,结果输出及内置的图形展示。

      一 bioconductor包安装

      #PS:安装是否顺利,看缘分了(用本站weinformatics则不用担心这些问题,包都装好了,解决了你入门的最大问题)

      #source("https://bioconductor.org/biocLite.R")
      #biocLite("DOSE") 
      #biocLite("topGO")
      #biocLite("clusterProfiler")
      #biocLite("pathview")

      ##加载需要的R包

      library(DOSE)
      library(org.Hs.eg.db)
      library(topGO)
      library(clusterProfiler)
      library(pathview)

      二、数据载入及转化

          需要将输入的基因格式转为clusterProfiler可分析的格式,通过功能函数bitr实现各种ID之间的转化。通过keytypes函数可查看所有支持及可转化类型

      keytypes(org.Hs.eg.db) 

      1.向量方式读入#适合少量基因

      x <- c("AASDH","ABCB11","ADAM12","ADAMTS16","ADAMTS18")
      test = bitr(x, #数据集
      fromType="SYMBOL", #输入为SYMBOL格式
      toType="ENTREZID",  # 转为ENTERZID格式
      OrgDb="org.Hs.eg.db") #人类 数据库
      head(test,2)
      
          SYMBOL ENTREZID
      1    AASDH   132949
      2   ABCB11     8647

      2.基因list文件读入

      data <- read.table("gene",header=FALSE) #单列基因名文件
      data$V1 <- as.character(data$V1) #需要character格式,然后进行ID转化
      #将SYMBOL格式转为ENSEMBL和ENTERZID格式 
      test1 = bitr(data$V1, fromType="SYMBOL", toType=c("ENSEMBL", "ENTREZID"), OrgDb="org.Hs.eg.db")
      head(test1,2)
      
          SYMBOL         ENSEMBL ENTREZID
      1    AASDH ENSG00000157426   132949
      2   ABCB11 ENSG00000073734     8647
      3. 内置示例数据
      data(geneList, package="DOSE") #富集分析的背景基因集
      gene <- names(geneList)[abs(geneList) > 2]
      gene.df <- bitr(gene, fromType = "ENTREZID", toType = c("ENSEMBL", "SYMBOL"), OrgDb = org.Hs.eg.db)
      head(gene.df,2)
      
        ENTREZID         ENSEMBL SYMBOL
      1     4312 ENSG00000196611   MMP1
      2     8318 ENSG00000093009  CDC45

      三、 GO分析

      3.1 groupGO 富集分析

      ggo <- groupGO(gene = test1$ENTREZID, OrgDb = org.Hs.eg.db, ont = "CC",level = 3,readable = TRUE)

      3.2 enrichGO 富集分析

      ego_ALL <- enrichGO(gene = test1$ENTREZID, 
                      universe = names(geneList), #背景基因集
                      OrgDb = org.Hs.eg.db, #没有organism="human",改为OrgDb=org.Hs.eg.db
      				#keytype = 'ENSEMBL',
                      ont = "ALL", #也可以是 CC  BP  MF中的一种
                      pAdjustMethod = "BH", #矫正方式 holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”, “fdr”, “none”中的一种
                      pvalueCutoff = 1, #P值会过滤掉很多,可以全部输出
                      qvalueCutoff = 1,
      				readable = TRUE) #Gene ID 转成gene Symbol ,易读
      head(ego_ALL,2)
           ONTOLOGY         ID                                                Description
      GO:0002887       BP GO:0002887 negative regulation of myeloid leukocyte mediated immunity
      GO:0033004       BP GO:0033004                negative regulation of mast cell activation
                 GeneRatio  BgRatio      pvalue  p.adjust    qvalue                   geneID Count
      GO:0002887     2/121 10/11461 0.004706555 0.7796682 0.7796682              CD300A/CD84     2
      GO:0033004     2/121 10/11461 0.004706555 0.7796682 0.7796682              CD300A/CD84     2
      其中:ONTOLOGY:CC  BP  MF 
      GO ID: Gene Ontology数据库中唯一的标号信息
      Description :Gene Ontology功能的描述信息
      GeneRatio:差异基因中与该Term相关的基因数与整个差异基因总数的比值
      BgRation:所有( bg)基因中与该Term相关的基因数与所有( bg)基因的比值
      pvalue: 富集分析统计学显著水平,一般情况下, P-value < 0.05 该功能为富集项
      p.adjust 矫正后的P-Value
      qvalue:对p值进行统计学检验的q值
      geneID:与该Term相关的基因
      Count:与该Term相关的基因数

      3.2.1 setReadable函数进行转化

      ego_MF <- enrichGO(gene = test1$ENTREZID, universe = names(geneList),OrgDb = org.Hs.eg.db,ont = "MF", pAdjustMethod = "BH",pvalueCutoff = 1,qvalueCutoff = 1,readable = FALSE)
      ego_MF1 <- setReadable(ego_MF, OrgDb = org.Hs.eg.db)
      3.3 GSEA分析 (暂略)

      3.4 输出结果及图像

      3.4.1 输出结果

      write.csv(summary(ego_ALL),"ALL-enrich.csv",row.names =FALSE)

      3.4.2 绘制图形

      ##可视化--点图
      dotplot(ego_MF,title="EnrichmentGO_MF_dot")#点图,按富集的数从大到小的
      ##可视化--条形图
      barplot(ego_MF, showCategory=20,title="EnrichmentGO_MF")#条状图,按p从小到大排,绘制前20个Term
      ##可视化--
      plotGOgraph(ego_MF)

      3.5 过滤

      3.5.1 利用内置cutoff设置阈值,by设定指标

      #可能会有问题,还不太清楚
      ego_MF.fil <- simplify(ego_MF, cutoff=0.05, by="pvalue", select_fun=min)

      3.5.2 存储结果然后过滤

      ego_ALL.sig <- ego_ALL[ego_ALL$pvalue <= 0.01]

      过滤后为数据框,不能用自带的参数直接绘制,可以使用ggplot2进行绘制。(暂略)

      四、KEGG分析

      4.1 候选基因进行通路分析

      kk <- enrichKEGG(gene = test1$ENTREZID,
                       organism = 'hsa', #KEGG可以用organism = 'hsa'
                       pvalueCutoff = 1)
      head(kk,2)
                    ID                                      Description GeneRatio  BgRatio
      hsa04750 hsa04750 Inflammatory mediator regulation of TRP channels      5/53  97/7387
      hsa04020 hsa04020                        Calcium signaling pathway      6/53 182/7387
                     pvalue   p.adjust    qvalue                              geneID Count
      hsa04750 0.0006135305 0.08589427 0.0807277             40/3556/3708/5608/79054     5
      hsa04020 0.0018078040 0.12654628 0.1189345        493/1129/2066/3707/3708/4842     6

      4.2 富集结果及图形展示

      4.2.1 结果输出文件

      write.csv(summary(kk),"KEGG-enrich.csv",row.names =FALSE)

      4.2.1 KEGG 气泡图

      dotplot(kk,title="Enrichment KEGG_dot")

      4.2.2 查看特定通路图

      hsa04750 <- pathview(gene.data  = geneList,
      
      pathway.id = "hsa04750", #上述结果中的hsa04750通路
      
                           species    = "hsa",
      
                           limit      = list(gene=max(abs(geneList)), cpd=1))
      
      

      五、注释文件、注释库

      如果clusterProfiler包没有所需要物种的内置数据库,可以通过自定义注释文件或者自建注释库的方法进行富集分析。

      5.1 自定义注释文件

      5.1.1 待富集的基因list

      data(geneList, package="DOSE")
      deg <- names(geneList)[abs(geneList)>2] 

      5.1.2 读入注释文件:

      ## downloaded from http://www.disgenet.org/ds/DisGeNET/results/all_gene_disease_associations.tar.gz
      gda <- read.delim("all_gene_disease_associations.tsv")
      head(gda,1)
      
       geneId geneSymbol diseaseId    diseaseName       score NofPmids NofSnps source
      1      1       A1BG  C0001418 Adenocarcinoma 0.002732912        1       0  LHGDN

      自定义的两个注释文件(重要)

      disease2gene=gda[, c("diseaseId", "geneId")]
      disease2name=gda[, c("diseaseId", "diseaseName")]

      5.1.3 enricher函数进行富集分析

      x = enricher(deg, TERM2GENE=disease2gene, TERM2NAME=disease2name)
      head(summary(x),1)
         ID                Description GeneRatio   BgRatio       pvalue     p.adjust
      C3642345 C3642345 Luminal A Breast Carcinoma    11/198 110/17074 6.057328e-08 0.0001488891
                     qvalue                                              geneID Count
      C3642345 0.0001287342 9833/55355/3620/891/9122/2167/3169/5304/2625/9/5241    11

      5.2 数据库

      5.2.1 查看当前支持的物种

      Bioconductor – BiocViews

      5.2.2 通过AnnotationHub自建OrgDb注释库并富集

      library(AnnotationHub)
      hub <- AnnotationHub() #较慢
      query(hub, "Cricetulus")
      Cgriseus <- hub[["AH48061"]]
      sample_gene <- sample(keys(Cgriseus), 100)
      str(sample_gene)
      sample_test <- enrichGO(sample_gene, OrgDb=Cgriseus, pvalueCutoff=1, qvalueCutoff=1)
      head(summary(sample_test))

      六、参考资料

      Statistical analysis and visualization of functional profiles for genes and gene clusters

      以上,clusterProfiler包进行GO、KEGG的富集分析方法,结果输出及内置的图形展示。待补充富集结果ggplot2绘制,GSEA等相关分析。

      请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!

      • 分享:
      作者头像
      一览

      上一篇文章

      ggplot2-plotly|让你的火山图“活”过来
      2020年3月8日

      下一篇文章

      R包里的帮助文档是怎么写的
      2020年3月8日

      你可能也喜欢

      2-1675088548
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      30 1月, 2023
      9-1675131201
      如何快速批量修改 Git 提交记录中的用户信息
      26 1月, 2023
      8-1678501786
      肿瘤细胞通过改变CD8+ T细胞中的丙酮酸利用和琥珀酸信号来调控抗肿瘤免疫应答。
      7 12月, 2022

      留言 取消回复

      要发表评论,您必须先登录。

      搜索

      分类

      • R语言
      • TCGA数据挖掘
      • 单细胞RNA-seq测序
      • 在线会议直播预告与回放
      • 数据分析那些事儿分类
      • 未分类
      • 生信星球
      • 老俊俊的生信笔记

      投稿培训

      免费

      alphafold2培训

      免费

      群晖配置培训

      免费

      最新博文

      Nature | 单细胞技术揭示衰老细胞与肌肉再生
      301月2023
      lncRNA和miRNA生信分析系列讲座免费视频课和课件资源包,干货满满
      301月2023
      如何快速批量修改 Git 提交记录中的用户信息
      261月2023
      logo-eduma-the-best-lms-wordpress-theme

      (00) 123 456 789

      weinfoadmin@weinformatics.cn

      恒诺新知

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      链接

      • 课程
      • 事件
      • 展示
      • 问答

      支持

      • 文档
      • 论坛
      • 语言包
      • 发行状态

      推荐

      • iHub汉语代码托管
      • iLAB耗材管理
      • WooCommerce
      • 丁香园论坛

      weinformatics 即 恒诺新知。ICP备案号:粤ICP备19129767号

      • 关于我们
      • 博客
      • 联系
      • 成为一名讲师

      要成为一名讲师吗?

      加入数以千计的演讲者获得100%课时费!

      现在开始

      用你的站点账户登录

      忘记密码?

      还不是会员? 现在注册

      注册新帐户

      已经拥有注册账户? 现在登录

      close
      会员购买 你还没有登录,请先登录
      • ¥99 VIP-1个月
      • ¥199 VIP-半年
      • ¥299 VIP-1年
      在线支付 激活码

      立即支付
      支付宝
      微信支付
      请使用 支付宝 或 微信 扫码支付
      登录
      注册|忘记密码?