ATAC + RNA 多组学：SCENIC+ 多组学调控网络分析

作者: SeekGene

时长: 60 分钟

字数: 14.0k 字

更新: 2026-07-07

阅读: 0 次

1. 教程简介

本教程基于 SCENIC+（Single-Cell Regulatory Network Inference and Clustering Plus）框架，整合 SeekArc 单细胞多组学数据（RNA + ATAC），系统性地构建"转录因子 → 顺式调控元件 → 靶基因"的完整调控网络。

SCENIC+ 的核心优势在于：利用染色质可及性信息约束 TF-靶基因关系的推断，相比传统 SCENIC 方法，显著降低假阳性连接，提升调控网络的可解释性。

完整分析流程（step0 → step4）：

step0：Seurat → AnnData 转换

功能：从 Seurat 对象中提取 RNA 和 ATAC assay 的 counts 矩阵，转换为 AnnData（h5ad）格式，供下游 pycisTopic 和 SCENIC+ pipeline 使用。同时支持外部元数据文件的合并与可选降采样。
输入：
- *.rds — Seurat 对象文件（需同时包含 RNA 和 ATAC assay）
- meta.tsv — 外部元数据文件（可选，TSV 格式，含细胞类型注释列）
输出：
- scRNA.h5ad — RNA counts 矩阵 + 细胞元数据（cells × genes）
- scATAC.h5ad — ATAC counts 矩阵 + 细胞元数据（cells × peaks，peak 名称标准化为 chr:start-end 格式）
关键参数：
- --input_rds：输入 Seurat RDS 文件路径
- --meta_path：外部元数据文件路径（可选）
- --celltype_col：细胞类型注释列名（默认 CellAnnotation）
- --output_dir：输出目录
- --downsample / --downsample_num：是否降采样及每类细胞保留数量
运行环境：R（Seurat + reticulate + anndata），Python 环境用于 h5ad 写入
运行脚本：0_seurat_to_anndata.R
调用示例：

bash

   Rscript scripts/0_seurat_to_anndata.R \
     --input_rds /path/to/seurat_object.rds \
     --meta_path /path/to/metadata.tsv \
     --celltype_col CellAnnotation \
     --output_dir /path/to/output

step1：pycisTopic 主题建模

功能：从 scATAC.h5ad 中加载 ATAC counts 矩阵，构建 CistopicObject，运行 Mallet LDA 主题建模，输出二值化主题和 consensus regions BED 文件。
输入：
- scATAC.h5ad — step0 输出的 ATAC AnnData
- *-blacklist.v2.bed — ENCODE 黑名单文件（按物种选择）
输出：
- cistopic_obj_with_models.pkl — 包含 LDA 模型的 CistopicObject
- binarized_topics.pkl — 二值化主题（每个 topic 的 DA regions）
- consensus_regions.bed — 所有区域的 BED 文件（供 step2 cisTarget 使用）
关键参数：
- --project_dir：项目根目录
- --species：物种名（mouse / human / fly / chicken / rat）
- --n_topics：LDA 主题数量（默认 40）
- --n_iter：Mallet 迭代次数（默认 500）
- --n_cpu：并行 CPU 数（默认 32）
运行环境：Python（pycisTopic + Mallet）
运行脚本：1_pycistopic_analysis.py
调用示例：

bash

python scripts/1_pycistopic_analysis.py \
     --project_dir /path/to/project \
     --scatac_h5ad /path/to/scATAC.h5ad \
     --blacklist /path/to/species-blacklist.v2.bed \
     --mallet_dir /path/to/mallet \
     --n_topics 40

step2：cisTarget 自定义数据库构建

功能：基于 step1 输出的 consensus_regions.bed，结合参考基因组 FASTA 和 Aertslab motif 集合，通过 Cluster Buster 扫描构建 cisTarget motif 打分数据库（rankings + scores），为 step3 的 motif 富集分析提供必需输入。
输入：
- consensus_regions.bed — step1 输出的区域 BED 文件
- genome.fa — 参考基因组 FASTA（按物种选择）
- *.chrom.sizes — 染色体大小文件（按物种选择）
- *.cb — Aertslab motif collection 文件
- cbust — Cluster Buster 可执行文件
输出：
- <species>_custom.regions_vs_motifs.rankings.feather — motif 富集 rankings 数据库
- <species>_custom.regions_vs_motifs.scores.feather — motif 富集 scores 数据库
关键参数：
- --project_dir：项目根目录
- --species：物种名（mouse / human / fly / chicken / rat）
运行环境：Bash + Python（cisTarget 数据库工具链）
运行脚本：2_create_cistarget_db.sh
调用示例：

bash

bash scripts/2_create_cistarget_db.sh \
     --project_dir /path/to/project \
     --species mouse \
     --genome_fasta /path/to/genome.fa \
     --chrom_sizes /path/to/species.chrom.sizes \
     --cbust_path /path/to/cbust \
     --motif_dir /path/to/motifs/singletons \
     --script_dir /path/to/create_cisTarget_databases \
     --python_exec /path/to/python

step3：SCENIC+ pipeline 核心推断

功能：整合 step0~2 的所有产出，初始化 Snakemake pipeline、生成 config.yaml、执行调控网络推断。完成 TF→调控元件→靶基因的关系推断、cisTarget motif 富集、差异可及性分析和 AUCell 活性评分。
输入：
- scRNA.h5ad — step0 输出
- cistopic_obj_with_models.pkl / binarized_topics.pkl — step1 输出
- <species>_custom.*.feather — step2 输出（rankings + scores）
输出：
- scplusmdata.h5mu — 最终 MuData 对象（含 direct/extended eRegulon 的 AUC 矩阵）
- ctx_results.hdf5 — cisTarget motif 富集结果
- dem_results.hdf5 — 差异可及性分析结果
- ACC_GEX.h5mu — ATAC+GEX 联合 MuData
- cistromes_direct.h5ad / cistromes_extended.h5ad — direct/extended cistrome
- AUCell_direct.h5mu / AUCell_extended.h5mu — direct/extended AUCell 活性矩阵
- 其他中间文件：genome_annotation.tsv、search_space.tsv、eRegulons 等
关键参数：
- --project_dir：项目根目录
- --species：物种名（mouse / human / fly / chicken / rat）
- --n_cpu：Snakemake 并行 CPU 数（默认 16）
- --run_snakemake：是否自动执行 Snakemake（默认仅生成 config）
运行环境：Python（SCENIC+ + Snakemake）
调用示例（仅生成 config，不运行 Snakemake）：
运行脚本：3_run_scenicplus_pipeline.py

bash

python scripts/3_run_scenicplus_pipeline.py \
     --project_dir /path/to/project \
     --scRNA_h5ad /path/to/scRNA.h5ad \
     --cistopic_obj /path/to/cistopic_obj_with_models.pkl \
     --binarized_topics /path/to/binarized_topics.pkl \
     --cistarget_dir /path/to/cisTarget \
     --motif_annotation_dir /path/to/motif_snapshots \
     --species mouse \
     --n_cpu 16

调用示例（生成 config + 自动运行 Snakemake）：

bash

python scripts/3_run_scenicplus_pipeline.py \
     --project_dir /path/to/project \
     --scRNA_h5ad /path/to/scRNA.h5ad \
     --cistopic_obj /path/to/cistopic_obj_with_models.pkl \
     --binarized_topics /path/to/binarized_topics.pkl \
     --cistarget_dir /path/to/cisTarget \
     --motif_annotation_dir /path/to/motif_snapshots \
     --species mouse \
     --n_cpu 16
     --run_snakemake

step4：SCENIC+ 数据后处理（本教程主体）

本 notebook 涵盖 step4 的完整后处理分析流程，对 step3 输出的 MuData 对象进行质量评估、过滤和可视化：

数据加载与元数据检查：读取 scplusmdata.h5mu，检查 direct/extended eRegulon 元数据。
SCENICplus 对象构建：将 MuData 转换为 SCENIC+ 内部对象，整合 Cistarget 和 DEM 结果。
eRegulon 质量评估与过滤：计算 TF 表达-AUC 相关性、Gene/Region AUC 一致性，筛选高质量 eRegulon。
eRegulon 相似性分析：构建 AUC 相关性热图与 Jaccard 交集热图，评估功能冗余性。
Regulon 特异性评分（RSS）：计算 eRegulon 在各细胞类型中的特异性得分，识别细胞类型特异调控因子。
降维与多模态可视化：基于 eRegulon AUC 矩阵进行 UMAP/t-SNE 降维，将 TF 表达、Gene AUC、Region AUC 映射到同一底图。
热图-点图联合展示：整合 RSS 打分与 TF 表达 Z-score，构建综合调控特征图谱。
分析结果持久化：将完整分析结果序列化为 pickle 文件，供后续复用。

💡 说明：文档中展示的图表仅为部分代表性的可视化结果。如需查看所有细胞类型或 eRegulon 的完整分析结果（包含详细的数据表格与所有高清绘图），请前往 ./result/ 目录进行查阅。

python

import os
import sys
import pandas as pd
import numpy as np
import scanpy as sc
import mudata
from plotnine import *
from sklearn.preprocessing import StandardScaler
from scenicplus.RSS import regulon_specificity_scores

python

# Parameters

#data_dir：SCENIC+ pipeline 输出目录路径，包含 scplusmdata.h5mu、ctx_results.hdf5、dem_results.hdf5 等文件。
data_dir = '/path/to/scenicplus_pipeline'

#celltype_col：细胞类型注释列名，用于分组分析和可视化着色。
celltype_col = 'Celltype'

#R2G_positive：是否仅保留 R→G 正向调控关系的 eRegulon（+/+ 或 -/+）。
R2G_positive = True

#AUC_Gene_based_Region_based_cor：Gene-based 与 Region-based AUC 相关性的最低阈值（默认 0.5）。
AUC_Gene_based_Region_based_cor = 0.5

#TF_eRegulon_cor：TF 表达与 eRegulon AUC 相关性的绝对值阈值（默认 0.5）。
TF_eRegulon_cor = 0.5

#min_target_genes：eRegulon 最少靶基因数，过滤低信息量的 regulon。
min_target_genes = 10

2. 数据加载与元数据检查

2.1 读取 MuData 对象

从 SCENIC+ pipeline 的输出中加载包含多组学数据的 MuData 对象。该对象整合了 RNA 表达矩阵、ATAC 可及性矩阵、direct/extended eRegulon 的 AUC 矩阵以及丰富的元数据信息。

python

# Cell 1: 只运行一次（初始化）
import warnings
from IPython.utils.io import capture_output
from itables import init_notebook_mode
import itables.options as opt

opt.warn_on_undocumented_option = False
warnings.filterwarnings("ignore", category=SyntaxWarning, module=r"itables\.typing")

