Down stream of RNA-seq 03 (GSVA)

最新推荐文章于 2024-07-17 20:06:24 发布

壹壹玖玖陆拾柒

最新推荐文章于 2024-07-17 20:06:24 发布

阅读量354

点赞数 9

文章标签：生物学 r语言

本文链接：https://blog.csdn.net/ymy2019/article/details/135383673

版权

GSVA 基因集变异分析

Code & Figures：

### 读取特定表格的表达矩阵（这里选用31个样本）

dim(exprSet_focus)
table(group_list)

library(limma)
library(edgeR)

dge <- DGEList(counts=exprSet_focus)
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
logCPM[1:4,1:4]

X <- as.data.frame(logCPM)
d <- 'msigdb_v2023.2.Mm_GMTs/'

gmts <- list.files(d,pattern = 'all')
gmts
library(ggplot2)
library(clusterProfiler)
library(org.Mm.eg.db)
library(GSVA) 
library(GSEABase)

keytypes(org.Mm.eg.db)
# 如果ENSG的ID具有小数点，用下面两句代码舍掉小数点
# library(stringr)
# rownames(DEG) <- str_sub(rownames(DEG), start = 1, end = 15)
X$ENSEMBL <- rownames(X)
# bitr in clusterProfiler
df <- bitr( rownames(X), fromType = "ENSEMBL", toType = c( "ENTREZID" ), OrgDb = org.Hs.eg.db )
head(df)
df <- unique(df, by = "ENTREZID")
X <- merge(X, df, by='ENSEMBL')
head(X)
dim(X)
X <- cbind(X[ncol(X)], X[-c(1,ncol(X))])
X <- avereps(X, ID = X$ENTREZID)

rownames(X) <- X[,1]
exp <- X[,2:ncol(X)]

f  <- './Step07-es_max.Rda'
if(!file.exists(f)){
  es_max <- lapply(gmts, function(gmtfile){ 
        #gmtfile=gmts[8];gmtfile
        geneset <- getGmt(file.path(d,gmtfile))  
        es.max <- gsva(exp, geneset, 
                       mx.diff=FALSE, verbose=FALSE, 
                       parallel.sz=1)
        return(es.max)
      })
  save(es_max, file = f)
}
adjPvalueCutoff <- 0.001
logFCcutoff <- log2(2)
f  <- './Step07-gsva_msigdb.Rda'
if(!file.exists(f)){
  es_deg <- lapply(es_max, function(es.max){
    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("Patch_Infarction-MI_Infarction",
                                   levels = design)
    
    contrast.matrix ##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较
    
    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)
  })
  
  
  gmts
  
  save(es_max,es_deg,file='./Step07-gsva_msigdb.Rda')
}


load(file = f)
  
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.3,]
    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,width=20,
                         annotation_col = ac,show_colnames = F,
                         filename = paste0('./Step07-gsva_',strsplit(gmts[i],'[.]')[[1]][1],'.png'))
      
  }
})
  
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 = './Step07-GSVA_DEG.csv')

Bugs：getGmt 调用的matrixStats

如果按照这个安装，需要在mac上下载Xcode,但是Xcode在软件商店安装会出现说电脑容量不够，看网上的经验说解压缩后是35g左右，Xcode15，所以这一步，getGmt做通路富集我放在学院集群跑去了，得到的Rda再下载到本地做后续的差异分析和可视化。
其中当不晓得getGmt为什么报错时，曾经以为是我读入的gen_list有问题，于是找了msigdbr的包，而不是从GSEA下载的数据库，教程不错，是生信技能树的头牌讲师小洁老师的推文【MsigDB的那些基因集合都R语言化，想做GSEA和GSVA，就可以不用去下载gmt文件啦！】https://www.jianshu.com/p/0098baf2df46
https://zhuanlan.zhihu.com/p/518145829?utm_id=0（下图的来源教程）

msigdbr_species()
这里详细写有如何从msigdbr包找特定的KEGG和GO对应的库

但是感觉还是上述Codes &Figures中的从本地读取gmts更有效
循环同时把每个database都富集和差异分析了。

但小洁老师这个教程是很好的初探索msigdbr包使用的比较全面的教程！

还有这个教程也还不错https://blog.csdn.net/qq_42090739/article/details/121788483
https://blog.csdn.net/weixin_49320263/article/details/131530673


注意在gsva函数中，cdf="Gaussian" ,#"Gaussian" for logCPM,logRPKM,logTPM,
"Poisson" for counts
http://www.manongjc.com/detail/63-dwgdiggecdgqxfv.html

这儿还有个画弦图的教程
https://blog.csdn.net/weixin_49320263/article/details/131530673