743. 单细胞基因调控网络SCENIC+¶

SCENIC（经典版）：
  scRNA-seq数据
  → 找基因共表达模块（GRNBoost2/GENIE3）
  → motif分析验证（cisTarget）
  → 确认regulon
  → 只知道：TF → 靶基因

SCENIC+（升级版）：
  scRNA-seq + scATAC-seq
  → 找候选增强子（开放染色质区域）
  → motif分析找TF结合位点
  → 关联增强子和靶基因
  → 知道：TF → 增强子(enhancer) → 靶基因
  → 更完整的调控图谱！

# 安装pySCENIC
conda create -n scenic python=3.8  # 创建环境
conda activate scenic
pip install pyscenic  # 安装pySCENIC

# 下载必需的数据库文件
# 1. TF列表
wget https://resources.aertslab.org/cistarget/tf_lists/allTFs_hg38.txt

# 2. motif排名数据库（约10GB）
wget https://resources.aertslab.org/cistarget/databases/homo_sapiens/hg38/refseq_r80/mc_v10_clust/gene_based/hg38_10kbp_up_10kbp_down_full_tx_v10_clust.genes_vs_motifs.rankings.feather

# 3. motif注释
wget https://resources.aertslab.org/cistarget/motif2tf/motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl

# ===== pySCENIC标准流程 =====
import loompy as lp  # 导入loom文件处理
import pandas as pd
import scanpy as sc

# ===== 第一步：GRN推断（找共表达模块）=====
# 用GRNBoost2算法找共表达的TF-基因对

# 命令行方式（推荐，更快）：
"""
pyscenic grn \
  expression.loom \           # 输入：表达矩阵（loom格式）
  allTFs_hg38.txt \           # TF列表
  -o adj.csv \                # 输出：TF-基因邻接矩阵
  --num_workers 8 \           # 使用8个线程
  --method grnboost2          # 使用GRNBoost2算法（比GENIE3快）
"""

# ===== 第二步：cisTarget motif富集（验证调控关系）=====
# 在共表达模块中，检查TF的结合motif是否真的富集在靶基因的启动子区

"""
pyscenic ctx \
  adj.csv \                   # 输入：第一步的邻接矩阵
  hg38_rankings.feather \     # motif排名数据库
  --annotations_fname motifs.tbl \  # motif注释文件
  --expression_mtx_fname expression.loom \  # 表达矩阵
  --output reg.csv \          # 输出：regulon文件
  --mask_dropouts \           # 忽略dropout（零值）
  --num_workers 8             # 线程数
"""

# ===== 第三步：AUCell评分（量化regulon活性）=====
# 对每个细胞计算每个regulon的活性得分

"""
pyscenic aucell \
  expression.loom \           # 输入：表达矩阵
  reg.csv \                   # 输入：regulon文件
  --output aucell.loom \      # 输出：含AUCell得分的loom文件
  --num_workers 8             # 线程数
"""

# ===== 分析SCENIC结果 =====
import scanpy as sc
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# 加载AUCell结果
adata = sc.read_loom("aucell.loom")  # 读取含AUCell得分的数据

# AUCell矩阵：行=细胞，列=regulon
# 值=该细胞中该regulon的活性得分

# 按细胞类型查看regulon活性
sc.tl.rank_genes_groups(
    adata,
    groupby="cell_type",  # 按细胞类型分组
    method="wilcoxon"     # Wilcoxon秩和检验
)

# 查看每种细胞类型的特征regulon
sc.pl.rank_genes_groups(
    adata,
    n_genes=10,  # 每组Top10
    save="regulon_markers.png"
)

# 绘制regulon活性热图
# 选择最有差异的regulon
top_regulons = []
for group in adata.obs["cell_type"].unique():
    result = sc.get.rank_genes_groups_df(adata, group=group)
    top_regulons.extend(result.head(5)["names"].tolist())
top_regulons = list(set(top_regulons))

sc.pl.heatmap(
    adata,
    var_names=top_regulons,  # 显示的regulon
    groupby="cell_type",     # 按细胞类型分组
    standard_scale="var",    # 标准化
    cmap="RdYlBu_r",
    save="regulon_heatmap.png"
)

# ===== SCENIC+流程概述 =====
# SCENIC+需要同时有scRNA-seq和scATAC-seq数据

# 安装
# pip install scenicplus

# ===== 流程步骤 =====
"""
Step 1: 预处理
  - scRNA-seq: Scanpy标准流程
  - scATAC-seq: ArchR/Signac处理，找开放染色质peaks

Step 2: 候选增强子区域
  - 从scATAC-seq的peaks中识别候选增强子
  - 过滤启动子区域（TSS±2kb通常不算增强子）

Step 3: Motif富集
  - 在候选增强子中搜索TF结合motif
  - 使用cisTarget数据库

Step 4: 增强子-基因关联
  - 计算增强子accessibility和邻近基因表达的相关性
  - 相关性高 → 该增强子可能调控该基因

Step 5: 构建eRegulon
  - TF → (通过motif结合)增强子 → (通过相关性)靶基因
  - 这就是增强子驱动的调控网络(eRegulon)

Step 6: eRegulon活性评分
  - 类似AUCell，为每个细胞评估eRegulon活性
"""

# 简化的SCENIC+运行代码
import scenicplus  # 导入SCENIC+

# 具体实现因版本差异较大，建议参考官方教程
# https://scenicplus.readthedocs.io/

# ===== SCENIC速查 =====

# 安装
pip install pyscenic

# 三步流程
pyscenic grn expression.loom TFs.txt -o adj.csv --method grnboost2
pyscenic ctx adj.csv rankings.feather --annotations motifs.tbl -o reg.csv
pyscenic aucell expression.loom reg.csv -o aucell.loom

# 必需数据库文件
# 1. TF列表: allTFs_hg38.txt
# 2. Rankings: hg38_rankings.feather (~10GB)
# 3. Motif注释: motifs.tbl

# 结果解读
# Regulon = TF + 靶基因集合
# AUCell score = 每个细胞中regulon的活性（0-1）
# 高AUCell → 该TF在该细胞中活跃

# SCENIC+ 额外需要
# scATAC-seq数据（peak矩阵）
# 输出：eRegulon = TF → 增强子 → 靶基因