with capture_output():
    init_notebook_mode(all_interactive=False)

python

scplus_mdata = mudata.read(os.path.join(data_dir, "scplusmdata.h5mu"))

python

from pathlib import Path

prefix = "scATAC_counts"
# 假设 celltype_col 已定义
# celltype_col = "cell_type"  # 示例
#data_dir=Path("./scenicplus_pipeline/")
celltype_col_full = f"{prefix}:{celltype_col}"
outdir = Path("./result")  # 将字符串转换为 Path 对象
outdir.mkdir(parents=True, exist_ok=True)

2.2 direct eRegulon 元数据预览

核心步骤解析：

提取元数据：从 MuData 的 uns 属性中提取 direct_e_regulon_metadata，包含 eRegulon 的基本信息（TF 名称、靶基因数、调控方向等）。
数值精度优化：将所有数值列保留 3 位有效数字，提升表格的可读性。
导出原始数据：将原始元数据导出为 CSV 文件，供后续离线分析或文献比对使用。

💡 说明： direct eRegulon 指 TF 通过直接结合染色质可及性区域来调控靶基因的调控关系。这类 eRegulon 的调控证据最强，假阳性率最低。

python

# ===== 1) 导入 =====
import pandas as pd
import numpy as np

# ===== 2) 准备数据 =====
df_metadata = scplus_mdata.uns["direct_e_regulon_metadata"].copy()

df_display = df_metadata.copy()
# 导出原始数据（确保 "./result/" 目录已存在）
df_display.to_csv("./result/direct_e_regulon_metadata_raw.csv", index=False)

# 数值列保留 3 位有效数字
numeric_cols = df_display.select_dtypes(include=[np.number]).columns
for col in numeric_cols:
    df_display[col] = df_display[col].apply(
        lambda x: float(f"{x:.3g}") if pd.notnull(x) else x
    )

# ===== 3) 简单查看前几行 =====
print("DataFrame shape:", df_display.shape)
print("\n--- 前 5 行 ---")
print(df_display.head())

output

DataFrame shape: (157905, 15)

--- 前 5 行 ---

Region Gene importance_R2G rho_R2G importance_x_rho \\
0 chr5:180815523-180816104 MGAT1 0.1440 0.2190 0.0316
1 chr7:73594104-73594910 TBL2 0.0310 0.0969 0.0030
2 chr6:158995255-158995808 TAGAP 0.0518 0.0907 0.0047
3 chr11:47390837-47391602 SPI1 0.0879 0.6070 0.0533
4 chr2:144412741-144415133 ZEB2 0.0593 0.2640 0.0157

importance_x_abs_rho TF is_extended eRegulon_name \\
0 0.0316 BACH1 False BACH1_direct_+/+
1 0.0030 BACH1 False BACH1_direct_+/+
2 0.0047 BACH1 False BACH1_direct_+/+
3 0.0533 BACH1 False BACH1_direct_+/+
4 0.0157 BACH1 False BACH1_direct_+/+

Gene_signature_name Region_signature_name importance_TF2G \\
0 BACH1_direct_+/+_(813g) BACH1_direct_+/+_(1732r) 3.230
1 BACH1_direct_+/+_(813g) BACH1_direct_+/+_(1732r) 0.915
2 BACH1_direct_+/+_(813g) BACH1_direct_+/+_(1732r) 1.250
3 BACH1_direct_+/+_(813g) BACH1_direct_+/+_(1732r) 1.260
4 BACH1_direct_+/+_(813g) BACH1_direct_+/+_(1732r) 2.480

regulation rho_TF2G triplet_rank
0 1.0 0.388 7740.0
1 1.0 0.181 57800.0
2 1.0 0.205 36000.0
3 1.0 0.524 24300.0
4 1.0 0.553 48800.0

2.3 extended eRegulon 元数据预览

与 direct eRegulon 类似，extended eRegulon 通过基因距离或染色质互作等间接信息推断 TF-靶基因关系。这类 eRegulon 覆盖范围更广，但证据强度略低于 direct 类型。两者的对比分析有助于区分直接与间接调控信号。

python

# ===== 提取与预处理 extended_e_regulon_metadata =====
df_metadata_ext = scplus_mdata.uns["extended_e_regulon_metadata"].copy() 

df_display_ext = df_metadata_ext.copy() 
df_display_ext.to_csv("./result/extended_e_regulon_metadata_raw.csv", index=False) 

# 找出所有的数值列，并保留 3 位有效数字
numeric_cols_ext = df_display_ext.select_dtypes(include=[np.number]).columns 
for col in numeric_cols_ext: 
    df_display_ext[col] = df_display_ext[col].apply(
        lambda x: float(f'{x:.3g}') if pd.notnull(x) else x
    ) 

# ===== 简单查看前几行 =====
print("Extended DataFrame shape:", df_display_ext.shape)
print("\n--- 前 5 行 ---")
print(df_display_ext.head())

output

Extended DataFrame shape: (48420, 15)

--- 前 5 行 ---

Region Gene importance_R2G rho_R2G \\
0 chr4:138241948-138242957 SLC7A11 0.0890 0.106
1 chr14:75301587-75302301 FOS 0.0414 0.429
2 chr2:27010539-27011332 TMEM214 0.0735 0.115
3 chr7:112509759-112511021 IFRD1 0.0585 0.152
4 chr14:68732843-68733669 ACTN1 0.0572 0.510

importance_x_rho importance_x_abs_rho TF is_extended \\
0 0.00942 0.00942 ATF4 True
1 0.01780 0.01780 ATF4 True
2 0.00849 0.00849 ATF4 True
3 0.00890 0.00890 ATF4 True
4 0.02920 0.02920 ATF4 True

eRegulon_name Gene_signature_name Region_signature_name \\
0 ATF4_extended_+/+ ATF4_extended_+/+_(21g) ATF4_extended_+/+_(22r)
1 ATF4_extended_+/+ ATF4_extended_+/+_(21g) ATF4_extended_+/+_(22r)
2 ATF4_extended_+/+ ATF4_extended_+/+_(21g) ATF4_extended_+/+_(22r)
3 ATF4_extended_+/+ ATF4_extended_+/+_(21g) ATF4_extended_+/+_(22r)
4 ATF4_extended_+/+ ATF4_extended_+/+_(21g) ATF4_extended_+/+_(22r)

importance_TF2G regulation rho_TF2G triplet_rank
0 1.430 1.0 0.129 4600.0
1 0.746 1.0 0.449 33600.0
2 0.823 1.0 0.156 22000.0
3 0.907 1.0 0.271 24500.0
4 0.891 1.0 0.232 23800.0

3. SCENICplus 对象构建

3.1 MuData 转 SCENICplus 对象

核心步骤解析：

Cistarget 结果整合：加载 CuB（Cis-regulatory element target binding）数据库的 motif 富集结果，将 TF motif 与染色质可及性区域关联。
DEM 结果整合：加载差异可及性分析（Differential Enrichment Matrix）结果，识别细胞类型特异的调控元件。
对象转换：将 MuData 格式的多组学数据转换为 SCENIC+ 内部对象结构，为后续 eRegulon 质量控制和过滤做准备。

⚠️ 注意： mudata_to_scenicplus 函数要求输入的 MuData 对象包含完整的 SCENIC+ pipeline 输出模态（scRNA_counts、scATAC_counts、各类 AUC 矩阵等）。如果 MuData 结构不完整，转换会抛出 KeyError。

python

from scenicplus.scenicplus_class import mudata_to_scenicplus

python

scplus_obj = mudata_to_scenicplus(
    mdata = scplus_mdata,
    path_to_cistarget_h5 = os.path.join(data_dir, "ctx_results.hdf5"),
    path_to_dem_h5 = os.path.join(data_dir, "dem_results.hdf5")
)

python

import numpy as np
import pandas as pd
from scenicplus.regulon_qc.quality_metrics import calculate_correlation
import pandas as pd
import re
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
def generate_pseudobulks_fixed(
    scplus_mudata,
    variable,
    modality,
    nr_cells_to_sample=10,
    nr_pseudobulks_to_generate=100,
    seed=555,
    normalize_data=False
):
    if variable not in scplus_mudata.obs.columns:
        raise ValueError(f"{variable} not in scplus_mudata.obs.columns")
    if modality not in scplus_mudata.mod.keys():
        raise ValueError(f"{modality} not in scplus_mudata.mod.keys()")

    np.random.seed(seed)
    data_matrix = scplus_mudata[modality].to_df()

    if normalize_data:
        # log1p(CPM), 保持二维 DataFrame，不要 sum(1)
        libsize = data_matrix.sum(axis=1).replace(0, np.nan)
        data_matrix = np.log1p(data_matrix.div(libsize, axis=0) * 1e6).fillna(0)

    variable_to_cells = (
        scplus_mudata.obs.groupby(variable).apply(lambda x: list(x.index)).to_dict()
    )

    variable_to_mean_data = {}
    for group, cells in variable_to_cells.items():
        num_to_sample = min(nr_cells_to_sample, len(cells))
        if len(cells) == 0:
            continue
        for i in range(nr_pseudobulks_to_generate):
            sampled_cells = np.random.choice(cells, size=num_to_sample, replace=False)
            variable_to_mean_data[f"{group}_{i}"] = data_matrix.loc[sampled_cells].mean(axis=0)

    return pd.DataFrame(variable_to_mean_data).T

（1）靶基因数提取函数 从 eRegulon 名称中提取靶基因数量（如 "TF_+/+(150g)" → 150），用于后续质量散点图的纵轴展示。

python

# 0) 复用你的提取函数
def extract_n_targets(reg_name):
    m = re.search(r"\((\d+)[gr]\)$", str(reg_name))
    return int(m.group(1)) if m else np.nan

（2）TF 表达量伪批量构建 对 TF 表达矩阵按细胞类型分组生成伪批量样本，并做 log1p(CPM) 标准化，消除测序深度差异。

python

# 1) TF 表达 pseudobulk（共用）
expr_pb = generate_pseudobulks_fixed(
    scplus_mudata=scplus_mdata,
    variable=celltype_col_full,
    modality='scRNA_counts',
    nr_cells_to_sample=10,
    nr_pseudobulks_to_generate=100,
    seed=555,
    normalize_data=True
)

# 2) 分别计算 direct / extended
df_direct = None
df_extended = None
modalities = ['direct_gene_based_AUC', 'extended_gene_based_AUC']

