下面是学徒写的《GEO数据挖掘课程》的配套笔记(第6篇)
除了GSEA之外,还有很多其他类型的分析方式,用到的统计学原理有差别。如GSVA,SSGSEA, PGSEA
GSVA与GSEA的差别在于,这种方法不需要对基因进行排序,因此也意味着不需要首先进行其他的统计学分析,如基因在样本之间的表达差异,如变化倍数,然后根据变化值从高到低进行排序。只需要样本内基因的排序,每个样本内部可以根据基因表达的count值来进行排序,从而在样本内部是否有基因富集。针对每个样本进行分析。
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转换结果
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对应
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得分矩阵可以用来做差异分析,看哪些通路在两个分组中存在差异,类似于基因表达差异分析。
这些分析,基本上读一下我五年前在生信技能树的表达芯片的公共数据库挖掘系列推文 就明白了;
视频观看方式
我把3年前的收费视频课程:3年前的GEO数据挖掘课程你可以听3小时或者3天甚至3个月,免费到B站:
然后马上就有了3千多学习量,而且有学员给出来了图文并茂版本万字笔记,让我非常感动!