Day 11 充分理解GSVA和GSEA

單個基因集的GSEA分析

http://software.broadinstitute.org/gsea/msigdb/index.jsp

GSEA官網(wǎng)列出關(guān)于msigdb基因集的詳細(xì)信息

C1是通過染色體位置進(jìn)行富集分析

load('./CORO2A/CORO2A_up_and_down_gene.Rdata')
load(file = './Rdata/DEGs_results_CORO2A_Si2_vs_NC.Rdata')

gene_diff<- c(up_gene,Down_gene)


# DEG<- DEG[rownames(DEG)%in%up_gene,]
# DEG<- DEG[rownames(DEG)%in%Down_gene,]

geneList=DEG$log2FoldChange

names(geneList)=rownames(DEG)

geneList=sort(geneList,decreasing = T)
#選擇gmt文件（MigDB中的全部基因集）
d='D:\\GSVA_and_GSEA_Gene_Sets/MsigDB/symbols/' ### 特別注意這里是兩杠，第一根\的意思就是轉(zhuǎn)字符？
gmts <- list.files(d,pattern = 'all')
gmts
#### 第二：關(guān)于GSEA的理解
library(GSEABase) # BiocManager::install('GSEABase')
## 下面使用lapply循環(huán)讀取每個gmt文件，并且進(jìn)行GSEA分析
gmtfile=gmts[1]
geneset <- read.gmt(file.path(d,gmtfile)) 
print(paste0('Now process the ',gmtfile))
egmt <- GSEA(geneList,
             TERM2GENE=geneset, 
             pvalueCutoff = 0.9,
             verbose=FALSE)
head(egmt)
gseaplot(egmt, geneSetID = rownames(egmt[1,]))
image.png
image.png

多個基因集和GSEA和GSVA分析

第一：要下載相應(yīng)的GSEA基因集

第二：關(guān)于GSEA 理解輸入輸出

1.需要輸入一個geneList選取差異表達(dá)的基因

2.需要輸入gmts數(shù)據(jù)集（從GSEA官網(wǎng)下載）：

d='./MsigDB/symbols'
gmts <- list.files(d,pattern = 'all')
gmts

3.lapply函數(shù)是對gmts基因集的循環(huán)，如果下載全基因集或者只關(guān)注一個基因集，那就不需要循環(huán)：如注釋信息 # gmtfile=gmts[2]:相當(dāng)于只取第二個基因集

  1. geneset <- read.gmt(file.path(d,gmtfile)) ：用read.gmt 這個函數(shù)讀取gmt文檔信息
  2. egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
    ? head(egmt)：相當(dāng)于還是主函數(shù)。主要的數(shù)據(jù)分析
  3. gseaplot(egmt, geneSetID = rownames(egmt[1,])) ：選取了一條基因集畫圖：選取的基因集的名字是rownames(egmt[1,])

7： gsea_results_list<- lapply(gsea_results, function(x){
? cat(paste(dim(x@result)),'\n')
? x@result
}) ##########這個的意思應(yīng)該是吧gsea_results進(jìn)行了操作列出 gsea_results_list

8：gsea_results_df <- do.call(rbind, gsea_results_list) ###########把gsea_results_list做轉(zhuǎn)化成data.frame。do.call神奇操作

9：gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE') 。選擇有差異的基因集進(jìn)行畫圖。此處要注意gsea_results[[2]]其中的數(shù)字2 要和KEGG_CELL_CYCLE是對應(yīng)的。所以畫圖前一定要用notepad++查看下載下來的基因集和數(shù)字的對應(yīng)關(guān)系。

load(file = 'anno_DEG.Rdata')
  geneList=DEG$logFC
  names(geneList)=DEG$symbol
  geneList=sort(geneList,decreasing = T)
  #選擇gmt文件（MigDB中的全部基因集）
  d='./MsigDB/symbols'
  gmts <- list.files(d,pattern = 'all')
  gmts
 #### 第二：關(guān)于GSEA的理解
  library(GSEABase) # BiocManager::install('GSEABase')
  ## 下面使用lapply循環(huán)讀取每個gmt文件，并且進(jìn)行GSEA分析
  gsea_results <- lapply(gmts, function(gmtfile){
    # gmtfile=gmts[2]
    geneset <- read.gmt(file.path(d,gmtfile)) 
    print(paste0('Now process the ',gmtfile))
    egmt <- GSEA(geneList, TERM2GENE=geneset, verbose=FALSE)
    head(egmt)
    # gseaplot(egmt, geneSetID = rownames(egmt[1,]))
    
    return(egmt)
  })
  # 上面的代碼耗時(shí)，所以保存結(jié)果到本地文件
  save(gsea_results,file = 'gsea_results.Rdata')
  #提取gsea結(jié)果，熟悉這個對象
  gsea_results_list<- lapply(gsea_results, function(x){
    cat(paste(dim(x@result)),'\n')
    x@result
  })
  gsea_results_df <- do.call(rbind, gsea_results_list)
  gseaplot(gsea_results[[2]],'KEGG_CELL_CYCLE') 
  gseaplot(gsea_results[[2]],'FARMER_BREAST_CANCER_CLUSTER_6') 
  gseaplot(gsea_results[[5]],'GO_CONDENSED_CHROMOSOME_OUTER_KINETOCHORE') 
  
}

