价值超200元 GSVA基因集变异分析--一个课题的开始
本篇文章价值超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.25ssgsea.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)
下面是这篇文章的示例文件和代码