27. 基因调控网络#
关键要点
我们用 R 准备了 RNA 和 ATAC 对象，以便用 FigR 和 cisTopic 处理。
我们用 FigR 计算了 DORC 分数，并把它们以散点图、热图和网络的形式可视化。
环境设置
安装 conda:
在创建环境之前，请确保 conda 已安装在您的系统中。
保存 yml 内容:
将 yml 选项卡中的内容复制到名为
environment.yml的文件中。
创建环境:
打开终端或命令提示符。
运行以下命令：
conda env create -f environment.yml
激活环境:
创建好环境后，使用以下命令激活它：
conda activate <environment_name>
替换
<environment_name>，名称就是在environment.yml文件中指定的那个。在 yml 文件里，它看起来像这样：name: <environment_name>
验证安装:
通过运行以下命令，检查环境是否创建成功：
conda env list
name: gene-regulatory-networks-atac
channels:
- conda-forge
- anaconda
- bioconda
dependencies:
- conda-forge::python=3.8 # this package version is required for figR
- conda-forge::r-base=4.2.2
- conda-forge::r-umap
- conda-forge::r-devtools
- conda-forge::r-cairo
- conda-forge::r-ggrastr
- conda-forge::r-patchwork
- conda-forge::r-gplots
- conda-forge::r-networkd3 # network viz
- conda-forge::r-r2d3 # network viz
- conda-forge::r-webshot # network viz
- conda-forge::phantomjs # network viz
- conda-forge::r-imager # network viz
- conda-forge::jupyterlab
- conda-forge::jupyter_server=1.18.1
- conda-forge::r-rmpfr
- conda-forge::r-irkernel
- conda-forge::r-dplyr
- bioconda::bioconductor-zellkonverter>=1.8
# - bioconda::bioconductor-bsgenome.hsapiens.ucsc.hg38=1.4.4 # this package will have to be installed while running
- bioconda::bioconductor-complexheatmap
- bioconda::bioconductor-motifmatchr
- bioconda::bioconductor-rcistarget
- bioconda::bioconductor-chromvar
27.1. 动机#
把染色质可及性和基因表达放在一起分析以理解基因调控，是很有帮助的，因为在基因调控的控制过程中，这两者之间存在机制性的关联，这种关联由转录因子（TF）和其他表观遗传调节因子介导 [Spitz and Furlong, 2012]。简而言之，在基因表达调控的早期阶段，被注释为启动子和局部/远端增强子的调控区域会被调用，而染色质可及性的增加或减少，可以用作其活性变化的替代指标。因此，在一定基因组邻近距离（如小于 200 Kbp）内，近端和远端可及元件（由 ATAC-seq 测量）与靶基因（由 RNA-seq 测量）之间的全局正相关或负相关，可用于在基因调控网络（GRN）推断中注释基因组层面的调控关系。利用描述基因（RNA）和峰（ATAC）特征的测序数据，在峰与基因之间构建相关性矩阵的工具，有助于总结出强的峰-基因相互作用。
27.1.1. 使用 RNA 和 ATAC 特征的基因调控网络推断#
在这个笔记本里，我们使用 FigR [Kartha et al., 2022] 来在 NeurIPS 数据集的一个供体上描述这种 GRN 构建步骤。这个笔记本的准备脚本也会调用 cisTopic [Bravo González-Blas et al., 2019] 生成峰聚类或 topics ，均来自 ATAC-seq 计数矩阵。在计算 RNA-ATAC 相关性时，FigR 使用 ChromVAR [Schep et al., 2017] 把转录因子的基序（motif）映射到比对上的峰。本笔记本中描述的主要处理步骤，是根据 FigR 针对 SHARE-seq 数据的核心教程改编的 教程
声明：在撰写本章时，已经有几种利用单细胞 RNA 和 ATAC 信息进行 GRN 推断的方法。出于与上一章所述相似的缺乏基准测试原因，我们不能说哪种方法在大多数场景下都表现最好。在给定可用数据的情况下，我们总体上建议把其中任何一种作为起点，用来推断初步的 GRN。
27.1.2. 安装并加载 FigR 包#
# the installation of this package is required for the proper execution of this notebook.
if(!suppressMessages(require("FigR"))){
suppressMessages(devtools::install_github("caleblareau/BuenColors")) # the package BuenColors is also a devtools dependency to install FigR.
suppressMessages(devtools::install_github("buenrostrolab/FigR"))
}
suppressMessages(library(FigR))
27.1.3. 使用 zellkonverter 把 h5ad 转换为 SingleCellExperiment#
library(zellkonverter)
library(SingleCellExperiment)
使用 zellkonverter 加载完整的 NeurIPS 数据集。然后，使用以下列取出一个供体（s1d1）的子集 batch。本教程将使用这个供体。
sce <- readH5AD("../../data/openproblems_bmmc_multiome_genes_filtered.h5ad")
sce <- sce[colData(sce)$batch == 's1d1',]
为了加快计算，我们定义了一个特征子集
ncells = -1 # if using all cells
nfeatures_rna = 10000 # -1 if using all rna features
nfeatures_atac = 10000 # -1 if using all atac features
取出供体 s1d1 的子集后，按 RNA 和 ATAC 把特征拆分到两个对象中（如果出现新的索引，可以用以下方式检查 rownames(sce))
RNA <- sce[1:13431,]
ATAC <- sce[13432:nrow(sce),]
下载 hg38 的 TF 注释列表，使用与“仅 RNA”教程相同的基因标签。其他注释，例如 HumanTFs,也可以使用。
if(!file.exists('allTFs_hg38.txt'))
download.file('https://raw.githubusercontent.com/aertslab/SCENICprotocol/master/example/allTFs_hg38.txt', 'allTFs_hg38.txt')
tf_names = rownames(read.table('allTFs_hg38.txt', row.names=1))
print(length(tf_names))
[1] 1797
# subset by cells:
if(ncells != -1){
RNA <- RNA[,1:ncells]
ATAC <- ATAC[,1:ncells]
}
if(nfeatures_rna != -1){
# select all TFs and a subset of non-TFs
is_tf <- rownames(RNA) %in% tf_names
index_tf <- which(is_tf)
index_not_tf <- which(!is_tf)
set.seed(nfeatures_rna)
RNA <- RNA[c(index_tf, sample(index_not_tf, nfeatures_rna - length(index_tf))),]
}
根据变异系数（coefficient of variation）准则选择峰，改编自 EpiScanpy 的实现
if(nfeatures_atac != -1){
frac_atac <- rowSums(assay(ATAC)) / ncol(ATAC)
acc_score <- abs(0.5 - frac_atac)
peak_mask <- rownames(ATAC) %in% names(sort(acc_score))[1:nfeatures_atac]
ATAC <- ATAC[peak_mask,]
}
print(c(dim(ATAC), dim(RNA)))
is_tf <- rownames(RNA) %in% tf_names
[1] 10000 6224 10000 6224
UMAP <- reducedDim(ATAC, 'GEX_X_umap')
# counts variable (for later functions)
assay(ATAC) <- as(assay(ATAC), 'sparseMatrix')
counts(ATAC) <- assay(ATAC)
使用细胞数对 SummarizedExperiment 对象做预处理
# Remove genes with zero expression across all cells
# RNA <- RNA[Matrix::rowSums(RNA) != 0,]
27.1.4. cisTopic 的准备与执行#
cisTopic 应用潜在狄利克雷分配（Latent Dirichlet Allocation，LDA）来定义若干个“主题（topic）”，这些主题能够概括 ATAC-seq 检索到的、所观察到的染色质可及性计数的变异。在这里，我们安装并计算 cisTopic 对象（使用 scATAC-seq），并在下游用作 FigR 的输入。
通过 devtools 安装 cisTopic（如有需要）。由于 cisTopic 有其特定的依赖项，这个安装大约需要十分钟。（可能需要重启内核。）
if(!suppressMessages(require("cisTopic")))
suppressMessages(devtools::install_github("aertslab/cisTopic"))
suppressMessages(library(cisTopic))
下面的脚本需要指定核心数（nCores). 根据现有计算资源增加
nCores = 2
cistopic_bkp_path <- "../../data/openproblems_bmmc_multiome_genes_filtered_atac_s1d1_counts_cisTopic.rds"
if(nfeatures_atac != -1)
cistopic_bkp_path <- paste0("../../data/openproblems_bmmc_multiome_genes_filtered_atac_s1d1_counts_cisTopic_npeaks", nfeatures_atac, ".rds")
print(c(file.exists(cistopic_bkp_path), cistopic_bkp_path))
[1] "TRUE"
[2] "../../data/openproblems_bmmc_multiome_genes_filtered_atac_s1d1_counts_cisTopic_npeaks10000.rds"
通过 cisTopic 进行的这一分析，使我们能够检索并解读两个特定的值，并以折线图呈现：
一条曲线，展示 cisTopic 计算的主对数似然随主题数的变化。这条曲线出现平台，表明相比之前的主题，再增加主题已不是解释数据的主要因素。
似然对所声明主题数的一阶导数。这种可视化有助于评估似然值的收敛，以及随着主题数增加、在零附近的随机振荡。
overwrite = TRUE # to visualize plots, this will always execute. Modify to FALSE to avoid replacing output from previous runs.
if(!file.exists(cistopic_bkp_path) || overwrite){
# the number of topics to test
n_topics = 1:25
atac <- as.matrix(counts(ATAC))
atac <- as.data.frame(atac)
# we need to work out the names of the rownames, and replace - into : to match the chromosome:start-end required format by cisTopic
chr <- sapply(strsplit(rownames(atac),"-"), `[`, 1)
start <- sapply(strsplit(rownames(atac),"-"), `[`, 2)
end <- sapply(strsplit(rownames(atac),"-"), `[`, 3)
rownames(atac) <- paste0(chr, ':', start, '-', end)
cisTopicObject <- createcisTopicObject(atac, project.name='neurips_s1d1')
cisTopicObject <- runCGSModels(cisTopicObject, topic=c(n_topics), # , 5:15, 20, 25), # topic=c(2, 5:15, 20, 25),
seed=987, nCores=nCores, burnin = 90,
iterations = 100, addModels=FALSE)
cisTopicObject <- selectModel(cisTopicObject, type='maximum')
# cisTopicObject
cisTopicObject <- runUmap(cisTopicObject, target='cell')
topic.mat <- modelMatSelection(cisTopicObject, 'cell', 'Probability')
topic.mat <- t(topic.mat)
topic.mat <- as.matrix(topic.mat)
saveRDS(topic.mat, cistopic_bkp_path)
}
[1] "Formatting data..."
[1] "Exporting data..."
[1] "Running models..."
Loading required package: umap
cisTopic 对象生成后，我们可以用它来提取相关特征
cisAssign <- readRDS(cistopic_bkp_path)
dim(cisAssign) # Cells x Topics
- 6224
- 25
使用 cisTopic 的主题矩阵计算 kNN 图
library(dplyr)
library(FNN)
set.seed(123)
cellkNN <- get.knn(cisAssign, k = 30)$nn.index
dim(cellkNN)
- 6224
- 30
使用 UMAP 生成的嵌入可视化细胞
colData(ATAC)$cellAnnot <- colData(ATAC)$cell_type
colData(ATAC)$UMAP1 <- UMAP[,1]
colData(ATAC)$UMAP2 <- UMAP[,2]
# Plot
library(ggplot2)
options(repr.plot.width=12, repr.plot.height=6)
colData(ATAC) %>% as.data.frame() %>% ggplot(aes(UMAP1,UMAP2,color=cellAnnot)) +
geom_point(size=0.5) + # scale_color_manual(values=annoCols)+
theme_classic() + guides(colour = guide_legend(override.aes = list(size=2))) + cowplot::theme_cowplot(font_family='Sans')
输入数据的准备工作已经完成。现在，开始执行 FigR 算法。
# if the hg38 genome is not installed successfully during environment building, it can be installed here.
# this might require a kernel restart
if(!suppressMessages(require('BSgenome.Hsapiens.UCSC.hg38'))){
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("BSgenome.Hsapiens.UCSC.hg38")
}
suppressMessages(library(BSgenome.Hsapiens.UCSC.hg38))
# check object dimensions
c(dim(ATAC), dim(RNA))
- 10000
- 6224
- 10000
- 6224
下面的脚本需要指定核心数（nCores). 根据现有计算资源增加
为 runGenePeakcorr 准备数据
library(Matrix)
RNAmat <- as.matrix(assay(RNA))
dim(RNAmat)
- 10000
- 6224
ATAC_df <- as.data.frame(as.matrix(counts(ATAC)))
# df <- as.data.frame(as.matrix(counts(ATAC)))
ATAC_df$seqnames <- sapply(strsplit(rownames(ATAC_df),"-"), `[`, 1)
ATAC_df$start <- sapply(strsplit(rownames(ATAC_df),"-"), `[`, 2)
ATAC_df$end <- sapply(strsplit(rownames(ATAC_df),"-"), `[`, 3)
ATAC_df <- subset(ATAC_df, grepl('chr', rownames(ATAC_df)))
dim(ATAC_df)
- 9997
- 6227
ATAC.se <- makeSummarizedExperimentFromDataFrame(ATAC_df)
counts(ATAC.se) <- assay(ATAC.se)
assay(ATAC.se) <- as(assay(ATAC.se), 'sparseMatrix')
用 2 个核心运行 FigR 的 runGenePeakcorr 函数
# This snippet can be run interactively, but it takes a long time.
nCores = 2
bkp_path_ciscorr <- '../../data/openproblems_bmmc_multiome_genes_filtered_s1d1_ciscorr.rds'
if(nfeatures_atac != -1)
cistopic_bkp_path <- paste0("../../data/openproblems_bmmc_multiome_genes_filtered_s1d1_ciscorr_npeaks", nfeatures_atac, ".rds")
if(!file.exists(bkp_path_ciscorr)){
cisCorr <- FigR::runGenePeakcorr(ATAC.se = ATAC.se,
RNAmat = RNAmat,
genome = "hg38", # One of hg19, mm10 or hg38
nCores = nCores,
p.cut = NULL, # Set this to NULL and we can filter later
n_bg = 250)
saveRDS(cisCorr, bkp_path_ciscorr)
}
cisCorr <- readRDS(bkp_path_ciscorr)
按 p 值过滤相关的峰-基因相关性
cisCorr.filt <- cisCorr %>% filter(pvalZ <= 0.05)
print(c('all associations', nrow(cisCorr)))
print(c('filtered associations', nrow(cisCorr.filt)))
[1] "all associations" "3914"
[1] "filtered associations" "674"
if(nrow(cisCorr.filt) == 0)
print('increase the number of cells/peaks/genes to discover more associations')
stopifnot(nrow(cisCorr.filt) > 0)
一旦计算出峰与基因的相关性，FigR 就会把 ATAC-seq 峰归入“调控染色质域”（Domains of Regulatory Chromatin，DORC）。这些分组有助于描述 TF 的 RNA 表达水平，与某个基因周围多个染色质可及元件的整体变化之间的关系。重要的是，如果某个 DORC 中的多个染色质可及元件也包含与所关注 TF 相关的 DNA 结合基序，那么这两条证据线对于推断 TF 激活因子（TF 基序富集 + 正的 DORC-TF RNA 相关）或 TF 抑制因子（TF 基序富集 + 负的 DORC-TF RNA 相关）都会很有用。
下面的可视化按“与某 TF 强相关的峰的数量”对 TF 进行排序，并提供哪些 TF 与 DORC 关系最密切的信息。
library(ggrepel)
# Determine DORC genes
options(repr.plot.width=12, repr.plot.height=6)
cowplot_mono <- cowplot::theme_cowplot(font_family='sans')
dorcGenes <- cisCorr.filt %>% dorcJPlot(cutoff=2, # Default
returnGeneList = TRUE, family='sans') # + cowplot_mono
theme_set(theme_gray(base_family = "sans"))
dorcGenes + theme(font='roboto') # font problem during visualization of y/x axes (Roboto)
NULL
在此修改截断值，以恢复至少 30 个基因。否则，runFigRGRN 将无法执行。
stopifnot(length(dorcGenes) > 30)
获取 DORC 分数
dorcMat <- getDORCScores(ATAC.se, dorcTab=cisCorr.filt, geneList=dorcGenes, nCores=nCores)
# Smooth DORC scores (using cell KNNs)
Running DORC scoring for 47 genes: MGAT4A
ADGRG1
APBB1
BEST1
C15orf39
CASP8
CCL5
CCR7
CD300LF
CD3E
CD86
COLQ
CR1L
CTSB
DNAJA1
DNPH1
EPN2
FPR1
GATA3
GM2A
KLHL36
LILRB2
LONRF1
MYADM
NCF2
NFIA
NOLC1
P2RX1
PRF1
PRKCQ-AS1
PTCH1
RFPL1S
RHOQ
RXRA
SESN3
SLC16A6
SLC2A9
SMAD7
STRN3
SWAP70
TEC
THEMIS
TLE1
TM9SF2
TNRC6B
ZNF471
ZNRF1
........
Normalizing scATAC counts ..
SummarizedExperiment object input detected .. Centering counts under assayCentering counts for cells sequentially in groups of size 5000 ..
Computing centered counts for cells: 1 to 5000 ..
Computing centered counts per cell using mean reads in features ..
Computing centered counts for cells: 5001 to 6224 ..
Computing centered counts per cell using mean reads in features ..
Merging results..
Done!
Computing DORC scores ..
Running in parallel using 2 cores ..
Time Elapsed: 0.89201807975769 secs
stopifnot(nrow(cellkNN) == ncol(dorcMat))
rownames(cellkNN) <- colnames(dorcMat)
要用多个核心执行 smoothScoresNN 函数，需要 doParallel。
library(doParallel)
# Smooth dorc scores using cell KNNs (k=30)
dorcMat.s <- smoothScoresNN(NNmat = cellkNN[,1:30], mat = dorcMat, nCores = nCores)
Number of cells in supplied matrix: 6224
Number of genes in supplied matrix: 47
Number of nearest neighbors being used per cell for smoothing: 30
| | 0%Running in parallel using 2 cores ..
|======================================================================| 100%
Merging results ..
Time Elapsed: 14.4018683433533 secs
stopifnot(nrow(cellkNN) == ncol(RNAmat))
rownames(cellkNN) <- colnames(RNAmat)
# Smooth RNA using cell KNNs
# This takes longer since it's all genes
RNAmat.s <- smoothScoresNN(NNmat = cellkNN[,1:30], mat = RNAmat, nCores = nCores)
Number of cells in supplied matrix: 6224
Number of genes in supplied matrix: 10000
Number of nearest neighbors being used per cell for smoothing: 30
| | 0%Running in parallel using 2 cores ..
|======================================================================| 100%
Merging results ..
Time Elapsed: 1.23741063674291 mins
library(ggplot2)
library(ggrastr)
# Visualize on pre-computed UMAP
umap.d <- as.data.frame(colData(ATAC)[,c("UMAP1","UMAP2")])
DORC 分数计算完成后，我们可以探索 TF 与靶基因之间的关联，或各细胞类型中 TF 的整体表达。这两个指标有助于理解各细胞类型的调控。
print(length(dorcGenes))
print(dorcGenes)
[1] 47
[1] "MGAT4A" "ADGRG1" "APBB1" "BEST1" "C15orf39" "CASP8"
[7] "CCL5" "CCR7" "CD300LF" "CD3E" "CD86" "COLQ"
[13] "CR1L" "CTSB" "DNAJA1" "DNPH1" "EPN2" "FPR1"
[19] "GATA3" "GM2A" "KLHL36" "LILRB2" "LONRF1" "MYADM"
[25] "NCF2" "NFIA" "NOLC1" "P2RX1" "PRF1" "PRKCQ-AS1"
[31] "PTCH1" "RFPL1S" "RHOQ" "RXRA" "SESN3" "SLC16A6"
[37] "SLC2A9" "SMAD7" "STRN3" "SWAP70" "TEC" "THEMIS"
[43] "TLE1" "TM9SF2" "TNRC6B" "ZNF471" "ZNRF1"
与 TF 表达强相关的基因之一是 核因子 IA（NFIA）。我们可以进一步考察这个案例：检查它的基因表达水平，以及与其基序在 ATAC 峰上相关联的 dorcGene。
marker_gene = 'NFIA'
dorcg <- plotMarker2D(umap.d,dorcMat.s,markers = c(marker_gene),maxCutoff = "q0.99",
colorPalette = "brewer_heat") + ggtitle(paste0(marker_gene, ' DORC'))
Plotting NFIA
rnag <- plotMarker2D(umap.d,RNAmat.s,markers = c(marker_gene),maxCutoff = "q0.99",
colorPalette = "brewer_purple") + ggtitle(paste0(marker_gene, ' RNA'))
Plotting NFIA
在这里，我们用 patchwork 可视化 dorcg 和 rnag 对象。这样可以逐细胞地对基因表达和相关的 DORC 分数进行直观比较。
options(repr.plot.width=12, repr.plot.height=6)
library(patchwork)
dorcg + cowplot_mono + rnag + cowplot_mono
通过目视比较，我们可以发现 NFIA 表达与 DORC 分数之间的吻合并非一一对应，也就是说，在某些细胞类型聚类中 NFIA 强烈表达，但这些表达水平与“假定被 NFIA 占据的 ATAC 峰”之间却没有强相关。
dim(dorcMat.s)
- 47
- 6224
figR.d <- runFigRGRN(ATAC.se = ATAC.se, # Must be the same input as used in runGenePeakcorr()
dorcTab = cisCorr.filt, # Filtered peak-gene associations
genome = "hg38",
dorcMat = dorcMat.s,
rnaMat = RNAmat.s,
nCores = nCores)
Assuming peak indices in Peak field
Removing genes with 0 expression across cells ..
Getting peak x motif matches ..
Determining background peaks ..
Using 50 iterations ..
Testing 626 TFs
Testing 47 DORCs
Running FigR using 2 cores ..
|======================================================================| 100%Finished!
Merging results ..
27.2. 结果可视化#
TF-DORC 调控分数（散点图）。y 轴表示某个 TF 在 DORC 中的基序富集，而 x 轴表示该 TF 表达与这些峰之间的相关性（相对于背景的 Z 检验）。直观地说，这种可视化让我们能够找出可能的 TF 激活因子（右上）和 TF 抑制因子（左上）。
require(ggplot2)
require(ggrastr)
require(BuenColors) # https://github.com/caleblareau/BuenColors
options(repr.plot.width=10, repr.plot.height=8)
figR.d %>%
ggplot(aes(Corr.log10P,Enrichment.log10P,color=Score)) +
ggrastr::geom_point_rast(size=0.01,shape=16) +
theme_classic() +
scale_color_gradientn(colours = jdb_palette("solar_extra"),limits=c(-3,3),oob = scales::squish,breaks=scales::breaks_pretty(n=3)) +
cowplot_mono + cowplot::theme_cowplot(font_family='sans' )
基于排名的驱动 TF 可视化，以获得一份直观的、正负向排列的相关候选 TF 列表。
drivers <- rankDrivers(figR.d, rankBy = "meanScore",interactive = FALSE)
options(repr.plot.width=15, repr.plot.height=8)
drivers + geom_text_repel(family='roboto') + geom_label_repel(family='roboto') + cowplot::theme_cowplot(font_family='roboto') + labs(font_family='roboto')
Ranking TFs by mean regulation score across all DORCs ..
Warning message:
“ggrepel: 39 unlabeled data points (too many overlaps). Consider increasing max.overlaps”
Warning message:
“ggrepel: 51 unlabeled data points (too many overlaps). Consider increasing max.overlaps”
Warning message:
“ggrepel: 55 unlabeled data points (too many overlaps). Consider increasing max.overlaps”
接下来的可视化，尝试按每个相关基序所关联的“被激活或被抑制的靶基因数量”来展示这些被注释的基序。这样可以观察到顶部基序在所研究样本中可能激活或抑制某些基因程序的覆盖范围。
options(repr.plot.width=10, repr.plot.height=10)
rankDrivers(figR.d,score.cut = 1, rankBy = "nTargets", interactive = FALSE, fontsiz)
Ranking TFs by total number of associated DORCs ..
Using absolute score cut-off of: 1 ..
下面的可视化是一张基于热图的可视化，它基于 TF-DORC-靶基因关联，展示候选 TF 及其可能调控的强相关基因的 DORC 分数。（行 = 靶基因，列 = TF）。
library(grid)
library(ComplexHeatmap)
options(repr.plot.width=10, repr.plot.height=10)
pushViewport(viewport(gp = gpar(fontfamily = "sans")));
heatmap <- plotfigRHeatmap(figR.d = figR.d,
score.cut = 1,
TFs = unique(figR.d$Motif),
column_names_gp = gpar(fontsize=10), # from ComplexHeatmap
show_row_dend = FALSE # from ComplexHeatmap
)
draw(heatmap, newpage = FALSE)
popViewport()
Using absolute score cut-off of: 1 ..
Using Score as value column: use value.var to override.
Plotting 45 DORCs x 46TFs
这种可视化使我们能够快速识别主要的 TF，以及它们如何通过其基因组邻域中的 DORC 来调控潜在的靶基因。
作为最后一个探索性可视化，还可以使用 networkD3 包检索出一个关联网络。这样可以探索“跨多个 TF 相互重叠的基因集”，以及可能由共同 TF 控制的、相互重叠的调节子。
library(networkD3)
library(r2d3)
library(imager)
# generate the network
d3 <- plotfigRNetwork(figR.d,
score.cut = 1,
TFs = unique(figR.d$Motif),
weight.edges = TRUE)
# in a local session, the network can be manipulated interactively
d3
27.2.1. 网络另存为 HTML/ PNG#
# If pandoc if not found by R, here the bin path in the environment has to be provided e.g. `envs/best_practices_regulons_rnanatac/bin`
# rmarkdown::find_pandoc(dir = "bin_path") # e.g. bin_path = '~/miniconda3/envs/best_practices_regulons_rnanatac/bin'
# is saved as an image and shown for exploratory purposes outside of this notebook
save_d3_html(
d3,
'network_tutorial_rna_n_atac.html',
)
# the network is saved as an image and shown for online purposes.
save_d3_png(
d3,
'network_tutorial_rna_n_atac.png',
width = 350,
height = 450,
delay = 4.0,
zoom = 1.6,
)
im <- load.image('network_tutorial_rna_atac.png')
plot(im)
用于执行此笔记本的软件包日志
sessionInfo()
R version 4.2.2 (2022-10-31)
Platform: x86_64-conda-linux-gnu (64-bit)
Running under: Ubuntu 20.04 LTS
Matrix products: default
BLAS/LAPACK: /home/rio/miniconda3/envs/best_practices_regulons_rnanatac/lib/libopenblasp-r0.3.21.so
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
[4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
[7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
[10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] grid parallel stats4 stats graphics grDevices utils
[8] datasets methods base
other attached packages:
[1] imager_0.42.19 magrittr_2.0.3
[3] r2d3_0.2.6 networkD3_0.4
[5] ComplexHeatmap_2.14.0 BuenColors_0.5.6
[7] MASS_7.3-58.3 patchwork_1.1.2
[9] ggrastr_1.0.1 doParallel_1.0.17
[11] iterators_1.0.14 foreach_1.5.2
[13] ggrepel_0.9.3 BSgenome.Hsapiens.UCSC.hg38_1.4.5
[15] BSgenome_1.66.3 rtracklayer_1.58.0
[17] Biostrings_2.66.0 XVector_0.38.0
[19] FNN_1.1.3.2 umap_0.2.10.0
[21] cisTopic_0.3.0 SingleCellExperiment_1.20.0
[23] zellkonverter_1.8.0 FigR_0.1.0
[25] motifmatchr_1.20.0 chromVAR_1.20.0
[27] cowplot_1.1.1 ggplot2_3.4.1
[29] dplyr_1.1.1 SummarizedExperiment_1.28.0
[31] Biobase_2.58.0 GenomicRanges_1.50.0
[33] GenomeInfoDb_1.34.9 IRanges_2.32.0
[35] S4Vectors_0.36.0 BiocGenerics_0.44.0
[37] MatrixGenerics_1.10.0 matrixStats_0.63.0
[39] Matrix_1.5-3
loaded via a namespace (and not attached):
[1] utf8_1.2.3 reticulate_1.26
[3] R.utils_2.12.2 tidyselect_1.2.0
[5] poweRlaw_0.70.6 RSQLite_2.3.0
[7] AnnotationDbi_1.60.2 htmlwidgets_1.6.2
[9] gmp_0.7-1 BiocParallel_1.32.5
[11] munsell_0.5.0 codetools_0.2-19
[13] DT_0.27 pbdZMQ_0.3-9
[15] miniUI_0.1.1.1 withr_2.5.0
[17] RcisTarget_1.18.2 colorspace_2.1-0
[19] filelock_1.0.2 knitr_1.42
[21] uuid_1.1-0 rstudioapi_0.14
[23] pbmcapply_1.5.1 labeling_0.4.2
[25] repr_1.1.6 GenomeInfoDbData_1.2.9
[27] lgr_0.4.4 bit64_4.0.5
[29] farver_2.1.1 rprojroot_2.0.3
[31] basilisk_1.9.12 vctrs_0.6.1
[33] generics_0.1.3 xfun_0.38
[35] float_0.3-1 R6_2.5.1
[37] clue_0.3-64 ggbeeswarm_0.7.1
[39] bitops_1.0-7 arrow_10.0.1
[41] cachem_1.0.7 DelayedArray_0.24.0
[43] assertthat_0.2.1 promises_1.2.0.1
[45] BiocIO_1.8.0 scales_1.2.1
[47] beeswarm_0.4.0 gtable_0.3.3
[49] Cairo_1.6-0 processx_3.8.0
[51] bmp_0.3 seqLogo_1.64.0
[53] rlang_1.1.0 GlobalOptions_0.1.2
[55] text2vec_0.6.3 splines_4.2.2
[57] lazyeval_0.2.2 yaml_2.3.7
[59] reshape2_1.4.4 httpuv_1.6.9
[61] tools_4.2.2 feather_0.3.5
[63] nabor_0.5.0 gplots_3.1.3
[65] ellipsis_0.3.2 RColorBrewer_1.1-3
[67] Rcpp_1.0.10 plyr_1.8.8
[69] base64enc_0.1-3 sparseMatrixStats_1.10.0
[71] zlibbioc_1.44.0 purrr_1.0.1
[73] RCurl_1.98-1.10 ps_1.7.3
[75] basilisk.utils_1.10.0 openssl_2.0.5
[77] GetoptLong_1.0.5 cluster_2.1.4
[79] here_1.0.1 data.table_1.14.8
[81] RSpectra_0.16-1 readbitmap_0.1.5
[83] circlize_0.4.15 mlapi_0.1.1
[85] fitdistrplus_1.1-8 hms_1.1.3
[87] mime_0.12 evaluate_0.20
[89] xtable_1.8-4 RhpcBLASctl_0.23-42
[91] XML_3.99-0.14 jpeg_0.1-10
[93] AUCell_1.20.2 shape_1.4.6
[95] compiler_4.2.2 tibble_3.2.1
[97] KernSmooth_2.23-20 crayon_1.5.2
[99] R.oo_1.25.0 htmltools_0.5.5
[101] tiff_0.1-11 later_1.3.0
[103] tzdb_0.3.0 TFBSTools_1.36.0
[105] snow_0.4-4 tidyr_1.3.0
[107] DBI_1.1.3 readr_2.1.4
[109] cli_3.6.1 R.methodsS3_1.8.2
[111] igraph_1.4.1 pkgconfig_2.0.3
[113] GenomicAlignments_1.34.0 rsparse_0.5.1
[115] dir.expiry_1.6.0 TFMPvalue_0.0.9
[117] IRdisplay_1.1 plotly_4.10.1
[119] annotate_1.76.0 lda_1.4.2
[121] vipor_0.4.5 DirichletMultinomial_1.40.0
[123] webshot_0.5.4 callr_3.7.3
[125] stringr_1.5.0 digest_0.6.31
[127] pracma_2.4.2 CNEr_1.34.0
[129] graph_1.76.0 rmarkdown_2.21
[131] DelayedMatrixStats_1.20.0 GSEABase_1.60.0
[133] restfulr_0.0.15 shiny_1.7.4
[135] Rsamtools_2.14.0 gtools_3.9.4
[137] rjson_0.2.21 lifecycle_1.0.3
[139] jsonlite_1.8.4 viridisLite_0.4.1
[141] askpass_1.1 fansi_1.0.4
[143] pillar_1.9.0 lattice_0.20-45
[145] KEGGREST_1.38.0 fastmap_1.1.1
[147] httr_1.4.5 survival_3.5-5
[149] GO.db_3.16.0 glue_1.6.2
[151] png_0.1-8 bit_4.0.5
[153] stringi_1.7.12 blob_1.2.4
[155] doSNOW_1.0.20 caTools_1.18.2
[157] memoise_2.0.1 Rmpfr_0.9-1
[159] IRkernel_1.3.2
27.3. Quiz#
27.3.1. 理论#
为什么峰-基因关联被认为在机制上是成立的？
峰-基因关联之所以被认为在机制上成立，是因为它们将可及的染色质区域（峰）与其靶基因联系起来，反映了调控元件影响基因表达的物理相互作用。
在处理 ATAC 和 RNA 数据时，染色质峰通常是否比基因更多？在构建 GRN 时，这在结构层面会带来什么后果？
一般来说，ATAC-seq 数据中的染色质峰多于基因，这导致基因调控网络（GRN）中存在多对一的关系，即多个调控元件可以控制单个基因。
在染色质可及性层面和基因调控层面，什么样的转录因子分别被视为激活因子／抑制因子？
在染色质可及性层面，转录因子（TF）激活因子通过促进开放的染色质状态来提高可及性，而抑制因子通过诱导关闭状态来降低可及性；在基因调控方面，激活因子增强转录，而抑制因子抑制转录。
在解读 ATAC+RNA 的 GRN 模型时，哪些额外的读出数据可与 scRNA-seq 和 scATAC-seq 互为补充？
与 scRNA-seq 和 scATAC-seq 互补的读出数据包括：用于刻画蛋白质-DNA 相互作用的 ChIP-seq，以及用于捕获染色质构象的 Hi-C，它们提供了额外层次的调控信息。
27.3.2. FigR#
什么是 DORC 分数？它如何有助于识别峰与基因之间的调控相互作用？
DORC（调控染色质域，Domain of Regulatory Chromatin）分数量化了与单个基因相关的多个染色质峰的总体可及性，通过将这些可及性域与基因表达水平相关联，有助于识别调控相互作用。
27.4. 参考文献#
[atacBGonzalezBMP+19]
Carmen Bravo González-Blas, Liesbeth Minnoye, Dafni Papasokrati, Sara Aibar, Gert Hulselmans, Valerie Christiaens, Kristofer Davie, Jasper Wouters, and Stein Aerts. Cistopic: cis-regulatory topic modeling on single-cell ATAC-seq data. Nat. Methods, 16(5):397–400, May 2019.
[atacKDH+22]
Vinay K Kartha, Fabiana M Duarte, Yan Hu, Sai Ma, Jennifer G Chew, Caleb A Lareau, Andrew Earl, Zach D Burkett, Andrew S Kohlway, Ronald Lebofsky, and Jason D Buenrostro. Functional inference of gene regulation using single-cell multi-omics. Cell Genom, September 2022.
[atacSWBG17]
Alicia N Schep, Beijing Wu, Jason D Buenrostro, and William J Greenleaf. chromVAR: inferring transcription-factor-associated accessibility from single-cell epigenomic data. Nat. Methods, 14(10):975–978, October 2017.
[atacSF12]
François Spitz and Eileen E M Furlong. Transcription factors: from enhancer binding to developmental control. Nat. Rev. Genet., 13(9):613–626, 2012.
27.5. 贡献者#
我们衷心感谢以下人员的贡献：
27.5.2. 审阅者#
Lukas Heumos
Anna Schaar