三阴性乳腺癌表达数据探索笔记之GSVA分析 / 开普饭

学徒和学员已经陆续出师，是时候把生信技能树的舞台交给后辈了！

下面是学徒写的《GEO数据挖掘课程》的配套笔记（第6篇）

除了GSEA之外，还有很多其他类型的分析方式，用到的统计学原理有差别。如GSVA，SSGSEA， PGSEA

GSVA与GSEA的差别在于，这种方法不需要对基因进行排序，因此也意味着不需要首先进行其他的统计学分析，如基因在样本之间的表达差异，如变化倍数，然后根据变化值从高到低进行排序。只需要样本内基因的排序，每个样本内部可以根据基因表达的count值来进行排序，从而在样本内部是否有基因富集。针对每个样本进行分析。

数据准备：

表达矩阵，需要进行ID转换，需要SYMBOL号,这根据下载的数据集类型，和GSEA用到的数据集,从MSigDB 下载
需要分组信息
基因集(gene_list)

第一步：表达矩阵的探针名转换为SYMBOL

load(file = 'step1-output.Rdata') # 每次都要检测数据 dat[1:4,1:4] library(hgu133plus2.db) ids=toTable(hgu133plus2SYMBOL) #通过看hgu133plus2.db这个包的说明书知道提取probe_id（探针名）和symbol（基因名）的对应关系的表达矩阵的函数为toTable head(ids) dat=dat[ids$probe_id,] dat[1:4,1:4] ids$median=apply(dat,1,median) #对dat这个矩阵按行操作，取每一行的中位数，将结果添加到ids矩阵median列 ids=ids[order(ids$symbol,ids$median,decreasing = T),] #对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids ids=ids[!duplicated(ids$symbol),] #将symbol这一列去除重复项 dat=dat[ids$probe_id,] rownames(dat)=ids$symbol dat[1:4,1:4]

第二步：进行GSVA分析，获得GSVA得分矩阵

X=dat table(group_list) ## Molecular Signatures Database (MSigDb) d='D:/bioinformation/曾老师任务/GSE76275-TNBC-master/MSigDB/symbols/' gmts=list.files(d,pattern = 'all') gmts library(ggplot2) library(clusterProfiler) library(org.Hs.eg.db) library(GSVA) # BiocManager::install('GSVA') es_max <- lapply(gmts, function(gmtfile){ #gmtfile=gmts[8];gmtfile geneset <- getGmt(file.path(d,gmtfile)) es.max <- gsva(X, geneset, mx.diff=FALSE, verbose=FALSE, parallel.sz=1) return(es.max) }) #es_max是一个长度为8的list,其中的每个元素都是基因集数据库对应的GSVA得分矩阵

第三步：对GSVA得分矩阵分别进行差异分析

adjPvalueCutoff <- 0.001 logFCcutoff <- log2(2) es_deg <- lapply(es_max, function(es.max){ #es.max <- es_max[[2]] table(group_list) dim(es.max) design <- model.matrix(~0+factor(group_list)) colnames(design)=levels(factor(group_list)) rownames(design)=colnames(es.max) design library(limma) contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"), levels = design) contrast.matrix<-makeContrasts("TNBC-noTNBC", levels = design) contrast.matrix ##这个矩阵声明，我们要把TNBC组跟noTNBC进行差异分析比较 deg = function(es.max,design,contrast.matrix){ ##step1 fit <- lmFit(es.max,design) ##step2 fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果 fit2 <- eBayes(fit2) ## default no trend !!! ##eBayes() with trend=TRUE ##step3 res <- decideTests(fit2, p.value=adjPvalueCutoff) summary(res) tempOutput = topTable(fit2, coef=1, n=Inf) nrDEG = na.omit(tempOutput) #write.csv(nrDEG2,"limma_notrend.results.csv",quote = F) head(nrDEG) return(nrDEG) } re = deg(es.max,design,contrast.matrix) nrDEG=re head(nrDEG) return(nrDEG) }) #es_deg是一个列表，包含全部GSVA得分矩阵的差异分析矩阵，与es_max对应

第四步：绘制热图展示差异显著的通路

library(pheatmap) lapply(1:length(es_deg), function(i){ # i=2 print(i) dat=es_max[[i]] df=es_deg[[i]] df=df[df$P.Value<0.01 & abs(df$logFC) > 0.5,] #筛选出差异比较显著的通路 print(dim(df)) if(nrow(df)>5){ n=rownames(df) dat=dat[match(n,rownames(dat)),] ac=data.frame(g=group_list) rownames(ac)=colnames(dat) rownames(dat)=substring(rownames(dat),1,50) pheatmap::pheatmap(dat, fontsize_row = 8,height = 11, annotation_col = ac,show_colnames = F, filename = paste0('gsva_',strsplit(gmts[i],'[.]')[[1]][1],'.pdf')) } }) #以上的代码会绘制出8个数据集的全部热图 adjPvalueCutoff <- 0.001 logFCcutoff <- log2(2) df=do.call(rbind ,es_deg) es_matrix=do.call(rbind ,es_max) df=df[df$P.Value<0.01 & abs(df$logFC) > 0.5,] write.csv(df,file = 'GSVA_DEG.csv') #将所有GSVA的得分差异显著的结果保存为一个csv,便于检查

结果解读：

GSVA对数据库中的每一个通路在每个样本中算了一个值，相当于GSEA的enrichment score, 如果得分越高，说明这个通路在该样本中被改变的越严重。

得到的GSVA得分矩阵可以用来做差异分析，看哪些通路在两个分组中存在差异，类似于基因表达差异分析。

这些分析，基本上读一下我五年前在生信技能树的表达芯片的公共数据库挖掘系列推文 就明白了；

视频观看方式

我把3年前的收费视频课程：3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月，免费到B站:

这个课程超级棒，B站免费学习咯：https://m.bilibili.com/video/BV1dy4y1C7jz
配套代码在GitHub哈：https://github.com/jmzeng1314/GSE76275-TNBC
TCGA数据库挖掘，代码在：https://github.com/jmzeng1314/TCGA_BRCA
GTEx数据库挖掘，代码在：https://github.com/jmzeng1314/gtex_BRCA
METABRIC数据库挖掘，代码在：https://github.com/jmzeng1314/METABRIC

然后马上就有了3千多学习量，而且有学员给出来了图文并茂版本万字笔记，让我非常感动！

扫描下面二维码马上就可以学习起来啦，笔记需要至少半个小时来阅读哦！

三阴性乳腺癌表达数据探索笔记之GSVA分析

相关推荐