GEO数据挖掘6

sunqi

2020/7/13

概述

使用SigDB（Molecular Signatures Database）基因集进行富集分析，包含8个系列

H: hallmark gene sets （癌症）特征基因集合，共50组，最常用；
C1: positional gene sets 位置基因集合，根据染色体位置，共326个，用的很少
C2: curated gene sets：（专家）校验基因集合，基于通路、文献等：
C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分
C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；
C5: GO gene sets：Gene Ontology 基因本体论，包括BP（生物学过程biological process，细胞原件cellular component和分子功能molecular function三部分）
C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO 发表芯片数据
C7: immunologic signatures: 免疫相关基因集合。

相较于KEGG，SigDB数据集包含的功能更多

GSEA分析

对 MigDB中的全部基因集做GSEA分析。

library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)


# 安装需要的包
# BiocManager::install("GSEABase")
# 导入kegg已经注释好的数据
load(file = 'anno_DEG.Rdata')
geneList=DEG$logFC
names(geneList)=DEG$symbol
geneList=sort(geneList,decreasing = T)

# 导入MigDB全部基因集
# 需要在 MigDB 官网下载 为gmt文件
# 下载好的路径
d='MsigDB/symbols'
gmts <- list.files(d,pattern = 'all')
# 导入8个序列的基因集名
#
library(GSEABase)
## 使用lapply循环读取每个gmt文件并进行GSEA分析
# 考验你计算机能力的时候到了
f='gsea_results.Rdata'
if(!file.exists(f)){
  #相比较apply，lapply较多的用于list的循环操作
  gsea_results <- lapply(gmts, function(gmtfile){
    #读取gmt文件
  geneset <- read.gmt(file.path(d,gmtfile))
  # 打印目前进行的基因库
  print(paste0('Now process the ',gmtfile))
  # 进行GSEA分析
  egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
  head(egmt)
  # 返回计算的结果
  return(egmt)
  })
  # 保存结果到本地文件
  save(gsea_results,file = f)
}
# 如果上述操作完成之后
# 再次运算就可以直接导入
load(file = f)
#提取gsea结果，
gsea_results_list<- lapply(gsea_results, function(x){
  cat(paste(dim(x@result)),'n')
  x@result
})

## 3 11
## 996 11
## 186 11
## 233 11
## 671 11
## 95 11
## 1591 11
## 27 11

# docall 函数能够对list使用dataframe结构的函数，下行为合并结果
gsea_results_df <- do.call(rbind, gsea_results_list)

# 选择有差异的基因集进行画图,第一个参数为基因集，第二个参数需要对应前面参数
# 具体可以查看gsea_results里面的行名
gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE')

gseaplot(gsea_results[[2]],'FARMER_BREAST_CANCER_CLUSTER_6')

gseaplot(gsea_results[[5]],'GO_CONDENSED_CHROMOSOME_OUTER_KINETOCHORE')