首先感谢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等等，需要什么下载什么
图片.png
-
第五步，读取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值出点图
默认p.ajust
p值
dotplot(KEGG,split=".sign")+facet_grid(~.sign) #出点图，并且分面激活和抑制
图片.png
dotplot(KEGG,split=".sign")+facet_wrap(~.sign,scales = "free") #换个显示方式
图片.png
library(enrichplot)
#特定通路作图
gseaplot2(KEGG,1,color="red",pvalue_table = T) # 按第一个做二维码图，并显示p值
图片.png
gseaplot2(KEGG,1:10,color="red") #按第一到第十个出图，不显示p值
图片.png
网友评论