从单细胞转录组推断转录因子活性？decoupleR (四)：我又双叒叕会了

[ixxmu]

https://mp.weixin.qq.com/s/oM1-3SZXgPufZXmC5Z8B2A

[ixxmu]

从单细胞转录组推断转录因子活性？decoupleR (四)：我又双叒叕会了 by 生信益站

生信益站，一点就有益！「祝友友们天天开心，月月发 CNS~」

各位友友记得将益站设为星标哈

这样才能更容易在每天早上八点收到我们的推送哈~

关于 decoupleR，站长已经写了 3 篇推文了：

1. 哪些基因激活/抑制了你的通路？decoupleR (一)：这个我会
2. 如何用转录组推断转录因子活性？decoupleR (二)：这个我也会！
3. 从单细胞转录组推断通路活性？decoupleR（三）：这个我还会

今天给大家带来 decoupleR 第 4 篇——从 scRNA-seq 推断转录因子活性。

加载包

首先，我们需要加载相关包，Seurat 处理 scRNA-seq 数据，decoupleR 使用统计方法。

## We load the required packages
library(Seurat)
library(decoupleR)

# Only needed for data handling and plotting
library(dplyr)
library(tibble)
library(tidyr)
library(patchwork)
library(ggplot2)
library(pheatmap)

加载数据集

Seurat 在这里，我们使用了小插图中使用的数据的下采样版本 。我们可以这样打开数据：

inputs_dir <- system.file("extdata", package = "decoupleR")
data <- readRDS(file.path(inputs_dir, "sc_data.rds"))

此时，可以观察到不同的细胞类型：

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

CollecTRI 网络

CollecTRI 是一个综合资源，包含由 12 个不同资源编译而成的 TF 及其转录目标的精选集合。与我们之前的 DoRothEA 网络和其他基于文献的 GRN 相比，该集合提供了更多的转录因子覆盖范围，并且在识别受干扰的 TF 方面具有卓越的性能。与 DoRothEA 类似，相互作用按其调节模式（激活或抑制）进行加权。

在本例中，我们将使用人类版本（也可以使用小鼠和大鼠）。我们可以使用 decoupleR 从 OmniPath 检索它。该参数 split_complexes 保留复合体或将它们分成子单元，默认情况下我们建议将复合体保持在一起。

net <- get_collectri(organism='human', split_complexes=FALSE)
net
#> # A tibble: 43,178 × 3
#>    source target   mor
#>    <chr>  <chr>  <dbl>
#>  1 MYC    TERT       1
#>  2 SPI1   BGLAP      1
#>  3 SMAD3  JUN        1
#>  4 SMAD4  JUN        1
#>  5 STAT5A IL2        1
#>  6 STAT5B IL2        1
#>  7 RELA   FAS        1
#>  8 WT1    NR0B1      1
#>  9 NR0B2  CASP1      1
#> 10 SP1    ALDOA      1
#> # ℹ 43,168 more rows

使用单变量线性模型 (ULM) 进行活性推断

为了推断 TF 富集分数，我们将运行单变量线性模型 ( ulm) 方法。对于数据集中的每个样本 ( mat) 和网络中的每个 TF net，它都拟合一个线性模型，该模型仅根据 TF 的 TF-Gene 相互作用权重来预测观察到的基因表达。拟合后，获得的斜率 t 值就是分数。如果为正，我们认为 TF 处于活跃状态；如果为负，我们认为 TF 处于非活跃状态。

要运行 decoupleR 方法，我们需要一个输入矩阵 ( mat)、一个输入先验知识网络/资源 ( 上一步得到的 net) 以及我们要使用的网络列名。

# Extract the normalized log-transformed counts
mat <- as.matrix(data@assays$RNA@data)

# Run ulm
acts <- run_ulm(mat=mat, net=net, .source='source', .target='target',
                .mor='mor', minsize = 5)