for modality in modalities:
    auc_pb = generate_pseudobulks_fixed(
        scplus_mudata=scplus_mdata,
        variable=celltype_col_full,
        modality=modality,
        nr_cells_to_sample=10,
        nr_pseudobulks_to_generate=100,
        seed=555,
        normalize_data=False
    )
    
    common_idx = expr_pb.index.intersection(auc_pb.index)
    expr_sub = expr_pb.loc[common_idx]
    auc_sub = auc_pb.loc[common_idx]
    
    reg_to_tf = {reg: reg.split('_(')[0].split('_')[0] for reg in auc_sub.columns}
    valid_regs = [reg for reg, tf in reg_to_tf.items() if tf in expr_sub.columns]
    
    A = pd.DataFrame({reg: expr_sub[reg_to_tf[reg]] for reg in valid_regs}, index=expr_sub.index)
    B = auc_sub[valid_regs]
    mapping = {reg: reg for reg in valid_regs}
    
    df_corr = calculate_correlation(A=A, B=B, mapping_A_to_B=mapping)
    df_corr['n_targets'] = df_corr['B'].map(extract_n_targets)
    
    if modality == 'direct_gene_based_AUC':
        df_direct = df_corr
    else:
        df_extended = df_corr

3.2 TF 表达-AUC 相关性可视化

将 direct 和 extended 两种 eRegulon 的质量评估结果放在同一页进行对比。

核心步骤解析：

统一阈值：以 TF_eRegulon_cor 为相关性阈值，在图中标记高质量 eRegulon 的边界。
显著性编码：以 -log10(adjusted p-value) 为散点颜色，同时展示相关性和统计显著性。
双面板对比：左侧为 direct eRegulon，右侧为 extended eRegulon，直观评估两种调控证据类型的质量分布差异。

python


# 3) 创建 1×2 子图，画在同一页
fig, axes = plt.subplots(1, 2, figsize=(11, 5))  # 宽14，高5
datasets = [
    (df_direct, 'Direct eRegulons'),
    (df_extended, 'Extended eRegulons')
]

# 统一阈值（可按需调整）
rho_thresh = TF_eRegulon_cor  # 你之前用的变量
thresholds = {"rho": [-rho_thresh, rho_thresh]}

for i, (df, title) in enumerate(datasets):
    ax = axes[i]
    if df is None or df.empty:
        ax.text(0.5, 0.5, f"No data for {title}", transform=ax.transAxes, ha='center')
        continue
    
    df_clean = df.dropna(subset=['rho', 'n_targets']).copy()
    df_clean['adj_pval'] = df_clean['pval_adj'].clip(lower=1e-300)
    
    sc = ax.scatter(
        df_clean['rho'],
        df_clean['n_targets'],
        c=-np.log10(df_clean['adj_pval']),
        s=8
    )
    
    ax.set_xlabel("Correlation coefficient")
    ax.set_ylabel("nr. targets")
    ax.set_title(title, fontweight='bold')
    
    # 竖线阈值
    ax.vlines(
        x=thresholds["rho"],
        ymin=0,
        ymax=df_clean['n_targets'].max(),
        colors='black',
        linestyles='dashed',
        linewidths=1
    )
    for thresh in thresholds["rho"]:
        ax.text(thresh, df_clean['n_targets'].max(), str(thresh), va='bottom', ha='center', fontsize=9)
    
    sns.despine(ax=ax)

# 共用 colorbar（放右侧中间）
#cbar = fig.colorbar(sc, ax=axes, shrink=0.95, pad=0.2)
cax = fig.add_axes([1, 0.15, 0.02, 0.7])
cbar = fig.colorbar(sc, cax=cax)
#cbar = fig.colorbar(sc, ax=axes, location='right', shrink=0.9, pad=0.02)
cbar.set_label('-log10(adjusted p-value)', rotation=270, labelpad=20)

plt.tight_layout()
save_path = "./result/TF_eRegulon_cor_Direct_vs_Extended.pdf"
plt.savefig(save_path, bbox_inches='tight')

plt.show()
plt.close(fig)

output

/tmp/ipykernel_1828/996563594.py:53: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
/tmp/ipykernel_1828/996563594.py:55: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
/tmp/ipykernel_1828/996563594.py:55: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
/PROJ2/FLOAT/shumeng/apps/miniconda3/envs/scenicplus/lib/python3.11/site-packages/IPython/core/pylabtools.py:170: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.

💡 解读指南
全局说明： 散点图展示了所有 eRegulon 的 TF 表达-AUC 相关性（横轴）和靶基因数量（纵轴），颜色表示统计显著性。
横轴（Correlation coefficient）： TF 表达量与 eRegulon AUC 值的 Pearson 相关系数。正值越高，说明 TF 表达升高伴随靶基因调控活性升高，调控关系越可信。
纵轴（nr. targets）： eRegulon 包含的靶基因数量。
虚线阈值： ±TF_eRegulon_cor 的边界线。落在虚线外侧的 eRegulon 为高质量调控关系，将被保留用于下游分析。
颜色梯度： 颜色越亮表示 -log10(p-value) 越大，即统计学显著性越高。

4. eRegulon 一致性过滤

4.1 Gene-based 与 Region-based AUC 一致性筛选

为获得高可信度的 eRegulon，程序对 Gene-based 和 Region-based 两种 AUC 计算方式进行一致性评估，并结合 R2G 调控方向、靶基因数量等多重标准进行严格过滤。

核心步骤解析：

base name 归一化：将 eRegulon 名称简化为基础 TF 名称（如 STAT1_+_direct_(150g) → STAT1），以便在 Gene/Region 两种计算方式间对齐。
一致性过滤：计算 Gene-based 与 Region-based AUC 在每个 eRegulon 上的 Pearson 相关性，仅保留 |correlation| > AUC_Gene_based_Region_based_cor 的 eRegulon。这确保两种独立的 AUC 计算方法得出一致的结果。
R2G 方向过滤：当 R2G_positive = True 时，仅保留正向调控关系（+/+ 或 -/+），排除反向调控（+/− 或 −/−），这些可能代表间接或抑制性调控。
Direct vs Extended 去冗余：如果同一 TF 的 direct 和 extended 版本同时存在，仅保留 direct 版本，避免重复。
靶基因数过滤：仅保留靶基因数 > min_target_genes 的 eRegulon，过滤低信息量的调控关系。
签名提取：基于过滤后的 eRegulon，生成 Gene-based 和 Region-based 的靶基因签名，用于后续富集分析。

💡 说明：
Gene-based AUC：基于靶基因表达矩阵计算的 AUC，反映 eRegulon 在基因表达层面的调控活性。
Region-based AUC：基于染色质可及性区域计算的 AUC，反映 eRegulon 在染色质层面的调控活性。
两者一致性越高，说明 eRegulon 的调控信号在多个组学层面都能被独立验证，可靠性更强。

（1）eRegulon 基础名称提取 将完整的 eRegulon 名称简化为基础 TF 名称，用于 Gene/Region 两种 AUC 矩阵之间的对齐。

python

def _eregulon_base_name(x):
    return str(x).split("_(")[0]

（2）靶基因数提取（Gene-based） 从 Gene-based eRegulon 列名中提取靶基因数量，用于过滤低信息量的调控关系。

python

def _gene_targets_from_eregulon_col(x):
    m = re.search(r"\((\d+)g\)$", str(x))
    return int(m.group(1)) if m else None

python

auc_container = scplus_obj.uns.get("eRegulon_AUC")
if auc_container is None:
    raise KeyError("scplus_obj.uns does not contain 'eRegulon_AUC'")
if "Gene_based" not in auc_container or "Region_based" not in auc_container:
    raise KeyError("scplus_obj.uns['eRegulon_AUC'] must contain 'Gene_based' and 'Region_based'")

df1 = auc_container["Gene_based"].copy()
df2 = auc_container["Region_based"].copy()

df1_base = df1.copy()
df2_base = df2.copy()
df1_base.columns = [_eregulon_base_name(c) for c in df1_base.columns]
df2_base.columns = [_eregulon_base_name(c) for c in df2_base.columns]
df1_base = df1_base.groupby(level=0, axis=1).mean()
df2_base = df2_base.groupby(level=0, axis=1).mean()

common = sorted(set(df1_base.columns).intersection(df2_base.columns))
if len(common) == 0:
    raise ValueError("No overlapping eRegulon base names between Gene_based and Region_based AUC matrices")

correlations = df1_base[common].corrwith(df2_base[common], axis=0)
correlations = correlations[abs(correlations) > AUC_Gene_based_Region_based_cor]

if R2G_positive:
    keep_base = [x for x in correlations.index if ("+/+" in x) or ("-/+" in x)]
    print(f"[Debug] 开启 R2G_positive 过滤 (+/+ 或 -/+) 后剩余: {len(keep_base)}")
else:
    keep_base = list(correlations.index)

extended = [x for x in keep_base if "extended" in str(x)]
direct = [x for x in keep_base if "extended" not in str(x)]
keep_extended = [x for x in extended if x.replace("_extended", "_direct") not in direct]
keep_base = direct + keep_extended

keep_gene = [
    x for x in auc_container["Gene_based"].columns
    if _eregulon_base_name(x) in keep_base
]
keep_gene = [
    x for x in keep_gene
    if (_gene_targets_from_eregulon_col(x) is not None) and (_gene_targets_from_eregulon_col(x) > min_target_genes)
]
keep_region = [
    x for x in auc_container["Region_based"].columns
    if _eregulon_base_name(x) in keep_base
]
# 保证 Gene/Region 两侧保留相同的 eRegulon 基础名称
keep_gene_base = {_eregulon_base_name(x) for x in keep_gene}
keep_region = [x for x in keep_region if _eregulon_base_name(x) in keep_gene_base]

scplus_obj.uns["selected_eRegulons"] = {"Gene_based": keep_gene, "Region_based": keep_region}

output

[Debug] 开启 R2G_positive 过滤 (+/+ 或 -/+) 后剩余: 108

4.2 eRegulon 提取

将过滤后的 eRegulon 转换为靶基因签名列表，用于后续的相似性分析（Jaccard 交集）。每个 eRegulon 对应一个靶基因集合，签名之间可以进行交集/并集运算，评估功能冗余性。

python

from scenicplus.eregulon_enrichment import get_eRegulons_as_signatures

scplus_obj.uns["eRegulon_signatures"] = get_eRegulons_as_signatures(
    eRegulons=scplus_obj.uns["eRegulon_metadata"]
)
gene_sig_pool = scplus_obj.uns["eRegulon_signatures"]["Gene_based"]
region_sig_pool = scplus_obj.uns["eRegulon_signatures"]["Region_based"]

selected_base = set(
    [_eregulon_base_name(x) for x in scplus_obj.uns["selected_eRegulons"]["Gene_based"]]
    + [_eregulon_base_name(x) for x in scplus_obj.uns["selected_eRegulons"]["Region_based"]]
)
final_gene_sig = [k for k in gene_sig_pool.keys() if _eregulon_base_name(k) in selected_base]
final_region_sig = [k for k in region_sig_pool.keys() if _eregulon_base_name(k) in selected_base]

