{ "cells": [ { "cell_type": "markdown", "id": "fb7e7642", "metadata": {}, "source": [ "---\n", "title: scATAC-seq 拟时序分析：染色质动态轨迹解析\n", "author: SeekGene\n", "date: 2026-01-29\n", "tags:\n", " - ATAC + RNA 双组学\n", " - Notebooks\n", " - 分析指南\n", " - 拟时序分析\n", "---\n" ] }, { "cell_type": "markdown", "id": "aca702d2-45f2-420e-83d1-ac0361d75f6a", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# scATAC-seq 拟时序分析：染色质动态轨迹解析" ] }, { "cell_type": "markdown", "id": "770f4be6", "metadata": {}, "source": [ "## 背景介绍\n", "\n", "ATAC_Monocle3 分析流程专为单细胞 ATAC-seq 数据的轨迹推断与发育动力学研究设计。该流程基于 Signac/Seurat 的 ATAC 数据预处理，结合 Monocle3 的伪时序分析能力，能够揭示细胞在发育、分化或应答过程中的染色质可及性变化轨迹。通过构建细胞状态转变路径，ATAC_Monocle3 有助于解析调控元件动态开放、关键转录因子活性变化及其对基因调控网络的影响。\n", "\n", "### 核心原理\n", "\n", "- 轨迹推断（Trajectory Inference）：利用 Monocle3 对高维 ATAC-seq 数据进行降维与聚类，重建细胞状态转变的伪时序路径。\n", "- 动态可及性分析：沿伪时序分析染色质开放区域的动态变化，识别与细胞命运决定相关的调控元件。\n", "- 调控网络解析：结合基因活力与 motif 富集，推断关键转录因子及其调控网络在细胞发育过程中的作用。\n", "\n", "### 目的与意义\n", "- 揭示细胞分化/发育过程中的染色质重塑规律\n", "- 识别驱动细胞命运转变的关键调控元件与转录因子\n", "- 为疾病机制、发育生物学等领域提供动态表观遗传学证据" ] }, { "cell_type": "markdown", "id": "51b4b2f1", "metadata": {}, "source": [ "## 参数设置\n", "\n", "请先挂载将要分析的云平台相关流程数据，挂载方式请参照 jupyter 使用教程\n", "\n", "具体参数填写：\n", "\n", "- rds：挂载的流程相关的 rds 数据，在/home/{UserName}/workspace/project/{UserName}/目录下，如下列项目数据/home/shumeng/workspace/data/AY1752167131383/\n", "- meta：与 rds 同目录下的 metadata 文件 \n", "- outdir: 分析结果保存路径\n", "- clusters_col：细胞类型列名，选择要分析的细胞类型或者聚类对应的标签。假如想对注释好的细胞类型进行分析，则选择其对应的标签，例如 CellAnnotation，与<细胞类型>配合使用。\n", "- celltypes：细胞类型，多选，选择要分析的细胞类型或者聚类结果，如 Monocyte 和 Macrophage 等。\n", "- root：做为发育起始位点的 cluster 或 celltype。\n", "- SAMple_col：样本分类标签，如 SAMple，是 meta.data 的列名\n", "- SAMples：用于分析的样本名称\n", "- use_partition：use_partition 是传给 Monocle3 learn_graph 的参数，控制是否按分区（partition）分别学习轨迹图。use_partition = TRUE（推荐，默认）：在每个分区内各自拟合轨迹，不跨分区连线，避免把本不相连的细胞硬连起来。\n", "- downSAMple：是否下采样\n", "- downSAMple_num：如果下采样，每个 cluster 或 celltype 下采样多少细胞" ] }, { "cell_type": "code", "execution_count": null, "id": "fa725331", "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "# Parameters\n", "outdir = \"/PROJ2/FLOAT/advanced-analysis/prod/17801-49246-86550-580868389543631482/_build/html/\"\n", "rds = \"/PROJ2/FLOAT/advanced-analysis/prod/17801-49246-86550-580868389543631482/input.rds\"\n", "meta = \"/PROJ2/FLOAT/advanced-analysis/prod/17801-49246-86550-580868389543631482/meta.txt\"\n", "root = \"9\"\n", "clusters_col = \"resolution.0.5_d30\"\n", "celltypes = \"9,2\"\n", "sample_col = \"Sample\"\n", "samples = \"NSC_0722_arc,TSC_0722_arc\"\n", "use_partition = \"TRUE\"\n", "downsample = \"TRUE\"\n", "downsample_num = \"1000\"\n" ] }, { "cell_type": "markdown", "id": "45228bef", "metadata": {}, "source": [ "## 环境准备\n", "\n", "**R 包加载**\n", "\n", "请选择 common_r 这个环境进行该整合教程的学习" ] }, { "cell_type": "code", "execution_count": 2, "id": "2dd37f51-7539-4df9-8968-7051a2961983", "metadata": { "editable": true, "execution": { "iopub.execute_input": "2025-08-21T10:15:40.449467Z", "iopub.status.busy": "2025-08-21T10:15:40.446953Z", "iopub.status.idle": "2025-08-21T10:15:40.659275Z", "shell.execute_reply": "2025-08-21T10:15:40.657077Z" }, "slideshow": { "slide_type": "" }, "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "suppressPackageStartupMessages({\n", " library(Signac)\n", " library(Seurat)\n", " library(SeuratWrappers)\n", " library(monocle3)\n", " #library(cicero)\n", " library(Matrix)\n", " library(ggplot2)\n", " library(patchwork)\n", " library(parallel)\n", "})\n" ] }, { "cell_type": "markdown", "id": "54213216", "metadata": {}, "source": [ "## Monocle3 轨迹分析\n", "\n", "\n", "Monocle3[1] 是由华盛顿大学 Cole Trapnell 实验室开发的单细胞组学分析框架，专用于重建细胞状态转变轨迹并计算伪时间（pseudotime）。在人为指定根细胞（root）的前提下，Monocle3 基于学习到的轨迹图（principal graph），以根细胞为起点计算到各细胞的测地距离，从而定义伪时间这一连续标度，用于刻画细胞在发育、分化或应答过程中的进程位置。借此可以将离散的细胞状态进行排序，重建连续的生物学变化过程。在 scATAC-seq 背景下，伪时间反映的是染色质可及性模式的渐进变化，有助于揭示关键调控元件与转录因子的时序性活动以及基因调控网络的动态重塑。" ] }, { "cell_type": "code", "execution_count": null, "id": "bda7628a-008b-47ee-90b0-6a77f0fa46fe", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "outdir <- paste0(outdir,'/output')\n", "dir.create(outdir, recursive = TRUE)\n", "\n", "celltypes <- strsplit(celltypes,\",\")[[1]]\n", "samples <- strsplit(samples,\",\")[[1]]\n", "\n", "use_partition <- as.logical(use_partition)\n", "downsample <- as.logical(downsample)\n", "downsample_num <- as.numeric(downsample_num)\n", "\n" ] }, { "cell_type": "markdown", "id": "cd15069c", "metadata": {}, "source": [ "### 数据读取\n", "\n", "读取输入对象与元数据，完成条形码对齐、细胞/样本筛选与默认 Assay 设置，为后续 Monocle3 轨迹分析准备数据。\n", "- 读取 RDS 与导入 metadata，并按条形码对齐\n", "- 按细胞类型与样本筛选对象\n", "- 设置默认 Assay 为 ATAC，可选下采样以控制规模" ] }, { "cell_type": "code", "execution_count": null, "id": "52cede64-3a7d-4c2c-be04-b87b43ae7898", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "obj <- readRDS(rds)\n", "if (meta != \"\"){\n", " meta <- read.table(meta,header=T,sep='\\t',check.names=F)\n", " rownames(meta) <- meta$barcodes\n", " obj <- AddMetaData(obj, meta)\n", "}\n", "obj@meta.data[[clusters_col]] <- as.character(obj@meta.data[[clusters_col]])\n", "obj <- subset(obj, subset = !!sym(clusters_col) %in% celltypes)\n", "obj <- subset(obj, subset = !!sym(sample_col) %in% samples)\n", "table(obj@meta.data[[clusters_col]])\n", "\n", "DefaultAssay(obj) <- \"ATAC\"\n", "Idents(obj) <- obj@meta.data[[clusters_col]]\n", "if (downsample){\n", " obj <-subset(obj,cells=WhichCells(obj,downsample=downsample_num))\n", "}\n", "table(obj@meta.data[[clusters_col]])\n" ] }, { "cell_type": "markdown", "id": "1a3e8de5", "metadata": {}, "source": [ "### CDS 构建与伪时间推断\n", "本节将构建适用于 Monocle3 的 `cell_data_set`，并完成轨迹学习与伪时间推断。请结合数据连通性选择是否启用 `use_partition`，在确定生物学起点后设置 `root` 以确定伪时间方向；当细胞数较大时，可启用 `downSAMple` 以提升运行效率而不改变整体趋势。输出包括学习到的轨迹图、细胞伪时间以及关键分支的可视化结果。```" ] }, { "cell_type": "code", "execution_count": 73, "id": "1e287255-f18a-4ffc-8ee8-00a729f7be42", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "\n", "\n", "\n", " |======================================================================| 100%" ] } ], "source": [ "cds <- as.cell_data_set(obj)\n", "cds <- detect_genes(cds)\n", "cds <- cluster_cells(cds = cds, reduction_method = \"UMAP\")\n", "cds <- learn_graph(cds, use_partition = use_partition)\n", "root_cells=rownames(obj@meta.data[obj@meta.data[[clusters_col]]==root,])\n", "cds <- order_cells(cds, reduction_method = \"UMAP\", root_cells = root_cells)\n", "saveRDS(cds, file.path(outdir, 'cds.rds'))" ] }, { "cell_type": "markdown", "id": "7bcaaca9", "metadata": {}, "source": [ "### Monocle3 结果可视化\n", "\n", "展示并解读轨迹推断的核心可视化，帮助理解细胞状态转换的时序与分支结构。\n", "\n", "- 轨迹图（UMAP + principal graph）：\n", " 在 UMAP 上叠加学习到的轨迹主干，观察细胞在低维空间中的连续过渡与分支走向。\n", "- 伪时间分布（Pseudotime on UMAP）：\n", " 以渐变色显示细胞伪时间，评估从根细胞出发的进程方向是否与生物学预期一致。\n", "- 特征趋势：\n", " 沿伪时间展示峰/基因活力或 motif 活性趋势，识别阶段性调控特征。\n", "\n", "解读要点：\n", "- 连续性与方向性是否合理（受 `root` 选择影响）。\n", "- 分支是否对应已知的细胞亚群或状态转变。\n", "- 分区之间应无硬连线（`use_partition=TRUE`），跨分区连接需谨慎解读。\n", "- 下采样仅改变点密度，不应改变整体轨迹形态与趋势。" ] }, { "cell_type": "code", "execution_count": null, "id": "6d1120ed-b9aa-40be-90db-cca964def890", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "p1=plot_cells(\n", " cds = cds,\n", " color_cells_by = \"pseudotime\",\n", " label_cell_groups = F,\n", " show_trajectory_graph = T\n", ")\n", "\n", "p2=plot_cells(\n", " cds = cds,\n", " color_cells_by = \"pseudotime\",\n", " label_cell_groups = F,\n", " show_trajectory_graph = F\n", ")\n", "p3=plot_cells(\n", " cds = cds,\n", " color_cells_by = clusters_col,\n", " show_trajectory_graph = T,\n", " label_cell_groups = F\n", ")\n", "p4=plot_cells(\n", " cds = cds,\n", " color_cells_by = clusters_col,\n", " show_trajectory_graph = F,\n", " label_cell_groups = F\n", ")\n", "ggsave(p1+p2, file = file.path(outdir, 'pseudotime.png'), width = 12, height = 5, dpi=300)\n", "ggsave(p3+p4, file = file.path(outdir, 'celltype.png'), width = 12, height = 5, dpi=300)\n", "ggsave(p1+p2, file = file.path(outdir, 'pseudotime.pdf'), width = 12, height = 5)\n", "ggsave(p3+p4, file = file.path(outdir, 'celltype.pdf'), width = 12, height = 5)" ] }, { "cell_type": "code", "execution_count": null, "id": "db5ffe9c-399f-48ef-ab1a-0b052b85c6f5", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "input_cds_lin <- cds[,is.finite(pseudotime(cds))]\n", "pData(input_cds_lin)$Pseudotime <- pseudotime(input_cds_lin)\n", "pData(input_cds_lin)$cell_subtype <- cut(pseudotime(input_cds_lin), 10)\n", "binned_input_lin <- aggregate_by_cell_bin(input_cds_lin, \"cell_subtype\")\n", "acc_fits <- fit_models(binned_input_lin, cores = ifelse(detectCores() > 32, 25, trunc(0.8*detectCores())), \n", " model_formula_str = \"~Pseudotime + num_genes_expressed\" )\n", "fit_coefs <- coefficient_table(acc_fits)\n", "\n", "significant_sites <- subset(fit_coefs, term == \"Pseudotime\" & q_value < .05)\n" ] }, { "cell_type": "code", "execution_count": 55, "id": "6905220b-3a77-4e85-92ef-005fca61df75", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "important_cols <- c(\"site_name\", \"estimate\", \"normalized_effect\", \"p_value\", \"q_value\")\n", "important_data <- significant_sites[, important_cols]\n", "\n", "write.csv(important_data, file.path(outdir, \"pseudotime_differentially_accessible_sites.csv\"), row.names = FALSE)" ] }, { "cell_type": "code", "execution_count": null, "id": "530cc7a0-46e0-4a7e-9436-83609cf797e0", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "top_increasing <- significant_sites[order(-significant_sites$estimate), ][1:5, ] # 随时间增加最显著\n", "top_decreasing <- significant_sites[order(significant_sites$estimate), ][1:5, ] # 随时间减少最显著\n", "top_increasing <- top_increasing[!is.na(top_increasing$site_name), ]\n", "top_decreasing <- top_decreasing[!is.na(top_decreasing$site_name), ]" ] }, { "cell_type": "code", "execution_count": 72, "id": "3164b8ec-b391-43f1-881d-f30f7245c1d6", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "p5=plot_accessibility_in_pseudotime(input_cds_lin[top_increasing$site_name])\n", "p6=plot_accessibility_in_pseudotime(input_cds_lin[top_decreasing$site_name])\n", "ggsave(p5, file = file.path(outdir, 'top_increasing.png'), width = 5, height = 8, dpi=300)\n", "ggsave(p6, file = file.path(outdir, 'top_decreasing.png'), width = 5, height = 8, dpi=300)\n", "ggsave(p5, file = file.path(outdir, 'top_increasing.pdf'), width = 5, height = 8)\n", "ggsave(p6, file = file.path(outdir, 'top_decreasing.pdf'), width = 5, height = 8)" ] }, { "cell_type": "markdown", "id": "6dd553ed", "metadata": {}, "source": [ "\n", "**4.3.1 拟时序轨迹图解读**\n", "\n", "- 每个点代表一个细胞；黑色线段为学习到的轨迹主干与分支路径。\n", "- 颜色越暗表示越接近发育起始端，颜色越亮表示越接近发育末端（伪时间由深→浅递增）。\n", "- 白底圆圈中的编号为根节点（root）的标记。\n", "- 不与根节点处于同一 partition 的细胞不计算伪时间，统一以灰色显示，通常为离群或不连通成分。" ] }, { "cell_type": "code", "execution_count": null, "id": "153d1bbf-196a-457b-b0e2-96c25c971a33", "metadata": {}, "outputs": [], "source": [ "options(repr.plot.width = 20, repr.plot.height = 8)\n", "print(p1+p2)" ] }, { "cell_type": "markdown", "id": "f94d4783", "metadata": {}, "source": [ "\n", "**4.3.2 细胞类型轨迹图解读**\n", "\n", "- **点与颜色**：每个点代表一个细胞；不同颜色对应不同细胞群/注释类别。\n", "- **轨迹主干**：黑色线为学习到的轨迹主干与分支路径，指示状态转换走向。\n", "- **分支节点**：黑底圆圈数字表示分支节点，连接不同的分化命运（类比枝丫）。\n", "- **终末状态**：灰底圆圈表示终末分化状态，是细胞命运的最终结果（类比树叶）。\n", "- **编号说明**：节点编号为随机标签，无先后含义；需结合伪时间（pseudotime）判定分化方向。" ] }, { "cell_type": "code", "execution_count": null, "id": "6ff6e266-ca88-4d70-8a6d-0ccd38e24504", "metadata": {}, "outputs": [], "source": [ "options(repr.plot.width = 20, repr.plot.height = 8)\n", "print(p3+p4)" ] }, { "cell_type": "markdown", "id": "11be8ab9", "metadata": {}, "source": [ "\n", "**4.3.3 可及性随伪时间变化趋势解读**\n", "\n", "本节展示峰/基因活力在伪时间轴上的动态变化，用于识别早期、晚期与瞬时调控模块。\n", "- 方法：按伪时间排序细胞，对所选特征进行平滑（如 LOESS），绘制热图/均值曲线。\n", "- 结果解读：上调→晚期开放；下调→早期开放；双相→阶段性调控；分支特异曲线指示命运分化差异。\n", "- 实用建议：特征按细胞内 z-score 标准化；对趋势聚类成模块并做 GO/TF 富集；跨分区请分别解读。" ] }, { "cell_type": "code", "execution_count": null, "id": "cc0e60c1-df05-452f-b536-9f5b9d72860f", "metadata": {}, "outputs": [], "source": [ "options(repr.plot.width = 20, repr.plot.height = 8)\n", "print(p5+p6)" ] }, { "cell_type": "markdown", "id": "60694d7d-d5d1-49b5-946a-e0571778f28a", "metadata": { "editable": true, "jp-MarkdownHeadingCollapsed": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## 结果文件\n", "\n", "- pseudotime.pdf/png - 拟时序轨迹图\n", "- celltype.pdf/png - 细胞类型轨迹图\n", "- top_increasing.pdf/png - 随时间可及性增加最显著的五个位点\n", "- top_decreasing.pdf/png - 随时间可及性减少最显著的五个位点\n", "\n", "结果目录：目录下包括此项分析所涉及的所有图片以及表格等结果文件\n", "" ] }, { "cell_type": "markdown", "id": "5a536d97-1ca8-4d46-bf92-75526e415226", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## 文献案例解析\n", "\n", "* **文献一：** \n", "文献《Single cell transcriptional and chromatin accessibility profiling redefine cellular heterogeneity in the adult human kidney》[2]中研究者使用 Monocle3 对近端小管(PT)细胞亚群进行了深入分析，发现了一个表达 VCAM1 的特殊 PT 亚群，这个亚群在肾脏损伤后具有再生潜能。通过伪时间分析，研究者追踪了从常规 PT 细胞向 VCAM1+PT 细胞转变的动态过程，并识别出 NF-κB 信号通路在这一转变过程中的关键作用。此外，他们发现了一个与 VCAM1 基因相关的共可及区域，这个区域富集了 NF-κB 结合位点，证实了 NF-κB 在调控 VCAM1 表达中的作用。\n", "" ] }, { "cell_type": "markdown", "id": "0e622a5c-6196-4f95-bb82-f298661eb656", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## 参考资料\n", "[1] Trapnell C, Cacchiarelli D, Grimsby J, Pokharel P, Li S, Morse M, Lennon NJ, Livak KJ, Mikkelsen TS, Rinn J L The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014 Apr;32(4):381-386. doi: 10.1038/nbt.2859. Epub 2014 Mar 23. PMID: 24658644; PMCID: PMC4122333.\n", "\n", "[2] Muto Y, Wilson P C, Ledru N et al. Single cell transcriptional and chromatin accessibility profiling redefine cellular heterogeneity in the adult human kidney[J]. Nat Commun, 2021, 12: 2190.\n" ] }, { "cell_type": "markdown", "id": "059455e6-b826-4582-b83d-aa42d6a4594c", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## 附录\n", "\n", "- .csv/.xls/.txt ：结果数据表格文件，文件以逗号或制表符（Tab）分隔。Linux/Mac 用户使用 less 或 more 命令查看；Windows 用户使用高级文本编辑器 Notepad++ 等查看，也可以用 Microsoft Excel 打开。\n", "\n", "- .png：结果图像文件，位图，无损压缩。\n", "\n", "- .pdf：结果图像文件，矢量图，可以放大和缩小而不失真，方便用户查看和编辑处理，可使用 Adobe Illustrator 进行图片编辑，用于文章发表等。" ] }, { "cell_type": "code", "execution_count": null, "id": "eead1d28-722f-448e-af69-e4767713ee22", "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "common_r", "language": "R", "name": "common_r" }, "language_info": { "codemirror_mode": "r", "file_extension": ".r", "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", "version": "4.3.3" } }, "nbformat": 4, "nbformat_minor": 5 }