单细胞富集分析系列：

• 单细胞之富集分析-1：单细胞GSEA分析流程
• 单细胞之富集分析-2：批量GSEA和GSVA分析
• 单细胞之富集分析-3：GO和KEGG富集分析及绘图
• 单细胞之富集分析-4：分组水平GSVA
• 单细胞之富集分析-5：一张图展示多组富集分析结果

0. 简介

异常细胞信号会引起癌症等其他疾病，并且是常见的治疗的靶点。常可以通过基因的表达来推断某个信号通路的活性。然而，只考虑基因表达对通路的作用往往忽略了翻译后修饰的作用，并且下游信号代表非常特定的实验条件。在这里，作者提出介绍PROGENy，这是一种通过利用大量公开可用的扰动实验，来克服了这两个局限性的方法。与现有方法不同，PROGENy可以（i）恢复已知驱动基因突变的作用，（ii）提供或改善药物的marker，以及（iii）区分致癌和肿瘤抑制途径，以确保患者生存。

PROGENy可以从基因表达数据中推断14种信号通路（雄激素，雌激素，EGFR，低氧，JAK-STAT，MAPK，NFkB，PI3K，p53，TGFb，TNFa，Trail，VEGF和WNT）的通路活性。默认情况下，途径活动推断是基于相应的途径扰动后前100个响应性最高的基因的基因集，我们将其称为途径的足迹基因。为每个足迹基因分配一个权重，该权重表示对路径扰动进行调节的强度和方向。途径得分是通过表达和足迹基因权重乘积的加权总和计算得出的。
例如：在pbmc单细胞数据中识别通路活性。

1. 分析Bulk数据

1.1 导入演示数据(芯片数据)

library(airway)
library(DESeq2)
data(airway)

# import data to DESeq2 and variance stabilize
dset = DESeqDataSetFromMatrix(assay(airway),
                              colData=as.data.frame(colData(airway)), design=~dex)
dset = estimateSizeFactors(dset)
dset = estimateDispersions(dset)
gene_expr = getVarianceStabilizedData(dset)

# annotate matrix with HGNC symbols
library(biomaRt)
mart = useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
genes = getBM(attributes = c("ensembl_gene_id","hgnc_symbol"),
              values=rownames(gene_expr), mart=mart)
matched = match(rownames(gene_expr), genes$ensembl_gene_id)
rownames(gene_expr) = genes$hgnc_symbol[matched]
dim(gene_expr)
# [1] 64102     8
gene_expr[1:5,1:5]
#          SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
# TSPAN6     9.531699   9.170317   9.674546   9.421197  10.027595
# TNMD       5.302948   5.302948   5.302948   5.302948   5.302948
# DPM1       9.055928   9.347026   9.235703   9.278863   9.166625
# SCYL3      8.358430   8.273162   8.212763   8.316611   8.136760
# C1orf112   6.967075   7.001886   6.596359   6.882393   7.060850

1.2 分析和绘图

library(progeny)
pathways = progeny(gene_expr, scale=FALSE)
controls = airway$dex == "untrt"
ctl_mean = apply(pathways[controls,], 2, mean)
ctl_sd = apply(pathways[controls,], 2, sd)
pathways = t(apply(pathways, 1, function(x) x - ctl_mean))
pathways = apply(pathways, 1, function(x) x / ctl_sd)

library(pheatmap)
myColor = colorRampPalette(c("Darkblue", "white","red"))(100)
pheatmap(t(pathways),fontsize=14, show_rownames = F,
         color=myColor, main = "PROGENy", angle_col = 45, treeheight_col = 0,  
         border_color = NA)

还可以加上分组标签，看起来会更直观

2. 分析单细胞数据

2.1 导入数据

library(Seurat)
library(progeny)
library(tidyr)
library(tibble)
pbmc <- readRDS("pbmc.rds") #注释好的单细胞数据集

2.2 提取细胞类型标签

查看每个细胞类型这14个通路的活性

## We create a data frame with the specification of the cells that belong to 
## each cluster to match with the Progeny scores.
CellsClusters <- data.frame(Cell = names(Idents(pbmc)), 
                            CellType = as.character(Idents(pbmc)),
                            stringsAsFactors = FALSE)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

2.3 计算

## We compute the Progeny activity scores and add them to our Seurat object
## as a new assay called Progeny. 
pbmc <- progeny(pbmc, scale=FALSE, organism="Human", top=500, perm=1, 
                return_assay = TRUE)
pbmc@assays$progeny
# Assay data with 14 features for 2638 cells
# First 10 features:
# Androgen, EGFR, Estrogen, Hypoxia, JAK-STAT, MAPK, NFkB, p53, PI3K, TGFb 

pbmc对象下面的assay槽下面多了progeny

pbmc@assays[["progeny"]]@data下储存的是14个通路在2638个细胞中的评分

pbmc@assays[["progeny"]]@data[1:6,1:4]
#          AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1
# Androgen         86.88406         57.92163        94.287118        11.135108
# EGFR             28.34279        -15.00708        -2.652041         1.738702
# Estrogen        -51.13204        -67.55945       -64.452782       -37.103230
# Hypoxia          96.56731        122.76491       135.750523        96.936620
# JAK-STAT        163.11016        269.55272       315.408767       609.408980
# MAPK            -10.39436        -59.37868       -17.161151       -61.763826

2.4 整合同一细胞类型的评分

## We can now directly apply Seurat functions in our Progeny scores. 
## For instance, we scale the pathway activity scores. 
pbmc <- Seurat::ScaleData(pbmc, assay = "progeny") 

## We transform Progeny scores into a data frame to better handling the results
progeny_scores_df <- 
  as.data.frame(t(GetAssayData(pbmc, slot = "scale.data", 
                               assay = "progeny"))) %>%
  rownames_to_column("Cell") %>%
  gather(Pathway, Activity, -Cell) 
dim(progeny_scores_df)
# [1] 36932     3

## We match Progeny scores with the cell clusters.
progeny_scores_df <- inner_join(progeny_scores_df, CellsClusters)

## We summarize the Progeny scores by cellpopulation
summarized_progeny_scores <- progeny_scores_df %>% 
  group_by(Pathway, CellType) %>%
  summarise(avg = mean(Activity), std = sd(Activity))
dim(summarized_progeny_scores)
# [1] 126   4

## We prepare the data for the plot
summarized_progeny_scores_df <- summarized_progeny_scores %>%
  dplyr::select(-std) %>%   
  spread(Pathway, avg) %>%
  data.frame(row.names = 1, check.names = FALSE, stringsAsFactors = FALSE) 

2.5 绘图

paletteLength = 100
myColor = colorRampPalette(c("Darkblue", "white","red"))(paletteLength)

progenyBreaks = c(seq(min(summarized_progeny_scores_df), 0, 
                      length.out=ceiling(paletteLength/2) + 1),
                  seq(max(summarized_progeny_scores_df)/paletteLength, 
                      max(summarized_progeny_scores_df), 
                      length.out=floor(paletteLength/2)))
progeny_hmap = pheatmap(t(summarized_progeny_scores_df),fontsize=12, 
                        fontsize_row = 10, 
                        color=myColor, breaks = progenyBreaks, 
                        main = "PROGENy (500)", angle_col = 45,
                        treeheight_col = 0,  border_color = NA)

2.6 绘制单个通路的FeaturePlot

DefaultAssay(pbmc) <- 'progeny'
p1= FeaturePlot(pbmc,features = "NFkB", coord.fixed = T, order = T, cols = viridis(10))
p2=FeaturePlot(pbmc,features = "MAPK", coord.fixed = T, order = T, cols = viridis(10))
p1|p2