scplus_obj.uns["selected_eRegulon"] = {"Gene_based": final_gene_sig, "Region_based": final_region_sig}
#print(f"[selected_eRegulon] Gene_based: {len(final_gene_sig)}, Region_based: {len(final_region_sig)}")

#print(sorted(scplus_obj.uns.keys()))

5. eRegulon 相似性分析

5.1 AUC 相关性热图构建

在 eRegulon 过滤后，为进一步评估各 eRegulon 之间的功能冗余性，程序计算基于 AUC 矩阵的 Pearson 相关性，并结合层次聚类生成热图。

核心步骤解析：

矩阵合并与标准化：将指定 signature_keys（如 Gene-based）的 AUC 矩阵按列合并。scale=True 时对数据进行 Z-score 标准化，消除量纲差异。
全零列过滤：剔除 AUC 总和为 0 的 eRegulon（在目标细胞类型中无活性），避免干扰聚类。
距离矩阵与层次聚类：以 1 - correlation 为距离度量，使用 average 链接方法进行层次聚类。通过 fcluster 在指定阈值处切割树状图，将相似的 eRegulon 归为同一簇。
聚类重排热图：按聚类标签对 eRegulon 排序，使功能相似的 eRegulon 在热图中相邻展示，便于识别调控模块。

5.2 Jaccard 交集热图构建

除了基于 AUC 值的相关性外，程序还从靶基因集合层面评估 eRegulon 之间的相似性。Jaccard 系数衡量两个 eRegulon 的靶基因交集占并集的比例，直接反映调控目标的重叠程度。

核心步骤解析：

Jaccard 系数计算：对于每对 eRegulon，计算其靶基因集合的 Jaccard 指数 = |A∩B| / |A∪B|。值域 [0, 1]，值越高说明两个 eRegulon 的靶基因越相似。
交集归一化（可选）：method='intersect' 时，以第一个 eRegulon 的靶基因为分母计算交集比例，生成非对称矩阵。
对称化处理：对 intersect 方法进行 (M + M^T) / 2 对称化，确保聚类结果合理。
层次聚类与热图：与相关性热图相同的聚类流程，但距离度量基于 1 - Jaccard。

python

from itertools import combinations
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import sklearn

from scipy.cluster.hierarchy import linkage, fcluster
from scipy.spatial.distance import squareform


def custom_correlation_heatmap(
    scplus_obj,
    auc_key='eRegulon_AUC',
    signature_keys=['Gene_based'],
    scale=False,
    linkage_method='average',
    fcluster_threshold=0.1,
    selected_regulons=None,
    cmap='viridis',
    fontsize=10,
    ax=None
):
    if scale:
        data_mat = pd.concat([
            pd.DataFrame(
                sklearn.preprocessing.StandardScaler().fit_transform(
                    scplus_obj.uns[auc_key][x].T
                ),
                index=scplus_obj.uns[auc_key][x].T.index.to_list(),
                columns=scplus_obj.uns[auc_key][x].T.columns
            )
            for x in signature_keys
        ]).T
    else:
        data_mat = pd.concat(
            [scplus_obj.uns[auc_key][x] for x in signature_keys],
            axis=1
        )

    if selected_regulons is not None:
        subset = [x for x in selected_regulons if x in data_mat.columns]
        data_mat = data_mat[subset]

    all_zero_eregs = data_mat.columns[np.where(data_mat.sum(axis=0) == 0)]
    if len(all_zero_eregs) > 0:
        data_mat = data_mat.drop(columns=all_zero_eregs)

    correlations = data_mat.corr()
    similarity = 1 - correlations

    # 保证对称且对角线为0，避免数值误差影响 squareform
    similarity = (similarity + similarity.T) / 2
    np.fill_diagonal(similarity.values, 0)

    Z = linkage(
        np.clip(squareform(similarity), 0, similarity.to_numpy().max()),
        method=linkage_method
    )
    labels = fcluster(Z, fcluster_threshold, criterion='distance')
    labels_order = np.argsort(labels)

    clustered = data_mat.iloc[:, labels_order]
    correlations = clustered.corr()

    if ax is None:
        fig, ax = plt.subplots(figsize=(10, 10))

    sns.heatmap(
        data=correlations,
        cmap=cmap,
        square=True,
        ax=ax,
        robust=True,
        cbar=True,
        cbar_kws={"shrink": 0.75, "label": "Correlation"},
        xticklabels=True,
        yticklabels=True
    )

    ax.tick_params(axis='x', labelsize=fontsize, rotation=90)
    ax.tick_params(axis='y', labelsize=fontsize, rotation=0)
    ax.set_title("Correlation heatmap", fontsize=fontsize + 2, fontweight='bold')

    return ax


def _jaccard(signature1, signature2):
    s_signature1 = set(signature1)
    s_signature2 = set(signature2)
    intersect = len(s_signature1 & s_signature2)
    union = len(s_signature1) + len(s_signature2) - intersect
    return intersect / union if union != 0 else 0


def _intersect_norm_by_one(signature1, signature2):
    s_signature1 = set(signature1)
    s_signature2 = set(signature2)
    intersect = len(s_signature1 & s_signature2)
    return intersect / len(s_signature1) if len(s_signature1) != 0 else 0


def custom_jaccard_heatmap(
    scplus_obj,
    method='jaccard',
    gene_or_region_based='Gene_based',
    signature_key='eRegulon_signatures',
    selected_regulons=None,
    linkage_method='average',
    fcluster_threshold=0.1,
    cmap='viridis',
    fontsize=10,
    vmin=None,
    vmax=None,
    ax=None
):
    signatures = scplus_obj.uns[signature_key][gene_or_region_based]

    if selected_regulons is not None:
        signatures = {k: signatures[k] for k in signatures.keys() if k in selected_regulons}

    signatures_names = list(signatures.keys())
    n_signatures = len(signatures_names)

    jaccards = np.zeros((n_signatures, n_signatures), dtype=float)

    for i, signature_1 in enumerate(signatures_names):
        for j, signature_2 in enumerate(signatures_names):
            if i == j:
                jaccards[i, j] = 1.0
            else:
                if method == 'jaccard':
                    jaccards[i, j] = _jaccard(signatures[signature_1], signatures[signature_2])
                elif method == 'intersect':
                    jaccards[i, j] = _intersect_norm_by_one(
                        signatures[signature_1], signatures[signature_2]
                    )
                else:
                    raise ValueError("method must be 'jaccard' or 'intersect'")

    # 用于聚类的矩阵必须是对称的
    if method == 'intersect':
        cluster_mat = (jaccards + jaccards.T) / 2
    else:
        cluster_mat = jaccards.copy()

    similarity = 1 - cluster_mat
    similarity = (similarity + similarity.T) / 2
    np.fill_diagonal(similarity, 0)

    Z = linkage(squareform(similarity), method=linkage_method)
    labels = fcluster(Z, fcluster_threshold, criterion='distance')
    labels_order = np.argsort(labels)

    # 热图显示原始矩阵；如果 method='intersect'，它允许非对称
    clustered_df = pd.DataFrame(
        jaccards,
        index=signatures_names,
        columns=signatures_names
    ).iloc[labels_order, labels_order]

    if ax is None:
        fig, ax = plt.subplots(figsize=(10, 10))

    sns.heatmap(
        data=clustered_df,
        cmap=cmap,
        square=True,
        ax=ax,
        robust=True,
        cbar=True,
        cbar_kws={
            "shrink": 0.75,
            "label": "Jaccard" if method == 'jaccard' else 'Intersection'
        },
        xticklabels=True,
        yticklabels=True,
        vmin=vmin,
        vmax=vmax
    )

    ax.tick_params(axis='x', labelsize=fontsize, rotation=90)
    ax.tick_params(axis='y', labelsize=fontsize, rotation=0)
    ax.set_title(
        "Jaccard heatmap" if method == 'jaccard' else "Intersection heatmap",
        fontsize=fontsize + 2,
        fontweight='bold'
    )

    return ax

5.3 双面板热图输出

将 AUC 相关性热图与 Jaccard 交集热图并排放置，从连续值和集合两个维度全面评估 eRegulon 的相似性。

python

# =========================
# 主绘图部分：两个热图放在同一个 PDF 中
# =========================
sns.set_style("white")

fig, axes = plt.subplots(
    1, 2,
    figsize=(26, 13)
    #constrained_layout=True
)

selected_regs = scplus_obj.uns['selected_eRegulons']['Gene_based']

custom_correlation_heatmap(
    scplus_obj,
    auc_key='eRegulon_AUC',
    signature_keys=['Gene_based'],
    scale=False,
    selected_regulons=selected_regs,
    linkage_method='average',
    fcluster_threshold=0.1,
    cmap='viridis',
    fontsize=10,
    ax=axes[0]
)

custom_jaccard_heatmap(
    scplus_obj,
    method='intersect',   # 这里可改成 'jaccard'
    gene_or_region_based='Gene_based',
    signature_key='eRegulon_signatures',
    selected_regulons=selected_regs,
    linkage_method='average',
    fcluster_threshold=0.1,
    cmap='viridis',
    fontsize=10,
    vmin=None,
    vmax=None,
    ax=axes[1]
)

save_path = "./result/correlation_and_intersection_heatmaps.pdf"
plt.savefig(save_path, format="pdf", bbox_inches="tight")

plt.show()
plt.close(fig)

💡 解读指南
全局说明： 左侧为 AUC 相关性热图，右侧为交集比例热图。两个热图共享相同的 eRegulon 排序（层次聚类结果）。
AUC 相关性热图： 反映 eRegulon 在细胞间的活性模式是否相似。高相关性的 eRegulon 可能调控同一通路或处于同一调控层级。
交集热图（Intersect）： 反映 eRegulon 的靶基因重叠程度。高交集的 eRegulon 可能具有功能冗余性或互为备份。
联合解读： 如果两个 eRegulon 在两种热图中都高度相似，说明它们可能代表同一调控轴的不同方面；如果仅在 AUC 热图中相似而靶基因不重叠，可能暗示间接调控关系。

6. Regulon 特异性评分（RSS）

6.1 MuData 过滤与一致性对齐

在计算 Regulon 特异性评分之前，需要将前面步骤得到的 eRegulon 过滤结果回传到 MuData 对象，确保 RSS 计算只基于高质量的 eRegulon。

核心步骤解析：

一致性交集：取 Gene-based 和 Region-based 两侧都保留的 eRegulon base 名称的交集，作为统一的过滤标准。
四模态过滤：对 direct/extended × gene/region 四种 AUC 模态分别进行变量过滤，确保下游分析的一致性。
AUC 矩阵合并：将 direct 和 extended 的 Gene-based/Region-based AUC 矩阵分别合并为统一的 AnnData 对象，为后续降维可视化做准备。

