6.1 Gene regulatory networks¶

动机¶

一旦处理了单细胞基因组数据，就可以在其基因组上下文中解析观察到的特征之间的重要关系。在我们的基因组中，基因的激活由细胞核中的RNA转录机制控制，该机制激活局部（启动子）或远端的顺式调控元件（增强子），以控制每个基因产生的RNA量。

概念上，基因调控网络（GRN）是指图形表示某些控制转录的基因，即“转录因子”（TF），直接控制其目标基因的转录率（顺式调控）。同时，这些目标基因一旦被激活，可以负责控制其他下游目标基因（反式调控）。在计算上，推断GRN的方法考虑基因和染色质可及性特征的协同变化，以识别可以同时与少数TF相关联的模块。由同一TF的活动控制的一组基因被定义为一个调控子。

除了协同变化外，几种方法还认识到收集和注入先验知识数据的重要性，例如TF在基因组中结合的位置或先前报告的TF-目标基因关联，以预定义基因-基因边缘，从而指导支持这种证据类型的GRN推断。迄今为止，GRN的基准测试不同于图像推断等其他机器学习任务，因为标记数据稀少且难以验证。出于这个原因，几种方法开发了自己的基准，主要使用真实数据将新方法与其他方法进行比较。GRN推断方法的泛化是调控基因组学中的一个活跃讨论话题，与前几章不同，需要一个目前难以达成一致并被社区一致使用的金标准。

本章的目的是展示如何生成GRN的通用流程，使用我们已知的方式呈现的方法，这些方法允许以较少的软件依赖性进行验证。由于GRN工具的一般基准测试限制，我们推荐这些工具，但不能说它们在所有可能的情况下表现最好。我们更愿意推荐这些工具作为起点，鉴于可用数据，然后以最低的计算成本探索它们生成的GRN表示。

从公共数据中收集TF调控子¶

TF-调控子已在学术研究、编译这些数据的数据库以及如ENCODE的联合项目中进行了注释。作为一般建议，可以检查TTRUST {cite}grn:Han2015-yb、DoRothEA {cite}grn:Garcia-Alonso2019-ly、KnockTF {cite}grn:Feng2020-nl等数据库以获取真核生物的TF-调控子。对于原核生物，RegulonDB {cite}grn:Santos-Zavaleta2018-be是一个公认的数据库。

TF调控子的局限性¶

TF-调控子的来源和可信度存在一些限制，例如数据来源和实验读数。如果调控子的数据来源与感兴趣的细胞类型不匹配，则结果无法放入特定生物系统的上下文中。或者，如果实验读数没有测量由感兴趣的TF直接引起的顺式调控事件，我们的TF-调控子可能包含反式调控事件。强烈建议使用在等效生物系统中收集的TF调控子。然而，如果这些不可用，由于强制使用先验知识，结果的解释可能会有偏差。

使用RNA数据生成GRN¶

我们将使用工具SCENIC探索scRNA-seq数据和预测的TF调控子。具体来说，我们将在NeurIPS 2021数据集的一个供体上执行SCENIC {cite}grn:Aibar2017-tp并解释结果。本笔记本中描述的主要处理步骤是基于SCENIC核心教程教程改编的。

分析结果¶

本笔记本允许通过推断基因调控网络来解释单细胞RNA-seq数据，以获得在细胞分化、转变和扰动等背景下解释基因表达的潜在转录因子和目标基因之间的关联。

Dataset description¶

Cell-type clusters of 100,000 human PBMCs and the NeurIPS dataset, which contain healthy donors as well as COVID-19 patients{cite}grn:Schulte-Schrepping2020.

Environment setup¶

In [1]:

import warnings

warnings.filterwarnings("ignore")

import pyscenic
import loompy as lp
import scanpy as sc
import numpy as np
import seaborn as sns
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib import rcParams
import os

import warnings warnings.filterwarnings("ignore") import pyscenic import loompy as lp import scanpy as sc import numpy as np import seaborn as sns import pandas as pd import matplotlib.pyplot as plt from matplotlib import rcParams import os

NeurIPs 数据集的准备¶