acts
#> # A tibble: 80,640 × 5
#>    statistic source condition        score  p_value
#>    <chr>     <chr>  <chr>            <dbl>    <dbl>
#>  1 ulm       MYC    AAACATACAACCAC-1 13.5  3.54e-41
#>  2 ulm       MYC    AAACGCTGTTTCTG-1  8.09 6.78e-16
#>  3 ulm       MYC    AACCTTTGGACGGA-1 12.0  1.04e-32
#>  4 ulm       MYC    AACGCCCTCGTACA-1 10.4  3.11e-25
#>  5 ulm       MYC    AACGTCGAGTATCG-1 10.8  4.20e-27
#>  6 ulm       MYC    AACTCACTCAAGCT-1 10.1  5.15e-24
#>  7 ulm       MYC    AAGATGGAAAACAG-1 11.4  5.79e-30
#>  8 ulm       MYC    AAGATTACCGCCTT-1 13.4  1.64e-40
#>  9 ulm       MYC    AAGCCATGAACTGC-1 11.8  5.50e-32
#> 10 ulm       MYC    AAGGTCTGCAGATC-1 12.3  1.56e-34
#> # ℹ 80,630 more rows

可视化

根据获得的结果，我们将它们存储在我们的对象中作为一个新的 assay，称为 tfsulm：

# Extract ulm and store it in tfsulm in pbmc
data[['tfsulm']] <- acts %>%
  pivot_wider(id_cols = 'source', names_from = 'condition',
              values_from = 'score') %>%
  column_to_rownames('source') %>%
  Seurat::CreateAssayObject(.)

# Change assay
DefaultAssay(object = data) <- "tfsulm"

# Scale the data
data <- ScaleData(data)
data@assays$tfsulm@data <- data@assays$tfsulm@scale.data

这种新的 assay 可用于绘制活性图。在这里，我们观察了推断的 PAX5 跨细胞活性，它在 B 细胞中特别活跃。有趣的是，PAX5 是一种已知的对 B 细胞身份和功能至关重要的 TF。从目标基因的“足迹”推断活性比仅仅查看给定 TF 的分子读数更能提供信息，这里的一个例子是 PAX5 的基因表达，其本身的信息量并不大：

p1 <- DimPlot(data, reduction = "umap", label = TRUE, pt.size = 0.5) +
  NoLegend() + ggtitle('Cell types')
p2 <- (FeaturePlot(data, features = c("PAX5")) &
  scale_colour_gradient2(low = 'blue', mid = 'white', high = 'red')) +
  ggtitle('PAX5 activity')
DefaultAssay(object = data) <- "RNA"
p3 <- FeaturePlot(data, features = c("PAX5")) + ggtitle('PAX5 expression')
DefaultAssay(object = data) <- "tfsulm"
p1 | p2 | p3

左中右以此为Umap图、PAX5活性图、PAX5表达图

探索

我们还可以查看 top 25 个较大变异的 TF 每组的平均活性是多少：

n_tfs <- 25
# Extract activities from object as a long dataframe
df <- t(as.matrix(data@assays$tfsulm@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))

# Get top tfs with more variable means across clusters
tfs <- df %>%
  group_by(source) %>%
  summarise(std = sd(mean)) %>%
  arrange(-abs(std)) %>%
  head(n_tfs) %>%
  pull(source)

# Subset long data frame to top tfs and transform to wide matrix
top_acts_mat <- df %>%
  filter(source %in% tfs) %>%
  pivot_wider(id_cols = 'cluster', names_from = 'source',
              values_from = 'mean') %>%
  column_to_rownames('cluster') %>%
  as.matrix()

# Choose color palette
palette_length = 100
my_color = colorRampPalette(c("Darkblue", "white","red"))(palette_length)

my_breaks <- c(seq(-3, 0, length.out=ceiling(palette_length/2) + 1),
               seq(0.05, 3, length.out=floor(palette_length/2)))

# Plot
pheatmap(top_acts_mat, border_color = NA, color=my_color, breaks = my_breaks)

在这里，我们可以观察到其他已知标记 TF 的出现，EBF1 代表 B 细胞，RFX5 代表骨髓谱系，EOMES 代表淋巴系。

---

OK，今天的分享到此为止。咱们明天见~

联系站长

❝

对本篇文章有疑问，或者有科研服务需求的友友可以在益站发消息留言，也欢迎各位童鞋扫下面的二维码加入我们的 QQ 交流群。

❞

科研服务

最新QQ群

您的关注、点赞、在看、转发是对益站最大的鼓励和支持哈