使用decoupleR一次性实现11种基因集的活性打分（R与Python我都要）

就是因为考虑到绝大部分小伙伴是Python和R编程语言的二选一，所以为了自己的工具使用更广泛，很多开发者会特意分发不同版本的软件。其中对基因集等进行活性打分推断的包实在是太多了，但是同时开发了R代码和python代码的应该不多见，且集合了11种打分方法，我们推荐decoupleR：

官网github链接：https://github.com/saezlab/decoupleR

文献信息：

Badia-i-Mompel P.,。。。。2022. decoupleR: ensemble of computational methods to infer biological activities from omics data. Bioinformatics Advances. https://doi.org/10.1093/bioadv/vbac016

11个被囊括的包如下：

decoupleR包的功能：

• bulk RNA-seq 中通路活性推断
• scRNA-seq 中通路活性推断
• bulk RNA-seq 中转录因子活性推断
• scRNA-seq 中转录因子活性推断
• 激酶和转录因子活性估计

安装此包：

options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
options("repos"=c(CRAN="https://mirrors.westlake.edu.cn/CRAN/"))

install.packages('remotes')
remotes::install_github('saezlab/decoupleR')

下面来简单看一下其中一个功能：Pathway activity inference in bulk RNA-seq

https://saezlab.github.io/decoupleR/articles/pw_bk.html

小试牛刀：bulk RNA-seq 中通路活性推断

在这个笔记本中，我们展示了如何使用decoupleR对一个bulk RNA测序数据集进行通路活性推断，该数据集中胰腺癌细胞系中的转录因子FOXA2被敲除。

数据包括3个野生型（WT）样本和3个敲除型（KO）样本，GEO编号为：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE119931，是一个表观调控的项目，其中转录组数据部分样品信息如下所示：

1、首先加载示例数据

作者已经对此数据集进行了预处理，并放到了包中：

# 加载
inputs_dir "extdata", package = "decoupleR")
data "bk_data.rds"))

# 简单探索
class(data)
str(data)

可以看到数据为一个list对象，包含三个数据：counts值（为the normalized log-transformed counts）、design分组矩阵、limma_ttop差异基因结果

List of 3
 $ counts    : spc_tbl_ [11,093 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)



    
  ..$ gene              : chr [1:11093] "NOC2L" "PLEKHN1" "PERM1" "ISG15" ...
  ..$ PANC1.WT.Rep1     : num [1:11093] 10.05 7.54 6.28 10.94 6.96 ...
  ..$ PANC1.WT.Rep2     : num [1:11093] 11.95 8.13 6.42 11.47 7.2 ...
  ..$ PANC1.WT.Rep3     : num [1:11093] 12.06 8.71 6.59 11.43 7.52 ...
  ..$ PANC1.FOXA2KO.Rep1: num [1:11093] 12.31 8.05 6.29 11.55 7.06 ...
  ..$ PANC1.FOXA2KO.Rep2: num [1:11093] 12.14 8.29 6.49 11.37 7.49 ...
  ..$ PANC1.FOXA2KO.Rep3: num [1:11093] 11.49 8.62 6.78 11.18 7.07 ...
  ..- attr(*, "spec")=
  .. .. cols(
  .. ..   gene = col_character(),
  .. ..   PANC1.WT.Rep1 = col_double(),
  .. ..   PANC1.WT.Rep2 = col_double(),
  .. ..   PANC1.WT.Rep3 = col_double(),
  .. ..   PANC1.FOXA2KO.Rep1 = col_double(),
  .. ..   PANC1.FOXA2KO.Rep2 = col_double(),
  .. ..   PANC1.FOXA2KO.Rep3 = col_double()
  .. .. )
 $ design    : spc_tbl_ [6 × 2] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
  ..$ sample   : chr [1:6] "PANC1.WT.Rep1" "PANC1.WT.Rep2" "PANC1.WT.Rep3" "PANC1.FOXA2KO.Rep1" ...
  ..$ condition: chr [1:6] "PANC1.WT" "PANC1.WT" "PANC1.WT" "PANC1.FOXA2KO" ...
  ..- attr(*, "spec")=
  .. .. cols(
  .. ..   sample = col_character(),
  .. ..   condition = col_character()
  .. .. )
 $ limma_ttop: spc_tbl_ [11,093 × 7] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
  ..$ ID       : chr [1:11093] "RHBDL2" "PLEKHH2" "HEG1" "CLU" ...
  ..$ logFC    : num [1:11093] -1.82 -1.57 -1.73 -1.79 2.09 ...
  ..$ AveExpr  : num [1:11093] 8.6 7.7 8.64 12.22 9.61 ...
  ..$ t        : num [1:11093] -12.81 -10.79 -9.79 -9.76 8.95 ...
  ..$ P.Value  : num [1:11093] 0.00000303 0.00000993 0.0000194 0.00001976 0.00003552 ...
  ..$ adj.P.Val: num [1:11093] 0.0336 0.0548 0.0548 0.0548 0.0788 ...
  ..$ B        : num [1:11093] 3.95 3.27 2.84 2.83 2.43 ...
  ..- attr(*, "spec")=
  .. .. cols(
  .. ..   ID = col_character(),
  .. ..   logFC = col_double(),
  .. ..   AveExpr = col_double(),
  .. ..   t = col_double(),
  .. ..   P.Value = col_double(),
  .. ..   adj.P.Val = col_double(),
  .. ..   B = col_double()
  .. .. )

