事先声明，这是随笔，着重记录思路、踩坑以及解决过程，以记录自己的学习和之后可能需要复现。
首先推荐阅读几篇教程，以了解需要做什么：
1.GSVA的安装
我使用Bioconductor安装，非常方便。
2.获取选定通路的基因集
（这里不确定自己使用的最简单的方法，方法来自如下：
这个代码可以放心跑，虽然看起来很长，但是用自己的电脑也可以跑通。有一步需要的时间比较长，耐心等待即可。解析如下：
2.1 首先需要下载两个R包，但通常情况下可能已经装载在R里面了，因此直接library一手先。没有的话去装包，都很好装。
library(KEGGREST, quietly = TRUE)
library(tidyverse, quietly = TRUE)
2.2 之后加载推文里的函数
# 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])
}
2.3 正式获取通路部分
建议从这里开始读脚本。建议自己在交互模式下试一下用到的KEGGREST函数，看看返回数据的结构。
1.这是第一步，取得所有的KEGG通路列表
hsaList <- keggList("pathway", "hsa")
IDList <- names(hsaList) %% map_chr(pathwayID)
2. 将通路ID和通路名放在一个表格(tibble)里
hsaPathway <- tibble::tibble(pathway_id=IDList, pathway_name=hsaList)
head(hsaPathway, n=3) %% print()
3. 用前面定义函数，获得每个通路的基因，然后也放在表格里
pathwayFull <- hsaPathway %% dplyr::mutate(entrezgene_id=map_chr(pathway_id, geneEntrez), hgnc_symbol=map_chr(pathway_id, geneSymbol))
4.查看数据维度
dim(pathwayFull) %% print()
5.会有通路没有基因，我的话只需要有基因的，所以把无基因的移除
pathwayWithGene <- dplyr::filter(pathwayFull, !is.na(entrezgene_id) & hgnc_symbol != "")
write_tsv(pathwayWithGene, path="KEGGREST_WithGene.tsv")
dim(pathwayWithGene) %% print()
所以最后得到的应该是pathwayWithGene这个变量！
3. 选定需要的通路，输入关键字，记录下pathway_id
3.1 View(pathwayWithGene) 内容如下
image
3.2 选取自己需要的通路成为新的变量
例如：c("hsa04350", "hsa04310", "hsa04014", "hsa04390", "hsa04330","hsa04110","hsa04151","hsa04115")，输入代码：
pathwayWithGene_select <- pathwayWithGene[c("hsa04350", "hsa04310",
"hsa04014", "hsa04390", "hsa04330","hsa04110","hsa04151","hsa04115"),]
image
3.3 将通路写成gsva函数所需要的列表形式
genelist <- list()
for (i in seq(length(unique(pathwayWithGene_select$pathway_id)))){
genelist[i] <- pathwayWithGene_select$hgnc_symbol[i]
}
names(genelist) <- paste0("gs",c(1:length(unique(pathwayWithGene_select$pathway_id))))
save(genelist, file = "./genelist.RData")
结果如下：
image
到这里，已经完成了一半的工作，但是这个还是不能用，因为参照正常的数据集是这样的：
image
因此需要对我们的基因集进行改造：
一行函数轻松搞定，比较吃基本功!
genelist_2 <- lapply(genelist,function(x) unlist(strsplit(x,split = ",")))
image
4. 进行GSVA分析
4.1 首先取单细胞的表达量数据集
一定要是matrix格式的！
gene.expr <- as.matrix(sce_cnv[["RNA"]]@data)
进行分析
gsva.result <- GSVA::gsva(gene.expr, genelist_2,kcdf='Gaussian')
首先转换为数据框格式，之后保存结果
gsva.result_df <- as.data.frame(gsva.result)
save(gsva.result_df,file = "./gsva.result_df.RData")
4.2 可以对GSVA Score进行做热图等相关分析
……
