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

生信技能树

发布于 2020-10-26 10:59:36

4.1K0

发布于 2020-10-26 10:59:36

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

下面是学徒写的《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,便于检查