生信分析工具包单细胞学习生信绘图

单细胞之富集分析-4:分组水平GSVA

2021-07-17  本文已影响0人  Hayley笔记

单细胞富集分析系列:


1. 数据和基因集准备
library(Seurat)
library(msigdbr)
library(GSVA)
library(tidyverse)
library(clusterProfiler)
library(patchwork)
library(limma)
rm(list=ls())
pbmc <- readRDS("pbmc.rds")
#DefaultAssay(scRNA) <- "RNA"
#scRNA <- NormalizeData(scRNA)
genesets <- msigdbr(species = "Homo sapiens", category = "H") 
genesets <- subset(genesets, select = c("gs_name","gene_symbol")) %>% as.data.frame()
genesets <- split(genesets$gene_symbol, genesets$gs_name)
Idents(pbmc) <- "cell_type" 
expr <- AverageExpression(pbmc, assays = "RNA", slot = "data")[[1]]
expr <- expr[rowSums(expr)>0,]  #选取非零基因
expr <- as.matrix(expr)
head(expr)
#               Naive CD4 T Memory CD4 T  CD14+ Mono          B      CD8 T FCGR3A+ Mono         NK         DC Platelet
# AL627309.1    0.006007987   0.04854338 0.006065400 0.00000000 0.01995673   0.00000000 0.00000000 0.00000000        0
# AP006222.2    0.000000000   0.01088471 0.008397321 0.00000000 0.01157323   0.00000000 0.00000000 0.00000000        0
# RP11-206L10.2 0.007306336   0.00000000 0.000000000 0.02065031 0.00000000   0.00000000 0.00000000 0.08462847        0
# RP11-206L10.9 0.000000000   0.01050116 0.000000000 0.00000000 0.00000000   0.01200008 0.00000000 0.00000000        0
# LINC00115     0.014872162   0.03753737 0.031095957 0.03888541 0.01892413   0.01469374 0.06302423 0.00000000        0
# NOC2L         0.501822970   0.27253750 0.346874253 0.58653489 0.59129394   0.50026775 0.65705381 0.32861136        0
2. 分组GSVA分析
# gsva默认开启全部线程计算
gsva.res <- gsva(expr, genesets, method="gsva") 
saveRDS(gsva.res, "gsva.res.rds")
gsva.df <- data.frame(Genesets=rownames(gsva.res), gsva.res, check.names = F)
write.csv(gsva.df, "gsva_res.csv", row.names = F)
pheatmap::pheatmap(gsva.res, show_colnames = T, 
                   scale = "row",angle_col = "45",
                   color = colorRampPalette(c("navy", "white", "firebrick3"))(50))
3. 计算通路在case和control之间的显著性并绘图(limma做的)

pbmc3k的数据是有一个样本,在这里我们把9种细胞当成4个con5个病例组来计算case和control之间某通路的显著性。(仅供演示,实际操作时候前面计算每个sample的基因表达均值,这里再换成case和control就可以了)

group_list <- data.frame(sample = colnames(gsva.df)[-1], group = c(rep("con", 4), rep("case", 5)))
group_list
#         sample group
# 1  Naive CD4 T   con
# 2 Memory CD4 T   con
# 3   CD14+ Mono   con
# 4            B   con
# 5        CD8 T  case
# 6 FCGR3A+ Mono  case
# 7           NK  case
# 8           DC  case
# 9     Platelet  case
design <- model.matrix(~ 0 + factor(group_list$group))
colnames(design) <- levels(factor(group_list$group))
rownames(design) <- colnames(gsva.res)
design
#              case con
# Naive CD4 T     0   1
# Memory CD4 T    0   1
# CD14+ Mono      0   1
# B               0   1
# CD8 T           1   0
# FCGR3A+ Mono    1   0
# NK              1   0
# DC              1   0
# Platelet        1   0
# attr(,"assign")
# [1] 1 1
# attr(,"contrasts")
# attr(,"contrasts")$`factor(group_list$group)`
# [1] "contr.treatment"
# 构建差异比较矩阵
contrast.matrix <- makeContrasts(con-case, levels = design)

