2021-10-17-GSVA分析

准备gmt文件

1、在GSEA网站上下载相应的通路的gmt文件，[GSEA官网]（http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp#C2），在上面可以下载MSigDB所收录的基因集。下下来以后跑代码分析
比如我们先下载一个上面收录的kegg的通路集合,SYMBOL和NCBI ID都可以
2、自己构建数据集
下载了KEGG官网的数据集，两列数据，一列通路，一列基因集

library(KEGGREST, quietly = TRUE)
library(tidyverse, quietly = TRUE)
setwd("D:/SXFX/SXDT1/GSEA/")
# 返回信息很长，只取基因symbol.根据自己需要调整
symbolOnly <- function(x){
  items <- strsplit(x, ";", fixed = TRUE) %>% unlist()
  return(items[1])
}

# keggGet(x)[[1]]$GENE 数据基因名是个向量，其中奇数位置是 entrezgene_id 偶数位置是 symbol 
geneEntrez <- function(x){
  geneList <- keggGet(x)[[1]]$GENE
  if(!is.null(geneList)){
    listLength <- length(geneList)
    entrezList <- geneList[seq.int(from = 1, by = 2, length.out = listLength/2)]
    entrez <- stringr::str_c(entrezList, collapse = ",")
    return(entrez)
  }else{
    return(NA)
  }
  
}

# keggGet(x)[[1]]$GENE 数据基因名是个向量，其中奇数位置是 entrezgene_id 偶数位置是 symbol 
geneSymbol <- function(x){
  geneList <- keggGet(x)[[1]]$GENE
  if(!is.null(geneList)){
    listLength <- length(geneList)
    symbolList <- geneList[seq.int(from = 2, by = 2, length.out = listLength/2)] %>% map_chr(symbolOnly)
    symbol <- stringr::str_c(symbolList, collapse = ",")
    return(symbol)
  }else{
    return("")
  }
  
}

# 取得 hsaxxxxx 这种通路ID
pathwayID <- function(x){
  items <- strsplit(x, ":", fixed = TRUE) %>% unlist()
  return(items[2])
}

# 建议从这里开始读脚本。建议自己在交互模式下试一下用到的KEGGREST函数，看看返回数据的结构。
# 这是第一步，取得所有的KEGG通路列表
hsaList <- keggList("pathway", "hsa")
IDList <- names(hsaList) %>% map_chr(pathwayID)

# 将通路ID和通路名放在一个表格(tibble)里
hsaPathway <- tibble::tibble(pathway_id=IDList, pathway_name=hsaList)
head(hsaPathway, n=3) %>% print()

# 用前面定义函数，获得每个通路的基因，然后也放在表格里
pathwayFull <- hsaPathway %>% dplyr::mutate(entrezgene_id=map_chr(pathway_id, geneEntrez), hgnc_symbol=map_chr(pathway_id, geneSymbol))

# 保存数据
write_tsv(pathwayFull, path="KEGGREST.tsv")
dim(pathwayFull) %>% print()

# 会有通路没有基因，我的话只需要有基因的，所以把无基因的移除
pathwayWithGene <- dplyr::filter(pathwayFull, !is.na(entrezgene_id) & hgnc_symbol != "")
write_tsv(pathwayWithGene, path="KEGGREST_WithGene.tsv")
dim(pathwayWithGene) %>% print()

构建gmt格式的文件

gmt<-data.frame()
for (i in 1:345){
  dt<-data.frame(unlist(str_split(pathwayWithGene$hgnc_symbol[i],',')))
  colnames(dt)[1]='gene'
  dt$term=unlist(str_split(pathwayWithGene$pathway_name[i],'-'))[1]
  gmt<-rbind(gmt,dt)
}
gmt<-gmt[,c('term','gene')]

GeneID2kegg_list<<- tapply(gmt[,1],as.factor(gmt[,2]),function(x) x)
kegg2GeneID_list<<- tapply(gmt[,2],as.factor(gmt[,1]),function(x) x)

write.gmt <- function(geneSet=kegg2symbol_list,gmt_file='kegg2symbol.gmt'){
  sink( gmt_file )
  for (i in 1:length(geneSet)){
    cat(names(geneSet)[i])
    cat('\tNA\t')
    cat(paste(geneSet[[i]],collapse = '\t'))
    cat('\n')
  }
  sink()
}
write.gmt(kegg2GeneID_list,'kegg2symbol.gmt')

GSVA 分析

rm(list = ls())
options(stringsAsFactors = F)
setwd("D:/SXFX/Analysis_data/WH/21-08-07/")
exp<-read.csv('expression_matrix.2.csv',row.names = 1)
suppressMessages(library(GSVA))
suppressMessages(library(GSVAdata))
suppressMessages(library(GSEABase))
suppressMessages(library(limma))
setwd("D:/SXFX/SXDT1/GSEA/")
geneset <- getGmt('kegg2symbol.gmt')  
es <- gsva(as.matrix(exp), geneset,
           min.sz=10, max.sz=500, verbose=TRUE)