378_单细胞转录因子活性推断¶

library(viper)       # TF活性推断
library(dorothea)    # TF-靶基因调控网络数据库
library(Seurat)      # 单细胞分析
library(dplyr)

# 步骤1：加载 DoRothEA 调控网络（取A-C级，高可信度）
data(dorothea_hs)    # 人类TF网络（hs=human sapiens）
# 小鼠用：data(dorothea_mm)

# 只保留高可信度的A-C级调控关系
regulons <- dorothea_hs %>%
  filter(confidence %in% c("A", "B", "C"))    # 过滤低可信度（D、E级证据弱）

cat("保留的TF-靶基因对数:", nrow(regulons), "\n")
cat("涉及TF数:", length(unique(regulons$tf)), "\n")

# 步骤2：准备表达矩阵（需要Z-score标准化的数据）
# VIPER需要normalized + scaled的表达矩阵
expr_matrix <- GetAssayData(seurat_obj, layer = "scale.data")  # 使用scaled data

# 步骤3：将 DoRothEA 格式转换为 VIPER 的 regulon 格式
regulon_list <- df2regulon(regulons)    # dorothea提供的转换函数

# 步骤4：用 VIPER 计算 TF 活性
tf_activity <- viper(
  eset = expr_matrix,          # 基因表达矩阵（基因×细胞）
  regulon = regulon_list,      # TF调控网络
  nes.method = "none",         # NES计算方法
  method = "scale",            # 标准化方法
  minsize = 4,                 # TF至少需要4个靶基因才计算（过滤孤立TF）
  eset.filter = FALSE,
  cores = 4,                   # 多线程
  verbose = FALSE
)
# tf_activity: TF × 细胞矩阵，值为NES（正值=活跃，负值=抑制）

# 步骤5：将TF活性整合为 Seurat 对象的新Assay
seurat_obj[["dorothea"]] <- CreateAssayObject(data = tf_activity)
DefaultAssay(seurat_obj) <- "dorothea"  # 切换到TF活性assay

# 可视化：UMAP上某个TF的活性
FeaturePlot(seurat_obj,
  features = "TP53",      # 查看TP53的活性
  reduction = "umap",
  order = TRUE            # 高活性细胞显示在前景
)

# 热图：各细胞类型的TF活性差异
# 先切回RNA assay找到差异TF
DefaultAssay(seurat_obj) <- "dorothea"
viper_diff_tf <- FindMarkers(
  seurat_obj,
  ident.1 = "T cells",      # 组1细胞类型
  ident.2 = "B cells",      # 组2细胞类型
  min.pct = 0.1
)
# 查看Top差异TF
head(viper_diff_tf[order(viper_diff_tf$avg_log2FC, decreasing = TRUE), ])

import os
import pandas as pd
import numpy as np
import scanpy as sc
from pyscenic.utils import load_motifs
from pyscenic.transform import df2regulons
from pyscenic.aucell import aucell

# ============================================================
# pySCENIC 三步法（完整流程）
# 前两步建议用命令行运行（计算量大）
# ============================================================

# 步骤1：推断基因调控网络（GRN）
# 命令行运行：
os.system("""
pyscenic grn \
    --num_workers 8 \
    --output adj.csv \
    --method grnboost2 \
    filtered_loom.loom \
    human_TFs.txt
""")
# adj.csv: TF-靶基因候选对及重要性分数

# 步骤2：用顺式调控元件筛选靶基因（pruning）
os.system("""
pyscenic ctx \
    adj.csv \
    hg38_10kbp_up_10kbp_dw_full_tx_v10_clust.genes_vs_motifs.rankings.feather \
    --annotations_fname motifs-v10nr_clust-nr.hgnc-m0.001-o0.0.tbl \
    --expression_mtx_fname filtered_loom.loom \
    --mode "dask_multiprocessing" \
    --output regulons.csv \
    --num_workers 8
""")
# regulons.csv: 经过motif验证的高可信TF-靶基因调控单元

# 步骤3：计算每个细胞的 regulon 活性（AUCell）
adata = sc.read_h5ad("preprocessed.h5ad")
regulons_df = pd.read_csv("regulons.csv", index_col=0)

# 转换为regulon格式
regulons = df2regulons(regulons_df)

# AUCell：对每个细胞计算每个regulon的AUC分数
auc_mtx = aucell(
    adata.X,                     # 表达矩阵（细胞×基因）
    regulons,                    # TF调控单元
    num_workers=4,               # 多线程
    seed=42                      # 随机种子
)
# auc_mtx: 细胞×TF矩阵，值为AUC（0-1，越高TF越活跃）

# 整合结果到 adata
adata.obsm["X_aucell"] = auc_mtx.values     # 存储AUC矩阵
adata.uns["regulon_names"] = auc_mtx.columns.tolist()   # 存储TF名称

print(f"计算完成：{auc_mtx.shape[0]} 细胞，{auc_mtx.shape[1]} 个regulon")

# 可视化：用 heatmap 展示细胞类型和TF活性
sc.pl.heatmap(
    adata,
    var_names=auc_mtx.columns[:20].tolist(),  # Top20 TF
    groupby="cell_type",                        # 按细胞类型分组
    use_raw=False,
    layer=None,
    save="_scenic_tf_activity.pdf"
)