DoRothEA|CollecTRI—转录因子活性评分

DoRothEA|CollecTRI—转录因子活性评分 by 生信小博士

写在开头

jimmy老师在他的生信游民系列提出了很多学徒任务，仔细去看每一篇文章的动机，都是科研过程中经常会经历的。

今天我们分享的任务是来自：生物学功能注释三板斧。本次任务我计划分为四个部分来写，今天写第四部分。

正文部分

通常我们都是利用基因集合去做单细胞评分，那么是否可以通过有权重的基因集来计算细胞的转录活性？这就是DoRothEA开发出来的意义。

CollecTRI数据集简介

DoRothEA是一种包含转录因子（TF）与其靶标相互作用的基因集合。这些相互作用形成了被定义为调节子（regulons）的TF及其对应靶点的集合。

CollecTRI 来源的调控子总结了从 12 种不同转录因子（TF）-目标基因相互作用的数据库。该集合扩大了转录因子的覆盖范围，并与其他已知的 GRN 进行了比较。结果表明，在使用 knockTF 数据集根据基因表达数据识别受干扰的 TF 方面，该集合具有卓越的性能。

简单来说，CollecTRI 是第二代DoRothEA，对DoRothEA做了更新。

加载所需数据

本教程中，我就使用示例数据了。如果你使用自己的数据，把数据的格式调整为一样即可！


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

# Remove NAs and set row namescounts <- data$counts %>%  dplyr::mutate_if(~ any(is.na(.x)), ~ if_else(is.na(.x),0,.x)) %>%   column_to_rownames(var = "gene") %>%   as.matrix()head(counts);dim(counts)
design <- data$designdesign

后续计算评分所需要的表达矩阵，行为基因名，列为样本名称

加载权重基因集


net <- get_collectri(organism='human', split_complexes=FALSE)net#如果你是小鼠数据，改为mouse即可，参考之前的推文

如果你是小鼠数据，改为mouse即可，参考之前的推文：DoRothEA--单细胞转录因子评分

通路活性推断

为了推断通路活性，这里将运行加权平均法（wmean）。它通过首先将每个转录因子的下游分子乘以其相关的权重来推断该转录因子的通路活动，然后将它们求和得到一个富集分数 wmean。此外，可以对随机目标特征进行排列组合，以获得一个用于计算 z 分数 norm_wmean 的概率分布，或者通过将 wmean 乘以所得经验性 p 值的负对数来得到一个校正估计 corr_wmean。

在这个例子中，我们这里使用了 wmean，但我们也可以使用其他任何方法。要查看可用的方法，可以使用 show_methods()。

代码非常简单：

# Run wmeansample_acts <- run_wmean(mat=counts, net=net, .source='source', .target='target',                         .mor='mor', times = 100, minsize = 5)sample_acts

这样我们就得到了每个样本的通路活性，这里的通路活性是以长数据的形式提供的。

而我们之前分享的内容得到的通路都是矩阵的形式：

如果你你想把长数据改为我们常见的宽数据格式，也就是通路矩阵，这里直接给出代码：

 # Transform to wide matrix  sample_acts_mat <- sample_acts %>%    filter(statistic == 'norm_wmean') %>%    pivot_wider(id_cols = 'condition', names_from = 'source',                values_from = 'score') %>%    column_to_rownames('condition') %>%    as.matrix();sample_acts_mat[1:6,1:6]

这里要注意通路矩阵的行、列各代表什么

得到这个通路矩阵之后，可视化又是各显神通啦！

如果你想了解更多长宽数据：

可视化一

这里我们选择变化较大的前25各tf进行可视化

 n_tfs <- 25    # Get top tfs with more variable means across clusters  tfs <- sample_acts %>%    group_by(source) %>%    summarise(std = sd(score)) %>%    arrange(-abs(std)) %>%    head(n_tfs) %>%    pull(source)  sample_acts_mat <- sample_acts_mat[,tfs] ; sample_acts_mat[,1:9]

可视化代码

 # Scale per sample  sample_acts_mat <- scale(sample_acts_mat)    # 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(sample_acts_mat, border_color = NA, color=my_color, breaks = my_breaks)

可视化二

所有的可视化都离不开前面得到的这个矩阵：

我们直接使用var来获取变异度最大的前50个tf进行可视化

  # Transform to wide matrix  sample_acts_mat <- sample_acts %>%    filter(statistic == 'norm_wmean') %>%    pivot_wider(id_cols = 'condition', names_from = 'source',                values_from = 'score') %>%    column_to_rownames('condition') %>%    as.matrix();sample_acts_mat[1:6,1:6]    
  # Transpose the matrix  transposed_mat <- t(sample_acts_mat)    # Calculate the variance of each gene  gene_variance <- apply(transposed_mat, 1, var)    # Sort the genes by variance and select the top 50  top_50_genes <- names(sort(gene_variance, decreasing = TRUE)[1:50])    # Extract the top 50 genes from the transposed matrix  top_50_genes_mat <- transposed_mat[top_50_genes, ]    # Now, install and load the pheatmap package if you haven't already  # install.packages("pheatmap")  library(pheatmap)    pheatmap:: pheatmap(top_50_genes_mat,scale = 'row'   )

加上组别信息

    design=as.data.frame(design)   rownames(design)=design$sample ;design   # Visualize the top 50 genes using pheatmappheatmap:: pheatmap(top_50_genes_mat,scale = 'row',cluster_cols = FALSE,show_colnames = FALSE,                    annotation_col = design                      )

可视化三

请思考：为什么上述的这些可视化结果不太好看呢？

原因就在于我们用的是变化最大的TF。如果我们使用差异TF来画热图就会很好看，这个问题就留给读者尝试吧！

参考：https://github.com/saezlab/CollecTRIhttps://github.com/saezlab/CollecTRI/blob/main/scripts/casestudy/02.2_pbmc_TFactivity.Rhttps://bioconductor.org/packages/release/bioc/vignettes/decoupleR/inst/doc/tf_bk.html

ixxmu / mp_duty