三阴性乳腺癌表达数据探索笔记之GSVA分析
时间:2022-07-28
本文章向大家介绍三阴性乳腺癌表达数据探索笔记之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]
ID转换结果
- 第二步:进行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得分矩阵
- 第三步:对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对应
GSVA得分差异分析矩阵
- 第四步:绘制热图展示差异显著的通路
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得分差异矩阵热图
- 结果解读:
GSVA对数据库中的每一个通路在每个样本中算了一个值,相当于GSEA的enrichment score, 如果得分越高,说明这个通路在该样本中被改变的越严重。
得到的GSVA得分矩阵可以用来做差异分析,看哪些通路在两个分组中存在差异,类似于基因表达差异分析。
这些分析,基本上读一下我五年前在生信技能树的表达芯片的公共数据库挖掘系列推文 就明白了;
- 解读GEO数据存放规律及下载,一文就够
- 解读SRA数据库规律一文就够
- 从GEO数据库下载得到表达矩阵 一文就够
- GSEA分析一文就够(单机版+R语言版)
- 根据分组信息做差异分析- 这个一文不够的
- 差异分析得到的结果注释一文就够
视频观看方式
我把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千多学习量,而且有学员给出来了图文并茂版本万字笔记,让我非常感动!
- Android+struts2+json方式模拟手机登录功能
- iOS 获取通讯录里边的电话号码AddressBook
- InvocationTargetException异常解析
- Spring Cloud构建微服务架构:服务注册与发现(Eureka、Consul)【Dalston版】
- java基础多线程之共享数据
- Spring Boot自动化配置的利弊及解决之道
- Java四种引用解析以及在Android的应用
- java基础之泛型
- java基础之反射
- 第四章 正则表达式回溯法原理
- Spring Boot属性配置文件详解
- StickyListHeaders的使用
- Android自绘跑马灯控件,你会了吗?
- HashMap原理解析
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法