昨天的推送原来Scenic转录因子分析升级到Scenic+了，这篇Nat Cancer做了一个示范应用中，钱博士和大家推送介绍了SCENIC+的文献应用。钱博士前几天和我聊了一下这个工具，刚好最近要用到这一部分的分析，简单分享下我的学习笔记。
对于有单细胞转录组和ATAC-seq测序数据的同学，使用SCENIC+进行转录因子分析再合适不过了。Nat Cancer的脚本是这两种测序技术进行联合分析非常好的示例代码，有兴趣的同学可以自行访问github学习（SC_TALL/SCENIC+ at main · tanlabcode/SC_TALL）。
github提供的脚本包含以下几个主要模块：
准备 FASTA 文件并基于 ATAC 峰创建峰值文件
Prepare_FASTA.sh：准备包含感兴趣区域的 FASTA 文件。
Prepare_Data_BMP_TSpec.R：使用 ATAC-seq 峰值生成 cistarget 数据库所需的bed文件。
Create_Cistarget_Motif.sh：以正确的格式准备数据库
结合ATAC-seq数据库结果，构建本地数据库以获得最佳分析结果。
T_ALL_Setup_SCENICPlus_BMP_TSpec.ipynb：设置 SCENIC+ 的 Jupyter Notebook 环境，包括 pycistopic 预处理
配置 Jupyter Notebook 环境并安装 SCENIC+ 所需的依赖项。
使用 pycistopic 对输入数据进行分析和预处理，例如构建染色质可及性特征矩阵。
基于 config.yaml 参数运行 SCENIC+ 管道
创建 config.yaml 配置文件，设定物种、参考基因组和分析参数。
启动 SCENIC+ 管道进行转录调控网络推断和调控因子识别。
后处理和结果检查（Jupyter Notebook）
整理 SCENIC+ 输出结果，例如调控网络和调控因子目标基因集。
T_ALL_Review_SCENICPlus_BMP_TSpec.ipynb：在Jupyter Notebook 中可视化调控网络，分析关键调控因子和目标基因。
这几个脚本主要用来提取ATAC中的坐标信息，然后构建本地数据库，最后用于scenic分析。尽管scenic+提供了预构建好的数据库，不过还是推荐自己本地构建数据库，获得最佳的分析效果。
首先学习下ATAC坐标信息的提取——Prepare_Data_BMP_TSpec.R
# 加载必要的R包，用于分析Seurat对象、处理基因组数据、绘图等。library(Seurat)library(SeuratDisk)library(SeuratData)library(anndata)library(Signac)library(Matrix)library(ggplot2)# 设置工作目录到SCENICPlus项目所在路径。setwd("/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus")
RNA-seq 数据的初步处理
# 加载RNA数据并查看特定列的分组统计信息rna = readRDS("etp.25.g1.downsample.1711.each.patient.rds")table(rna$comparison.bmp.vs.t.specified) # 统计BMP与T的比较信息table(rna$ETP) # 统计ETP的分组情况# 再次加载一个新的RNA-seq数据集，并检查其分组信息rna = readRDS("../GEO_Upload/Seurat/Full_CITE_seq_40_TALL.rds")table(rna$ETP) # 检查ETP分类table(rna$orig.ident) # 检查原始身份信息table(rna$is.blast.viscello) # 检查特定条件# 筛选数据集中满足`is.blast.viscello`为TRUE的细胞rna = subset(rna, subset = is.blast.viscello == TRUE)# 将身份标记为原始身份Idents(rna) = "orig.ident"# 随机下采样RNA数据，每组1146个细胞rna.subtype = subset(rna, downsample = 1146)table(rna.subtype$ETP) # 检查下采样后的ETP分类table(rna.subtype$orig.ident) # 检查下采样后的身份table(rna.subtype$is.blast.viscello) # 检查特定条件
ATAC-seq 数据的处理
# 加载ATAC-seq数据并检查其分组信息atac.blasts = readRDS("t.all.40.atac.blasts.1146.each.with.v4.anno.updated.RDS")table(atac.blasts$comparison.bmp.vs.t.specified) # 统计BMP与T的比较信息table(atac.blasts$ETP) # 统计ETP分类# 筛选数据中`comparison.bmp.vs.t.specified`列中不为“Other”的样本atac.blasts = subset(atac.blasts, comparison.bmp.vs.t.specified == "Other", invert = T)table(atac.blasts$comparison.bmp.vs.t.specified) # 重新检查分组table(atac.blasts$ETP) # 重新检查ETP分类# 绘制按分组着色的降维图DimPlot(atac.blasts, group.by = "comparison.bmp.vs.t.specified") + coord_fixed()
SCENIC+ 数据保存
RNA-seq 数据保存为 Loom 格式
# 设置默认的分析对象为RNA，并保存为Loom格式供SCENIC+使用DefaultAssay(rna) <- "RNA"rna.loom <- as.loom(rna, filename = "Data_SCENICplus/TALL_BMP_TSpec_RNA.loom", verbose = TRUE)rna.loom$close_all()DefaultAssay(rna.subtype) <- "RNA"rna.loom <- as.loom(rna.subtype, filename = "Data_SCENICplus/TALL_ETP_Subtype_RNA.loom", verbose = TRUE)rna.loom$close_all()
ATAC-seq 数据保存为稀疏矩阵格式
# 将ATAC-seq数据的计数矩阵转换为稀疏矩阵格式并保存counts_matrix <- as.matrix(atac.blasts@assays$ATAC@counts)counts_sparse <- Matrix::Matrix(counts_matrix, sparse = TRUE)writeMM(counts_sparse, file = "Data_SCENICplus/ATAC_Peaks_Sparse.mtx")# 保存细胞名称和区域名称cell_names <- colnames(counts_sparse)region_names <- rownames(counts_sparse)write.table(cell_names, file = "Data_SCENICplus/ATAC_Cell_Names.txt", col.names = FALSE, row.names = FALSE, quote = FALSE)write.table(region_names, file = "Data_SCENICplus/ATAC_Region_Names.txt", col.names = FALSE, row.names = FALSE, quote = FALSE)# 保存元数据metadata_frame <- atac.blasts@meta.datawrite.table(metadata_frame, file = "Data_SCENICplus/ATAC_Metadata.txt", sep = "\t", quote = FALSE, row.names = FALSE)
区域名称转换为 BED 格式
# 定义函数将区域名称转换为BED格式（染色体、起始位置、终止位置）convert_to_bed <- function(region) {parts <- unlist(strsplit(region, "-")) # 分割染色体和位置chr <- parts[1]start <- as.integer(parts[2])end <- as.integer(parts[3])bed <- c(chr, start, end)return(bed)}# 应用转换函数并保存为BED文件bed_data <- apply(as.matrix(region_names), 1, convert_to_bed)write.table(t(bed_data), "Data_SCENICplus/ATAC_Region_Names.bed", sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE)
上述代码就准备好了scenic+数据库构建和下游分析的文件，接下预处理fasta文件——Prepare_FASTA.sh。
# ------------------------------# 环境变量设置# ------------------------------# BED文件，包含ATAC-seq区域信息REGION_BED="/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/Data_SCENICplus/ATAC_Region_Names.bed"# 基因组FASTA文件GENOME_FASTA="/mnt/isilon/tan_lab/sussmanj/Single_Cell_Tools/ScenicPlus/Genome_Files/hg38.fa"# 染色体大小文件CHROMSIZES="/mnt/isilon/tan_lab/sussmanj/Single_Cell_Tools/ScenicPlus/Genome_Files/hg38.chrom.sizes"# 输出文件前缀DATABASE_PREFIX="T_ALL_40"# 脚本目录路径SCRIPT_DIR="/mnt/isilon/tan_lab/sussmanj/Single_Cell_Tools/ScenicPlus/create_cisTarget_databases"# ------------------------------# 运行create_fasta_with_padded_bg_from_bed.sh脚本# ------------------------------"${SCRIPT_DIR}/create_fasta_with_padded_bg_from_bed.sh" \${GENOME_FASTA} \ # 输入基因组FASTA文件${CHROMSIZES} \ # 输入染色体大小文件${REGION_BED} \ # 输入区域的BED文件hg38_T_ALL_1kb_bg_padding.fa \ # 输出带背景填充的FASTA文件名1000 \ # 背景填充大小（bp）yes # 是否允许基因组边界填充
构建数据库——Create_Cistarget_Motif.sh
# ------------------------------# 环境变量设置# ------------------------------# 单个调控元件（motif）文件的目录CBDIR="/mnt/isilon/tan_lab/sussmanj/Single_Cell_Tools/ScenicPlus/Motif_Collection/v10nr_clust_public/singletons"# 输入的FASTA文件，包含基因组序列及1kb背景填充FASTA_FILE="/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/Database_Files/hg38_T_ALL_1kb_bg_padding.fa"# motif列表文件，包含所有用于分析的motif名称MOTIF_LIST="/mnt/isilon/tan_lab/sussmanj/Single_Cell_Tools/ScenicPlus/Motif_Collection/v10nr_clust_public/motifs.txt"# 脚本路径，指向create_cistarget_motif_databases.py脚本SCRIPT_DIR="/mnt/isilon/tan_lab/sussmanj/Single_Cell_Tools/ScenicPlus/create_cisTarget_databases"# 输出数据库文件的前缀DATABASE_PREFIX="T_ALL"# ------------------------------# 运行create_cistarget_motif_databases.py脚本# ------------------------------"${SCRIPT_DIR}/create_cistarget_motif_databases.py" \-f ${FASTA_FILE} \ # 指定FASTA输入文件-M ${CBDIR} \ # 指定motif集合所在目录-m ${MOTIF_LIST} \ # 指定motif列表文件-o ${DATABASE_PREFIX} \ # 指定输出文件的前缀--bgpadding 1000 \ # 为背景区域提供的填充大小（以bp为单位）-t 40 # 指定使用的线程数
配置yaml文件，使用snakemake运行scenic+
input_data:# cisTopic对象文件，包含Cistopic分析模型cisTopic_obj_fname: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/Data_SCENICplus/ATAC_cistopic_obj_with_model.pkl"# GEX数据（基因表达数据）的anndata文件路径GEX_anndata_fname: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/Data_SCENICplus/TALL_BMP_TSpec_RNA.h5ad"# region_sets文件夹路径，包含基因组区域集合region_set_folder: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/Data_SCENICplus/region_sets"# ctx数据库文件，包含区域与转录因子（TFs）之间的关系ctx_db_fname: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/Database_Files/T_ALL.regions_vs_motifs.rankings.feather"# dem数据库文件，包含区域与转录因子（TFs）的富集得分dem_db_fname: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/Database_Files/T_ALL.regions_vs_motifs.scores.feather"# 转录因子注释文件路径path_to_motif_annotations: "/mnt/isilon/tan_lab/sussmanj/Single_Cell_Tools/ScenicPlus/Motif_Collection/v10nr_clust_public/snapshots/motifs-v10-nr.hgnc-m0.00001-o0.0.tbl"output_data:# 为prepare_GEX_ACC输出的combined GEX ACC数据文件路径（.h5mu格式）combined_GEX_ACC_mudata: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/ACC_GEX.h5mu"# motif富集结果输出文件路径（.hdf5格式）dem_result_fname: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/dem_results.hdf5"ctx_result_fname: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/ctx_results.hdf5"# motif富集结果的HTML输出文件路径output_fname_dem_html: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/dem_results.html"output_fname_ctx_html: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/ctx_results.html"# prepare_menr的输出文件路径（.h5ad格式）cistromes_direct: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/cistromes_direct.h5ad"cistromes_extended: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/cistromes_extended.h5ad"# 转录因子名称的输出文件路径（.txt格式）tf_names: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/tf_names.txt"# 基因组注释文件路径（.tsv格式）genome_annotation: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/genome_annotation.tsv"# 染色体大小的输出文件路径（.tsv格式）chromsizes: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/chromsizes.tsv"# 搜索空间的输出文件路径（.tsv格式）search_space: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/search_space.tsv"# 转录因子与基因关系的输出文件路径（.tsv格式）tf_to_gene_adjacencies: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/tf_to_gene_adj.tsv"# 区域与基因关系的输出文件路径（.tsv格式）region_to_gene_adjacencies: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/region_to_gene_adj.tsv"# 电子调控网络（eGRN）输出文件路径（.tsv格式）eRegulons_direct: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/eRegulon_direct.tsv"eRegulons_extended: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/eRegulons_extended.tsv"# AUCell分析结果的输出文件路径（.h5mu格式）AUCell_direct: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/AUCell_direct.h5mu"AUCell_extended: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/AUCell_extended.h5mu"# SCPlus分析结果的输出文件路径（.h5mu格式）scplus_mdata: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Output/scplusmdata.h5mu"params_general:# 临时目录路径temp_dir: "/mnt/isilon/tan_lab/sussmanj/Temp/ETP_ALL/SCENICPlus/SCENICPlus_Pipeline/Temporary"# 使用的CPU数量n_cpu: 72# 随机种子值seed: 69420params_data_preparation:# prepare_GEX_ACC的参数bc_transform_func: "\"lambda x: f'{x}-nonmultiome'\""is_multiome: False # 是否为多组学数据key_to_group_by: "comparison.bmp.vs.t.specified" # 用于分组的键nr_cells_per_metacells: 5 # 每个meta细胞的细胞数量# prepare_menr的参数direct_annotation: "Direct_annot"extended_annotation: "Orthology_annot"# 下载基因组注释的参数species: "hsapiens" # 物种biomart_host: "http://www.ensembl.org" # Biomart服务器地址# 搜索空间的参数search_space_upstream: "0 500000"search_space_downstream: "0 500000"search_space_extend_tss: "10 10" # 目标基因的TSS扩展范围params_motif_enrichment:# motif富集分析的物种和注释版本species: "homo_sapiens"annotation_version: "v10nr_clust"# motif相似性FDR阈值motif_similarity_fdr: 0.001orthologous_identity_threshold: 0.0 # 同源基因身份阈值annotations_to_use: "Direct_annot Orthology_annot" # 使用的注释类型# motif富集分析的相关参数fraction_overlap_w_dem_database: 0.4 # 在DEM数据库中与区域的重叠比例dem_max_bg_regions: 500 # DEM数据库中最多使用的背景区域数量dem_balance_number_of_promoters: True # 是否平衡启动子数量dem_promoter_space: 1_000 # 启动子的空间大小dem_adj_pval_thr: 0.05 # 调整后的p值阈值dem_log2fc_thr: 1.0 # log2倍数变化阈值dem_mean_fg_thr: 0.0 # 背景均值阈值dem_motif_hit_thr: 3.0 # motif命中阈值fraction_overlap_w_ctx_database: 0.4 # 在CTX数据库中与区域的重叠比例ctx_auc_threshold: 0.005 # AUC阈值ctx_nes_threshold: 3.0 # NES阈值ctx_min_regions: 10 # 最小区域数
input_data: 输入文件路径，包括染色质可及性数据、基因表达数据和基序数据库。
output_data: 输出文件路径，用于保存分析过程中的各类结果。
params_general: 通用参数，如临时目录、CPU 数量、随机种子等。
params_data_preparation: 输入数据处理和调控搜索空间的参数。
params_motif_enrichment: 基序富集分析的特定参数设置。
params_inference: 转录因子（TF）推断与 eRegulon 生成的参数.