我们将要查看的当前数据集已经过预处理（见前几章），包含69,249个细胞，并注释为22种细胞类型。由于批量整合的考虑，我们将仅展示GRN方法在一个标签为s1d1（n=6,224个细胞）上的使用。

Load the full dataset

In [2]:

adata = sc.read_h5ad("../../data/openproblems_bmmc_multiome_genes_filtered.h5ad")
adata.shape

adata = sc.read_h5ad("../../data/openproblems_bmmc_multiome_genes_filtered.h5ad") adata.shape

Out[2]:

(69249, 129921)

Subset RNA only features using the label GEX

In [3]:

rna = adata[:, adata.var.feature_types == "GEX"]
del adata
sc.pp.log1p(rna)
rna.obs.batch.value_counts()

rna = adata[:, adata.var.feature_types == "GEX"] del adata sc.pp.log1p(rna) rna.obs.batch.value_counts()

Out[3]:

s4d8     9876
s4d1     8023
s3d10    6781
s1d2     6740
s1d1     6224
s2d4     6111
s2d5     4895
s3d3     4325
s4d9     4325
s1d3     4279
s2d1     4220
s3d7     1771
s3d6     1679
Name: batch, dtype: int64

In [4]:

rna.shape

rna.shape

Out[4]:

(69249, 13431)

在这里，我们选择高变异基因（HVGs），以减少后续步骤的计算需求。不定义HVGs也可以执行，但需要额外的计算内存和时间。

In [5]:

sc.pp.highly_variable_genes(rna, batch_key="batch", flavor="seurat")

sc.pp.highly_variable_genes(rna, batch_key="batch", flavor="seurat")

In [6]:

sc.set_figure_params(facecolor="white")

sc.set_figure_params(facecolor="white")

这是所有供体的所有基因表达数据细胞的嵌入。可以看出，在群组方面，供体对细胞类型具有主导作用，并且尚未执行批量校正方法。

In [7]:

sc.pl.embedding(rna, "GEX_X_umap", color=["cell_type", "batch"])

sc.pl.embedding(rna, "GEX_X_umap", color=["cell_type", "batch"])

这里我们只观察供体 s1d1 的细胞。观察单一供体的细胞时，可以根据其提供的注释识别细胞簇。

In [8]:

adata_batch = rna[rna.obs.batch == "s1d1", :]
sc.pl.embedding(adata_batch, "GEX_X_umap", color=["cell_type", "batch"])

adata_batch = rna[rna.obs.batch == "s1d1", :] sc.pl.embedding(adata_batch, "GEX_X_umap", color=["cell_type", "batch"])

Preparation of SCENIC¶

使用 loompy，我们将基因表达值转换为 loom 文件。pyscenic 需要这种文件格式作为输入。此外，文件 allTFs_hg38.txt 定义了与转录因子相关的基因符号列表，在衡量这些基因与其他基因之间的关联时需要考虑这些符号。

In [9]:

## this file has to be downloaded if not found
!wget -nc https://raw.githubusercontent.com/aertslab/SCENICprotocol/master/example/allTFs_hg38.txt

## this file has to be downloaded if not found !wget -nc https://raw.githubusercontent.com/aertslab/SCENICprotocol/master/example/allTFs_hg38.txt

File ‘allTFs_hg38.txt’ already there; not retrieving.

In [10]:

tfs_path = "allTFs_hg38.txt"

tfs_path = "allTFs_hg38.txt"

In [11]:

loom_path = "data/neurips_processed_input.loom"
loom_path_output = "data/neurips_processed_output.loom"
tfs = [tf.strip() for tf in open(tfs_path)]

loom_path = "data/neurips_processed_input.loom" loom_path_output = "data/neurips_processed_output.loom" tfs = [tf.strip() for tf in open(tfs_path)]

作为验证，建议检查标注为转录因子（TFs）的基因是否包含在提供的输入数据中。如果未观察到高覆盖率的基因符号（例如50%或更多），可能是由于.var中声明的特征名称存在问题，例如使用了Ensembl ID而不是基因符号，或者如果基因组装是鼠类或使用不同大小写风格的其他物种，则可能存在大小写差异。

