参考下面链接作者的方法：
https://www.jianshu.com/p/bb7442bb8cad
自己程序整理的代码：用自己的标志将训练集的数据由表达谱分为高低风险组，然后DEGS，最后进行GSEA分析，从Msigdb数据库下载的数据为c2.cp.reactome.v7.4.entrez.gmt
setwd("E:/1.RStudio/geo_ovarian_GPL96_570/geo_array_ovarian/OVCP3")
Exp_low_train<-read.table("./Exp_low_train.txt",header=T,stringsAsFactors = F,sep = "\t")
Exp_High_train<-read.table("./Exp_High_train.txt",header=T,stringsAsFactors = F,sep = "\t")
DEG<-cbind(Exp_low_train,Exp_High_train)
group <-rep(c('low','High'), c(297,264) )#C1是i类组，C2是2类组
#3. 加载limma包，构建表达矩阵
library(limma)
design <- model.matrix(~0+factor(group))
colnames(design)=levels(factor(group))
rownames(design)=colnames(DEG)
#4. 规定哪一组数据与哪一组数据比较
contrast.matrix<-makeContrasts("low-High",levels=design) #比较design矩阵中HC与MDD组
#5. 获取差异表达数据
##step1
fit <- lmFit(DEG,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#加change列,标记上下调基因
logFC_t=0#设置logFC的阈值
P.Value_t = 0.05#设置P.Val的阈值
k1 = (nrDEG$adj.P.Val < P.Value_t)&(nrDEG$logFC < -logFC_t)#下调基因
k2 = (nrDEG$adj.P.Val < P.Value_t)&(nrDEG$logFC > logFC_t)#上调基因
change = ifelse(k1,
"down",
ifelse(k2,"up","stable"))#标记上下调基因
table(change)
deg <- mutate(nrDEG,change)#添加到deg表格中
#过滤出上下调的基因作为背景基因
library(dplyr)
deg %>% filter(!change=='stable') -> deg1
#6. 导出差异表达数据
write.table(deg1,"Train_limma_deg.txt",sep="\t",quote=F,row.names=T)
######################基于Msigdb数据库的GSEA分析
gene<-read.table("./Msigdb.txt",header=T,stringsAsFactors = F,sep = "\t")
library(clusterProfiler)
library(org.Hs.eg.db)
library(stringr)
load('all_intersect_data\\tcga_data.RData')
entrezid<-tcga_GID0[match(gene[,1],row.names(tcga_Exp0))]
logfc<-gene[,2]
gene_df<-as.data.frame(cbind(entrezid,logfc))
geneList<-gene_df $logfc #第二列可以是folodchange，也可以是logFC
names(geneList)=gene_df $entrezid #使用转换好的ID
geneList=sort(geneList,decreasing = T) #从高到低排序
reacmt<-read.gmt("c2.cp.reactome.v7.4.entrez.gmt") #读gmt文件
reac<-GSEA(geneList,TERM2GENE = reacmt) #GSEA分析
library(ggplot2)
dotplot(KEGG) #出点图
dotplot(KEGG,color="pvalue") #按p值出点图
dotplot(reac,split=".sign")+facet_grid(~.sign) #出点图，并且分面激活和抑制
dotplot(reac,split=".sign")+facet_wrap(~.sign,scales = "free") #换个显示方式
library(enrichplot)
#特定通路作图
gseaplot2(reac,1,color="red",pvalue_table = T) # 按第一个做二维码图，并显示p值
gseaplot2(reac,1:10,color="red") #按第一到第十个出图，不显示p值
MSigDB:基因集数据库的介绍：
https://www.jianshu.com/p/f2929d2e3f2e
查询某一通路中的所有基因：
https://www.keyangou.com/topic/1167