第三：GSVA的輸入輸出

計(jì)算每個樣本在基因集中的富集程度，然后計(jì)算Fold change 和P value。做類似于熱圖的圖。所以后面的很多算法和制作熱圖一致。

  1. 關(guān)注input：函數(shù)需要輸入X=dat這個expression data，row為基因名，col為sample 。另外還需要gmt這個數(shù)據(jù)集。

  2. d='./MSigDB/symbols/'
    gmts=list.files(d,pattern = 'all') 通過list.files這個可以列出當(dāng)前路徑下的所有g(shù)mt文件

  3. es_max <- lapply(gmts, function(gmtfile) ；使用lapply函數(shù)進(jìn)行循環(huán)，對內(nèi)一個下載的基因集進(jìn)行操作；如果我們只關(guān)心一個基因集，就不需要這個循環(huán)了比如 #gmtfile=gmts[8]

  4. geneset <- getGmt(file.path(d,gmtfile))
    ? es.max <- gsva(X, geneset,
    ? mx.diff=FALSE, verbose=FALSE,
    ? parallel.sz=1)

#######這一步很重要，是函數(shù)的最重要步驟。geneset <- getGmt(file.path(d,gmtfile)) 通過這個getGmt這個函數(shù)可以從gmt 文件中獲得geneset；gsva這個主函數(shù)通過輸入X這個表達(dá)矩陣和geneset數(shù)據(jù)集進(jìn)行計(jì)算Fold change 和P value。

5adjPvalueCutoff <- 0.001
? logFCcutoff <- log2(2)####進(jìn)行閾值調(diào)整

6.帥選差異表達(dá)的基因集##類似做熱圖和火山圖（因?yàn)閷?shí)驗(yàn)設(shè)計(jì)是本身存在不同的分組信息）；

6.1

es_deg <- lapply(es_max, function(es.max)這個函數(shù)對所有的基因集得到的es.max進(jìn)行循環(huán)操作

6.2

分組內(nèi)容如下：design <- model.matrix(~0+factor(group_list))
? colnames(design)=levels(factor(group_list))
? rownames(design)=colnames(es.max)

6.3 deg = function(es.max,design,contrast.matrix) 這個函數(shù)對es.max按照分組信息進(jìn)行基因集富集情況進(jìn)行差異分析（類似于尋找差異表達(dá)基因）

7.最終會輸出ex.max 和ex.deg這兩個數(shù)據(jù)

8.最后作圖，做差異富集的圖，（類似差異表達(dá)的基因）作圖的時(shí)候回報(bào)錯，我標(biāo)注在代碼上：詳見代碼注釋

### 對 MigDB中的全部基因集 做GSVA分析。
## 還有ssGSEA, PGSEA
{
  load(file = 'step1-output.Rdata')
  # 每次都要檢測數(shù)據(jù)
  dat[1:4,1:4]  
  library(hgu133plus2.db)
  ids=toTable(hgu133plus2SYMBOL)
  head(ids)
  dat=dat[ids$probe_id,]
  dat[1:4,1:4] 
  ids$median=apply(dat,1,median)
  ids=ids[order(ids$symbol,ids$median,decreasing = T),]
  ids=ids[!duplicated(ids$symbol),]
  dat=dat[ids$probe_id,]
  rownames(dat)=ids$symbol
  dat[1:4,1:4]  
  
  X=dat
  table(group_list)
  ## Molecular Signatures Database (MSigDb) 
  d='./MSigDB/symbols/'
  gmts=list.files(d,pattern = 'all')
  gmts
  library(ggplot2)
  library(clusterProfiler)
  library(org.Hs.eg.db)
  library(GSVA) # BiocManager::install('GSVA')
  
  if(T){
    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)
    })
    adjPvalueCutoff <- 0.001
    logFCcutoff <- log2(2)
    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("TNBC-noTNBC",
                                     levels = design)
      
      contrast.matrix ##這個矩陣聲明，我們要把progres.組跟stable進(jìn)行差異分析比較
      
      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='gsva_msigdb.Rdata')
  
  load(file='gsva_msigdb.Rdata')
  
  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,]   #######在這里調(diào)閾值，保證輸出
    print(dim(df))
    if(nrow(df)>5){    #######nrow(df)>5:只有差異富集的基因集大于五個以上才畫圖
      n=rownames(df)
      dat=dat[match(n,rownames(dat)),]  ######n=rownames(df),篩選表達(dá)矩陣數(shù)據(jù)
      ac=data.frame(g=group_list)   ########對應(yīng)后面的圖中的annotation_col
      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'))
      
    }
})
  
  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')
}
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容