In [12]:

# as a general QC. We inspect that our object has transcription factors listed in our main annotations.
print(
    f"%{np.sum(adata_batch.var.index.isin(tfs))} out of {len(tfs)} TFs are found in the object"
)

# as a general QC. We inspect that our object has transcription factors listed in our main annotations. print( f"%{np.sum(adata_batch.var.index.isin(tfs))} out of {len(tfs)} TFs are found in the object" )

1160 out of 1797 TFs are found in the object

是否使用HVG或所有特征可以通过将标志use_hvg更改为False来定义。

In [13]:

use_hvg = True
if use_hvg:
    mask = (adata_batch.var["highly_variable"] == True) | adata_batch.var.index.isin(
        tfs
    )
    adata_batch = adata_batch[:, mask]

use_hvg = True if use_hvg: mask = (adata_batch.var["highly_variable"] == True) | adata_batch.var.index.isin( tfs ) adata_batch = adata_batch[:, mask]

为NeurIPS供体创建一个loom文件。如果在此步骤中出现错误，请确认基因、细胞ID等标签是否正确定义。

In [14]:

row_attributes = {
    "Gene": np.array(adata_batch.var.index),
}
col_attributes = {
    "CellID": np.array(adata_batch.obs.index),
    "nGene": np.array(np.sum(adata_batch.X.transpose() > 0, axis=0)).flatten(),
    "nUMI": np.array(np.sum(adata_batch.X.transpose(), axis=0)).flatten(),
}

lp.create(loom_path, adata_batch.X.transpose(), row_attributes, col_attributes)

row_attributes = { "Gene": np.array(adata_batch.var.index), } col_attributes = { "CellID": np.array(adata_batch.obs.index), "nGene": np.array(np.sum(adata_batch.X.transpose() > 0, axis=0)).flatten(), "nUMI": np.array(np.sum(adata_batch.X.transpose(), axis=0)).flatten(), } lp.create(loom_path, adata_batch.X.transpose(), row_attributes, col_attributes)

一旦生成了loom文件，我们执行pyscenic以生成转录因子（TFs）和基因之间的关联。TF-基因关联由GRNBoost推断，并通过TF和目标基因之间的方向性权重进行总结。该分析的输出是一个表格，汇总了所有报告的关联及其重要性权重。

Some steps below require indicating a number of cores (num_workers). Increase according to computing resources available

In [15]:

num_workers = 3

num_workers = 3

In [16]:

outpath_adj = "adj.csv"
if not os.path.exists(outpath_adj):
    !pyscenic grn {loom_path} {tfs_path} -o $outpath_adj --num_workers {num_workers}

outpath_adj = "adj.csv" if not os.path.exists(outpath_adj): !pyscenic grn {loom_path} {tfs_path} -o $outpath_adj --num_workers {num_workers}

Show the top of TF-target associations

In [17]:

results_adjacencies = pd.read_csv("adj.csv", index_col=False, sep=",")
print(f"Number of associations: {results_adjacencies.shape[0]}")
results_adjacencies.head()

results_adjacencies = pd.read_csv("adj.csv", index_col=False, sep=",") print(f"Number of associations: {results_adjacencies.shape[0]}") results_adjacencies.head()

# associations: 683484

Out[17]:

	TF	target	importance
0	SOX6	SLC4A1	186.835278
1	SOX6	ANK1	162.917264
2	SOX6	HBA2	143.750583
3	SOX6	SLC25A37	127.603415
4	HMGB2	UBE2S	123.732197

可视化权重的分布，以便一般检查从pyscenic获得的分位数和阈值。根据pyscenic的GRN步骤，重要性评分呈现单峰分布，负值/正值分别表示TF-基因关联的重要性较低/较高。从该分布的右尾部，我们可以恢复TFs和潜在目标基因之间最相关的相互作用，这些相互作用由基因表达值和pyscenic分析支持。

In [18]:

plt.hist(np.log10(results_adjacencies["importance"]), bins=50)
plt.xlim([-10, 10])

