GSEA分析也是功能分析的一种,GSEA的结果图想必大家也不陌生,接下来就让小编带大家画一下炫酷的基因集富集分析图吧~GSEA富集分析可以用GO的基因集也可以用KEGG的基因集,今天来看一下GSEA-KEGG分析 ###########加载 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) 因为entrez ID 进行下一步分析比较准所以利用bitr函数转换一下ID,这里会转换###########ID转换 AML_genelist_up_l<-bitr(unique(row.names(AML_genelist_up)), fromType='ENSEMBL',toType='ENTREZID', OrgDb='org.Hs.eg.db',drop = TRUE)#转换ID 由于表达谱的基因ID一般是ensemblID,而 gseKEGG接受的ID是entrezID所以构建genelist要将entrezID、ensemblID和logFC合并
###########信息合并 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')
|