事先聲明弧呐，這是隨筆，著重記錄思路冒签、踩坑以及解決過程，以記錄自己的學習和之后可能需要復現(xiàn)钟病。
首先推薦閱讀幾篇教程萧恕，以了解需要做什么：
1.GSVA的安裝
我使用Bioconductor安裝票唆，非常方便。
2.獲取選定通路的基因集
（這里不確定自己使用的最簡單的方法屹徘，方法來自如下：
KEGG通路數(shù)據(jù)庫下載(人種)-BeeBee生信
這個代碼可以放心跑走趋，雖然看起來很長，但是用自己的電腦也可以跑通噪伊。有一步需要的時間比較長簿煌，耐心等待即可。解析如下：
2.1 首先需要下載兩個R包鉴吹，但通常情況下可能已經(jīng)裝載在R里面了啦吧，因此直接library一手先。沒有的話去裝包拙寡，都很好裝授滓。
library(KEGGREST, quietly = TRUE)
library(tidyverse, quietly = TRUE)
2.2 之后加載推文里的函數(shù)
# keggGet(x)[[1]]$GENE 數(shù)據(jù)基因名是個向量，其中奇數(shù)位置是 entrezgene_id 偶數(shù)位置是 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 數(shù)據(jù)基因名是個向量肆糕，其中奇數(shù)位置是 entrezgene_id 偶數(shù)位置是 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函數(shù)，看看返回數(shù)據(jù)的結(jié)構诚啃。
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. 用前面定義函數(shù)，獲得每個通路的基因始赎，然后也放在表格里
pathwayFull <- hsaPathway %% dplyr::mutate(entrezgene_id=map_chr(pathway_id, geneEntrez), hgnc_symbol=map_chr(pathway_id, geneSymbol))
4.查看數(shù)據(jù)維度
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) 內(nèi)容如下
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函數(shù)所需要的列表形式
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")
結(jié)果如下：
image
到這里办斑，已經(jīng)完成了一半的工作外恕，但是這個還是不能用，因為參照正常的數(shù)據(jù)集是這樣的：
image
因此需要對我們的基因集進行改造：
一行函數(shù)輕松搞定乡翅，比較吃基本功!
genelist_2 <- lapply(genelist,function(x) unlist(strsplit(x,split = ",")))
image
4. 進行GSVA分析
4.1 首先取單細胞的表達量數(shù)據(jù)集
一定要是matrix格式的鳞疲！
gene.expr <- as.matrix(sce_cnv[["RNA"]]@data)
進行分析
gsva.result <- GSVA::gsva(gene.expr, genelist_2,kcdf='Gaussian')
首先轉(zhuǎn)換為數(shù)據(jù)框格式，之后保存結(jié)果
gsva.result_df <- as.data.frame(gsva.result)
save(gsva.result_df,file = "./gsva.result_df.RData")
4.2 可以對GSVA Score進行做熱圖等相關分析
……
