decoupleR推断单细胞相关通路及转录因子活性（内含代码）

 

组学数据往往会有着高维度大样本量的特点，很难在机制上用一个合适的作用方式去解读，利用算法在先验知识的基础上去推断相关的生物活性是一种常用的研究手段。decoupleR是一种利用组学定量数据结合先验知识推断生物活性的集成算法，其综合了不同的统计方法，其中包括AUCell、GSEA、GSVA、VIPER等算法，在实际应用中相较于其他方法表现突出。decoupleR有R和Python两个版本，下面将讲述decoupleR基于R对单细胞定量数据进行通路及转录因子活性推断分析。

 

 

 

decoupleR对单细胞数据进行通路活性推断

一、先验知识：PROGENy模型简介

PROGENy (Pathway RespOnsive GENes for activity )：从大量公开的扰动实验中获取关键通路的核心共同基因，并精选出与癌症、疾病相关信号通路形成的资源库，最初PROGENy只收集了11个人相关的通路，随着数据的更新，现阶段一共有14个人和小鼠相关的通路及核心基因，通路如下：

 

 

二、decoupleR安装

if (!require("BiocManager", quietly = TRUE))

    install.packages("BiocManager")
BiocManager::install("decoupleR")

 

三、使用decoupleR获取PROGENy通路

library(decoupleR) 

net <- get_progeny(organism = 'human', top = 500)

注释：

第一列：PROGENy通路。

第二列：通路中核心基因。

第三列：通路与基因间的权重。

第四列：通路与基因间权重显著性Pvalue。

 

四、多元线性模型推断通路活性

i. 算法简介：

基因表达定量矩阵和通路与基因集间的权重做为输入数据，使用多元线性模型去推断每个细胞（观测值）与通路间的活性，在模型拟合中获取的线性模型斜率的t-values做为细胞（观测值）和通路间活性打分（enrichment score），活性打分为正理解为通路路在细胞中是活跃的，反之不活跃。

 

ii. 使用10X PBMCs单细胞抽样数据进行通路活性推断：

library(Seurat)

inputs_dir <- system.file("extdata", package = "decoupleR")

data <- readRDS(file.path(inputs_dir, "sc_data.rds"))

mat <- as.matrix(data@assays$RNA@data)

acts <- run_mlm(mat=mat, net=net, .source='source', .target='target',

                .mor='weight', minsize = 5)

acts

 

注释：

第一列：多元线性模型标签。

第二列：PROGENy通路。

第三列：细胞barcode。

第四列：通路在细胞中enrichment score。

第五列：通路在细胞中enrichment score 显著性pvalue。

 

iii. 结果展示：

1. 将enrichment score 嵌入Seurat对象中。

library(dplyr)

library(tibble)

library(tidyr)

library(patchwork)

library(ggplot2)

data[['pathwaysmlm']] <- acts %>%

  pivot_wider(id_cols = 'source', names_from = 'condition',

              values_from = 'score') %>%

  column_to_rownames('source') %>%

  Seurat::CreateAssayObject(.)

DefaultAssay(object = data) <- "pathwaysmlm"

data <- ScaleData(data)

data@assays$pathwaysmlm@data <- data@assays$pathwaysmlm@scale.data

 

2. 筛选Trail通路活性使用Seurat展示。

p1 <- DimPlot(data, reduction = "umap", label = TRUE, pt.size = 0.5) + 

  NoLegend() + ggtitle('Cell types')

p2 <- (FeaturePlot(data, features = c("Trail")) & 

  scale_colour_gradient2(low = 'blue', mid = 'white', high = 'red')) +

  ggtitle(' Trail activity')

从上图可以看出，Trail通路主要在B细胞和FCGR3A+Mono细胞中处于活性较高。

 

3. 展示每个通路在细胞类型中平均活性。

df <- t(as.matrix(data@assays$pathwaysmlm@data)) %>%

  as.data.frame() %>%

  mutate(cluster = Idents(data)) %>%

  pivot_longer(cols = -cluster, names_to = "source", values_to = "score") %>%

  group_by(cluster, source) %>%

  summarise(mean = mean(score))