plt.hist(np.log10(results_adjacencies["importance"]), bins=50) plt.xlim([-10, 10])

Out[18]:

(-10.0, 10.0)

由于目标基因在启动子区域具有DNA基序（序列特异性DNA基序），这些基序可以用来将TFs与目标基因连接起来。接下来，我们使用TF关联到转录起始位点（TSSs）的注释来细化这一注释。

Download TSS annotations precalculated by Aerts's lab

In [19]:

!wget -nc https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc9nr/gene_based/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather

!wget -nc https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc9nr/gene_based/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather

File ‘hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather’ already there; not retrieving.

In [20]:

import glob
import os

# ranking databases
db_glob = "*feather"
db_names = " ".join(glob.glob(db_glob))

import glob import os # ranking databases db_glob = "*feather" db_names = " ".join(glob.glob(db_glob))

Download a catalog of motif-to-TF associations

In [21]:

!wget -nc https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl

!wget -nc https://resources.aertslab.org/cistarget/motif2tf/motifs-v9-nr.hgnc-m0.001-o0.0.tbl

File ‘motifs-v9-nr.hgnc-m0.001-o0.0.tbl’ already there; not retrieving.

In [22]:

# motif databases
motif_path = "motifs-v9-nr.hgnc-m0.001-o0.0.tbl"

# motif databases motif_path = "motifs-v9-nr.hgnc-m0.001-o0.0.tbl"

使用启动子区域基序及其基因关联的目录，通过修剪先前的邻接关系来检索子集邻接关系。此步骤在普通硬件上可能需要几分钟时间。

In [23]:

if not os.path.exists('reg.csv'):
    !pyscenic ctx adj.csv \
        {db_names} \
        --annotations_fname {motif_path} \
        --expression_mtx_fname {loom_path} \
        --output reg.csv \
        --mask_dropouts \
        --num_workers {num_workers} > pyscenic_ctx_stdout.txt

if not os.path.exists('reg.csv'): !pyscenic ctx adj.csv \ {db_names} \ --annotations_fname {motif_path} \ --expression_mtx_fname {loom_path} \ --output reg.csv \ --mask_dropouts \ --num_workers {num_workers} > pyscenic_ctx_stdout.txt

为了探索报告的候选基因，建议根据相对贡献的排名或通过视觉定义的顶分位数阈值来探索输出，以获得高信噪比。

Define custom quantiles for further exploration

In [24]:

import numpy as np

n_genes_detected_per_cell = np.sum(adata_batch.X > 0, axis=1)
percentiles = pd.Series(n_genes_detected_per_cell.flatten().A.flatten()).quantile(
    [0.01, 0.05, 0.10, 0.50, 1]
)
print(percentiles)

import numpy as np n_genes_detected_per_cell = np.sum(adata_batch.X > 0, axis=1) percentiles = pd.Series(n_genes_detected_per_cell.flatten().A.flatten()).quantile( [0.01, 0.05, 0.10, 0.50, 1] ) print(percentiles)

0.01     101.0
0.05     144.0
0.10     171.0
0.50     259.0
1.00    1387.0
dtype: float64

下图显示了每个细胞检测到的基因分布。此可视化便于在下一步中定义参数--auc_threshold。具体来说，默认参数--auc_threshold为0.05，在此图中将导致选择144个基因，作为每个细胞的参考，用于AUCell计算。此参数的修改会影响AUCell计算的AUC值估计。

In [25]:

fig, ax = plt.subplots(1, 1, figsize=(8, 5), dpi=100)
sns.distplot(n_genes_detected_per_cell, norm_hist=False, kde=False, bins="fd")
for i, x in enumerate(percentiles):
    fig.gca().axvline(x=x, ymin=0, ymax=1, color="red")
    ax.text(
        x=x,
        y=ax.get_ylim()[1],
        s=f"{int(x)} ({percentiles.index.values[i]*100}%)",
        color="red",
        rotation=30,
        size="x-small",
        rotation_mode="anchor",
    )
ax.set_xlabel("# of genes")
ax.set_ylabel("# of cells")
fig.tight_layout()

