GSVA基因集打分/多通路limma差异分析流程及可视化学习和整理

最新推荐文章于 2025-02-27 12:24:38 发布

凑齐六个字吧

最新推荐文章于 2025-02-27 12:24:38 发布

阅读量9.8k

点赞数 28

分类专栏：分析工具文章标签：人工智能数据分析数据库

本文链接：https://blog.csdn.net/zfyyzhys/article/details/141576596

版权

分析工具专栏收录该内容

43 篇文章

订阅专栏

基因集变异分析

GSVA（Gene Set Variation Analysis，基因集变异分析） 是一种无监督的基因集富集分析方法，用于评估样本层面上基因集（如 KEGG 路径、GO 术语）的表现。与传统的基因集富集分析不同，GSVA 不是直接检测不同基因集在不同实验条件下的富集程度，而是计算每个样本中预定义的基因集的富集得分，从而评估基因集在不同样本或条件下的变异性。

应用场景：

评估基因集的表现：GSVA 通过计算样本中基因集的富集得分，评估基因集在不同样本中的表现。这使得研究者可以在无监督的背景下评估不同样本中的基因集活动水平，而无需先定义实验组和对照组。

样本分类：GSVA 可以用于对样本进行分类或聚类分析，帮助识别样本间基因集的差异，从而揭示潜在的生物学机制。

医学科研：GSVA 在医学科研中具有重要应用，可以根据患者的基因表达谱评估关键基因集的活性水平，进而指导个性化治疗方案的选择。

适用范围/条件：

单细胞/bulk RNA数据：GSVA 适用于处理来自单细胞或bulkRNA 测序的基因表达数据。

基因集分析：适用于任何预定义的基因集，如 KEGG 路径、GO 术语、TF 调控网络等。

无监督分析：GSVA 主要用于无监督分析背景下的基因集富集分析，不需要预先定义样本的分组信息。

GSVA 与其他算法的差别

与 GSEA（Gene Set Enrichment Analysis）的差别：

分析方式：GSEA 主要用于检测不同实验组（如实验组和对照组）之间基因集富集差异，而 GSVA 是在每个样本层面计算基因集的富集得分。

输入数据：GSEA 需要预先定义样本分组，而 GSVA 则不需要分组信息，可以直接在无监督的情况下计算基因集在每个样本中的表现。

分析结果：GSEA 的结果是基因集在不同条件下的富集显著性（如 p 值），而 GSVA 的结果是每个样本中基因集的富集得分（一个连续值）。

与 ssGSEA（Single-sample Gene Set Enrichment Analysis） 的差别：

分析思路：ssGSEA 是 GSEA 的扩展，用于计算每个样本中基因集的富集得分。它与 GSVA 相似，都是针对单个样本进行的分析。

算法实现：GSVA 和 ssGSEA 在计算基因集富集得分时采用不同的算法。GSVA 使用的是非参数核密度估计方法，而 ssGSEA 是基于样本排序的积分方法。

应用场景：虽然二者应用场景类似，但 GSVA 通常被认为在处理批量样本时更为稳定，而 ssGSEA 更直接用于单样本分析。

与其他基因集富集算法（如 CAMERA, ROAST）的差别：

统计假设：CAMERA 和 ROAST 等方法基于线性模型和统计假设，用于评估预定义基因集在不同实验条件下的富集显著性。这些方法通常要求有明确的样本分组信息。

适用范围：GSVA 不依赖于分组信息，因此在无监督分析或复杂实验设计中更为灵活，而 CAMERA 和 ROAST 则更适合在明确分组的实验设计中使用。

GSVA基因集打分分析步骤及热图

1、加载R包、导入数据

准备了乳酸化和缺氧两个基因集

rm(list = ls())  
setwd("./2-gsva/")
library(readxl)
library(GSVA)

load(file = "combat.Rdata")
exp[1:4,1:4]
#           GSM3106326 GSM3106327 GSM3106328 GSM3106329
# LINC01128   8.000049   8.044479   7.311868   7.610713
# SAMD11      7.499611   7.562370   7.412649   7.239141
# KLHL17      7.662076   7.672331   7.627330   7.564212
# PLEKHN1     7.666884   7.889537   7.640606   7.746440
head(pd)[1:4,1:4]
#                                                title                  samples   GSE_num group
# GSM3106326 Pulmonary arterial hypertension patient 1   idiopathic PAH patient GSE113439   PAH
# GSM3106327 Pulmonary arterial hypertension patient 2 patient with PAH and CHD GSE113439   PAH
# GSM3106328 Pulmonary arterial hypertension patient 3 patient with PAH and CTD GSE113439   PAH
# GSM3106329 Pulmonary arterial hypertension patient 4 patient with PAH and CHD GSE113439   PAH
lactate <- read_excel("Hypoxia Related Genes.xls", sheet = 1)
names(lactate) <- "lactate"
hypoxia <- read_excel("Hypoxia Related Genes.xls", sheet = 2)
names(hypoxia) <- "hypoxia"

