GSVA：pathway级别的差异分析

GSVA其实就是pathway级别的差异分析
标准差异分析通常是不够的，定位到成百上千个有统计学显著变化的差异表达基因后，同样是有成百上千个生物学功能注释（GO功能和KEGG通路），普通的超几何分布检验已经不能满足大家多元化的分析了，所以就有了大家耳熟能详的GSEA分析，以及绝大部分人比较陌生的GSVA分析。
GSVA分析的文章发表于2013年，GSVA: gene set variation analysis for microarray and RNA-Seq data 同样是broad 研究生出品，其在2005年PNAS发表的gsea已经高达1.4万的引用了，不过这个GSVA才不到300。去年我就介绍过一波它的分析流程，在：使用GSVA方法计算某基因集在各个样本的表现非常简单的代码，所以各个培训机构，公司人员都开始学习和二次创作进而分享。考虑到最近邮箱接收的GSVA提问比较多，我这里还是得再次归纳总结一波，这次我准备从GSVA其实就是pathway级别的差异分析的角度来分享。
一般人做数据挖掘，到差异基因的生物学功能注释（GO功能和KEGG通路）就结束了，进而也就是去使用一些网页工具，比如string，出一些花花绿绿的图表，比如PPI网络图。实际上，使用了GSVA，可以把成百上千个生物学功能注释（GO功能和KEGG通路）转换为新的表达矩阵，就是具体的每个通路在各个样本的基因集变异分析（Gene Set Variation Analysis，GSVA）值，我们把它当作一般的矩阵文件，进行差异表达分析，热图绘制，火山图绘制。
下面我们以文献 Metabolic remodeling contributes towards an immune‐suppressive phenotype in glioblastoma 为例，欣赏它的两个图表，文章发表在Cancer Immunology, Immunotherapy (2019)
https://doi.org/10.1007/s00262-019-02347-3
基于PATHWAY的热图
虽然作者这里使用的代谢组学数据：

Global metabolomic profiling was performed on patient-derived glioblastoma (GBM; n=80) and LGA (n=28) tumor samples using LG/GC–MS.
本质上仍然是记录表达量。
![image-20190928112733036]http://www.bio-info-trainee.com/wp-content/uploads/2019/12/image-20190928112733036.png)

基于PATHWAY的火山图
![image-20190928112802174]http://www.bio-info-trainee.com/wp-content/uploads/2019/12/image-20190928112802174.png)

PATHWAY的具体含义
pathway在我这里是基因集的别名，其中msigdb有着丰富的基因集，MSigDB（Molecular Signatures Database）数据库中定义了已知的基因集合：http://software.broadinstitute.org/gsea/msigdb 包括H和C1-C7八个系列（Collection），每个系列分别是：

H: hallmark gene sets （癌症）特征基因集合，共50组，最常用；
C1: positional gene sets 位置基因集合，根据染色体位置，共326个，用的很少；
C2: curated gene sets：（专家）校验基因集合，基于通路、文献等：
C3: motif gene sets：模式基因集合，主要包括microRNA和转录因子靶基因两部分
C4: computational gene sets：计算基因集合，通过挖掘癌症相关芯片数据定义的基因集合；
C5: GO gene sets：Gene Ontology 基因本体论，包括BP（生物学过程biological process，细胞原件cellular component和分子功能molecular function三部分）；
C6: oncogenic signatures：癌症特征基因集合，大部分来源于NCBI GEO 发表芯片数据；
C7: immunologic signatures: 免疫相关基因集合。