提取count：

# Remove NAs and set row names
counts $counts %>%
  dplyr::mutate_if(~ any(is.na(.x)), 
                   ~ dplyr::if_else(is.na(.x), 0, .x)) %>% 
  tibble::column_to_rownames(var = "gene") %>% 
  as.matrix()

head(counts)

提取The design meta-data:

design $design
design

提取limma结果：会用到结果中的t值

# Extract t-values per gene
deg $limma_ttop %>%
  dplyr::select(ID, t) %>% 
  dplyr::filter(!is.na(t)) %>% 
  tibble::column_to_rownames(var = "ID") %>%
  as.matrix()

head(deg)

2、评分使用基因集PROGENy

访问网址：https://saezlab.github.io/progeny/

PROGENy是一个包含了经过整理的通路及其靶基因的集合，并且为每个相互作用提供了权重。在这个例子中，我们将使用人类权重（也提供了其他生物体的权重），并且我们将使用按p值排名的前500个responsive genes。以下是每个通路的简要描述：

• 雄激素（Androgen）：参与男性生殖器官的生长和发育。
• 表皮生长因子受体（EGFR）：在哺乳动物细胞中调节生长、存活、迁移、凋亡、增殖和分化。
• 雌激素（Estrogen）：促进女性生殖器官的生长和发育。
• 缺氧（Hypoxia）：在氧气水平低时促进血管生成和代谢重编程。
• JAK-STAT 信号通路：涉及免疫、细胞分裂、细胞死亡和肿瘤形成。
• 丝裂原活化蛋白激酶（MAPK）：整合外部信号并促进细胞生长和增殖。
• 核因子κB（NFkB）：调节免疫反应、细胞因子产生和细胞存活。
• 肿瘤蛋白p53（p53）：调节细胞周期、凋亡、DNA修复和肿瘤抑制。
• 磷脂酰肌醇3激酶（PI3K）：促进生长和增殖。
• 转化生长因子β（TGFb）：涉及大多数组织的发展、稳态和修复。
• 肿瘤坏死因子α（TNFa）：介导造血、免疫监视、肿瘤退化和防止感染。
• TRAIL：诱导凋亡。
• 血管内皮生长因子（VEGF）：介导血管生成、血管通透性和细胞迁移。
• WNT 信号通路：在发育过程中调节器官形态发生和组织修复。

获取：

    
# 这里需要安装一下OmnipathR包
# BiocManager::install("OmnipathR")

net 'human', top = 500)
net
head(net)
table(net$source)

# 保存一下，这里不是很好下载，后面看看是不是有什么地方可以加速这个下载
saveRDS(net,file = "net.rds")

3、使用decoupleR评估PROGENy基因集在此数据中的活性

活性推断采用的模型为：Multivariate Linear Model (MLM)，MLM在文章中方法的可靠性排前三。推断的打分如果为positive，表示 the pathway is active；如果是negative，则inactive。

# Run mlm
sample_acts                                   net = net, 
                                  .source = 'source', 
                                  .target = 'target',
                                  .mor = 'weight', 
                                  minsize = 5)
sample_acts

结果如下：

4、打分结果可视化

首先提取每个样本中每条通路的score：

# Transform to wide matrix
sample_acts_mat %
                   tidyr::pivot_wider(id_cols = 'condition', 
                                      names_from = 'source',
                                      values_from = 'score') %>%
                   tibble::column_to_rownames('condition') %>%
                   as.matrix()