2、GSVA分析

# 制作基因集list
# 其中要把基因集变成matrix哦~ 
dataset <- list("hypoxia"= as.matrix(hypoxia),"lactate"= as.matrix(lactate))

# GSVA
gsvaPar <- gsvaParam(exp, dataset,maxDiff = TRUE)
gsvaPar 
# A GSVA::gsvaParam object
# expression data:
#   matrix [18837, 132]
#     rows: LINC01128, SAMD11, ..., SPRR2G (18837 total)
#     cols: GSM3106326, GSM3106327, ..., GSM1291010 (132 total)
# using assay: none
# gene sets:
#   list
#     names: 493 hypoxia-related genes (1 total)
#     unique identifiers: ADM, ABCB1, ..., VPS11 (493 total)
# gene set size: [1, Inf]
# kcdf: Gaussian
# tau: 1
# maxDiff: TRUE
# absRanking: FALSE
expr_geneset <- gsva(gsvaPar)
dim(expr_geneset)
head(expr_geneset)[1:2,1:4]
#          GSM3106326 GSM3106327  GSM3106328 GSM3106329
# hypoxia  0.08122476  0.1611271 -0.04422083 -0.1121407
# lactate -0.17011732  0.1384359  0.33231938  0.0823520

3、可视化-热图

library(ComplexHeatmap)
data_total <- as.data.frame(t(expr_geneset))
identical(rownames(data_total),rownames(pd))

# 检查数据类型是否正确转换为数值型
library(dplyr)
data_total <- data_total %>% mutate_all(as.numeric)
# 标准化-方法3
# 范围限定在-2至2
data_total <- scale(data_total)
range(data_total)
data_total[data_total > 2] <- 2
data_total[data_total + 2 < 0] <- -2

# 也可以用circlize创建颜色映射
# 定义颜色渐变
library(circlize)
color_fun <- colorRamp2(c(-2, 0, 2), c("#336699", "white", "tomato"))

# pd$type 和 pd$samples 是注释数据怎么上色
type_colors <- c("control" = "green", "PAH" = "red")
samples_colors <- c("all PAH" = "blue", "CTEPH patient" = "orange", "FD" = "purple",
                    "idiopathic PAH patient" = "brown", "normal control" = "yellow",
                    "patient with PAH and CHD" = "pink", "patient with PAH and CTD" = "magenta",
                    "pulmonary arterial hypertension (PAH) patient" = "cyan")
GSE_colors <- c("GSE113439" = "navy","GSE53408"="darkolivegreen","GSE117261"="darkgoldenrod")

# 注意这里是行注释！如果是顶部注释需要用HeatmapAnnotation
rowAnno <- rowAnnotation(type = pd$group,samples = pd$samples,GSE_num =pd$GSE_num,
                                col = list(type = type_colors, 
                                           samples = samples_colors,
                                           GSE_num = GSE_colors))

ComplexHeatmap::Heatmap(data_total, 
                        na_col = "white",
                        col = color_fun,  # 添加颜色映射函数
                        show_column_names = T, #是否显示列名
                        row_names_side = "left",
                        name = "fraction",
                        row_order = c(rownames(data_total)[c(grep("control",pd$group),
                                                                grep("PAH",pd$group))]),
                        row_split = pd$group, 
                        row_title = NULL,
                        cluster_columns = F,
                        left_annotation = rowAnno,
                        #top_annotation = columnAnno,
                        #heatmap_width = unit(20, "cm"),  # 调整热图宽度
                        #row_dend_width = unit(1, "cm"),  # 调整聚类树宽度
                        cluster_rows = F, # 行聚类
                        row_names_gp = gpar(fontsize = 0),  # 调整行名字体大小
                        column_names_gp = gpar(fontsize = 12)
                        #cluster_columns = FALSE # 列聚类。                     
                        )

多条通路GSVA分析及limma差异分析步骤及热图

1、GSVA分析

数据跟上面的一样

###GSVA分析
library(msigdbr)
library(GSVA)
library(tidyverse)
library(clusterProfiler)
library(patchwork)
library(limma)
library(data.table)

