最近帮师弟做一个GSVA分析，发现已经忘的差不多了，顺便整理一下吧。
基因集变异分析（Gene Set Variation Analysis, GSVA）可以将一个基因表达矩阵转换成基因集表达矩阵，并对基因集进行差异分析；并且进一步对使用表达矩阵中目的基因的表达数据与基因集表达举证做相关性分析，初步探究目的基因的可能相关通路。
1.曾老师标准操作-，-！
rm(list = ls())
options(stringsAsFactors = F)
2.加载需要用到的包
library(devtools)
library(AnnoProbe)
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
library(GSVA)
library(GSEABase)
3.数据预处理
load(file = 'countdata_grouplist.Rdata') ##加载数据，包括表达矩阵与样本分组信息
exprSet<-countdata
dim(exprSet)
exprSet=exprSet[apply(exprSet,1, function(x) sum(x>1) > 5),]
dat=log2(edgeR::cpm(exprSet)+1)
4.通路基因集的获取
首先需要下载基因集Gmt格式文件，我的基因集是通过msigdb官网下载的，这里使用了C2目录下的基因集，共2871个基因集。
geneSets<-getGmt("c2.cp.v7.2.symbols.gmt")#下载好的gmt文件放在R语言当前工作路径下即可通过getGmt函数直接加载
5.GSVA分析
Result=gsva(dat,
geneSets,
min.sz=10,
max.sz=500,
verbose=TRUE,
mx.diff=F,
parallel.sz=1)
dim(Result)
###函数中的参数都代表什么意义可以通过help函数自行查看
###环境变量中得到的Result即为基因集表达矩阵
6.对基因集表达矩阵进行差异分析，这里使用limma包完成
library(limma)
design<-model.matrix(~0+factor(group_list))
colnames(design)<-levels(factor(group_list))
rownames(design)<-colnames(Result)
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix<-makeContrasts("Normal-OA",levels = design)
fit<-lmFit(Result,design)
fit2<-contrasts.fit(fit,contrast.matrix)
fit2<-eBayes(fit2)
tempOutput<-topTable(fit2,coef = 1,n=Inf)
nrDEG<-na.omit(tempOutput)##nrDEG即为基因集差异表达结果
View(nrDEG)
7.对差异表达结果完成火山图绘制
#数据处理
df<-nrDEG
df$v<- -log10(df$P.Value)
df$g<-ifelse(df$P.Value>0.01,'stable',ifelse(df$logFC>0.3,'up',ifelse(df$logFC<(-0.3),'down','stable')))##对上调、下调及无变化基因集进行标注
table(df$g)
#画图
library(ggplot2)
g = ggplot(data=df,
aes(x=logFC, y=-log10(P.Value),
color=g)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab("log2 fold change") + ylab("-log10 p-value") +
theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red'))
print(g)
8.目的基因与通路的相关性分析
#这里使用表达矩阵中的基因TNFSF15进行演示，看看TNFSF15的表达与那些通路存在相关性。
KDM3A<-as.data.frame(dat["TNFSF15",])
KDM3A<-t(TNFSF15)
rownames(TNFSF15)<-"TNFSF15"
cordat<-rbind(Result,TNFSF15)#将TNFSF15的表达数据与基因集的整合在一起
##包装batch_cor函数
batch_cor <- function(gene){
y <- as.numeric(cordat[gene,])
rownames <- rownames(cordat)
do.call(rbind,future_lapply(rownames, function(x){
dd <- cor.test(as.numeric(cordat[x,]),y,type="spearman")
data.frame(gene=gene,mRNAs=x,cor=dd$estimate,p.value=dd$p.value )
}))
}
#开线程，完成相关性分析
library(future.apply)
plan(multiprocess)
system.time(dd <- batch_cor("TNFSF15"))
##保存结果dd,用于画图分析
save(dd,file="TNFSF15_cordata.Rdata")
9.使用棒棒图展示相关性分析结果
dd$abscor<-abs(dd$cor)
dd = dd[order(dd[,3]),]##按相关性数值排序
##提取作图数据，正相关前十与负相关前十
tmp_1<-dd[1:10,]
tmp_2<-dd[2313:2323,]
tmp<-rbind(tmp_1,tmp_2)
##画图
library(ggpubr)
p<-ggdotchart(tmp, x = "mRNAs", y = "cor",
color = "p.value",
sorting = "descending",
add = "segments",
add.params = list(color = "lightgray", size = 1.5),
rotate = TRUE,
#group = "cyl",
dot.size = "abscor",
ggtheme = theme_pubr(),
xlab=""
)
print(p)
p<-p + gradient_color(c("blue","red")) + grids(linetype = "dashed")
p + grids(linetype = "dashed")
网友评论