专栏名称: 生信技能树
生物信息学学习资料分析,常见数据格式及公共数据库资料分享。常见分析软件及流程,基因检测及癌症相关动态。
目录
相关文章推荐
51好读  ›  专栏  ›  生信技能树

单细胞数据的GSVA

生信技能树  · 公众号  ·  · 2024-06-25 12:32

正文


单细胞数据的GSVA和芯片、bulk转录组的GSVA没有本质区别,就使用AverageExpression获取平均表达量得到新的表达矩阵再计算即可。

1.加载数据和R包

获得每种细胞的平均表达量

rm(list = ls())
library(Seurat)
library(GSVA)
library(clusterProfiler)
load("seu.obj.Rdata")
table(Idents(seu.obj))
#
##  Naive CD4 T   CD14+ Mono            B        CD8 T           NK FCGR3A+ Mono 
##         1675         1206          598          406          337          125 
##     Platelet           DC 
##           48           88

exp  =  AverageExpression(seu.obj)[[1]]
#exp =  AggregateExpression(seu.obj)[[1]]
exp  =  as.matrix(exp)
exp  =  exp[rowSums(exp)>0,] 
exp[1:4,1:4]

#
#          Naive CD4 T  CD14+ Mono          B      CD8 T
## TSPAN6    0.01890007 0.000000000 0.00000000 0.00446691
## DPM1      0.50764534 0.398461857 0.52602493 0.49951298
## SCYL3     0.10701976 0.049771894 0.10397003 0.12101561
## C1orf112  0.02653607 0.005093801 0.05426134 0.02747031

Seurat v5 提示建议用AggregateExpression做伪bulk转录组分析, 那个是用来求和的,目前查到的文献和教程都是使用平均值,这里就木有改动.

2.做GSVA

gmt文件下载自GSEA-msigdb官网

h_df = read.gmt("h.all.v2023.2.Hs.symbols.gmt")[,c(2,1)]
h_list = unstack(h_df)
ES  =  gsva(exp, h_list)
ES[1:4,1:4]

3.热图可视化

library(pheatmap)
pheatmap(ES, scale = "row",angle_col = "45",
         color = colorRampPalette(c("navy""white""firebrick3"))(50))

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI,多一点数据认知,让他们的科研上一个台阶:







请到「今天看啥」查看全文