##### GSVA个性化分析 #####
library(tidyverse)
library(ggplot2)
library(clusterProfiler)
library(GSEABase)
library(GSVA)
# setwd('/public/igenebook/GRCh38_p13_AJRS2200903005_20210414/personal_analysis/heatmap')##创建gmt文件转list对象函数
gmt2list <- function(gmtfile){sets <- as.list(read_lines(gmtfile))for(i in 1:length(sets)){tmp = str_split(sets[[i]], '\t')n = length(tmp[[1]])names(sets)[i] = tmp[[1]][1]sets[[i]] = tmp[[1]][3:n]rm(tmp, n)}return(sets)
}
#读取基因集数据库（gmts数据集:从GSEA官网MsigDB下载）
#读入gmt文件，这个可以从MSigDB上下载，这边选的上gene symbol根据自己的data来选择
gmt_file="/public/database/MSigDB/msigdb_v7.4_GMTs/c2.cp.kegg.v7.4.symbols.gmt"
geneset = gmt2list(gmt_file)
geneset <- getGmt(gmt_file)
geneset2 <- read.gmt(gmt_file)#读入exp文件
exp <- read.table('/public/igenebook/GRCh38_p13_AJRS2200903005_20210414/AJRS2200903005/Result/03.Exp/All_sample_counts.xls',row.names = 1,header = T,sep = '\t')# gsva分析
es <- gsva(as.matrix(exp), geneset,verbose=TRUE)
# es <- gsva(as.matrix(exp), TERM2GENE=geneset2,verbose=TRUE)
head(es)
dim(es)library(limma)
adjPvalueCutoff <- 0.1
logFCcutoff <- 1
group_list = c(rep("U251NC",3),rep("U251SH1",3))
design <- model.matrix(~ factor(group_list)) # ctr在前,treatment在后
colnames(design)=levels(factor(group_list))
row.names(design)<-colnames(es)contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix<-makeContrasts("Trt-ck",levels = design)
# contrast.matrix ##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较##step1
fit <- lmFit(es, design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)fit <- eBayes(fit2)
allGeneSets <- topTable(fit, coef="ctr_vs_treatment", number=Inf)
DEgeneSets <- topTable(fit, coef="ctr_vs_treatment", number=Inf,p.value=adjPvalueCutoff, adjust="BH")
nrDEG = na.omit(DEgeneSets)
es_heatmap <- es[row.names(es) %in% row.names(nrDEG),]library(pheatmap)
library(patchwork)
#绘制热图
annotation_col = data.frame(Type = factor(rep(c("U251NC", "U251SH1"), each=3))#Time = 1:3
)
rownames(annotation_col) = c(paste("U251NC",1:3, sep = "_"), paste("U251SH1", 1:3, sep = "_"))p <- pheatmap(es_heatmap, show_rownames=1, show_colnames = F,annotation_col = annotation_col,fontsize=9, width=15, height=12)
png('GSVA.png')
p
dev.new()
#也可绘制火山图

GSVA：pathway级别的差异分析相关推荐

单细胞测序GSVA及下游limma差异分析
1.GSVA/基因集变异分析定义:将分析的功能单元从基因向基因集进行改变,进行基因集(通路)级别的差异分析. 2.分析原理:将基因在不同样本间的表达矩阵(列为样本,行为基因名)转化成基因集在样本间的表 ...
在线作图|在线做Metastats组间差异分析
Metastats组间差异分析 Metastats组间差异分析与LEfSe类似,也是多用于寻找不同区组中微生物群落差异物种的一种工具.Metastats实际上是非参数多重检验和p值校正的整合,它对所有 ...
游戏全球化美术风格差异分析
日前,在2021年GDC期间,欧美出海整合营销公司NewCode联合阿里云,举办了一年一度的游戏出海市场营销峰会,邀请了业内包括IGN.Sensor Tower.Share Creators等各领域的 ...
多组差异分析的可视化，这样做最省心！
欢迎关注"生信修炼手册"! 下图是来自文献中的一个示例图,展示了不同亚型之间risk score值的差异分析结果从图形的主题可以看出,这是一幅基于ggplot2绘制的图片,但是如 ...
生信文献 | KLF7:高级别浆液性卵巢癌新的候选生物标志物和治疗靶点
这是一篇biomarker的生信文章,,也有很多实验验证. 癌症类型:高级别浆液性卵巢癌(HGSOC) 问题:发病机制不清楚,高异质性和化疗后复发耐药解决问题:需要进一步研究参与肿瘤进展的分子机制, ...
GSVA+limma差异通路分析+发散条形图
Nat Med 的图是真的好看啊 GSVA+limma 差异通路的分析,主要参考:跟着 Nat Med. 学作图 | GSVA+limma差异通路分析+发散条形图,过程略去不谈.主要想说一下,其中之前 ...
【生信分析】Pathway Network Visualizer (PaNeV)
PANEV:一个基于路径的网络可视化的R包 PANEV: an R package for a pathway-based network visualization 背景实现包安装和功能数据准 ...
跟着Cell学单细胞转录组分析(十三):单细胞GSVA分析|这个包涵盖大多数物种
之前我们发过GSVA分析(有了这个包,猪的GSEA和GSVA分析也不在话下(第一集),[后续来了]有了这个包,猪的GSEA和GSVA分析也不在话下(第二集)),接着单细胞系列,重新说一下GSVA分析. ...
有了这个包，猪的GSEA和GSVA分析也不在话下（第一集）
救救孩子吧,GSVA分析都是做人的,有现成的人的数据集,可是其他物种的就惨了,很难下手! 今天我们就说说小鼠,也是常见物种的GSVA分析,结合单细胞的数据! GSVA的作用不用多说了,大家都熟悉,至少 ...

GSVA：pathway级别的差异分析

GSVA：pathway级别的差异分析相关推荐

最新文章

热门文章