首先感謝Y叔的clusterprofiler神包轩猩，做富集分析優(yōu)點(diǎn)是在線爬取數(shù)據(jù)弥奸，結(jié)果很可信拳亿，但是缺點(diǎn)也是網(wǎng)絡(luò)問(wèn)題博杖，網(wǎng)絡(luò)差點(diǎn)就要等很久，不過(guò)GSEA有自帶GMT文件土铺，因此下載好離線數(shù)據(jù)雷激，這些就可以擺脫在線的問(wèn)題睡汹，單機(jī)就可以操作GSEA了
GSEA也就是基因集富集分析理郑，不需要分析是上調(diào)基因還是下調(diào)基因蹄溉，分析的是所有基因，所以結(jié)果應(yīng)該更可靠您炉，根據(jù)官方操作說(shuō)明柒爵， 需要有兩組數(shù)據(jù)，第一組是表達(dá)矩陣赚爵，第二組是分組棉胀，然后用軟件操作，但是至始至終出的圖奇丑無(wú)比冀膝，而且還只能是png格式唁奢，遠(yuǎn)達(dá)不到300dpi的發(fā)表級(jí)，雖然目前有很多再次作圖方法畸写，但依然比較曲折，Y叔的包解決了這個(gè)問(wèn)題氓扛，下面分享一下教程
- 第一步：數(shù)據(jù)格式
表格需要兩列枯芬， 第一列是基因名论笔，第二列是Foldchange或者Log2FC
以前我一直不明顯第二列是如何來(lái)的，明明是矩陣和分組千所，怎么算Fc狂魔，但是也怪自己學(xué)藝不精，其實(shí)有列矩陣和分組就可以算出來(lái)了淫痰，懂R的可以用各種差異分析包如“DESq”最楷，“l(fā)imma”和“edgeR”算出來(lái)，不懂R的用Excel也可以待错，可以看我之前的教程籽孙，其實(shí)Foldchange就是兩組數(shù)據(jù)的均值之比，也就是先算出各組的平均值火俄，然后拿比較組除以對(duì)照組就算出來(lái)了 (沒(méi)事還是要多讀書(shū)啊) - 第二步犯建，如基因名是symbol，需要基因ID轉(zhuǎn)換為entrezid格式
gene<-read.csv("你的文件.csv") #csv可讀性較好瓜客，帶表頭适瓦，需要有一列是symbol
library(clusterProfiler)
library(org.Hs.eg.db)
library(stringr)
gene<-str_trim(d$symbol,"both") #定義gene
#開(kāi)始ID轉(zhuǎn)換
gene=bitr(gene,fromType="SYMBOL",toType="ENTREZID",OrgDb="org.Hs.eg.db") #會(huì)有部分基因數(shù)據(jù)丟失，或者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")
- 第三步谱仪，定義基因列表玻熙，對(duì)logFC進(jìn)行從高到低排序
geneList<-gene_df $logFC #第二列可以是folodchange，也可以是logFC
names(geneList)=gene_df $ENTREZID #使用轉(zhuǎn)換好的ID
geneList=sort(geneList,decreasing = T) #從高到低排序
-
第四步疯攒，下載離線GMT文件
https://www.gsea-msigdb.org/gsea/downloads.jsp （需要先注冊(cè)嗦随，不過(guò)一個(gè)郵箱地址即可）
包含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) #出點(diǎn)圖
dotplot(KEGG,color="pvalue") #按p值出點(diǎn)圖
默認(rèn)p.ajust
p值
dotplot(KEGG,split=".sign")+facet_grid(~.sign) #出點(diǎn)圖呜舒，并且分面激活和抑制
圖片.png
dotplot(KEGG,split=".sign")+facet_wrap(~.sign,scales = "free") #換個(gè)顯示方式
圖片.png
library(enrichplot)
#特定通路作圖
gseaplot2(KEGG,1,color="red",pvalue_table = T) # 按第一個(gè)做二維碼圖锭汛，并顯示p值
圖片.png
gseaplot2(KEGG,1:10,color="red") #按第一到第十個(gè)出圖，不顯示p值
圖片.png
