DecoupleR/CollecTRI network-单细胞转录因子活性分析

DecoupleR包是一个可以从组学数据中分析内部生物学活性的计算方法集合。该R包内部收纳了11种分析方法，包括GSVA,GSEA,univariate linear model (ULM),VIPER等。

PMID：36699385

根据包含的计算方法，我们也能大概的知道decoupleR包能够做哪些分析了，包括基本分析(AUCell，Fast GSEA, GSVA和viper)，bulk/scRNA-seq的通路活性、转录因子分析，转录因子和激酶活性分析。这次我们主要学习一下如何使用该R包进行scRNA-seq的转录因子活性分析。

github链接：https://saezlab.github.io/decoupleR/index.html

那么在进行分析之前，首先需要理解一下"转录因子在特定细胞中的活性水平"这个概念。粗略来说可以理解为转录因子(TFs)对其靶基因发挥调控潜能的程度，所以这个调控潜能的程度高低并不一定意味着跟转录因子的表达量有密切的正相关，而是TFs跟这个细胞中所表达的基因中的TFs靶基因的表达情况的“相关性打分”。

PMID：33135076

从代码分析的角度来看，为了推断TFs的活性分数，需要采用单变量线性模型(ulm)方法。对于数据集中的每个样本(每一个单细胞)和TFs数据集中的每个TFs拟合一个线性模型，该模型是基于TFs的TFs-基因的相互作用权重来预测观察到的基因表达。拟合后，得到的斜率t值就是得分。如果它是阳性的，就解释为TFs是活跃的，如果它是阴性的，就认为它是非活跃的。

此外，我们还需要获得TFs的调控子(可以认为是TFs-靶基因集群)的数据集。目前推荐使用CollecTRI(比DoRothEA更全面),这个资源整合了多个来源的数据集，目前纳入了1186个TFs和43175个有一定证据的TF-基因相互作用数据。

PMID：37843125

具体分析流程

1.数据读取和加载R包

我这里采用的GSE148190-SM4455931-10X样本

Step1.final.Rdata已经是经过了降维、聚类、分群后的结果。

rm(list=ls())
library(Seurat)
library(decoupleR)
library(dplyr)
library(tibble)
library(tidyr)
library(patchwork)
library(ggplot2)
library(pheatmap)

load("step1.final.Rdata")
sce <- step1.final
table(Idents(sce))

2.先看一下聚类图

DimPlot(sce,label = T,repel = T) 

3.转录因子分析

#多线程模式
plan("multisession", workers = 4)
#获取转录因子数据
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
# ℹ Use `print(n = ...)` to see more rows

#提取数据
mat <- as.matrix(sce@assays$RNA@data)
#下面的运算会耗费一点时间，取决于线程数量和CPU算力
acts <- run_ulm(mat, net,minsize = 5,
                .source='source', .target='target',.mor='mor')
acts
# A tibble: 1,933,792 × 5
#   statistic source condition             score p_value
#   <chr>     <chr>  <chr>                 <dbl>   <dbl>
# 1 ulm       ABL1   AAACCTGAGGACGAAA-1  0.274   0.784  
# 2 ulm       ABL1   AAACCTGAGGGTGTTG-1  0.724   0.469  
# 3 ulm       ABL1   AAACCTGCAAGGACTG-1  2.29    0.0221 
# 4 ulm       ABL1   AAACCTGCAATCTGCA-1 -0.423   0.673  
# 5 ulm       ABL1   AAACCTGCACTGTTAG-1 -0.00305 0.998  
# 6 ulm       ABL1   AAACCTGGTCTAGTGT-1  3.01    0.00262
# 7 ulm       ABL1   AAACCTGTCGTACCGG-1  0.550   0.582  
# 8 ulm       ABL1   AAACGGGAGGCAATTA-1  1.51    0.131  
# 9 ulm       ABL1   AAACGGGAGGTGCTTT-1  0.312   0.755  
#10 ulm       ABL1   AAACGGGCAGCCTATA-1  0.896   0.370  
# ℹ 1,933,782 more rows
# ℹ Use `print(n = ...)` to see more rows

#提取获得的数据,修改格式并保存
sce[['tfsulm']] <- acts %>%
  pivot_wider(id_cols = 'source', 
              names_from = 'condition',
              values_from = 'score') %>%
  column_to_rownames('source') %>%
  Seurat::CreateAssayObject(.)
#更改默认的分析对象
DefaultAssay(object = sce) <- "tfsulm"
#缩放数据并重新保存
sce <- ScaleData(sce)
sce@assays$tfsulm@data <- sce@[email protected]

4.可视化

p1 <- DimPlot(sce, reduction = "umap", 
              label = TRUE, pt.size = 0.5) + 
      NoLegend() + ggtitle('Cell types')
p1

p2 <- (FeaturePlot(sce, features = c("E2F1")) & 
  scale_colour_gradient2(low = 'blue', mid = 'white', high = 'red')) +
  ggtitle('E2F1 activity')
p2

DefaultAssay(object = sce) <- "RNA"
p3 <- FeaturePlot(sce, features = c("E2F1")) + 
          ggtitle('PAX5 expression')
p3

DefaultAssay(object = sce) <- "tfsulm"
p1 + p2 + p3

挑选了一个很常见的转录因子E2F1，从表达量分析看它似乎”价值不大“，但从TFs活性分析来看其在CD4 T,CD8 T和B细胞中都显示出了较大的活性。

5.top25的TFs

n_tfs <- 25
#提取活性数据并修改数据格式
df <- t(as.matrix(sce@assays$tfsulm@data)) %>%
  as.data.frame() %>%
  mutate(cluster = Idents(sce)) %>%
  pivot_longer(cols = -cluster, #指定除cluster列以外的都需要转换
               names_to = "source", #列名放到source列中
               values_to = "score") %>% #值放到score列中
  group_by(cluster, source) %>% #将数据按cluster和source两列进行分组
  summarise(mean = mean(score))

#获得不同簇中活性差异较大的TFs
tfs <- df %>%
  group_by(source) %>%
  summarise(std = sd(mean)) %>%
  arrange(-abs(std)) %>%
  head(n_tfs) %>%
  pull(source)

#修改数据格式
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()

#绘图
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)))

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

参考文献：

1.Pau Badia-I-Mompel, Jesús Vélez Santiago, Jana Braunger et al. decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinform Adv. 2022 Mar 8;2(1):vbac016. doi: 10.1093/bioadv/vbac016 IF: NA NA NA.

2.Cynthia Z Ma, Michael R Brent et al. Inferring TF activities and activity regulators from gene expression data with constraints from TF perturbation data. Bioinformatics. 2021 Jun 9;37(9):1234-1245.

3.Sophia Müller-Dott, Eirini Tsirvouli, Miguel Vazquez et al. Expanding the coverage of regulons from high-confidence prior knowledge for accurate estimation of transcription factor activities. Nucleic Acids Res. 2023 Nov 10;51(20):10934-10949.

注：若对内容有疑惑或者发现有明确错误的朋友，请联系后台(希望多多交流)。更多内容可关注公众号：生信方舟

- END -