当前位置：首页 > news >正文

单细胞分析（21）——SCENIC 分析流程（singularity容器版）

news 2025/3/11 0:24:35

SCENIC 分析流程笔记

SCENIC (Single-Cell rEgulatory Network Inference and Clustering) 是一种基于单细胞 RNA 测序数据的调控网络分析方法，主要用于识别调控因子（TFs）及其靶基因（Regulons），并评估这些调控因子的活性。

1. 数据预处理及 LOOM 文件构建

1.1 读取单细胞数据

scRNA1 = readRDS("./tmp/scRNA_harmony_Mcell.rds")

1.2 提取基因表达矩阵

归一化表达矩阵 (normalized_data_matrix)
原始计数矩阵 (raw_count_matrix)

normalized_data_matrix <- GetAssayData(scRNA1, assay = "RNA", slot = "data")
raw_count_matrix <- GetAssayData(scRNA1, assay = "RNA", slot = "counts")

1.3 基因过滤

保留至少在 1% 细胞中表达 的基因。
仅保留 总体表达量高于 3% 细胞总数的基因。

gene_expression_sums <- Matrix::colSums(raw_count_matrix)
threshold_sum <- ncol(raw_count_matrix) * 0.03
genes_over_threshold <- gene_expression_sums > threshold_sumgenes_expressed_in_cells <- Matrix::rowSums(raw_count_matrix > 0) > (ncol(raw_count_matrix) * 0.01)
filtered_genes <- genes_over_threshold & genes_expressed_in_cells
exprMat_filtered <- raw_count_matrix[filtered_genes, ]

1.4 生成 LOOM 文件

library(SCopeLoomR)
loom <- SCopeLoomR::build_loom(file.name = "./tmp/Mcell.loom",dgem = exprMat_filtered,default.embedding = NULL
)
loom$close()

2. 共表达网络推断（GRN 构建）

Singularity 安装说明

SCENIC 分析流程依赖于 singularity 容器来运行 pyscenic，请确保系统已安装 singularity。

如果尚未安装，可以使用以下命令进行安装（以 Ubuntu 为例）：

sudo apt update && sudo apt install -y singularity-container

或者参考官方安装指南：Singularity 官方文档

2.1 运行 GRNBoost2 计算 TF-靶基因网络

time singularity run \
/home/user/container/aertslab-pyscenic-0.9.18.sif \
pyscenic grn \
--seed 777 \
--num_workers 20 \
--output /home/user/yard/project_analysis/output/Mcell_recluster/SCENIC/adj.tsv \
--method grnboost2 \
/home/user/yard/project_analysis/tmp/Mcell.loom \
/home/user/reference/scenic_resource/allTFs_hg38.txt

2.2 输出

adj.tsv：推测的 TF-靶基因共表达关系

3. Motif 富集分析（筛选 TF）

3.1 运行 motif 富集分析

time singularity run \
/home/user/container/aertslab-pyscenic-0.12.1.sif \
pyscenic ctx \
/home/user/yard/project_analysis/output/Mcell_recluster/SCENIC/adj.tsv \
/home/user/reference/scenic_resource/hg38__refseq-r80__10kb_up_and_down_tss.mc9nr.genes_vs_motifs.rankings.feather \
--annotations_fname /home/user/reference/scenic_resource/motifs-v9-nr.hgnc-m0.001-o0.0.tbl \
--expression_mtx_fname /home/user/yard/project_analysis/tmp/Mcell.loom \
--output /home/user/yard/project_analysis/output/Mcell_recluster/SCENIC/reg.csv \
--num_workers 20

3.2 输出

reg.csv：经过 motif 证据支持的 TF-靶基因调控网络

4. 计算调控因子活性（AUC 评分）

4.1 将 Regulon 转换为 GMT 文件

方法 1：使用`pyscenic to-gmt`命令

time singularity run \
/home/user/container/aertslab-pyscenic-0.12.1.sif \
pyscenic to-gmt \
/home/user/yard/project_analysis/output/Mcell_recluster/SCENIC/reg.csv \
--output /home/user/yard/project_analysis/output/Mcell_recluster/SCENIC/Mcell_regulon.regulons.gmt

方法 2：使用 Python 脚本转换

import sys
import pandas as pd
from pyscenic.cli.utils import load_signatures
# 输入和输出路径
regulon_file = "/home/user/yard/project_analysis/output/Mcell_recluster/SCENIC/reg.csv"
project = "/home/user/yard/project_analysis/output/Mcell_recluster/SCENIC/Mcell_regulon"
min_regulon_size = 10  # 过滤小的 regulon
def get_motif_logo(regulon):  base_url = "http://motifcollections.aertslab.org/v10nr_clust/logos/"  for elem in regulon.context:  if elem.endswith('.png'):  return base_url + elem
# 转换 regulons 为 GMT 格式
print("Transforming regulons to GMT file...")
sys.stdout.flush()
regulons = load_signatures(regulon_file)
select_cols = [i.name for i in regulons if len(i.genes) >= min_regulon_size]
gmt_file = project + ".regulons.gmt"
txt_file = project + ".regulons.txt"
with open(gmt_file, 'w') as fo1, open(txt_file, 'w') as fo2:for i in regulons:if i.name in select_cols:motif = get_motif_logo(i)genes = "  ".join(i.genes)tf = "%s(%sg)" % (i.transcription_factor, len(i.genes))fo1.write("%s  %s  %s
" % (tf, motif, genes))fo2.write("%s  %s  %s
" % (tf, motif, genes.replace("  ",",")))

4.2 计算 AUCell 评分

regulons <- clusterProfiler::read.gmt("./output/Mcell_recluster/SCENIC/Mcell_regulon.regulons.gmt")
seu <- ComputeModuleScore.Seurat(scRNA1, gene.sets = regulon.list, min.size = 10, cores = 16)
DefaultAssay(seu) <- "AUCell"

4.3 计算 cDC 细胞亚群 TF 活性

seu_sub = subset(seu, subset = MF_cluster_annotation %in% c('cDC_CD1C', 'cDC_LAMP3', 'cDC_CLEC9A'))
marker_aucell = FindAllMarkers(seu_sub, logfc.threshold = 0.1, only.pos = F)

4.4 可视化（热图）

library(pheatmap)
phData <- MinMax(scale(avgData), -2, 2)
pdf('./output/Mcell_recluster/SCENIC/cDC_regulon/cDC_heatmap.pdf', width = 4, height = 5)
pheatmap(phData, cluster_cols = F, cluster_rows = F,color = colorRampPalette(c("#118ab2", "#fdffb6", "#e63946"))(100))
dev.off()

5. 结果总结

基因共表达网络 (adj.tsv)：基于 GRNBoost2 计算 TF-靶基因关系。
Motif 富集分析 (reg.csv)：筛选出有 motif 证据支持的 TF-靶基因关系。
Regulon 活性评分 (AUCell 数据)：
- 计算不同细胞群体的 TF 活性。
- 通过热图可视化关键 TF 在不同细胞亚群中的调控模式。