扫码联系客服
###########加载
library(topGO)
library(enrichplot)
library(ggplot2)
library(org.Hs.eg.db)#人类基因组注释相关的包
library(DO.db)
library(clusterProfiler)###########差异分析结果
AML_diff<-read.table("edgerAML.xls",sep="\t",quote="",header= T,row.names = 1)
AML_genelist_up<-AML_diff[(which(AML_diff$logFC>0)),]
head(AML_diff)
head(AML_genelist_up)
###########ID转换
AML_genelist_up_l<-bitr(unique(row.names(AML_genelist_up)),
fromType="ENSEMBL",toType="ENTREZID",
OrgDb="org.Hs.eg.db",drop = TRUE)#转换ID
###########信息合并
AML_genelist_up<-as.data.frame(cbind(row.names(AML_genelist_up),AML_genelist_up$logFC))
names(AML_genelist_up)<-c("ENSEMBL","logFC")
AML_up<-merge(AML_genelist_up_l,AML_genelist_up,all=F)#将entrezID、ensemblID和logFC合并
AML_up<-AML_up[,-1]#去掉ensemblID只保留entrezID
gseKEGG的输入必须是排序后的geneList；需要两列：命名(每一个数字都有一个对应的名字，就是相应的基因ID了)；排序(是一串数字，数字是从高到低排序的)
###########排序
geneList<-as.numeric(as.character(AML_up[,2]))
names(geneList) = as.character(AML_up[,1])
geneList= sort(geneList, decreasing = TRUE)#构建geneList,并根据logFC由高到低排列
###########GSEA分析
AML_GSEA_KEGG_up<-gseKEGG(
geneList =geneList,
nPerm = 1000,#
keyType = 'kegg',#可以选择"kegg",'ncbi-geneid', 'ncib-proteinid' and 'uniprot'
organism = 'hsa'#定义物种,
#pvalueCutoff = 0.05, #自定义pvalue的范围
#pAdjustMethod = "BH" #校正p值的方法
)AML_GSEA_KEGG_up$Description#富集到那些基因集上
AML_GSEA_KEGG_up$enrichmentScore#富集得分
#根据enrichmentScore对GSEA的结果进行排序
sortAML_GSEA_KEGG_up<-AML_GSEA_KEGG_up[order(AML_GSEA_KEGG_up$enrichmentScore,decreasing=T)]
###########画图
gseaplot2(AML_GSEA_KEGG_up,row.names(sortAML_GSEA_KEGG_up))
###########美化
gseaplot2(AML_GSEA_KEGG_up,row.names(sortAML_GSEA_KEGG_up),#只显示前三个GSEA的结果
title="AML_GSEA_KEGG_up",#标题
color = c("#626262","#8989FF","#FF0404"),#颜色
pvalue_table = FALSE,
ES_geom = "line"#enrichment scored的展现方式 'line' or 'dot')
欢迎关注生信人
TCGA | 小工具 | 数据库 |组装| 注释 | 基因家族 | Pvalue
基因预测 |bestorf | sci | NAR | 在线工具 | 生存分析 | 热图
生信不死 | 初学者 | circRNA | 一箭画心| 十二生肖 | circos
舞台|基因组 | 黄金测序 | 套路 | 杂谈组装 | 进化 | 测序简史
