最近在微信群看到了一个交流，如何使用最少的代码完成GSEA分析，并且绘制美图！目前得分最高的是4句话，如下所示：
载入测试数据做GSEA
需要3个包，分别是：'clusterProfiler','enrichplot','patchwork'，然后是DOSE包里面有一个geneList的向量，它是排序好的基因列表，而且是entrezID形式，使用 gseKEGG 函数即可做gsea分析啦 ：
lapply(c('clusterProfiler','enrichplot','patchwork'),
function(x) {library(x, character.only = T)})
# Please go to https://yulab-smu.github.io/clusterProfiler-book/ for the full vignette.
data(geneList, package="DOSE")
#4312 8318 10874 55143 55388 991
#4.572613 4.514594 4.418218 4.144075 3.876258 3.677857
class(geneList)
#[1] "numeric"
kk2 <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 10000,
minGSSize = 10,
maxGSSize = 200,
pvalueCutoff = 0.05,
pAdjustMethod = "none" )
gseKEGG输入形式：将基因按照logFC进行从高到低排序，只需要基因列和logFC organism：物种，http://www.genome.jp/kegg/catalog/org_list.html nPerm：permutation numbers minGSSize：通路最小基因数 maxGSSize：通路最大基因数（一般可通过改变minGSSize，maxGSSize数目调整通路大小） pvalueCutoff：最小p值 pAdjustMethod：p值校正方法，"BH"
输入形式非常关键！通常差异分析后会形成如下所示 数据框：
lapply(c('org.Hs.eg.db','stringr','dplyr'),
function(x) {library(x, character.only = T)})
load(file = 'step2-output.Rdata')
head(deg)
# logFC AveExpr t P.Value adj.P.Val
# COL11A1 7.235897 9.241369 13.66359 1.499441e-14 3.500895e-10
# ZIC2 5.658879 7.917121 13.26803 3.244100e-14 3.787162e-10
# PGM5-AS1 -3.626344 7.758838 -12.66304 1.090802e-13 6.546434e-10
# PGM5 -3.020465 10.639885 -12.59561 1.251776e-13 6.546434e-10
# ANGPTL1 -4.141375 9.283651 -12.53414 1.419781e-13 6.546434e-10
# SHOX2 -5.304161 8.434325 -12.45164 1.682311e-13 6.546434e-10
# B
# COL11A1 22.87130
# ZIC2 22.15323
# PGM5-AS1 21.01867
# PGM5 20.88943
# ANGPTL1 20.77110
# SHOX2 20.61157
简单的差异分析看我六年前的表达芯片的公共数据库挖掘系列推文即可哈 ：
解读GEO数据存放规律及下载，一文就够 解读SRA数据库规律一文就够 从GEO数据库下载得到表达矩阵 一文就够 GSEA分析一文就够（单机版+R语言版） 根据分组信息做差异分析- 这个一文不够的 差异分析得到的结果注释一文就够
很容易把差异分析结果转为类似的 DOSE包里面有一个geneList的向量 ，代码如下所示：
## 不同的阈值，筛选到的差异基因数量就不一样，后面的超几何分布检验结果就大相径庭。
logFC_t=1.5
deg$g=ifelse(deg$P.Value>0.05,'stable',
ifelse( deg$logFC > logFC_t,'UP',
ifelse( deg$logFC < -logFC_t,'DOWN','stable') )
)
table(deg$g)
#DOWN stable UP
# 762 21771 815
# 转成ENTREZID
deg$symbol=rownames(deg)
df <- bitr(unique(deg$symbol), fromType = "SYMBOL",
toType = c( "ENTREZID"),
OrgDb = org.Hs.eg.db)
DEG=deg
DEG=merge(DEG,df,by.y='SYMBOL',by.x='symbol')
data_all_sort <- DEG %>%
arrange(desc(logFC))
geneList = data_all_sort$logFC #把foldchange按照从大到小提取出来
names(geneList) <- data_all_sort$ENTREZID #给上面提取的foldchange加上对应上ENTREZID
head(geneList)
#1301 57214 7546 92196 389336 23657
#7.235897 6.273463 5.658879 5.341371 4.958484 4.957802
这个geneList就是符合要求的，可以直接去走 gseKEGG 函数即可做gsea分析。
结果解读
class(kk2)
#[1] "gseaResult"
#attr(,"package")
#[1] "DOSE"
# 结果保存在kk2@result
colnames(kk2@result)
# [1] "ID" "Description" "setSize"
# [4] "enrichmentScore" "NES" "pvalue"
# [7] "p.adjust" "qvalues" "rank"
# [10] "leading_edge" "core_enrichment"
kegg_result <- as.data.frame(kk2)
ID ：信号通路 Description ：信号通路的描述 setSize ：富集到该信号通路的基因个数 enrichmentScore ：富集分数，也就是ES NES ：标准化以后的ES，全称normalized enrichment score pvalue：富集的P值 p.adjust ：校正后的P值 qvalues ：FDR （false discovery rate）错误发现率 rank ：排名 core_enrichment：富集到该通路的基因列表
rownames(kk2@result)[head(order(kk2@result$enrichmentScore))]
#[1] "hsa00360" "hsa04710" "hsa00350" "hsa00650" "hsa00982" "hsa04974"
# geneSetID对应表格的id，排序后取前6个和后六个
gseaplot2(kk2, geneSetID = rownames(kk2@result)[head(order(kk2@result$enrichmentScore))])+
gseaplot2(kk2, geneSetID = rownames(kk2@result)[tail(order(kk2@result$enrichmentScore))])
# 单独通路
gseaplot2(kk2,
title = "Focal adhesion", #设置title
"hsa04510", #绘制hsa04510通路的结果
color="red", #线条颜色
base_size = 20, #基础字体的大小
subplots = 1:2, #展示上2部分
pvalue_table = T) # 显示p值
GSEA结果解读：
第一步我们需要根据基因的logFC对基因进行排序 研究的这个数据集中是否包含我们的目的基因，计算Enrich score的原则就是，从前到后依次检查基因是否是我们当前研究的数据集所包含的，如果包含就加一个正值，如果不包含就加一个负值 横坐标表示基因列表的数量 黑色的竖线代表的是我们的目的基因，已经被排好序，如果竖线聚集在头部，称为头部效应，如果在尾部，称为尾部效应 GSEA也可以进行GO，KEGG，Reactome分析，找到对应的数据集即可
可视化
还有不同的展示方式哦
山脊图
ridgeplot(kk2, 10)
气泡图
dotplot(kk2)
可以看到Y叔的clusterProfiler包基本上承包了富集分析结果的可视化， Chapter 12 Visualization of Functional Enrichment Result ，链接是：http://yulab-smu.top/clusterProfiler-book/chapter12.html
1 Bar Plot
2 Dot plot
3 Gene-Concept Network
4 Heatmap-like functional classification
5 Enrichment Map
6 UpSet Plot
7 ridgeline plot for expression distribution of GSEA result
8 running score and preranked list of GSEA result
9 pubmed trend of enriched terms
10 goplot
11 browseKEGG
12 pathview from pathview package
你最喜欢哪一种可视化方法呢？
试试看把这个gsea分析扩充到msigdb的各个通路（学徒作业）
上面我们演示的是默认的kegg数据库的gsea分析，但是生物学功能数据库还有很多。这个时候可以看看msigdbr 这个 R 语言包，它提供了对 MSigDB（Molecular Signatures Database）的直接访问。MSigDB 是一个广泛使用的基因集合注释数据库，它包含了大量关于基因集的注释信息，这些信息可以用于各种基因表达分析，尤其是在癌症生物学、免疫学和其他基因组学研究领域。
我们的学徒作业就是使用上面的geneList变量里面的基因排序结果，再结合msigdbr 这个包里面的通路，循环做一下gsea分析。
完成学徒作业就可以免费参加生信进阶培训（线下）
需要线下，必须是来广东省中山市，我提供成体系的全套生物信息学培养资源，不可能说见一面都没有就给出去的。 不限时，你可以待三五天或者一两个月，甚至更长时间 生信进阶培训内容主要是ngs多组学和单细胞，在肿瘤相关科研领域的实战
因为免费所以需要你首先有生信基础能力
我们提供的是永久免费的生信进阶培训，并不适合于从零开始的那些跨专业或者说转化的小伙伴呢，需要自己提前进行生物信息学数据分析学习过程的计算机基础知识的打磨 ，我把它粗略的分成基于R语言的统计可视化，以及基于Linux的NGS数据处理：
有生信基础并且完成了考核题就可以预约参加我们的永久免费的生信进阶培训啦，麻烦先发简历到我邮箱（jmzeng1314@163.com） 哈！
下面的练习题可以配合生物信息学入门使用，也可以作为面试考核题哦！理论上完成了下面的任务就代表你是一个合格的生信游民啦！
介绍一下线下地址
在广东省中山市的市中心（中国唯一一个用伟人名字命名的城市），一生总得来一次吧！ 房价（买房租房住酒店费用）相当于隔壁珠海的一半，广东省会城市广州的四分之一，深圳的八分之一，香港的16分之一，年轻人的友好型城市！ 目前一楼的铺位以及一个专属的走廊和大阳台还没有装修，但是欢迎大家过来帮忙规划，后面可以拿到专属的纪念册！ 办公室是早九晚七开门 办公室WiFi是1000M电信宽带，还有服务器可以使用 办公室很多书籍可以随便看！ 我每天不定时出现在办公室，有机会跟我合照！！！ 附赠大湾区旅游指南一份！ 如果是需要跟我一起吃饭，也可以半价优惠，暂定999人民币一次！ 办公室一角实拍如下所示：
还等什么呢？
赶快发简历到我邮箱（jmzeng1314@163.com）吧，或者直接添加我微信哈！在前面提到的的生信共享办公室出租，有我的联系方式哈！也可以直接过来请我吃饭，999一次，绝对是你这辈子距离明星最近的一次，而且是历史最低价！！！
