四句话代码GSEA

生信技能树

发布于 2024-11-27 15:11:36

15900

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

最近在微信群看到了一个交流，如何使用最少的代码完成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

简单的差异分析看我六年前的表达芯片的公共数据库挖掘系列推文即可哈：

很容易把差异分析结果转为类似的 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 是一个广泛使用的基因集合注释数据库，它包含了大量关于基因集的注释信息，这些信息可以用于各种基因表达分析，尤其是在癌症生物学、免疫学和其他基因组学研究领域。