做了小鼠的RNA-seq，目的是想看下实验组和对照组小鼠在特定通路下差异情况，方法用GSVA
小鼠的分子特征数据库用msigdbr包提取，载入的包有些事用不到其实，不过怕万一报错需要别的方法所以都载入了，注释个人觉得非常详细的了哈

library(clusterProfiler)
library(tidyverse)
library(GSVA)
library(GSEABase)
library(org.Mm.eg.db)
library(msigdbr)
library(limma)
library(org.Hs.eg.db)
library(enrichplot)
library(tinyarray)


exp <- read.csv("seq_combine.csv",#读入文件
                header = T,#第一行为列明
                stringsAsFactors = F)#取消万恶的factor

rownames(exp) <- exp$Symbol#设定行名为基因名
exp <- exp[2:ncol(exp)]#删除第一列
exp <- as.matrix(exp)#转换为矩阵格式
Group <- c("Control","Control","Control",#构建分组信息
           "test","test","test")


genesets <- msigdbr(species = "Mus musculus", #设定为小鼠基因集
                    category = "C2") %>%#设定为C2基因集
  dplyr::select(.,"gs_name","gene_symbol")%>% #仅筛选保留通路名称和基因名
  as.data.frame()#设定为列表格式

#以下按照基因通路及基因名重新分组
genesets <- split(genesets$gene_symbol, genesets$gs_name)

#以下进行GSEA分析，由于是未进行标准化的数据，因此选择泊松分布
res<- GSVA::gsva(exp, genesets, 
                 method="ssgsea", 
                 kcdf="Poisson", 
                 parallel.sz=6)

write.csv(res,"all_result.csv")#导出结果

res <- as.data.frame(res)#转换结果

##下步简单但是比较关键，就是根据关键词筛选你要的通路，注意关键词为大写，筛选的结果保存在target_res 里面
target_res <- res[str_detect(rownames(res),"这里写你要通路的关键词"),]


#以下进行差异分析
design = model.matrix(~Group)#构建分组对象
fit = lmFit(target_res , design)#构建差异分析统计学模型
fit = eBayes(fit)#进行差异分析
DEG = topTable(fit, coef = 2, number = Inf)#进行差异分析结果取出
head(DEG)#查看差异结果

#以下绘制差异分析热图
draw_heatmap(target_res ,#目标通路结果
             Group,#分组信息
             cluster_cols=F, #列是否聚类
             show_rownames = T,#是否显示行名
             show_column_title = T,#显示title
             legend = T)#图注
#以下根据差异结果绘制火山图
draw_volcano(DEG,pkg = 4,logFC_cutoff = 0.5)
#导出差异结果
write.csv(DEG,"DEG.csv")