top_acts_mat <- df %>%

  pivot_wider(id_cols = 'cluster', names_from = 'source',

              values_from = 'mean') %>%

  column_to_rownames('cluster') %>%

  as.matrix()

write.csv(top_acts_mat,"PathwayScore.csv")

 

PathwayScore.csv为通路在每个细胞类型中平均的enrichment score 打分结果，下面使用SBCToolBox云平台对数据进行结果展示，将此文件上传至HeatMap App：

https://www4.shbiochip.com/V2023/HtM/HTM/

将数据上传至SBCToolBox热图程序，提交分析即可得到可视化结果：

 

从热图上可以直观的看出哪些通路在哪些细胞类型中活性高低。decoupleR对通路的活性推断依赖于先验知识：PROGENy数据库，对于癌症和疾病相关通路的研究是相对友好的，对于其他类型通路的研究则存在一定的局限性。

 

 

decoupleR对单细胞数据进行转录因子活性推断

 

decoupleR转录因子活性推断与通路活性推断类似，同样基于先验知识，这里使用的是CollecTRI（Collection of Transcriptional Regulatory Interactions）基因调控网络数据库，其收集了公共数据库以及文本挖掘出的高质量且可靠的TF regulons（regulon：受同一个转录因子调控的靶基因集），然后根据先验知识、转录因子角色及其他方式，对每一对TF-TargetGene调控关系进行权衡并划分为激活和抑制两种模式。

 

 

单变量线性模型推测转录因子活性

 

一、使用decoupleR调用CollecTRI调控网络数据

 

net <- get_collectri(organism='human', split_complexes=FALSE)

 

注释：

第一列：转录因子。

第二列：转录因子靶基因。

第三列：转录因子调控模式。

 

二、使用10X PBMCs单细胞抽样数据推测转录因子活性

i. 模型简介

每一个样品（细胞）表达数据与调控网络中每一个TF-TargetGene互作权重做为输出数据，使用单变量线性模型对转录因子活性进行推测，模型拟合得到的斜率对应的t-value做为转录因子富集得分，t-value为正值为激活，反之抑制。

 

ii. 单细胞转录因子活性预测

acts <- run_ulm(mat=mat, net=net, .source='source', .target='target',

                .mor='mor', minsize = 5)

注释：

第一列：为单变量线性模型标签。

第二列：转录因子。

第三例：细胞。

第四列：转录因子enrichment score。

第五列：转录因子enrichment score显著性pvalue。

 

iii. 转录因子活性展示

转录因子活性预测输出结果与通路活性输出结果结构一致，同样的逻辑，可以将转录因子活性打分内嵌到Seurat对象中，也可以计算出其在每个细胞类型中平均的转录活性。

展示转录因子PAX5在细胞群中的转录活性及转录表达情况。

 

 

筛选出25个转录因子，计算其在细胞类型中平均转录活性，使用SBCToolBox进行展示。

 

上图横向为转录因子，纵向为细胞类型，从上图可以清晰的看出转录因子在细胞类型中的活性高低。

 

参考：

1. Pau Badia-i-Mompel, Jesús Vélez Santiago, Jana Braunger, Celina Geiss, Daniel Dimitrov, Sophia Müller-Dott, Petr Taus, Aurelien Dugourd, Christian H Holland, Ricardo O Ramirez Flores, Julio Saez-Rodriguez, decoupleR: ensemble of computational methods to infer biological activities from omics data, Bioinformatics Advances, Volume 2, Issue 1, 2022, vbac016, https://doi.org/10.1093/bioadv/vbac016

2. Sophia Müller-Dott, Eirini Tsirvouli, Miguel Vazquez, Ricardo O Ramirez Flores, Pau Badia-i-Mompel, Robin Fallegger, Dénes Türei, Astrid Lægreid, Julio Saez-Rodriguez, Expanding the coverage of regulons from high-confidence prior knowledge for accurate estimation of transcription factor activities, Nucleic Acids Research, Volume 51, Issue 20, 10 November 2023, Pages 10934–10949, https://doi.org/10.1093/nar/gkad841

3. Schubert, M., Klinger, B., Klünemann, M. et al. Perturbation-response genes reveal signaling footprints in cancer gene expression. Nat Commun 9, 20 (2018). https://doi.org/10.1038/s41467-017-02391-6