# 差异分析,case vs. con
fit <- lmFit(gsva.res, design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
x <- topTable(fit2, coef = 1, n = Inf, adjust.method = "BH", sort.by = "P")
head(x)
#                                       logFC     AveExpr         t     P.Value adj.P.Val          B
# HALLMARK_MTORC1_SIGNALING        -0.3291987 -0.08593222 -4.477744 0.000747704 0.0373852 -0.2823204
# HALLMARK_PI3K_AKT_MTOR_SIGNALING -0.2988238 -0.16659423 -3.391304 0.005325363 0.1331341 -2.1754501
# HALLMARK_PEROXISOME              -0.1862704 -0.08599750 -2.611784 0.022660475 0.2801731 -3.5499137
# HALLMARK_PROTEIN_SECRETION       -0.2411741 -0.06878214 -2.523986 0.026641622 0.2801731 -3.7008368
# HALLMARK_CHOLESTEROL_HOMEOSTASIS -0.2317843 -0.08879774 -2.496598 0.028017314 0.2801731 -3.7476296
# HALLMARK_FATTY_ACID_METABOLISM   -0.2112693 -0.10325132 -2.242946 0.044472826 0.3706069 -4.1730598

#把通路的limma分析结果保存到文件
write.csv(x, "gsva_limma.csv", quote = F)

#输出t值,用做绘图的输入数据
pathway <- str_replace(row.names(x), "HALLMARK_", "")
df <- data.frame(ID = pathway, score = x$t)
head(df)
#                        ID    score
# 1        MTORC1_SIGNALING 4.477744
# 2 PI3K_AKT_MTOR_SIGNALING 3.391304
# 3              PEROXISOME 2.611784
# 4       PROTEIN_SECRETION 2.523986
# 5 CHOLESTEROL_HOMEOSTASIS 2.496598
# 6   FATTY_ACID_METABOLISM 2.242946
write.csv(df, "enrich_bar.csv", quote = F, row.names = F)
#按照score的值分组
cutoff <- 0
df$group <- cut(df$score, breaks = c(-Inf, cutoff, Inf),labels = c(1,2))
head(df)
#                        ID    score group
# 1        MTORC1_SIGNALING -4.477744     2
# 2 PI3K_AKT_MTOR_SIGNALING -3.391304     2
# 3              PEROXISOME -2.611784     2
# 4       PROTEIN_SECRETION -2.523986     2
# 5 CHOLESTEROL_HOMEOSTASIS -2.496598     2
# 6   FATTY_ACID_METABOLISM -2.242946     2

#按照score排序
sortdf <- df[order(df$score),]
sortdf$ID <- factor(sortdf$ID, levels = sortdf$ID)
head(sortdf)
#                        ID     score group
# 1        MTORC1_SIGNALING -4.477744     1
# 2 PI3K_AKT_MTOR_SIGNALING -3.391304     1
# 3              PEROXISOME -2.611784     1
# 4       PROTEIN_SECRETION -2.523986     1
# 5 CHOLESTEROL_HOMEOSTASIS -2.496598     1
# 6   FATTY_ACID_METABOLISM -2.242946     1
ggplot(sortdf, aes(ID, score, fill = group)) + geom_bar(stat = 'identity') + 
  coord_flip() + 
  scale_fill_manual(values = c('palegreen3', 'dodgerblue4'), guide = FALSE) + 
  #画2条虚线
  geom_hline(yintercept = c(-1,1), 
             color="white",
             linetype = 2, #画虚线
             size = 0.3) + #线的粗细
  #写label
  geom_text(data = subset(df, score < 0),
            aes(x=ID, y= 0.1, label=ID, color = group),
            size = 3, #字的大小
            hjust = "outward" ) +  #字的对齐方式
  geom_text(data = subset(df, score > 0),
            aes(x=ID, y= -0.1, label= paste0(" ", ID), color = group),#bar跟坐标轴间留出间隙
            size = 3, hjust = "inward") +  
  scale_colour_manual(values = c("black","black"), guide = FALSE) +
  
  xlab("") +ylab("t value of GSVA score")+
  theme_bw() + #去除背景色
  theme(panel.grid =element_blank()) + #去除网格线
  theme(panel.border = element_rect(size = 0.6)) + #边框粗细
  theme(axis.line.y = element_blank(), axis.ticks.y = element_blank(), axis.text.y = element_blank()) #去除y轴

ggsave("gsva.pdf", width = 6, height = 8)
上一篇下一篇

猜你喜欢

热点阅读