fig, ax = plt.subplots(1, 1, figsize=(8, 5), dpi=100) sns.distplot(n_genes_detected_per_cell, norm_hist=False, kde=False, bins="fd") for i, x in enumerate(percentiles): fig.gca().axvline(x=x, ymin=0, ymax=1, color="red") ax.text( x=x, y=ax.get_ylim()[1], s=f"{int(x)} ({percentiles.index.values[i]*100}%)", color="red", rotation=30, size="x-small", rotation_mode="anchor", ) ax.set_xlabel("# of genes") ax.set_ylabel("# of cells") fig.tight_layout()

此步骤将使用转录因子（TFs）计算曲线下面积（AUC）评分，总结每个细胞中观察到的基因表达如何通过上述TFs调控的目标基因调控关联起来。

Using the above-generated matrix of cell x TFs and those scores, we can calculate a new embedding using only those.

In [27]:

if not os.path.exists(loom_path_output):
    !pyscenic aucell $loom_path \
        reg.csv \
        --output {loom_path_output} \
        --num_workers {num_workers} > pyscenic_aucell_stdout.txt

if not os.path.exists(loom_path_output): !pyscenic aucell $loom_path \ reg.csv \ --output {loom_path_output} \ --num_workers {num_workers} > pyscenic_aucell_stdout.txt

In [ ]:

import json
import zlib
import base64

# collect SCENIC AUCell output
lf = lp.connect(loom_path_output, mode="r+", validate=False)
auc_mtx = pd.DataFrame(lf.ca.RegulonsAUC, index=lf.ca.CellID)
lf.close()

import json import zlib import base64 # collect SCENIC AUCell output lf = lp.connect(loom_path_output, mode="r+", validate=False) auc_mtx = pd.DataFrame(lf.ca.RegulonsAUC, index=lf.ca.CellID) lf.close()

In [28]:

import anndata as ad

ad_auc_mtx = ad.AnnData(auc_mtx)
sc.pp.neighbors(ad_auc_mtx, n_neighbors=10, metric="correlation")
sc.tl.umap(ad_auc_mtx)
sc.tl.tsne(ad_auc_mtx)

import anndata as ad ad_auc_mtx = ad.AnnData(auc_mtx) sc.pp.neighbors(ad_auc_mtx, n_neighbors=10, metric="correlation") sc.tl.umap(ad_auc_mtx) sc.tl.tsne(ad_auc_mtx)

WARNING: You’re trying to run this on 247 dimensions of `.X`, if you really want this, set `use_rep='X'`.
         Falling back to preprocessing with `sc.pp.pca` and default params.

Visualize the data base on the TF-regulons and auc_mtx generated.

In [29]:

adata_batch.obsm["X_umap_aucell"] = ad_auc_mtx.obsm["X_umap"]
adata_batch.obsm["X_tsne_aucell"] = ad_auc_mtx.obsm["X_tsne"]

adata_batch.obsm["X_umap_aucell"] = ad_auc_mtx.obsm["X_umap"] adata_batch.obsm["X_tsne_aucell"] = ad_auc_mtx.obsm["X_tsne"]

该UMAP可视化确认了SCENIC的信号能够捕捉到将大多数细胞群体继续细分为亚群的调控子。因此，TF调控子的信息能够实现细胞类型的识别。

In [30]:

sc.pl.embedding(adata_batch, basis="X_umap_aucell", color="cell_type")

sc.pl.embedding(adata_batch, basis="X_umap_aucell", color="cell_type")

A visualization of the tSNE values generated by SCENIC also confimrs this cell-type separation, for the majority of cell-types

In [31]:

sc.pl.embedding(adata_batch, basis="X_tsne_aucell", color="cell_type")

sc.pl.embedding(adata_batch, basis="X_tsne_aucell", color="cell_type")

Interpretation of results¶

In [32]:

import seaborn as sns

import seaborn as sns

In [33]:

auc_mtx["cell_type"] = adata_batch.obs["cell_type"]
mean_auc_by_cell_type = auc_mtx.groupby("cell_type").mean()

auc_mtx["cell_type"] = adata_batch.obs["cell_type"] mean_auc_by_cell_type = auc_mtx.groupby("cell_type").mean()