python

from scenicplus.RSS import (regulon_specificity_scores, plot_rss)
import anndata

# 将 scplus_mdata 也进行同样的 eRegulon 过滤，保证计算 RSS 时只使用保留下来的 eRegulon
selected_gene_regulons = scplus_obj.uns['selected_eRegulons']['Gene_based']
selected_region_regulons = scplus_obj.uns['selected_eRegulons']['Region_based']

# 统一口径：仅保留在 Gene/Region 两侧都存在的 eRegulon base 名称
selected_gene_base = {_eregulon_base_name(x) for x in selected_gene_regulons}
selected_region_base = {_eregulon_base_name(x) for x in selected_region_regulons}
selected_base_consensus = selected_gene_base.intersection(selected_region_base)

keep_gene_vars_direct = [v for v in scplus_mdata["direct_gene_based_AUC"].var_names 
                         if _eregulon_base_name(v) in selected_base_consensus]
keep_gene_vars_extended = [v for v in scplus_mdata["extended_gene_based_AUC"].var_names 
                           if _eregulon_base_name(v) in selected_base_consensus]
keep_region_vars_direct = [v for v in scplus_mdata["direct_region_based_AUC"].var_names 
                           if _eregulon_base_name(v) in selected_base_consensus]
keep_region_vars_extended = [v for v in scplus_mdata["extended_region_based_AUC"].var_names 
                             if _eregulon_base_name(v) in selected_base_consensus]

scplus_mdata_filt = scplus_mdata.copy()
scplus_mdata_filt.mod["direct_gene_based_AUC"] = scplus_mdata_filt.mod["direct_gene_based_AUC"][:, keep_gene_vars_direct]
scplus_mdata_filt.mod["extended_gene_based_AUC"] = scplus_mdata_filt.mod["extended_gene_based_AUC"][:, keep_gene_vars_extended]
scplus_mdata_filt.mod["direct_region_based_AUC"] = scplus_mdata_filt.mod["direct_region_based_AUC"][:, keep_region_vars_direct]
scplus_mdata_filt.mod["extended_region_based_AUC"] = scplus_mdata_filt.mod["extended_region_based_AUC"][:, keep_region_vars_extended]
scplus_mdata_filt.update()

# 为后续联合 UMAP 可视化显式构建过滤后的 Gene/Region AUC 矩阵
eRegulon_gene_AUC = anndata.concat(
    [scplus_mdata_filt["direct_gene_based_AUC"], scplus_mdata_filt["extended_gene_based_AUC"]],
    axis=1,
)
eRegulon_gene_AUC.obs = scplus_mdata_filt.obs.loc[eRegulon_gene_AUC.obs_names].copy()

eRegulon_region_AUC = anndata.concat(
    [scplus_mdata_filt["direct_region_based_AUC"], scplus_mdata_filt["extended_region_based_AUC"]],
    axis=1,
)
eRegulon_region_AUC.obs = scplus_mdata_filt.obs.loc[eRegulon_region_AUC.obs_names].copy()

6.2 RSS 计算与可视化

Regulon 特异性评分（RSS）衡量每个 eRegulon 在不同细胞类型中的特异性富集程度。高 RSS 值的 eRegulon 在特定细胞类型中活性最强，是识别细胞类型特异性调控因子的关键指标。

核心步骤解析：

RSS 计算：基于过滤后的 MuData 对象，计算每个 eRegulon 在各细胞类型中的特异性得分。得分越高，说明该 eRegulon 在该细胞类型中的活性越特异。
结果转置：将 RSS 矩阵转置为"eRegulon × 细胞类型"格式，便于阅读和导出。

python

rss = regulon_specificity_scores(
    scplus_mudata = scplus_mdata_filt,
    variable = celltype_col_full,
    modalities = ["direct_gene_based_AUC", "extended_gene_based_AUC"]
)
# 3. 展示 RSS 打分矩阵（前几行）
# 需求2：倒置展示 (行列转置：列为 CellType，行为 eRegulon)
df_rss_display = rss.copy().T.reset_index()

# 将第一列（原本的 regulon name 列）重命名为 'eRegulon'
df_rss_display.rename(columns={df_rss_display.columns[0]: 'eRegulon'}, inplace=True)

# 将所有数值列格式化，保留 3 位有效数字，方便阅读
numeric_cols = df_rss_display.select_dtypes(include=['float64', 'float32']).columns
for col in numeric_cols:
    df_rss_display[col] = df_rss_display[col].apply(lambda x: float(f'{x:.3g}') if pd.notnull(x) else x)

# 简单查看前几行
print("RSS matrix shape:", df_rss_display.shape)
print("\n--- 前 5 行 ---")
print(df_rss_display.head())

output

RSS matrix shape: (71, 11)

--- 前 5 行 ---

eRegulon Monocytes Pro B cells B cells T cells \\
0 BACH1_direct_+/+_(813g) 0.430 0.198 0.289 0.433
1 BACH2_direct_-/+_(69g) 0.301 0.178 0.205 0.588
2 BCL11B_direct_+/+_(100g) 0.216 0.187 0.236 0.630
3 BPTF_direct_+/+_(194g) 0.258 0.223 0.339 0.479
4 BRCA1_direct_+/+_(176g) 0.250 0.236 0.324 0.391

Erythroblast CMP NK cells Dividing B cells pDC Plasma Cells
0 0.213 0.188 0.276 0.211 0.191 0.178
1 0.184 0.178 0.333 0.182 0.179 0.174
2 0.180 0.176 0.331 0.197 0.177 0.174
3 0.236 0.193 0.288 0.247 0.192 0.176
4 0.335 0.199 0.255 0.295 0.190 0.174

python

import os
import matplotlib.pyplot as plt

out_dir = "./result/Regulon_Specificity_Score/"
os.makedirs(out_dir, exist_ok=True)

original_show = plt.show
plt.show = lambda *args, **kwargs: None

try:
    for cell_type in rss.index:
        safe_cell_type = str(cell_type).replace('/', '_').replace('\\', '_')
        
        # 提取单行数据
        rss_subset = rss.loc[[cell_type]]
        
        # 如果你的 plot_rss 真的可以接受 data_matrix 参数，就这么调
        plot_rss(
            data_matrix=rss_subset,
            top_n=5,
            num_columns=1
        )
        
        # 暴力抓取内存里的图保存
        fig = plt.gcf()
        save_path = os.path.join(out_dir, f"{safe_cell_type}.pdf")
        fig.savefig(save_path, format="pdf", bbox_inches='tight')
        plt.close(fig)

finally:
    plt.show = original_show

print(f"所有细胞类型的 RSS 图已成功拆分，并且 PDF 保存完整，路径: {out_dir}")

output

所有细胞类型的 RSS 图已成功拆分，并且 PDF 保存完整，路径: ./result/Regulon_Specificity_Score/

python

# ===== 重新组织：将前两张图合并到一个画布中展示 =====
if len(figures) >= 2:
    # 创建一个新的 1行2列 的画布，设置合适的尺寸
    fig_combined, axes = plt.subplots(1, 2, figsize=(16, 6)) 
    
    # 将第一张图的内容复制到左边坐标轴
    axes[0].imshow(figures[0].canvas.renderer.buffer_rgba())
    axes[0].set_title(figures[0].get_axes()[0].get_title(), fontsize=12)
    axes[0].axis('off') # 隐藏坐标轴
    
    # 将第二张图的内容复制到右边坐标轴
    axes[1].imshow(figures[1].canvas.renderer.buffer_rgba())
    axes[1].set_title(figures[1].get_axes()[0].get_title(), fontsize=12)
    axes[1].axis('off') # 隐藏坐标轴
    
    plt.tight_layout()
    plt.show() # 统一显示，避免遮挡
    
elif len(figures) == 1:
    display(figures[0])
else:
    print("没有生成任何图片")

💡 说明：
RSS 值解读：RSS 值越高，说明该 eRegulon 在该细胞类型中的特异性越强。通常 Top 5 RSS 的 eRegulon 代表该细胞类型的核心调控因子。
多模态整合：同时使用 direct 和 extended 两种 gene-based AUC 矩阵计算 RSS，综合利用直接与间接调控信号。

💡 解读指南
全局说明： RSS 图展示了每个细胞类型中活性最强的 Top eRegulon。条形图从高到低排列，最左侧为该细胞类型最具特异性的调控因子。
条形长度： 代表 RSS 得分。值越大，特异性越强。
跨细胞类型比较： 同一个 eRegulon 如果只在某个细胞类型中出现，说明它是该类型的高度特异调控因子；如果在多个细胞类型中都出现，则可能是广谱调控因子。
应用场景： 通过 RSS 排序，可以快速锁定每个细胞亚群的核心转录因子，为下游功能验证提供候选目标。

7. eRegulon 降维可视化

7.1 降维函数定义

基于过滤后的 eRegulon AUC 矩阵进行 UMAP/t-SNE 降维，展示细胞在调控空间中的分布模式。与基于基因表达的降维不同，eRegulon 降维直接反映调控活性的细胞异质性。

核心步骤解析：

AUC 矩阵合并：将指定 signature_keys（Gene-based / Region-based / 两者组合）的 AUC 矩阵按列拼接，构建统一的特征矩阵。
元数据对齐：从原始 MuData 对象中复制细胞元数据（如细胞类型注释）到临时 AnnData 对象，确保降维结果可以按细胞类型着色。
PCA → 邻接图 → UMAP/t-SNE：先对 AUC 矩阵进行 Z-score 标准化，然后 PCA 降维（n_pcs），基于 PCA 结果构建 KNN 图，最后生成 UMAP 和/或 t-SNE 嵌入。
结果保存：将降维坐标保存回 scplus_obj.dr_cell 字典，供后续可视化函数调用。

python

import os 
import pandas as pd 
import matplotlib.pyplot as plt 
import scanpy as sc 
import anndata 

