本篇文章价值超200元，为什么这么说呢，我在网上查询GSVA资料的时候看到一个课程，标价是200，23人学过。作为一个在线课程，录制完以后就可以持续的收益，不说远了，5年后还是会有人买这个课程，所以说价值远远不止200。而且我看了这个课程，里面是GEO的芯片数据作为示例，RNA-seq数据没有提到，而GSVA这个包对不同数据要求不一样。本篇文章我们分别从模拟数据，GEO芯片数据，TCGA的RNA-seq数据进行示范。所以说本篇文章价值超200应该没问题。
GSVA简介
基因集变异分析（Gene Set Variation Analysis， GSVA）是一种以非监督方式对一个简单群体评估通路活性变异的GSE方法。GSE基因集富集（gene set enrichment）方法通常被认为是一个生物信息学分析的终点，但是GSVA构成了一个构建以通路为中心的生物学模型的起点。
GSVA的目的
将一个基因表达矩阵转换成基因集表达矩阵
GSEA与GSVA的区别
相同点
无需寻找差异基因
异同点：
GSEA:计算的基本原理是扫描排序序列，当出现一个功能基因集中的基因时，就增加 ES 值， 反之， 就减少 ES 值， 所以在整个扫描过程中， ES是一个动态的值。
GSVA :基因矩阵转换成基因集矩阵。
GSVA及相关R包的安装
if(length(getOption("CRAN"))==0) options(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/")if(!require("GSVA")) BiocManager::install("GSVA")if(!require("GSEABase")) BiocManager::install("GSEABase")if(!require("GSVAdata")) BiocManager::install("GSVAdata")if(!require("pheatmap")) BiocManager::install("pheatmap")if(!require("ggplot2")) BiocManager::install("ggplot2")if(!require("tidyr")) BiocManager::install("tidyr")if(!require("dplyr")) BiocManager::install("dplyr")if(!require("limma")) BiocManager::install("limma")GSVA的运行注意事项
如果是芯片的数据是需要对数据进行过滤的，可以参考官方的文档
GSVA本身提供了四种算法，一般使用默认的算法就可以了
对于RNA-seq的数据，如果是read count可以选择Poisson分布，如果是均一化后的表达量值，可以选择默认参数高斯分布就可以了
GSVA参数
gsva(expr, gset.idx.list, annotation, method=c("gsva", "ssgsea", "zscore", "plage"), kcdf=c("Gaussian", "Poisson", "none"), abs.ranking=FALSE, min.sz=1, max.sz=Inf, parallel.sz=0, parallel.type="SOCK", mx.diff=TRUE, tau=switch(method, gsva=1, ssgsea=0.25, NA), ssgsea.norm=TRUE, verbose=TRUE)
expr
基因表达矩阵，行是基因，列是样本。
gset.idx.list
基因集列表，GMT文件
annotation
注释参数，默认情况下gsva函数会将表达矩阵和基因集的基因标识符匹配起来。
method
用于估计每个样本的基因集富集分数的方法。默认情况下，该值设置为gsva。还有ssgsea、zscore、plage。
kcdf
如果是RNA-seq的原始整数的read count 在使用gsva时需要设置kcdf="Possion",如果是取过log的RPKM,TPM等结果可以使用默认的值，所以如果我们在使用fpkm进行分析的时候使用默认参数即可。
abs.ranking
仅当mx.diff=TRUE时使用。默认abs.rank =FALSE，此时使用 modified Kuiper statistic 计算富集分数，当abs.rank =TRUE时，使用 original Kuiper statistic。
min.sz
设置基因集的最小数目
max.sz
设置基因集的最大数目
parallel.sz
设置并行计算时要使用的处理器数量。
mx.diff
提供了2种计算ES值（也称为GSVA score）的方法：mx.diff=TRUE时，ES值近似正态分布，mx.diff=FALSE时，ES值是一个双峰的分布。
tau
当 method="gsva" 时，tau=1；当 method="ssgsea" 时，tau=0.25
ssgsea.norm
当 method="ssgsea" 时，设置为 ssgsea.norm=TRUE，会对分数进行标准化，若ssgsea.norm=FALSE，则不进行标准化。
verbose
给出有关每个计算步骤的信息。默认是FALSE。
模拟数据的使用
#构造一个 30个样本，2万个基因的表达矩阵， 加上 100 个假定的基因集p <- 20000 ## number of genesn <- 30 ## number of samplesnGS <- 100 ## number of gene setsmin.sz <- 10 ## minimum gene set sizemax.sz <- 100 ## maximum gene set sizeX <- matrix(rnorm(p*n), nrow=p, dimnames=list(1:p, 1:n))dim(X)gs <- as.list(sample(min.sz:max.sz, size=nGS, replace=TRUE)) ## sample gene set sizesgs <- lapply(gs, function(n, p) sample(1:p, size=n, replace=FALSE), p) ## sample gene setses.max <- gsva(X, gs, mx.diff=FALSE, verbose=FALSE, parallel.sz=1)es.dif <- gsva(X, gs, mx.diff=TRUE, verbose=FALSE, parallel.sz=1)pheatmap::pheatmap(es.max)pheatmap::pheatmap(es.dif)
mx.diff=TRUE时，ES值近似正态分布，mx.diff=FALSE时，ES值是一个双峰的分布。
par(mfrow=c(1,2), mar=c(4, 4, 4, 1))
#绘制密度分布图plot(density(as.vector(es.max)), main="Maximum deviation from zero",xlab="GSVA score", lwd=2, las=1, xaxt="n", xlim=c(-0.75, 0.75), cex.axis=0.8)axis(1, at=seq(-0.75, 0.75, by=0.25), labels=seq(-0.75, 0.75, by=0.25), cex.axis=0.8)双峰分布
正态分布
plot(density(as.vector(es.dif)), main="Difference between largest\npositive and negative deviations",xlab="GSVA score", lwd=2, las=1, xaxt="n", xlim=c(-0.75, 0.75), cex.axis=0.8)axis(1, at=seq(-0.75, 0.75, by=0.25), labels=seq(-0.75, 0.75, by=0.25), cex.axis=0.8)dev.off()
一般会默认选择正态分布
芯片数据实操
导入数据（文末有数据获取方式）
GSE<-read.table('GSE.txt',header = T,sep = '\t')GSE[1:5,1:5]dim(GSE)
芯片数据预处理
GSE=as.matrix(GSE)rownames(GSE)=GSE[,1]exp=GSE[,2:ncol(GSE)]dimnames=list(rownames(exp),colnames(exp))mat=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)mat=avereps(mat)#Condense a microarray data object so that values for within-array replicate probes are replaced with their average.mat[1:5,1:5]dim(mat)
芯片数据标准化
mat=normalizeBetweenArrays(mat)gmt格式的基因集获取方式
登录GSEA网址，https://www.gsea-msigdb.org/gsea/index.jsp
点击箭头所示
这里需要邮箱注册，文件下载后这样的，共有2万多行（文末有获取方式）
加载gmt格式的基因集
geneSets <- getGmt("msigdb.v7.0.symbols.gmt")
GSVA分析
Result=gsva(mat, geneSets, min.sz=10, max.sz=500, verbose=TRUE, parallel.sz=1
查看，保存结果
Result[1:5,1:5]write.table(Result,file="gsvaresult.txt",sep="\t",quote=F,col.names=F)
用LIMMA包进行差异分析
logFCcutoff=0.3adjPvalueCutoff=0.05type=c( rep("con",25),rep("treat",25) )design=model.matrix(~ type)colnames(design)=c("con", "treat")fit=lmFit(Result, design)fit=eBayes(fit)dif=topTable(fit, coef="con", number=Inf,adjust.method="holm")查看，保存结果
dif[1:5,1:5]write.table(dif,file="dif.txt",sep="\t",quote=F)
火山图绘制
详见
TCGA数据分析系列之火山图colnames(dif)dif$type[(dif$adj.P.Val > 0.05|dif$adj.P.Val=="NA")|(dif$logFC < 0.3)& dif$logFC > -0.3] <- "none significant"dif$type[dif$adj.P.Val <= 0.05 & dif$logFC >= 0.3] <- "up-regulated"dif$type[dif$adj.P.Val <= 0.05 & dif$logFC <= -0.3] <- "down-regulated"
library(ggplot2)p = ggplot(dif,aes(logFC,-1*log10(adj.P.Val),color=type)) p + geom_point()
火山图美化（or 丑化）
x_lim <- max(dif$logFC,-dif$logFC) gg=p + geom_point( size=1) + xlim(-x_lim,x_lim) +labs(x="log2(FC)",y="-log10(adjP)")+ scale_color_manual(values =c("#A52A2A","grey","#f8766d"))+ geom_hline(aes(yintercept=-1*log10(0.05)),colour="black", linetype="dashed") + geom_vline(xintercept=c(-0.3,0.3),colour="black", linetype="dashed")print(gg)
热图绘制
详见
R语言学习系列之“多变的热图”筛选差异基因集矩阵
dif2 <- topTable(fit, coef="con", number=Inf, p.value=adjPvalueCutoff, adjust="holm", lfc=logFCcutoff)pheatmap_data <-Result[rownames(dif2),]
作图
annotation=c(rep("con",25),rep("treat",25))names(annotation)=colnames(Result)annotation=as.data.frame(annotation)
pheatmap(pheatmap_data, annotation=annotation, cluster_cols =F,show_rownames = F)
RNA-seq数据实操
导入数据
rm(list = ls())mRNAdata<-read.table('LIHC FPKM.txt',header = T, sep = '\t')
处理数据
row.names(mRNAdata)<-mRNAdata[,1]mRNAdata<-mRNAdata[,-1]#将数据框转换成矩阵mRNAdata = as.matrix(mRNAdata)mRNAdata[1:4,1:4]
加载GMT文件
geneSets <- getGmt("msigdb.v7.0.symbols.gmt")GSVA算分（由于数据量大，这一步比较耗时）
Result=gsva(mRNAdata, geneSets, min.sz=10, max.sz=500, verbose=T, parallel.sz=1)
Result[1:5,1:5]
根据TCGA样本名设计分组
（具体方法可见
TCGA差异分析及ggplot作图验证）
group_list=ifelse(as.numeric(substr(colnames(mRNAdata),14,15)) < 10,'tumor','normal')design <- model.matrix(~0+factor(group_list))colnames(design)=levels(factor(group_list))rownames(design)=colnames(mRNAdata)design用LIMMA包进行差异分析
fit <- lmFit(Result, design)cont.matrix=makeContrasts(contrasts=c('tumor-normal'),levels = design)fit2=contrasts.fit(fit,cont.matrix)fit2=eBayes(fit2)tempOutput = topTable(fit2, coef='tumor-normal', n=Inf)dif = na.omit(tempOutput)查看，保存结果
head(dif)write.table(dif,file="dif_RNAseq",sep="\t",quote=F)
绘制火山图
#火山图绘制colnames(dif)dif$type[(dif$adj.P.Val > 0.05|dif$adj.P.Val=="NA")|(dif$logFC < 0.5)& dif$logFC > -0.5] <- "none significant"dif$type[dif$adj.P.Val <= 0.05 & dif$logFC >= 0.5] <- "up-regulated"dif$type[dif$adj.P.Val <= 0.05 & dif$logFC <= -0.5] <- "down-regulated"
library(ggplot2)p = ggplot(dif,aes(logFC,-1*log10(adj.P.Val),color=type)) p + geom_point()
美化(or 丑化)火山图
x_lim <- max(dif$logFC,-dif$logFC) gg=p + geom_point( size=1) + xlim(-x_lim,x_lim) +labs(x="log2(FC)",y="-log10(adjP)")+ scale_color_manual(values =c("#A52A2A","grey","#f8766d"))+ geom_hline(aes(yintercept=-1*log10(0.05)),colour="black", linetype="dashed") + geom_vline(xintercept=c(-0.5,0.5),colour="black", linetype="dashed")print(gg)
绘制热图
数据准备
logFCcutoff=0.5adjPvalueCutoff=0.05dif2<-as.data.frame(dif)dif2$name<-row.names(dif2)dif2 <- filter(dif2, logFC>0.5|logFC< -0.5, adj.P.Val <= 0.05)row.names(dif2)<-dif2[,'name']
pheatmap_data <-Result[rownames(dif2),]
作图
names(group_list)=colnames(Result)group_list=as.data.frame(group_list)
pheatmap(pheatmap_data, annotation=group_list, cluster_cols =F,show_rownames=F)
下面是这篇文章的示例文件和代码