show the top N TF regulons by a color

In [34]:

top_n = 50
top_tfs = mean_auc_by_cell_type.max(axis=0).sort_values(ascending=False).head(top_n)
mean_auc_by_cell_type_top_n = mean_auc_by_cell_type[
    [c for c in mean_auc_by_cell_type.columns if c in top_tfs]
]

top_n = 50 top_tfs = mean_auc_by_cell_type.max(axis=0).sort_values(ascending=False).head(top_n) mean_auc_by_cell_type_top_n = mean_auc_by_cell_type[ [c for c in mean_auc_by_cell_type.columns if c in top_tfs] ]

一旦我们确定了研究的生物系统中涉及的主要转录因子调控子，我们可以基于每个细胞的评分或每种细胞类型的总体AUC（如下方蓝色热图所示），检查每个转录因子的活动。

In [35]:

sns.clustermap(
    mean_auc_by_cell_type_top_n,
    figsize=[15, 6.5],
    cmap="Blues",
    xticklabels=True,
    yticklabels=True,
)

sns.clustermap( mean_auc_by_cell_type_top_n, figsize=[15, 6.5], cmap="Blues", xticklabels=True, yticklabels=True, )

Out[35]:

<seaborn.matrix.ClusterGrid at 0x7fc377ccce10>

由于红色热图显示某些转录因子与特定细胞类型有很强的关联，我们可以通过验证它们的表达水平来进一步确认。这可以通过匹配我们想要突出的转录因子名称，并使用Scanpy的红色热图进行可视化来完成。

In [36]:

tf_names = top_tfs.index.str.replace("\(\+\)", "")
adata_batch_top_tfs = adata_batch[:, adata_batch.var_names.isin(tf_names)]

tf_names = top_tfs.index.str.replace("\(\+\)", "") adata_batch_top_tfs = adata_batch[:, adata_batch.var_names.isin(tf_names)]

In [37]:

sc.pl.matrixplot(
    adata_batch,
    tf_names,
    groupby="cell_type",
    cmap="Reds",
    dendrogram=True,
    figsize=[15, 5.5],
    standard_scale="group",
)

sc.pl.matrixplot( adata_batch, tf_names, groupby="cell_type", cmap="Reds", dendrogram=True, figsize=[15, 5.5], standard_scale="group", )

WARNING: dendrogram data not found (using key=dendrogram_cell_type). Running `sc.tl.dendrogram` with default parameters. For fine tuning it is recommended to run `sc.tl.dendrogram` independently.
WARNING: You’re trying to run this on 2785 dimensions of `.X`, if you really want this, set `use_rep='X'`.
         Falling back to preprocessing with `sc.pp.pca` and default params.

通过对TF调控子AUC评分和细胞类型的TF基因表达的视觉检查和比较，我们可以验证在某些情况下，某个TF的细胞类型特异性表达与特定细胞类型有关，其中相应的TF调控子也处于活跃状态。例如，RUNX2在pDCs中活跃，TCF7L2在CD16+ Mono中活跃，LEF1在CD4+ T激活细胞中活跃。额外的检查可以用于验证先前的见解和/或其他关联。

Takeaways¶

In this notebook, we have:

1. Setup a basic pipeline to prepare, execute, and verify SCENIC results using the NeurIPs 2021 dataset.
2. Detected TF-regulons that allow to explain the observed transcriptomes by the regulatory potential of a handful of TFs.

Quiz¶

Theory¶

1. How many genes in the human/mouse genome are considered TFs?
2. What is a TF-regulon?
3. How many DNA-binding motifs are there to each TFs? 3.1. Currently, are there more known DNA-binding motifs or TFs that bind those? 3.2. How can one reconcile the redundancy of these during the analysis?
4. Additional reading: Describe the futility theorem.

SCENIC¶

1. Describe the major steps in the SCENIC pipeline (three or more)?
2. What could happen if the number of TFs reported in your own data has a low overlap to the one externally retrieved?
3. What is the AUC score from pyscenic and how it can be helpful for interpretation of results per cell cluster?