首先感谢Y叔的clusterprofiler
神包,做富集分析优点是在线爬取数据,结果很可信,但是缺点也是网络问题,网络差点就要等很久,不过GSEA有自带GMT文件,因此下载好离线数据,这些就可以摆脱在线的问题,单机就可以操作GSEA了
GSEA
也就是基因集富集分析,不需要分析是上调基因还是下调基因,分析的是所有基因,所以结果应该更可靠,根据官方操作说明, 需要有两组数据,第一组是表达矩阵,第二组是分组,然后用软件操作,但是至始至终出的图奇丑无比,而且还只能是png格式,远达不到300dpi的发表级,虽然目前有很多再次作图方法,但依然比较曲折,Y叔的包解决了这个问题,下面分享一下教程
- 第一步:数据格式
表格需要两列, 第一列是基因名
,第二列是Foldchange
或者Log2FC
以前我一直不明显第二列是如何来的,明明是矩阵和分组,怎么算Fc,但是也怪自己学艺不精,其实有列矩阵和分组就可以算出来了,懂R的可以用各种差异分析包如“DESq”,“limma”和“edgeR”算出来,不懂R的用Excel也可以,可以看我之前的教程,其实Foldchange就是两组数据的均值之比,也就是先算出各组的平均值,然后拿比较组除以对照组就算出来了 (没事还是要多读书啊) - 第二步,如基因名是symbol,需要基因ID转换为entrezid格式
gene<-read.csv("你的文件.csv") #csv可读性较好,带表头,需要有一列是symbol
library(clusterProfiler)
library(org.Hs.eg.db)
library(stringr)
gene<-str_trim(d$symbol,"both") #定义gene
#开始ID转换
gene=bitr(gene,fromType="SYMBOL",toType="ENTREZID",OrgDb="org.Hs.eg.db") #会有部分基因数据丢失,或者ENSEMBL
## 去重
gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)
gene_df <- data.frame(logFC=gene$logFC, #可以是foldchange
SYMBOL = gene$symbol) #记住你的基因表头名字
gene_df <- merge(gene_df,gene,by="SYMBOL")
- 第三步,定义基因列表,对logFC进行从高到低排序
geneList<-gene_df $logFC #第二列可以是folodchange,也可以是logFC
names(geneList)=gene_df $ENTREZID #使用转换好的ID
geneList=sort(geneList,decreasing = T) #从高到低排序
-
第四步,下载离线GMT文件
https://www.gsea-msigdb.org/gsea/downloads.jsp (需要先注册,不过一个邮箱地址即可)
包含GO
,KEGG
,REACTOME
,HALLMARKS
等等,需要什么下载什么
第五步,读取GMT格式,直接GSEA分析并出图,以Kegg为例
kegmt<-read.gmt("c2.cp.kegg.v7.1.entrez.gmt") #读gmt文件
KEGG<-GSEA(geneList,TERM2GENE = kegmt) #GSEA分析
library(ggplot2)
dotplot(KEGG) #出点图
dotplot(KEGG,color="pvalue") #按p值出点图
dotplot(KEGG,split=".sign")+facet_grid(~.sign) #出点图,并且分面激活和抑制
dotplot(KEGG,split=".sign")+facet_wrap(~.sign,scales = "free") #换个显示方式
library(enrichplot)
#特定通路作图
gseaplot2(KEGG,1,color="red",pvalue_table = T) # 按第一个做二维码图,并显示p值
gseaplot2(KEGG,1:10,color="red") #按第一到第十个出图,不显示p值