• GEO生信数据挖掘（八）富集分析（GO 、KEGG、 GSEA 打包带走）

    第六节，我们使用结核病基因数据，做了一个数据预处理的实操案例。例子中结核类型，包括结核，潜隐进展，对照和潜隐，四个类别。第七节延续上个数据，进行了差异分析。 本节对差异基因进行富集分析。

    目录

    数据展示

    GO富集分析 -对基因名称映射基因ID

    GO富集分析 -从org.Hs.eg.db库中去匹配基因

    KEGG富集分析 （不详细讲了看注释）

    GSEA 富集分析

    更多复杂的图（关联网络图、八卦图 、弦图）

    数据展示

    差异基因计算完毕的指标如下图所示

    差异基因筛选后表达矩阵

    GO富集分析 -对基因名称映射基因ID

    加载数据

    1.  
    2. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    3. #+&&&&&&&&&&&&&&&&&&加载数据&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    4. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    5. load( "DEG_TB_LTBI_step13.Rdata")
    6. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    7. #+&&&&&&&&&&&&&&&&&&加载数据&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    8. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    9.  
    10.  
    11. library(clusterProfiler)
    12. library(org.Hs.eg.db)
    13.  
    14. #增加基因名
    15. all_diff$SYMBOL=rownames(all_diff)
    16.  
    17.  
    18. #基因名称转换注释
    19. gene_ids_DEG_TB_LTBI = bitr(geneID = rownames(dataset_TB_LTBI_DEG),fromType="SYMBOL",toType = c("ENTREZID","ENSEMBL","SYMBOL"),OrgDb = 'org.Hs.eg.db',drop =TRUE)
    20.  
    21. #合并 增加logFC 为后续GSEA富集分析所需数据准备
    22. gene_ids_DEG_TB_LTBI <- merge(gene_ids_DEG_TB_LTBI,all_diff,by="SYMBOL")
    1. #观察
    2. dim(gene_ids_DEG_TB_LTBI)
    3. head(gene_ids_DEG_TB_LTBI)
    4.  
    5. #获取基因ID ENSEMBL
    6. gene_ENSEMBL <- gene_ids_DEG_TB_LTBI$ENSEMBL
    7. gene_ENTREZID <- gene_ids_DEG_TB_LTBI$ENTREZID
    8. gene_SYMBOL<- gene_ids_DEG_TB_LTBI$SYMBOL

    经过映射，2813个差异基因得到2551个基因ID,下图为三种不同形式的基因名称，富集分析时，按需进行转换。

    GO富集分析 -从org.Hs.eg.db库中去匹配基因

    1. #Go富集分析，从库中去匹配
    2. go <- enrichGO(gene_SYMBOL,OrgDb = org.Hs.eg.db, ont='ALL',pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'SYMBOL')#进行GO富集，确定P值与Q值得卡值并使用BH方法对值进行调整。
    3. #查看富集结果
    4. dim(go)
    5. #导出GO富集的结果
    6. write.csv(go,file="go1.csv")

    绘制气泡图

    1. #绘制气泡图
    2. pdf(file="15aGO富集分析step15.pdf", width = 9, height = 6)
    3. dotplot(go,showCategory=20,label_format = 80)#气泡图
    4.  dev.off()

    三种不同类别的合并的气泡图（#CC细胞组件，MF分子功能，BP生物学过程）

    1. pdf(file="15bGO富集分析三组step15.pdf", width = 9, height = 6)
    2.  
    3. #CC细胞组件，MF分子功能，BP生物学过程
    4. goCC <- enrichGO(gene = gene_ENTREZID,  #基因列表(转换的ID)
    5.                keyType = "ENTREZID",  #指定的基因ID类型，默认为ENTREZID
    6.                OrgDb=org.Hs.eg.db,  #物种对应的org包
    7.                ont = "CC",   #CC细胞组件，MF分子功能，BP生物学过程
    8.                pvalueCutoff = 0.05,  #p值阈值
    9.                pAdjustMethod = "fdr",  #多重假设检验校正方式
    10.                minGSSize = 1,   #注释的最小基因集，默认为10
    11.                maxGSSize = 500,  #注释的最大基因集，默认为500
    12.                qvalueCutoff = 0.2,  #q值阈值
    13.                readable = TRUE)  #基因ID转换为基因名
    14.  
    15. goBP <- enrichGO(gene_ENTREZID,OrgDb = org.Hs.eg.db, ont='BP',pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENTREZID')
    16.  
    17. goMF <- enrichGO(gene_ENTREZID,OrgDb = org.Hs.eg.db, ont='MF',pAdjustMethod = 'BH',pvalueCutoff = 0.05, qvalueCutoff = 0.2,keyType = 'ENTREZID')
    18. #通过ggplot2将BP、MF、CC途径的富集结果挑选前8条绘制在一张图上
    19. barplot(go,label_format=100, split="ONTOLOGY")+ facet_grid(ONTOLOGY~.,scale="free")
    20.  
    21. dev.off()

    KEGG富集分析 （不详细讲了看注释）

    1.  
    2.  
    3. #+=================================================================
    4. #============================================================
    5. #+========KEGG富集分析 气泡图step16===================
    6. #+==========================================
    7. #+================================
    8.  
    9.  
    10. #KEGG富集分析
    11. pdf(file="16KEGG富集分析step16.pdf", width = 9, height = 6)
    12. kegg<- enrichKEGG(gene = gene_ENTREZID,   #基因列表(ENTREZID ID: 54490,51144,31,3906) 
    13.                   organism = "hsa",  #物种
    14.                   keyType = "kegg",  #指定的基因ID类型，默认为kegg
    15.                   minGSSize = 3, 
    16.                   maxGSSize = 500,
    17.                   pvalueCutoff = 0.05,  
    18.                   pAdjustMethod = "fdr", # pAdjustMethod = 'BH'
    19.                   qvalueCutoff = 0.02)
    20. #观察
    21. dim(kegg)
    22. #绘制气泡图
    23. dotplot(kegg)
    24. dev.off()
    25.  
    26. #kegg 增加可读性，对基因ID 转基因名
    27. kegg_enrich_results <- DOSE::setReadable(kegg, 
    28.                                          OrgDb="org.Hs.eg.db", 
    29.                                          keyType='ENTREZID') #ENTREZID to gene Symbol
    30.  
    31. #保存kegg结果
    32. write.csv(kegg_enrich_results@result,'KEGG_gene_up_enrichresults.csv') 
    33. #save(kegg_enrich_results, file ='KEGG_gene_up_enrichresults.Rdata')
    34.  
    35.  
    36. ##查看与选择所需通路
    37. kegg_enrich_results@result$Description[1:10] #查看前10通路
    38.  
    39. ###选择所需通路的ID号
    40. i=1
    41. select_pathway <- kegg_enrich_results@result$ID[i] #选择所需通路的ID号
    42.  
    43.  
    44.  
    45.  
    46. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    47. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    48. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    49. save(gene_ids_DEG_TB_LTBI,go,keggfile ="15_gene_ids_DEG_TB_LTBI.Rdata")
    50. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    51. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    52. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

    GSEA 富集分析

    1.  
    2. #+=================================================================
    3. #============================================================
    4. #+========GSEA 富集分析 气泡图step17===================
    5. #+==========================================
    6. #+================================
    7.  
    8. # GSEA 分析
    9. #需要把多个方法取并集
    10. #该方法的输入需要基因和 logFC 排序后的结果
    11. #不同方法 相同基因的的logFC值不一样，直接保留第一个重复基因
    12.  
    13.  
    14. library(stringr)
    15.  
    16.  
    17.  
    18. ## 去重  #去除NA值
    19. dim(gene_ids_DEG_TB_LTBI)
    20. colnames(gene_ids_DEG_TB_LTBI)
    21.  
    22. gene_list_df = gene_ids_DEG_TB_LTBI[,c('ENTREZID','logFC')]
    23. gene_list_df_na <- na.omit(gene_list_df)
    24.  
    25. gene_ids_TB_LTBI_distinct <- dplyr::distinct(gene_list_df_na,ENTREZID,.keep_all=TRUE) 
    26.  
    27. dim(gene_ids_TB_LTBI_distinct)
    28. gene_list=gene_ids_TB_LTBI_distinct$logFC   #提取logFC列
    29. names(gene_list)=gene_ids_TB_LTBI_distinct$ENTREZID             #加上ENTREZID
    30. gene_list_gsea = sort(gene_list, decreasing = T)  #降序排列
    31.  
    32.  
    33.  
    34.  
    35.  
    36. gsea_KEGG <- gseKEGG(gene_list_gsea,
    37.                      organism = "hsa",
    38.                      keyType = "kegg")
    39.  
    40. gsea_KEGG_d <- as.data.frame(gsea_KEGG)
    41.  
    42. gsea_KEGG_d
    43. #path 为需要展示的pathway id，这里展示的是enrichment score最高的4条通路
    44.  
    45.  
    46. t_index=order(gsea_KEGG_d$enrichmentScore,decreasing = T)
    47. path=rownames(gsea_KEGG[t_index,]) #选择展示的 pathwayrownames(gsea_KEGG[t_index,]) [1:4]
    48.  
    49. #作图 
    50. pdf(file="17GSEA富集分析step17.pdf", width = 9, height = 6)
    51. gseaplot2(gsea_KEGG,
    52.           path, 
    53.           subplots = 1:2, #展示前2个图
    54.           pvalue_table = T, #显示p值
    55.           title = "Olfactory transduction",  #设置title
    56.           base_size = 10, #字体大小
    57.           color="red") #线条颜色可选
    58. dev.off()
    59.  
    60.  
    61. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    62. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    63. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    64. save(gene_ids_DEG_TB_LTBI,go,kegg,gene_list_gsea,gsea_KEGG,file ="17_gene_ids_DEG_TB_LTBI.Rdata")
    65. #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    66. #+&&&&&&&&&&&&&&&&&&数据保存&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
    67. #&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&

    更多复杂的图（关联网络图、八卦图 、弦图）

    参考 最全的GO, KEGG, GSEA分析教程(R),你要的高端可视化都在这啦！[包含富集圈图] - 糖糖家的老张的文章 - 知乎 https://zhuanlan.zhihu.com/p/377356510

    原文链接：https://blog.csdn.net/qq_50898257/article/details/120588222

    1. #+=================================================================
    2. #============================================================
    3. #+========富集分析 更多的图step18===================
    4. #+==========================================
    5. #+================================
    6. library(clusterProfiler)
    7. library(enrichplot)
    8. #+富集基因与所在功能集/通路集的关联网络图：
    9. enrichplot::cnetplot(go,circular=FALSE,colorEdge = TRUE)#基因-通路关联网络图
    10. enrichplot::cnetplot(kegg,circular=FALSE,colorEdge = TRUE)#circluar为指定是否环化，基因过多时建议设置为FALSE
    11.  
    12. GO2 <- pairwise_termsim(go)
    13. KEGG2 <- pairwise_termsim(kegg)
    14. enrichplot::emapplot(GO2,showCategory = 50, color = "p.adjust", layout = "kk")#通路间关联网络图
    15. enrichplot::emapplot(KEGG2,showCategory =50, color = "p.adjust", layout = "kk")
    16.  
    17.  
    18. write.table(kegg$ID, file = "KEGG_IDs.txt", #将所有KEGG富集到的通路写入本地文件查看
    19.             append = FALSE, quote = TRUE, sep = " ",
    20.             eol = "\n", na = "NA", dec = ".", row.names = TRUE,
    21.             col.names = TRUE, qmethod = c("escape", "double"),
    22.             fileEncoding = "")
    23. #打印几条通路名称看看
    24. kegg$ID[1:3]
    25. #打开浏览器观察通路
    26. browseKEGG(kegg,"hsa04660")#选择其中的hsa05166通路进行展示
    27.  
    28.  
    29.  
    30.  
    31. #富集弦图
    32. genedata<-data.frame(ID=gene_ids_DEG_TB_LTBI$SYMBOL  ,logFC=gene_ids_DEG_TB_LTBI$logFC)
    33. write.table(go$ONTOLOGY, file = "GO_ONTOLOGYs.txt", #将所有GO富集到的基因集所对应的类型写入本地文件从而得到BP/CC/MF各自的起始位置如我的数据里是1，2103，2410
    34.             append = FALSE, quote = TRUE, sep = " ",
    35.             eol = "\n", na = "NA", dec = ".", row.names = TRUE,
    36.             col.names = TRUE, qmethod = c("escape", "double"),
    37.             fileEncoding = "")
    38.  
    39. '''
    40. 根据计算出的go 文件数量，调整
    41. '''
    42. GOplotIn_BP<-go[1:178,c(2,3,7,9)] #提取GO富集BP的前10行,提取ID,Description,p.adjust,GeneID四列
    43. GOplotIn_CC<-go[179:194,c(2,3,7,9)]#提取GO富集CC的前10行,提取ID,Description,p.adjust,GeneID四列
    44. GOplotIn_MF<-go[195:209,c(2,3,7,9)]#提取GO富集MF的前10行,提取ID,Description,p.adjust,GeneID四列
    45.  
    46. library(stringr)
    47. GOplotIn_BP$geneID <-str_replace_all(GOplotIn_BP$geneID,'/',',') #把GeneID列中的’/’替换成‘,’
    48. GOplotIn_CC$geneID <-str_replace_all(GOplotIn_CC$geneID,'/',',')
    49. GOplotIn_MF$geneID <-str_replace_all(GOplotIn_MF$geneID,'/',',')
    50. names(GOplotIn_BP)<-c('ID','Term','adj_pval','Genes')#修改列名,后面弦图绘制的时候需要这样的格式
    51. names(GOplotIn_CC)<-c('ID','Term','adj_pval','Genes')
    52. names(GOplotIn_MF)<-c('ID','Term','adj_pval','Genes')
    53. GOplotIn_BP$Category = "BP"#分类信息
    54. GOplotIn_CC$Category = "CC"
    55. GOplotIn_MF$Category = "MF"
    56.  
    57. BiocManager::install('GOplot')
    58. library(GOplot)
    59. circ_BP<-GOplot::circle_dat(GOplotIn_BP,genedata) #GOplot导入数据格式整理
    60. circ_CC<-GOplot::circle_dat(GOplotIn_CC,genedata) 
    61. circ_MF<-GOplot::circle_dat(GOplotIn_MF,genedata) 
    62.  
    63. chord_BP<-chord_dat(data = circ_BP,genes = genedata) #生成含有选定基因的数据框
    64. chord_CC<-chord_dat(data = circ_CC,genes = genedata) 
    65. chord_MF<-chord_dat(data = circ_MF,genes = genedata) 
    66. '''
    67. > chord_CC<-chord_dat(data = circ_CC,genes = genedata) 
    68. Error in `[<-`(`*tmp*`, g, p, value = ifelse(M[g] %in% sub2$genes, 1,  : 
    69.   subscript out of bounds
    70.   我去检查了go和genelist的数据结构发现，genelist里的gene用的是gene名，而go里的基因用的是基因ID，不一样了，所以跑不出结果，所以我把genelist的gene换成了基因ID，就能跑出来了。
    72. 作者：ff的小世界勿扰
    73. 链接：https://www.jianshu.com/p/ee4012fd253f
    74. 来源：简书
    75. 著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。
    76. '''
    77.  
    78. #可以画 数量太多了
    79. GOChord(data = chord_BP,#弦图
    80.         title = 'GO-Biological Process',space = 0.01,#GO Term间距
    81.         limit = c(1,1),gene.order = 'logFC',gene.space = 0.25,gene.size = 5,
    82.         lfc.col = c('red','white','blue'), #上下调基因颜色
    83.         process.label = 10) #GO Term字体大小
    84. GOChord(data = chord_CC,title = 'GO-Cellular Component',space = 0.01,
    85.         limit = c(1,1),gene.order = 'logFC',gene.space = 0.25,gene.size = 5,
    86.         lfc.col = c('red','white','blue'), 
    87.         process.label = 10) 
    88. GOChord(data = chord_MF,title = 'GO-Mollecular Function',space = 0.01,
    89.         limit = c(1,1),gene.order = 'logFC',gene.space = 0.25,gene.size = 5,
    90.         lfc.col = c('red','white','blue'), 
    91.         process.label = 10)
    92. '''
    93. Warning messages:
    94. 1: Using size for a discrete variable is not advised. 
    95. 2: Removed 15 rows containing missing values (`geom_point()`). 
    96. '''

    富集分析完毕！

    回顾我们用到方法，差异分析后进行富集分析，理论基础实际上就是简单的找不同，分析。

    实际应用种，由于基因之间存在关联，另一套分析理论考虑的是基因之间的相互作用，下节，我们来看非常火的WGCNA 共表达加权网络进行基因分析。

  • 相关阅读:
    三、循环语句
    浙江大学数据结构MOOC-课后习题-第十讲-排序4 统计工龄
    在 EMR 上使用 Delta Lake
    Mac代码文本编辑器Sublime Text 4
    R语言作业--第六章判别分析
    Git三剑客之基础部分
    WorldPop2000年至2020年的全中国的人口统计数据
    ubuntu bind9 主从配置
    毕业设计 基于大数据的股票量化分析与股价预测系统
    运动装备什么牌子好？运动装备品牌排行榜推荐
  • 原文地址：https://blog.csdn.net/zzh1464501547/article/details/126898776