sample_acts_mat

# Scale per feature
sample_acts_mat 

打分矩阵：

绘图：

# Color scale
colors "RdBu"))
colors.use 
my_breaks                seq(0.05,2, length.out = floor(100 / 2)))

# Plot
pheatmap::pheatmap(mat = sample_acts_mat,
                   color = colors.use,
                   border_color = "white",
                   breaks = my_breaks,
                   cellwidth = 20,
                   cellheight = 20,
                   treeheight_row = 20,
                   treeheight_col = 20)

打分结果：

5、也可以从差异表达基因（DEGs）的t值推断通路活性

# Run mlm
contrast_acts                                     net = net, 
                                    .source = 'source', 
                                    .target = 'target',
                                    .mor = 'weight', 
                                    minsize = 5)

contrast_acts

# 可视化
colors "RdBu")[c(2, 10)])

p                      mapping = ggplot2::aes(x = stats::reorder(source, score), 
                                            y = score)) + 
     ggplot2::geom_bar(mapping = ggplot2::aes(fill = score),
                       color = "black",
                       stat = "identity") +
     ggplot2::scale_fill_gradient2(low = colors[1], 
                                   mid = "whitesmoke", 
                                   high = colors[2], 



    
                                   midpoint = 0) + 
     ggplot2::theme_minimal() +
     ggplot2::theme(axis.title = element_text(face = "bold", size = 12),
              axis.text.x = ggplot2::element_text(angle = 45, 
                                                  hjust = 1, 
                                                  size = 10, 
                                                  face = "bold"),
              axis.text.y = ggplot2::element_text(size = 10, 
                                                  face = "bold"),
              panel.grid.major = element_blank(), 
              panel.grid.minor = element_blank()) +
     ggplot2::xlab("Pathways")

p

结果：与野生型（WT）相比，敲除型（KO）中的p53和TRAIL通路被抑制，而MAPKK和JAK-STAT通路似乎被激活，与上面的每个样本中活性打分热图基本一致。

我们可以进一步沿着它们的t值可视化每个通路中最有反应的基因，以解释结果。例如，让我们来看看属于MAPK通路的基因：

pathway 'MAPK'

df %
      dplyr::filter(source == pathway) %>%
      dplyr::arrange(target) %>%
      dplyr::mutate(ID = target, 
                    color = "3") %>%
      tibble::column_to_rownames('target')

inter 
df 
df['t_value'] 
df %
      dplyr::mutate(color = dplyr::if_else(weight > 0 & t_value > 0, '1', color)) %>%
      dplyr::mutate(color = dplyr::if_else(weight > 0 & t_value '2', color)) %>%
      dplyr::mutate(color = dplyr::if_else(weight  0, '2', color)) %>%
      dplyr::mutate(color = dplyr::if_else(weight '1', color))

colors "RdBu")[c(2, 10)])

p                      mapping = ggplot2::aes(x = weight, 
                                            y = t_value, 
                                            color = color)) + 
     ggplot2::geom_point(size = 2.5, 
                         color = "black") + 
     ggplot2::geom_point(size = 1.5) +
     ggplot2::scale_colour_manual(values = c(colors[2], colors[1], "grey")) +
     ggrepel::geom_label_repel(mapping = ggplot2::aes(label = ID)) + 
     ggplot2::theme_minimal() +
     ggplot2::theme(legend.position = "none") +
     ggplot2::geom_vline(xintercept = 0, linetype = 'dotted') +
     ggplot2::geom_hline(yintercept = 0, linetype = 'dotted') +
     ggplot2::ggtitle(pathway)

p

结果如下：

可以去试试看~

GSM3387462 PANC1.WT.Rep1
GSM3387463 PANC1.WT.Rep2
GSM3387464 PANC1.WT.Rep3
GSM3387465 PANC1.FOXA2-KO.Rep1
GSM3387466 PANC1.FOXA2-KO.Rep2
GSM3387467 PANC1.FOXA2-KO.Rep3
GSM3387468 CFPAC1.WT.Rep1
GSM3387469 CFPAC1.WT.Rep2
GSM3387470 CFPAC1.WT.Rep3 
GSM3387474 CFPAC1.FOXA2-KD.Rep1
GSM3387475 CFPAC1.FOXA2-KD.Rep2
GSM3387476 CFPAC1.FOXA2-KD.Rep3