# ==========================================
# 1. 优化后的核心降维函数 (修复了 Metadata 丢失的问题)
# ==========================================
def _run_dr(scplus_obj, mdata=None, methods=['umap', 'tsne'], scale=True, 
            signature_keys=['Gene_based', 'Region_based'], 
            suffix='', selected_regulons=None, 
            n_pcs=50, random_state=0): 
    # 获取所有的 AUC 数据并合并 
    auc_mats = [] 
    for k in signature_keys: 
        df = scplus_obj.uns['eRegulon_AUC'][k].copy() 
        if selected_regulons is not None: 
            # 过滤掉不在 selected_regulons 中的列 
            keep_cols = [c for c in df.columns if _eregulon_base_name(c) in [_eregulon_base_name(x) for x in selected_regulons]] 
            df = df[keep_cols] 
        auc_mats.append(df) 
        
    merged_auc = pd.concat(auc_mats, axis=1) 
    
    # 构造临时的 AnnData 
    adata = anndata.AnnData(X=merged_auc.values) 
    adata.obs_names = merged_auc.index 
    adata.var_names = merged_auc.columns 
    
    # 【关键修复】：直接从原始的 scplus_mdata 中把 obs (metadata) 复制过来
    if mdata is not None:
        common_cells = adata.obs_names.intersection(mdata.obs_names)
        adata = adata[common_cells].copy()
        adata.obs = mdata.obs.loc[adata.obs_names].copy()
    elif hasattr(scplus_obj, 'metadata_cell'):
        common_cells = adata.obs_names.intersection(scplus_obj.metadata_cell.index)
        adata = adata[common_cells].copy()
        adata.obs = scplus_obj.metadata_cell.loc[adata.obs_names].copy()
    
    # 降维预处理
    if scale: 
        sc.pp.scale(adata) 
        
    sc.tl.pca(adata, svd_solver='arpack', random_state=random_state) 
    sc.pp.neighbors(adata, n_pcs=min(n_pcs, adata.obsm['X_pca'].shape[1]), random_state=random_state) 
    
    # 初始化 dr_cell
    if not hasattr(scplus_obj, 'dr_cell') or scplus_obj.dr_cell is None: 
        scplus_obj.dr_cell = {} 
        
    # 根据需求跑 UMAP / tSNE，并同步保存到 scplus_obj
    if 'umap' in methods:
        sc.tl.umap(adata, random_state=random_state) 
        scplus_obj.dr_cell[f'eRegulons_UMAP{suffix}'] = adata.obsm['X_umap']
        print(f"[DR] Finished UMAP for {signature_keys}. Stored as eRegulons_UMAP{suffix}")
        
    if 'tsne' in methods:
        sc.tl.tsne(adata, random_state=random_state, use_rep='X_pca') 
        scplus_obj.dr_cell[f'eRegulons_tSNE{suffix}'] = adata.obsm['X_tsne']
        print(f"[DR] Finished tSNE for {signature_keys}. Stored as eRegulons_tSNE{suffix}")
        
    return adata

7.2 降维执行与可视化

核心步骤解析：

过滤条件提取：从 scplus_obj.uns 中获取前面步骤筛选后的 eRegulon 列表。
联合降维：将 Gene-based 和 Region-based 的 eRegulon 合并，在统一的调控空间中进行降维。这种联合方式可以捕获跨模态的调控信号。
双面板可视化：左侧 UMAP 图（不显示图例，避免遮挡），右侧 t-SNE 图（保留图例），同一细胞类型着色方案贯穿两张图。

python

# ==========================================
# 2. 执行计算 (传入 scplus_mdata 以保留 metadata)
# ==========================================

# 确保输出目录存在 
out_dir = "./result/eRegulon_reduction/" 
os.makedirs(out_dir, exist_ok=True) 

# 提取过滤条件
selected_gene_regulons = scplus_obj.uns['selected_eRegulons']['Gene_based'] 
selected_region_regulons = scplus_obj.uns['selected_eRegulons']['Region_based'] 
selected_regulons_combined = sorted(set(selected_gene_regulons) | set(selected_region_regulons))


# --- A. Gene-based 降维 ---
#adata_gene = _run_dr(scplus_obj, 
#                     mdata=scplus_mdata,   # <--- 这里传入了 scplus_mdata
#                     methods=['umap', 'tsne'], 
#                     signature_keys=['Gene_based'], 
#                     suffix='_gb', 
#                     selected_regulons=selected_gene_regulons)

# --- B. Region-based 降维 ---
#adata_region = _run_dr(scplus_obj, 
#                       mdata=scplus_mdata, # <--- 这里传入了 scplus_mdata
#                       methods=['umap', 'tsne'], 
#                       signature_keys=['Region_based'], 
#                       suffix='_rb', 
#                       selected_regulons=selected_region_regulons)

# --- C. Combined (Gene + Region) 降维 ---

adata_combined = _run_dr(scplus_obj, 
                         mdata=scplus_mdata, # <--- 这里传入了 scplus_mdata
                         methods=['umap', 'tsne'], 
                         signature_keys=['Gene_based', 'Region_based'], 
                         suffix='_combined', 
                         selected_regulons=selected_regulons_combined)


# ==========================================
# 3. 统一封装画图与保存逻辑
# ==========================================
def plot_and_save_reduction(adata, title_prefix, filename):
    fig, axes = plt.subplots(1, 2, figsize=(13, 4.5), constrained_layout=True) 
    
    # 左侧 UMAP
    sc.pl.umap(adata, color=celltype_col_full, ax=axes[0], 
               title=f"{title_prefix} UMAP", show=False, legend_loc=None) 
    # 右侧 tSNE
    sc.pl.tsne(adata, color=celltype_col_full, ax=axes[1], 
               title=f"{title_prefix} tSNE", show=False) 
    
    fig.savefig(os.path.join(out_dir, filename), format="pdf", bbox_inches='tight') 
    plt.show()
    plt.close(fig) # 关闭画布，不打印输出 

# 依次画图并保存 PDF
#plot_and_save_reduction(adata_gene, "Filtered Gene-based", "eRegulon_gene_AUC_reduction.pdf")
#plot_and_save_reduction(adata_region, "Filtered Region-based", "eRegulon_region_AUC_reduction.pdf")
plot_and_save_reduction(adata_combined, "Combined", "eRegulon_combined_AUC_reduction.pdf")

print("降维与绘图全部完成，PDF 已保存！")

output

[DR] Finished UMAP for ['Gene_based', 'Region_based']. Stored as eRegulons_UMAP_combined
[DR] Finished tSNE for ['Gene_based', 'Region_based']. Stored as eRegulons_tSNE_combined

/PROJ2/FLOAT/shumeng/apps/miniconda3/envs/scenicplus/lib/python3.11/site-packages/scanpy/plotting/_tools/scatterplots.py:378: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
/PROJ2/FLOAT/shumeng/apps/miniconda3/envs/scenicplus/lib/python3.11/site-packages/scanpy/plotting/_tools/scatterplots.py:378: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored

降维与绘图全部完成，PDF 已保存！

8. TF 表达与 AUC 联合 UMAP 可视化

将 TF 的 RNA 表达量、Gene-based AUC 和 Region-based AUC 三种信号映射到同一张 UMAP 底图上，实现从"调控发起（TF 表达）→ 染色质变化（Region AUC）→ 基因表达（Gene AUC）"的完整逻辑链条可视化。

核心步骤解析：

轴方向检测：自动判断 X_EXP 矩阵的基因/细胞轴方向（行 vs 列），兼容不同数据源的存储格式。
名称映射：将 eRegulon 名称中的 TF 部分与 X_EXP 中的基因名精确匹配，同时建立 Gene/Region AUC 列名的映射表。
Top eRegulon 选取：对每个细胞类型，选取 RSS 最高的 Top 5 eRegulon 作为可视化目标。
表达量平滑：当 TF 表达量的最大值 > 50 时，进行 log1p 平滑处理，压缩极端高表达值对颜色分布的影响。
UMAP 映射：将 TF 表达、Gene-based AUC、Region-based AUC 三个信号分别映射到三格 UMAP 图上，使用 vmax="p99" 截断前 1% 的极端值，确保中等表达特异性的可视化效果。

💡 说明：
log1p 平滑：单细胞 RNA 表达数据常存在极端高表达值，log1p 变换可以压缩动态范围，使颜色梯度更均匀。
p99 截断：将 UMAP 图例上限设为第 99 百分位数，避免少数极端细胞将整个颜色尺度拉偏。

使用 Combined 降维（_run_dr 输出）作为统一底图，绘制：

TF expression, 2) TF 对应 Gene-based AUC, 3) TF 对应 Region-based AUC

python

import os
import numpy as np
import matplotlib.pyplot as plt
import scanpy as sc

# 使用 Combined 降维（_run_dr 输出）作为统一底图，绘制： 
# 1) TF expression, 2) TF 对应 Gene-based AUC, 3) TF 对应 Region-based AUC 
if not hasattr(scplus_obj, "X_EXP"): 
    raise AttributeError("scplus_obj 缺少 X_EXP，无法提取 TF 表达量。") 

exp_df = scplus_obj.X_EXP 
if not hasattr(exp_df, "index") or not hasattr(exp_df, "columns"): 
    raise TypeError("scplus_obj.X_EXP 不是带基因和细胞索引的表结构。") 

umap_out_dir = "./result/TF_AUC_Combined_UMAP/" 
os.makedirs(umap_out_dir, exist_ok=True) 

# 自动识别 X_EXP 的轴方向：基因在行(index) 还是在列(columns) 
adata_cells = set(adata_combined.obs_names.astype(str)) 
idx_names = exp_df.index.astype(str) 
col_names = exp_df.columns.astype(str) 
idx_cell_hits = len(adata_cells.intersection(set(idx_names))) 
col_cell_hits = len(adata_cells.intersection(set(col_names))) 

# 哪个轴更像“细胞轴”，另一个轴即“基因轴” 
if idx_cell_hits >= col_cell_hits: 
    exp_cells_axis = "index" 
    exp_genes_axis = "columns" 
    exp_cell_names = exp_df.index 
    exp_gene_names = exp_df.columns 
else: 
    exp_cells_axis = "columns" 
    exp_genes_axis = "index" 
    exp_cell_names = exp_df.columns 
    exp_gene_names = exp_df.index 

print(f"[EXP] detected cells on {exp_cells_axis}, genes on {exp_genes_axis}") 

exp_gene_map = {} 
for g in exp_gene_names.astype(str): 
    exp_gene_map.setdefault(g.lower(), g) 

gene_var_map = {} 
for v in eRegulon_gene_AUC.var_names.astype(str): 
    gene_var_map.setdefault(_eregulon_base_name(v), v) 

region_var_map = {} 
for v in eRegulon_region_AUC.var_names.astype(str): 
    region_var_map.setdefault(_eregulon_base_name(v), v) 

def _to_1d_array(x): 
    if hasattr(x, "toarray"): 
        return x.toarray().ravel() 
    if hasattr(x, "todense"): 
        return np.asarray(x.todense()).ravel() 
    return np.asarray(x).ravel() 

common_cells = ( 
    adata_combined.obs_names 
    .intersection(eRegulon_gene_AUC.obs_names) 
    .intersection(eRegulon_region_AUC.obs_names) 
    .intersection(exp_cell_names) 
) 