#genesets <- msigdbr(species = "Homo sapiens") 
genesets <- msigdbr(species = "Homo sapiens", category = "C2", subcategory = "CP:KEGG")
genesets <- subset(genesets, select = c("gs_name","gene_symbol")) %>% as.data.frame()
genesets <- split(genesets$gene_symbol, genesets$gs_name)

# gsva默认开启全部线程计算
gsvaPar <- gsvaParam(exp, genesets,maxDiff = TRUE)
gsvaPar 
gsva.res <- gsva(gsvaPar)
dim(gsva.res)
head(gsva.res)[1:5,1:5]
#                                                 GSM3106326  GSM3106327  GSM3106328  GSM3106329  GSM3106330
# KEGG_ABC_TRANSPORTERS                            0.07197716 -0.17846775  0.18517308 -0.15940843 -0.25379687
# KEGG_ACUTE_MYELOID_LEUKEMIA                     -0.05004563 -0.05015776  0.17257837 -0.33558735  0.09424981
# KEGG_ADHERENS_JUNCTION                          -0.14788713 -0.10257106 -0.18413147  0.05595831  0.25874815
# KEGG_ADIPOCYTOKINE_SIGNALING_PATHWAY            -0.12701380 -0.36331421 -0.27864321 -0.41398011 -0.10789732
# KEGG_ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM  0.04323269  0.29484579  0.02247517 -0.15112383 -0.16768683
write.csv(gsva.res,"gsva.res.csv")

2、limma差异分析

# 不同通路之间的limma差异分析
# 一定要设置分组信息！
identical(colnames(gsva.res),rownames(pd))
s <- intersect(colnames(gsva.res),rownames(pd))
gsva.res <- gsva.res[,s]
pd <- pd[s,]

Group <- factor(pd$cluster,levels = c("c1","c2"))
design = model.matrix(~Group)
fit = lmFit(gsva.res,design)
fit = eBayes(fit)
pathway_res = topTable(fit,coef = 2,number = Inf)
pathway_res[1:5,1:5]
#                                             logFC      AveExpr         t      P.Value    adj.P.Val
# KEGG_PROTEASOME                          0.4510362 -0.302452399  8.637118 2.410744e-13 4.483983e-11
# KEGG_FRUCTOSE_AND_MANNOSE_METABOLISM     0.3944620 -0.002586752  6.339069 9.818417e-09 9.131127e-07
# KEGG_RNA_POLYMERASE                      0.4186369 -0.053443215  6.074686 3.152025e-08 1.772463e-06
# KEGG_VASCULAR_SMOOTH_MUSCLE_CONTRACTION -0.2652879 -0.029894344 -6.011044 4.162805e-08 1.772463e-06
# KEGG_P53_SIGNALING_PATHWAY               0.2808839 -0.035073110  5.980057 4.764685e-08 1.772463e-06
write.csv(pathway_res,"pathway_res.csv")
# 标准自定啊！
range(pathway_res$logFC)
range(pathway_res$adj.P.Val)
pathway_sig <- pathway_res[pathway_res$adj.P.Val<0.001 & abs(pathway_res$logFC) > 0.2,]

3、可视化热图

###热图
library(pheatmap)
library(stringr)
# 这里又把排名靠前的通路提取出来，去跟gsva之后的数据取交集，得到不同样本下的通路数据
data <- gsva.res[match(rownames(pathway_sig),rownames(gsva.res)),]

#对数据进行排序
identical(rownames(pd),colnames(data))
Group <- factor(pd$cluster,levels = c("c1","c2"))
sort_data <- order(Group)
n <- data[,sort_data]
pd <- pd[sort_data,]

# 调整颜色梯度
range(data)
breaksList = seq(-0.5, 0.5, by = 0.1)
colors <- colorRampPalette(c("#336699", "white", "tomato"))(length(breaksList))

#创建列和行注释
annCol <- data.frame(group = pd$cluster,
                     row.names = colnames(data),
                     stringsAsFactors = FALSE)
#annRow <- data.frame(row.names = rownames(data))
pheatmap(data,
         annotation_col = annCol,
         #annotation_row = annRow,
         color = colors,
         breaks = breaksList,
         cluster_rows = T,
         cluster_cols = FALSE,
         show_rownames = TRUE,
         show_colnames = FALSE,
         #gaps_col = cumsum(table(annCol$Type)),  # 使用排序后的列分割点
         #gaps_row = cumsum(table(annRow$Methods)), # 行分割
         fontsize_row = 6,
         fontsize_col = 6,
         annotation_names_row = FALSE
         )