下面是学徒写的《GEO数据挖掘课程》的配套笔记(第6篇)
回顾:
limma
clusterProfiler
基于超几何分布检验的富集分析做KEGG数据库的时候,它总共只有七千多个基因,人类总的背景基因有两万多个,被KEGG记住的只有6500个(一直在增加),假设一条通路有117个基因参与,我们的差异基因中有10个与之重合,这已经是很多了,超几何分布检验会判定是统计学显著。
总和 | 部分 | 所占比例 |
---|---|---|
KEGG:6500 | 全部差异基因:162 | 0.02492308 |
某通路:117 | 差异基因重合:10 | 0.08547009 |
但是这样的超几何分布检验过于依赖我们的差异基因的选择。差异基因靠的的logFC和P value,一直是被人诟病的。绝大部分基因都没有被考虑进去,如果想考虑全部的基因, 就需要替换到gsea和gsva算法。今天我们先介绍gsea,如下:
GSEA :gene set enrichment analysis,即基因富集,是常见的富集分析之一
基本原理:先给定一个排序的基因表L和一个预先定义的基因集S,比如编码某个代谢通路的产物的基因集,基因组的物理位置相近的基因,或者同一个GO注释下的基因。GSEA的目的是判断S里面的成员s在L里面是随机分布还是主要聚集在L的顶部和底部。这些基因排序的依据是在其不同表型状态下的差异表达,如果研究的基因集S的成员显著聚集在L的顶部或者底部,则说明此基因集成员对表型的差异有显著,也是我们所关注的基因集。
GSEA分析文件准备:
load(file = 'anno_DEG.Rdata') #1.导入差异分析结果
geneList=DEG$logFC #2.取出logFC结果和基因的对应关系
names(geneList)=DEG$symbol
geneList=sort(geneList,decreasing = T) #3.对基因进行排序,用于统计学分析,这些基因相对于总的基因来说是随机挑选的,如果这些基因集中聚集在某个位置,这就破坏了随机性。
#选择gmt文件(MigDB中的全部基因集)
d='./MsigDB/symbols'
gmts <- list.files(d,pattern = 'all') #总共有7个基因集,这里把这些基因集全部合并到一起,一次完成
gmts
#GSEA分析
library(GSEABase) # BiocManager::install('GSEABase')
## 下面使用lapply循环读取每个gmt文件,并且进行GSEA分析,自动将每个基因集都检测一遍,最终得到enrich score,pvalue等值
gsea_results <- lapply(gmts, function(gmtfile){
# gmtfile=gmts[2]
# 如果不确定循环里面的代码,就尝试赋值,测试下面的代码。
geneset <- read.gmt(file.path(d,gmtfile))
print(paste0('Now process the ',gmtfile))
egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
head(egmt)
# gseaplot(egmt, geneSetID = rownames(egmt[1,])) #绘图
return(egmt)
})
save(gset_results,file="gset_results.Rdata") #gset_results中包含在所有的egmt
gsea_results_list<- lapply(gsea_results, function(x){
cat(paste(dim(x@result)),'\n')
x@result #S4对象可以用@符号取出值
})
gsea_results_df <- do.call(rbind, gsea_results_list) #将list变为dataframe,便于肉眼观察
gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE',title = 'KEGG_CELL_CYCLE') #挑选的"KEGG_CELL_CYCLE"在gsea_results的第二个元素中,一定要对应好,否则会出错
gseaplot(gsea_results[[2]],'FARMER_BREAST_CANCER_CLUSTER_6')
gseaplot(gsea_results[[5]],'GO_CONDENSED_CHROMOSOME_OUTER_KINETOCHORE')
GSEA结果解读:
length(geneList)
可以得到图一:Enrichment score 高, 基因显著富集到右边,绿色的曲线出现高峰,而且,黑色竖线的分布也比较密集,说明该通路受到影响。
Enrichment score为负数,绿色曲线开口向上。
gsea_kegg.Rplot01
图二:Enrichment score 低,以20000为分界线,左边和右边的黑色竖线(代表基因)分布的密集程度差不多,而且绿色的曲线没有单独的峰值,所以并没有显著富集,该通路没有改变。
gsea_均匀分布.Rplot01
下载地址:https://www.gsea-msigdb.org/gsea/downloads.jsp
根据自己系统需要下载合适的版本
GSEA_JAVA版
MSIGDB数据库下载
文献写的超级清楚,我这里就不截图演示啦!
这些分析,基本上读一下我五年前在生信技能树的表达芯片的公共数据库挖掘系列推文 就明白了;
视频观看方式
我把3年前的收费视频课程:3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月,免费到B站:
然后马上就有了3千多学习量,而且有学员给出来了图文并茂版本万字笔记,让我非常感动!
扫描下面二维码马上就可以学习起来啦,笔记需要至少半个小时来阅读哦!