for ct in rss.index: 
    top5_eregulons = rss.loc[ct].sort_values(ascending=False).head(5).index.tolist() 
    for ereg in top5_eregulons: 
        tf_name = str(ereg).split("_")[0] 
        matched_tf = exp_gene_map.get(tf_name.lower()) 
        gene_col = gene_var_map.get(_eregulon_base_name(str(ereg))) 
        region_col = region_var_map.get(_eregulon_base_name(str(ereg))) 

        if not matched_tf: 
            print(f"Warning: TF {tf_name} not found in RNA matrix.") 
            continue 
        if gene_col is None or region_col is None: 
            print(f"Warning: regulon {ereg} not found in gene/region AUC matrices.") 
            continue 

        plot_adata = adata_combined[common_cells].copy() 

        # 提取表达量
        if exp_genes_axis == "index": 
            expr_vals = _to_1d_array(exp_df.loc[matched_tf, common_cells].values) 
        else: 
            expr_vals = _to_1d_array(exp_df.loc[common_cells, matched_tf].values) 
            
        if np.max(expr_vals) > 50:
            expr_vals = np.log1p(expr_vals)

        gene_auc_vals = _to_1d_array(eRegulon_gene_AUC[common_cells, gene_col].X) 
        region_auc_vals = _to_1d_array(eRegulon_region_AUC[common_cells, region_col].X) 

        expr_key = f"{tf_name}__expr" 
        gene_key = f"{tf_name}__gene_auc" 
        region_key = f"{tf_name}__region_auc" 
        
        plot_adata.obs[expr_key] = expr_vals 
        plot_adata.obs[gene_key] = gene_auc_vals 
        plot_adata.obs[region_key] = region_auc_vals 

        fig, axes = plt.subplots(1, 3, figsize=(15, 4.5)) 
        sc.pl.umap(plot_adata, color=expr_key, cmap="viridis", vmax="p99", show=False, ax=axes[0], title=f"{tf_name} Expression") 
        sc.pl.umap(plot_adata, color=gene_key, cmap="viridis", vmax="p99", show=False, ax=axes[1], title=f"{tf_name} Gene-based AUC") 
        sc.pl.umap(plot_adata, color=region_key, cmap="viridis", vmax="p99", show=False, ax=axes[2], title=f"{tf_name} Region-based AUC") 

        safe_ereg = str(ereg).replace("/", "_").replace("+", "pos").replace("-", "neg") 
        safe_ct = str(ct).replace("/", "_").replace(" ", "_") 
        pdf_path = os.path.join(umap_out_dir, f"{safe_ct}_combined_umap_{safe_ereg}.pdf") 
        fig.savefig(pdf_path, bbox_inches="tight") 
        plt.close(fig) 

print(f"Combined UMAP plots saved to {os.path.abspath(umap_out_dir)}")

output

[EXP] detected cells on index, genes on columns
Combined UMAP plots saved to /home/shumeng/workspace/project/shumeng/result/TF_AUC_Combined_UMAP

python

# ===== 展示最后一张图 =====
if last_fig is not None:
    # 直接显示最后一个 figure
    display(last_fig)
    print(f"展示最后一张图: {pdf_path}")
else:
    print("没有生成任何图片")

展示最后一张图: ./result/TF_AUC_Combined_UMAP/Plasma_Cells_combined_umap_RELB_direct_pos_pos_(69g).pdf

💡 说明：
Gene-based vs Region-based：Gene-based 降维反映基于靶基因表达模式的细胞分组；Region-based 降维反映基于染色质可及性模式的分组。
Combined 降维：同时使用两种 AUC 矩阵，捕获最全面的调控异质性信号，推荐作为主要分析结果。
注释代码：单独的 Gene-based 和 Region-based 降维代码已被注释。如需启用，取消对应注释即可。

9. 热图-点图联合可视化

9.1 SCENIC+ 内置 heatmap_dotplot

利用 SCENIC+ 内置的 heatmap_dotplot 函数，以点大小表示 Region-based AUC、以颜色表示 Gene-based AUC，在同一张图中同时展示两种模态的调控活性。

核心步骤解析：

元数据过滤：仅保留在 Gene/Region 两侧过滤后都存在的 regulon，确保点图展示的均为高质量 eRegulon。
双色道编码：颜色通道（Gene-based AUC）反映靶基因表达层面的调控活性；点大小通道（Region-based AUC）反映染色质可及性层面的调控活性。
水平布局：eRegulon 按行排列，细胞类型按列排列，便于逐一查看每个 eRegulon 的细胞类型特异性模式。

python

from scenicplus.plotting.dotplot import heatmap_dotplot

python

# 过滤后可用的列名
avail_gene = set(scplus_mdata_filt.mod["direct_gene_based_AUC"].var_names.astype(str))
avail_region = set(scplus_mdata_filt.mod["direct_region_based_AUC"].var_names.astype(str))

meta = scplus_mdata_filt.uns["direct_e_regulon_metadata"].copy()
meta = meta[
    meta["Gene_signature_name"].astype(str).isin(avail_gene) &
    meta["Region_signature_name"].astype(str).isin(avail_region)
].copy()

subset_names = meta["eRegulon_name"].astype(str).unique().tolist()

plot = heatmap_dotplot(
    scplus_mudata=scplus_mdata_filt,
    color_modality="direct_gene_based_AUC",
    size_modality="direct_region_based_AUC",
    group_variable=celltype_col_full,
    eRegulon_metadata_key="direct_e_regulon_metadata",
    color_feature_key="Gene_signature_name",
    size_feature_key="Region_signature_name",
    feature_name_key="eRegulon_name",
    subset_feature_names=subset_names,   # 关键：只画过滤后还存在的 regulon
    sort_data_by="direct_gene_based_AUC",
    orientation="horizontal",
    figsize=(19, 8)
)
plot = plot + theme(
    axis_text_x=element_text(
     #   rotation=90,      # 旋转90度
    #    hjust=0.5,          # 水平对齐：1=右对齐
     #   vjust=0.5,        # 垂直对齐：0.5=居中
        angle=90,
        size=12           # 字体大小
      #  color='black',     # 字体颜色
     #   margin={'t': 10},
       # va='top'
    ),
    axis_text_y=element_text(size=12),  # 也可以调整y轴
    legend_position = 'bottom'
)

# 显示图形
dotplot_out_dir = "./result/Regulon_Specificity_Score/"
os.makedirs(dotplot_out_dir, exist_ok=True)
dotplot_pdf = os.path.join(dotplot_out_dir, "direct_eRegulon_heatmap_dotplot.pdf")
plot.save(dotplot_pdf, format="pdf", width=19, height=8, units="in", verbose=False)
print(f"Heatmap-dotplot 已保存: {dotplot_pdf}")

output

Heatmap-dotplot 已保存: ./result/Regulon_Specificity_Score/direct_eRegulon_heatmap_dotplot.pdf

python

plot

<Figure Size: (1900 x 800)>

9.2 自定义 RSS-表达量联合点图

在前述 SCENIC+ 内置点图基础上，进一步构建以颜色表示 TF 表达量 Z-score、以点大小表示 RSS 得分的自定义联合点图，实现调控活性（RSS）与表达水平（Z-score）的双维度可视化。

核心步骤解析：

Top TF 提取：从每个细胞类型的 RSS 矩阵中提取 Top N（默认 6）eRegulon 作为可视化目标。
表达量提取与标准化：从 RNA 模态中提取对应 TF 的表达量，按细胞类型分组计算均值，再进行 Z-score 标准化（clip 至 [-2, 2]），消除绝对量纲差异。
动态阈值：以 RSS 的 5% 分位数为阈值，过滤低特异性的 eRegulon-细胞类型组合，避免图面过于拥挤。
联合编码：背景色块表示 TF 表达 Z-score（红色 = 高表达，蓝色 = 低表达），黑色圆点大小表示 RSS 得分，实现双模态信号在同一图面中的精确对齐。

python


import pandas as pd
import numpy as np
from plotnine import *
from sklearn.preprocessing import StandardScaler

print("Building custom RSS and Expression dotplot...")

# 1. 设置参数
group_var = celltype_col_full
top_n = 6

# 2. 提取 RNA modality
rna_keys = [k for k in scplus_mdata.mod.keys() if 'RNA' in k.upper() or 'GEX' in k.upper()]
if not rna_keys:
   raise ValueError("Cannot find RNA/GEX modality in scplus_mdata")
rna_key = rna_keys[0]
adata_rna = scplus_mdata.mod[rna_key]

# 同步细胞元数据
adata_rna.obs[group_var] = scplus_mdata.obs[group_var].values

# 3. 提取 Top TFs
selected_regulons = []
for ct in rss.index:
   selected_regulons.extend(rss.loc[ct].sort_values(ascending=False).head(top_n).index.tolist())
selected_regulons = list(dict.fromkeys(selected_regulons)) # 去重

# 4. 构建作图数据框
plot_data = []
celltypes = rss.index.tolist()

for reg in selected_regulons:
   tf_name = reg.split('_')[0]
   
   if tf_name in adata_rna.var_names:
       if hasattr(adata_rna[:, tf_name].X, 'todense'):
           expr_array = np.ravel(adata_rna[:, tf_name].X.todense())
       else:
           expr_array = np.ravel(adata_rna[:, tf_name].X)
           
       df_tmp = pd.DataFrame({'expr': expr_array, 'group': adata_rna.obs[group_var].values})
       mean_expr = df_tmp.groupby('group')['expr'].mean().reindex(celltypes).fillna(0)
   else:
       mean_expr = pd.Series(0, index=celltypes)
       
   # Z-score 标准化并限制范围 [-2, 2]
   z_expr = StandardScaler().fit_transform(mean_expr.values.reshape(-1, 1)).flatten()
   z_expr = np.clip(z_expr, -2, 2)
   
   for i, ct in enumerate(celltypes):
       plot_data.append({
           'CellType': ct,
           'eRegulon': reg,
           'RSS': rss.loc[ct, reg],
           'Expression_Zscore': z_expr[i]
       })

df_plot = pd.DataFrame(plot_data)

# 5. 动态计算气泡阈值 (60% 分位数)
rss_60th = np.percentile(df_plot['RSS'], 5)
rss_max = df_plot['RSS'].max()
print(f"Calculated 60th percentile for RSS: {rss_60th:.4f}")
print(f"Calculated Maximum RSS: {rss_max:.4f}")

# 将小于 60% 分位数的 RSS 设置为 NaN，使绘图时不显示气泡
df_plot['RSS_plot'] = df_plot['RSS'].where(df_plot['RSS'] >= rss_60th, np.nan)

# 设置 Categorical 保证坐标轴系顺序
df_plot['CellType'] = pd.Categorical(df_plot['CellType'], categories=celltypes[::-1], ordered=True)
df_plot['eRegulon'] = pd.Categorical(df_plot['eRegulon'], categories=selected_regulons, ordered=True)

