参考小丫画图：61
BiocManager::install("GSVA")
library(clusterProfiler)
library(limma)
library(ggplot2)
library(data.table)
library(ggpubr)
library(GSVA)
Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
setwd("/data/ff/ref/GSVA")
#查看msigdbr包里自带的物种
msigdbr_show_species()
#H: hallmark gene sets;C1: positional gene sets; C2: curated gene sets.KEGG ;C3: motif gene sets;
#C4: computational gene sets
#C5: GO gene sets;C6: oncogenic signatures;C7: immunologic signatures
#"Homo sapiens" "Mus musculus"
###############制作参考基因集----
#########################H
h <- msigdbr(species = "Homo sapiens", # 物种拉丁名
category = "H") #此处以hallmark为例，你也可以选择MSigDB的其他注释
# 示例数据表达矩阵的基因名是gene symbol，这里就选gene_symbol。
# 如果你的表达矩阵以ENTREZ ID作为基因名，就把下面这段的gene_symbol换成entrez_gene
h <- dplyr::select(h, gs_name, gene_symbol) %>% #或entrez_gene
as.data.frame %>%
split(., .$gs_name) %>%
lapply(., function(x)(x$gene_symbol)) #或entrez_gene
# 在每个geneset里面去掉重复的基因
gs <- lapply(h, unique)
# 预览过滤后的结果
head(gs)
# 保存到文件，方便以后重复使用
save(gs, file = "hallmark.gs.RData")
#########################C1
C1 <- msigdbr(species = "Homo sapiens", # 物种拉丁名
category = "C1") #此处以C1allmark为例，你也可以选择MSigDB的其他注释
# 示例数据表达矩阵的基因名是gene symbol，这里就选gene_symbol。
# 如果你的表达矩阵以ENTREZ ID作为基因名，就把下面这段的gene_symbol换成entrez_gene
C1 <- dplyr::select(C1, gs_name, gene_symbol) %>% #或entrez_gene
as.data.frame %>%
split(., .$gs_name) %>%
lapply(., function(x)(x$gene_symbol)) #或entrez_gene
# 在每个geneset里面去掉重复的基因
gs <- lapply(C1, unique)
# 预览过滤后的结果
head(gs)
# 保存到文件，方便以后重复使用
save(gs, file = "C1.positional.gene.sets.gs.RData")
#########################C2
C2 <- msigdbr(species = "Homo sapiens", # 物种拉丁名
category = "C2") #此处以C2allmark为例，你也可以选择MSigDB的其他注释
# 示例数据表达矩阵的基因名是gene symbol，这里就选gene_symbol。
# 如果你的表达矩阵以ENTREZ ID作为基因名，就把下面这段的gene_symbol换成entrez_gene
C2 <- dplyr::select(C2, gs_name, gene_symbol) %>% #或entrez_gene
as.data.frame %>%
split(., .$gs_name) %>%
lapply(., function(x)(x$gene_symbol)) #或entrez_gene
# 在每个geneset里面去掉重复的基因
gs <- lapply(C2, unique)
# 预览过滤后的结果
head(gs)
# 保存到文件，方便以后重复使用
save(gs, file = "C2.KEGG.gs.RData")
#########################C3
C3 <- msigdbr(species = "Homo sapiens", # 物种拉丁名
category = "C3") #此处以C3allmark为例，你也可以选择MSigDB的其他注释
# 示例数据表达矩阵的基因名是gene symbol，这里就选gene_symbol。
# 如果你的表达矩阵以ENTREZ ID作为基因名，就把下面这段的gene_symbol换成entrez_gene
C3 <- dplyr::select(C3, gs_name, gene_symbol) %>% #或entrez_gene
as.data.frame %>%
split(., .$gs_name) %>%
lapply(., function(x)(x$gene_symbol)) #或entrez_gene
# 在每个geneset里面去掉重复的基因
gs <- lapply(C3, unique)
# 预览过滤后的结果
head(gs)
# 保存到文件，方便以后重复使用
save(gs, file = "C3.motif.gene.sets.gs.RData")
#########################C4
C4 <- msigdbr(species = "Homo sapiens", # 物种拉丁名
category = "C4") #此处以C4allmark为例，你也可以选择MSigDB的其他注释
# 示例数据表达矩阵的基因名是gene symbol，这里就选gene_symbol。
# 如果你的表达矩阵以ENTREZ ID作为基因名，就把下面这段的gene_symbol换成entrez_gene
C4 <- dplyr::select(C4, gs_name, gene_symbol) %>% #或entrez_gene
as.data.frame %>%
split(., .$gs_name) %>%
lapply(., function(x)(x$gene_symbol)) #或entrez_gene
# 在每个geneset里面去掉重复的基因
gs <- lapply(C4, unique)
# 预览过滤后的结果
head(gs)
# 保存到文件，方便以后重复使用
save(gs, file = "C4.computational gene sets.gs.RData")
#########################C5
C5 <- msigdbr(species = "Homo sapiens", # 物种拉丁名
category = "C5") #此处以C5allmark为例，你也可以选择MSigDB的其他注释
# 示例数据表达矩阵的基因名是gene symbol，这里就选gene_symbol。
# 如果你的表达矩阵以ENTREZ ID作为基因名，就把下面这段的gene_symbol换成entrez_gene
C5 <- dplyr::select(C5, gs_name, gene_symbol) %>% #或entrez_gene
as.data.frame %>%
split(., .$gs_name) %>%
lapply(., function(x)(x$gene_symbol)) #或entrez_gene
# 在每个geneset里面去掉重复的基因
gs <- lapply(C5, unique)
# 预览过滤后的结果
head(gs)
# 保存到文件，方便以后重复使用
save(gs, file = "C5.GO.gs.RData")
#########################C6
C6 <- msigdbr(species = "Homo sapiens", # 物种拉丁名
category = "C6") #此处以C6allmark为例，你也可以选择MSigDB的其他注释
# 示例数据表达矩阵的基因名是gene symbol，这里就选gene_symbol。
# 如果你的表达矩阵以ENTREZ ID作为基因名，就把下面这段的gene_symbol换成entrez_gene
C6 <- dplyr::select(C6, gs_name, gene_symbol) %>% #或entrez_gene
as.data.frame %>%
split(., .$gs_name) %>%
lapply(., function(x)(x$gene_symbol)) #或entrez_gene
# 在每个geneset里面去掉重复的基因
gs <- lapply(C6, unique)
# 预览过滤后的结果
head(gs)
# 保存到文件，方便以后重复使用
save(gs, file = "C6.oncogenic.gs.RData")
#########################C7
C7 <- msigdbr(species = "Homo sapiens", # 物种拉丁名
category = "C7") #此处以C7allmark为例，你也可以选择MSigDB的其他注释
# 示例数据表达矩阵的基因名是gene symbol，这里就选gene_symbol。
# 如果你的表达矩阵以ENTREZ ID作为基因名，就把下面这段的gene_symbol换成entrez_gene
C7 <- dplyr::select(C7, gs_name, gene_symbol) %>% #或entrez_gene
as.data.frame %>%
split(., .$gs_name) %>%
lapply(., function(x)(x$gene_symbol)) #或entrez_gene
# 在每个geneset里面去掉重复的基因
gs <- lapply(C7, unique)
# 预览过滤后的结果
head(gs)
# 保存到文件，方便以后重复使用
save(gs, file = "C7.immunologic.gs.RData")
###################酌情考虑
为了减少冗余信息的干扰，
在每个geneset中，不应该出现重复的基因；↑
在两个或更多个pathway中出现的基因应该被彻底剔除。↓
# 在每个geneset里面去掉重复的基因
gs <- lapply(h, unique)
# 接下来去掉那些在两个或更多个pathways里出现过的genes
count <- table(unlist(gs))
keep <- names(which(table(unlist(gs)) < 2))
gs <- lapply(gs, function(x) intersect(keep, x))
# 过滤之后，很多pathway一个gene都不剩了，去掉这些
gs <- gs[lapply(gs, length) > 0]
# 预览过滤后的结果
head(gs)
#######################开始GSVA
################C2.KEGG----
(load("/data/ff/ref/GSVA/C2.KEGG.gs.RData")) #保存在当前文件夹
# 这一句就完成了GSVA分析
gsva_KEGG <- gsva(as.matrix(gsym.expr), gs)
head(gsva_KEGG)
# 把通路的表达量保存到文件
write.csv(gsva_KEGG, "gsva.C2.KEGG.csv", quote = F)
################C5.GO----
(load("/data/ff/ref/GSVA/C5.GO.gs.RData")) #保存在当前文件夹
# 这一句就完成了GSVA分析
gsva_GO <- gsva(as.matrix(gsym.expr), gs)
head(gsva_GO)
# 把通路的表达量保存到文件
write.csv(gsva_GO, "gsva.C5.GO.csv", quote = F)
################C7.IMMUNOLOGIC----
(load("/data/ff/ref/GSVA/C7.immunologic.gs.RData")) #保存在当前文件夹
# 这一句就完成了GSVA分析
gsva_C7 <- gsva(as.matrix(gsym.expr), gs)
head(gsva_C7)
# 把通路的表达量保存到文件
write.csv(gsva_C7, "gsva.C7.immunologic.csv", quote = F)
################HALLMARK----
(load("/data/ff/ref/GSVA/hallmark.gs.RData")) #保存在当前文件夹
# 这一句就完成了GSVA分析
gsva_hallmark <- gsva(as.matrix(gsym.expr), gs)
head(gsva_hallmark)
# 把通路的表达量保存到文件
write.csv(gsva_hallmark, "gsva.HALLMARK.csv", quote = F)
###########################合并总文件
gsva.all <- rbind(gsva_KEGG,gsva_GO,gsva_C7,gsva_hallmark)
# 把通路的表达量保存到文件
write.csv(gsva.all, "gsva.all.csv", quote = F)
网友评论