三阴性乳腺癌表达矩阵探索笔记之GSEA

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是学徒写的《GEO数据挖掘课程》的配套笔记（第6篇）

回顾：

找差异基因时，挑出的上调和下调的基因在250-500之间都是可以接受的。如果基因数量太少可以对阈值进行调整
- 工具limma
得到差异基因可以用来做基于超几何分布检验的富集分析,富集分析需要ENTREZID
- 工具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分析文件准备:

GSEA需要一个包含基因名以及logFC值的表格
KEGG和GO注释的数据库被做成了一个可操作的R包，这里的GSEA分析需要用到的数据库是需要自己下载的，MSigDB
根据下载数据库需要的ID进行ID转换，SYMBOL或者ENTREZID

  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结果解读：

第一步我们需要根据基因的logFC对基因进行排序
我们研究的这个染色体区域所包含的所有基因映射到我们排列好顺序的基因列表上，计算Enrich score的原则就是，从前到后依次检查基因是否是我们当前研究的染色体区域所包含的，如果包含就加一个正值，如果不包含就加一个负值，如果这个染色体区域中的所有基因被隐射到基因列表的一个集中区域，那么就会出现一个高峰。
横坐标表示所有探针的数量length(geneList)可以得到
黑色的竖线代表的是我们的目的基因，已经被排好序，如果竖线聚集在头部，称为头部效应，如果在尾部，称为尾部效应
GSEA也可以进行GO和KEGG分析，找到对应的数据集即可

图一：Enrichment score 高, 基因显著富集到右边，绿色的曲线出现高峰，而且，黑色竖线的分布也比较密集，说明该通路受到影响。

Enrichment score为负数，绿色曲线开口向上。

gsea_kegg.Rplot01

图二：Enrichment score 低，以20000为分界线，左边和右边的黑色竖线（代表基因）分布的密集程度差不多，而且绿色的曲线没有单独的峰值，所以并没有显著富集，该通路没有改变。

gsea_均匀分布.Rplot01

Java版的GSEA

下载地址：https://www.gsea-msigdb.org/gsea/downloads.jsp

根据自己系统需要下载合适的版本

GSEA_JAVA版

MSIGDB数据库下载

文献写的超级清楚，我这里就不截图演示啦！

这些分析，基本上读一下我五年前在生信技能树的表达芯片的公共数据库挖掘系列推文 就明白了；

视频观看方式

我把3年前的收费视频课程：3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月，免费到B站:

这个课程超级棒，B站免费学习咯：https://m.bilibili.com/video/BV1dy4y1C7jz
配套代码在GitHub哈：https://github.com/jmzeng1314/GSE76275-TNBC
TCGA数据库挖掘，代码在：https://github.com/jmzeng1314/TCGA_BRCA
GTEx数据库挖掘，代码在：https://github.com/jmzeng1314/gtex_BRCA
METABRIC数据库挖掘，代码在：https://github.com/jmzeng1314/METABRIC

然后马上就有了3千多学习量，而且有学员给出来了图文并茂版本万字笔记，让我非常感动！

扫描下面二维码马上就可以学习起来啦，笔记需要至少半个小时来阅读哦！