# 6. 画图
custom_plot = (
   ggplot(df_plot, aes(x='eRegulon', y='CellType'))
   + geom_tile(aes(fill='Expression_Zscore'), color='white', size=0.5)
   + scale_fill_cmap(cmap_name='RdYlBu_r', limits=(-2, 2))
   + geom_point(aes(size='RSS_plot'), color='black', na_rm=True)
   + scale_size_continuous(range=(1.5, 6), limits=(rss_60th, rss_max))
   + theme_minimal()
   + theme(
       axis_text_x=element_text(angle=90, hjust=1, vjust=1, size=10),
       axis_text_y=element_text(size=10),
       panel_grid_major=element_blank(),
       figure_size=(max(8, len(selected_regulons)*0.3), max(6, len(celltypes)*0.5))
   )
   + labs(x='', y='', fill='Expression\n(Z-score)', size='RSS')
)

# 显示图形
dotplot_out_dir = "./result/Regulon_Specificity_Score/"
os.makedirs(dotplot_out_dir, exist_ok=True)
dotplot_pdf = os.path.join(dotplot_out_dir, "RSS_expersion_heatmap_dotplot.pdf")
custom_plot.save(dotplot_pdf, format="pdf", width=19, height=8, units="in", verbose=False)
print(f"Heatmap-dotplot 已保存: {dotplot_pdf}")

output

Building custom RSS and Expression dotplot...n Calculated 60th percentile for RSS: 0.1751
Calculated Maximum RSS: 0.6351
Heatmap-dotplot 已保存: ./result/Regulon_Specificity_Score/RSS_expersion_heatmap_dotplot.pdf

python

custom_plot

<Figure Size: (1230 x 600)>

💡 解读指南
全局说明： 自定义点图将 RSS 得分（点大小）与 TF 表达水平（颜色）整合到同一矩阵中。每行代表一个细胞类型，每列代表一个 eRegulon。
颜色梯度： 红色表示 TF 在该细胞类型中高表达，蓝色表示低表达。Z-score > 0 意味着表达高于全局平均。
点大小： 点越大说明 RSS 得分越高，即该 eRegulon 对该细胞类型的特异性越强。
空位（无点）： 表示该 eRegulon 在该细胞类型中的 RSS 低于分位数阈值，特异性不显著。
联合解读： 理想的细胞类型特异调控因子应同时满足"点大 + 红色"——即高 RSS 特异性和高表达水平。如果点大但颜色偏蓝，可能暗示转录后调控机制。

10. 分析结果保存

将包含所有分析结果（过滤后的 eRegulon、降维坐标、元数据等）的 scplus_obj 对象序列化为 pickle 文件，供后续分析复用。

💡 说明：
使用 dill 而非标准 pickle，因为 dill 可以序列化更复杂的 Python 对象（如 lambda 函数、嵌套类等）。
序列化后的对象可通过 dill.load() 重新加载，所有 uns、dr_cell、X_EXP 等属性均完整保留。

使用原始字符串，不需要f-string

python

import dill

# 使用原始字符串，不需要f-string
pkl_file = "./scplus_obj.pkl"
with open(pkl_file, 'wb') as f:
    dill.dump(scplus_obj, f)

11. 完整流程所需的前置文件清单

按依赖顺序排列：

Step 0：Seurat → AnnData 转换

文件	说明	来源
`*.rds`	Seurat 对象，必须同时包含 RNA 和 ATAC 两个 assay	用户自带数据
`meta.tsv`（可选）	外部元数据表，TSV 格式，第一列为 cell barcode，需含 `celltype_col` 指定的列、`orig.ident`、`Sample` 列	用户自备

Step 1：pycisTopic 主题建模

文件	说明	来源
`scATAC.h5ad`	Step 0 输出的 ATAC AnnData（自动读取，无需额外准备）	step0 产物
`{species}-blacklist.v2.bed`	ENCODE 黑名单 BED 文件（按物种选择：mouse/human/fly/chicken/rat）	需提前下载，Aertslab 提供

Step 2：cisTarget 自定义数据库构建

文件	说明	来源
`consensus_regions.bed`	Step 1 输出的区域 BED（自动读取）	step1 产物
`genome.fa` + `genome.fa.fai`	参考基因组 FASTA + 索引（按物种选择）	需提前下载
`{species}.chrom.sizes`	染色体大小文件（mm10/hg38/dm6/GRCg7b/rn7）	需提前下载
`*.cb` 文件集合	Aertslab motif collection（`v10nr_clust_public/singletons/` 目录下所有 .cb 文件）	需提前下载，来自 SCENIC/AERTSLAB
`cbust` 可执行文件	Cluster Buster motif 扫描工具	需提前编译/下载
`create_fasta_with_padded_bg_from_bed.sh`	辅助脚本（从 bin/create_cisTarget_databases/ 调用）	仓库自带
`create_cistarget_motif_databases.py`	辅助脚本（从 bin/create_cisTarget_databases/ 调用）	仓库自带

Step 3：SCENIC+ Pipeline 核心推断

文件	说明	来源
`scRNA.h5ad`	Step 0 输出的 RNA AnnData	step0 产物
`cistopic_obj_with_models.pkl`	Step 1 输出的 CistopicObject	step1 产物
`binarized_topics.pkl`	Step 1 输出的二值化主题	step1 产物
`{species}_custom.*.rankings.feather`	Step 2 输出的 rankings 数据库	step2 产物
`{species}_custom.*.scores.feather`	Step 2 输出的 scores 数据库	step2 产物
`motifs-v10-nr.*.tbl`	Motif 注释表（snapshots 目录下，按物种选择）	包含在 motif collection 中
`genes.gtf`	GTF 基因组注释文件（用于生成 genome_annotation.tsv）	需提前下载，按物种选择
`{assembly}.chrom.sizes`	同 step2，需要放在 `data/` 目录下	同 step2

Step 4：后处理（notebook）

文件	说明	来源
`scplusmdata.h5mu`	Step 3 pipeline 最终输出的 MuData	step3 产物
`ctx_results.hdf5`	cisTarget 富集结果	step3 产物
`dem_results.hdf5`	差异可及性结果	step3 产物

需要提前下载并准备好的数据库/参考文件（共 6 类）

#	文件类型	物种覆盖	下载来源
1	参考基因组 FASTA	mouse (mm10), human (hg38), fly, chicken, rat	UCSC / Ensembl
2	染色体大小文件 (.chrom.sizes)	同物种	UCSC Genome Browser
3	GTF 基因组注释文件	同物种（需要包含 gene_name/gene_id、gene_type 字段）	Ensembl / GENCODE / CellRanger ARC
4	ENCODE 黑名单 (.blacklist.v2.bed)	mouse, human, fly 等	ENCODE Project
5	Aertslab motif collection（含 singletons/.cb + snapshots/.tbl）	跨物种（v10nr_clust）	SCENIC AERTSLAB
6	cbust 可执行文件	通用	AERTSLAB Cluster Buster

其中 1-4 是按物种各一份；5-6 是全项目共用一份。

ATAC + RNA 多组学：SCENIC+ 多组学调控网络分析 ​

1. 教程简介 ​

step0：Seurat → AnnData 转换 ​

step1：pycisTopic 主题建模 ​

step2：cisTarget 自定义数据库构建 ​

step3：SCENIC+ pipeline 核心推断 ​

step4：SCENIC+ 数据后处理（本教程主体） ​

2. 数据加载与元数据检查 ​

2.1 读取 MuData 对象 ​

2.2 direct eRegulon 元数据预览 ​

2.3 extended eRegulon 元数据预览 ​

3. SCENICplus 对象构建 ​

3.1 MuData 转 SCENICplus 对象 ​

3.2 TF 表达-AUC 相关性可视化 ​

4. eRegulon 一致性过滤 ​

4.1 Gene-based 与 Region-based AUC 一致性筛选 ​

4.2 eRegulon 提取 ​

5. eRegulon 相似性分析 ​

5.1 AUC 相关性热图构建 ​

5.2 Jaccard 交集热图构建 ​

5.3 双面板热图输出 ​

6. Regulon 特异性评分（RSS） ​

6.1 MuData 过滤与一致性对齐 ​

6.2 RSS 计算与可视化 ​

7. eRegulon 降维可视化 ​

7.1 降维函数定义 ​

7.2 降维执行与可视化 ​

8. TF 表达与 AUC 联合 UMAP 可视化 ​

9. 热图-点图联合可视化 ​

9.1 SCENIC+ 内置 heatmap_dotplot ​

9.2 自定义 RSS-表达量联合点图 ​

10. 分析结果保存 ​

11. 完整流程所需的前置文件清单 ​

Step 0：Seurat → AnnData 转换 ​

Step 1：pycisTopic 主题建模 ​

Step 2：cisTarget 自定义数据库构建 ​

Step 3：SCENIC+ Pipeline 核心推断 ​

Step 4：后处理（notebook） ​

需要提前下载并准备好的数据库/参考文件（共 6 类） ​

ATAC + RNA 多组学：SCENIC+ 多组学调控网络分析

1. 教程简介

step0：Seurat → AnnData 转换

step1：pycisTopic 主题建模

step2：cisTarget 自定义数据库构建

step3：SCENIC+ pipeline 核心推断

step4：SCENIC+ 数据后处理（本教程主体）

2. 数据加载与元数据检查

2.1 读取 MuData 对象

2.2 direct eRegulon 元数据预览

2.3 extended eRegulon 元数据预览

3. SCENICplus 对象构建

3.1 MuData 转 SCENICplus 对象

3.2 TF 表达-AUC 相关性可视化

4. eRegulon 一致性过滤

4.1 Gene-based 与 Region-based AUC 一致性筛选

4.2 eRegulon 提取

5. eRegulon 相似性分析

5.1 AUC 相关性热图构建

5.2 Jaccard 交集热图构建

5.3 双面板热图输出

6. Regulon 特异性评分（RSS）

6.1 MuData 过滤与一致性对齐

6.2 RSS 计算与可视化

7. eRegulon 降维可视化

7.1 降维函数定义

7.2 降维执行与可视化

8. TF 表达与 AUC 联合 UMAP 可视化

9. 热图-点图联合可视化

9.1 SCENIC+ 内置 heatmap_dotplot

9.2 自定义 RSS-表达量联合点图

10. 分析结果保存

11. 完整流程所需的前置文件清单

Step 0：Seurat → AnnData 转换

Step 1：pycisTopic 主题建模

Step 2：cisTarget 自定义数据库构建

Step 3：SCENIC+ Pipeline 核心推断

Step 4：后处理（notebook）

需要提前下载并准备好的数据库/参考文件（共 6 类）