做了小鼠的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")
网友评论