CATCH-UP：一种用于批量ATAC-Seq和ChIP-Seq数据的高通量上游分析流程

摘要

ATAC-seq 和 ChIP-seq 可用于对基因调控进行详细研究；然而，这些数据类型的处理具有挑战性，且不同研究团队之间通常存在不一致性。我们提出 CATCH-UP：一个易于使用的计算分析流程，可对新生成或已发表的 ATAC/ChIP-seq 数据集进行标准化、可重复的数据处理与分析。

摘要

转座酶可及性染色质检测（ATAC）与染色质免疫沉淀（ChIP）联合高通量测序（NGS）技术，已彻底改变了基因调控研究领域。由于这些技术所产生的高维数据集缺乏标准化分析流程，导致研究结果难以重复，进而造成已发表的处理后数据之间存在差异。这一问题的部分原因在于，针对此类数据的生物信息学分析工具种类繁多。其次，从原始数据生成完整且可解释的结果，需依次使用多种不同的生物信息学工具，而这些工具对计算技能的要求各不相同。此外，在数据处理过程中，质量控制环节存在多种可选方案，但尚未被统一采用。为解决上述问题，我们开发了完整的转座酶可及性染色质测序（ATAC-seq）与染色质免疫沉淀测序（ChIP-seq）上游分析流程（CATCH-UP），这是一个基于 Python、易于使用的分析流程，可将批量的 ChIP-seq 和 ATAC-seq 原始 fastq 文件处理为可可视化的 bigwig 轨迹图和峰值识别结果。该流程安装与运行简便，几乎无需计算背景知识。该流程具有模块化、可扩展性和并行化能力，适用于多种计算基础设施，便于方法学报告，从而实现对新产生或已发表数据集的可重复分析。

引言

基因表达必须受到严密调控，细胞才能建立并维持其正确的生物学功能。众所周知，异常的基因表达是多种疾病发病机制的基础，因此，大量研究关注于理解基因调控的机制¹。基因表达由启动子和增强子等调控元件促进。这些元件的序列中含有转录因子（TF）结合位点，当其处于活性状态时，可为转录因子的结合提供平台。转录因子在这些位点上的结合会导致核小体发生位移，从而增加DNA的可及性，并进一步提高转录机器的可进入性。由于这种可及性的增加，DNA的这些区域对DNase和Tn5等核酸酶和转座酶更为敏感，这一生化特性已被研究转录调控的研究人员广泛利用^2,3。

DNase-seq 和 ATAC-seq 使研究人员能够绘制全基因组范围内开放染色质区域、转录因子（TF）结合位点以及核小体定位图谱。在这两种技术中，由于 ATAC-seq 操作简便，仅需两步流程，且细胞用量少（每份样本仅需 50,000 个细胞，而 DNase-seq 每个重复需 100 万个细胞），因此在过去十年中其应用日益广泛。尽管 ATAC-seq 能够提供细胞群体中染色质整体景观的概览，但它通常无法确定具体是哪些蛋白质与基因组结合^4,5。为了鉴定特定蛋白质与基因组相互作用的具体位置，目前的金标准技术是染色质免疫沉淀（ChIP）-seq。ChIP-seq 首先通过化学方法固定细胞内的蛋白质-DNA 相互作用，然后使用针对目标蛋白的特异性抗体进行免疫沉淀（"pull-down"），从而富集与目标蛋白（POI）结合的 DNA 片段。这些 DNA 片段可通过测序揭示特定蛋白质（如转录因子）的基因组结合位点，或含有特定组蛋白修饰的位点¹。通过整合 ATAC-seq 与 ChIP-seq 数据集，可以构建出细胞群体中基因调控景观的详细图谱。

分析的基本工作流程如下：原始测序读段在比对到参考基因组（"mapping"）之前必须进行质量控制。成功比对的读段可进一步过滤，以去除低质量读段和PCR重复序列。为了可视化这些已比对并经过过滤的读段，需要计算这些读段在整个基因组上的"覆盖度"。由此生成的文件可作为"轨道"^6,7上传至基因组浏览器，例如多基因座视图（MLV）或UCSC基因组浏览器。对这些覆盖度轨道进行峰区识别，即"峰区 calling"，通常使用LanceOtron或MACS2^8,9等工具实现。最后，通过分析峰区的位置、形状和大小，可在样本或不同生物学条件之间进行比较。这些数据集的分析与整合是一个复杂的多步骤过程，可采用不同组合的生物信息学工具来完成。不同版本的工具之间可能存在不兼容性，从而影响数据处理的结果。此外，如nf-core¹⁰、panpipes¹¹、genpipes¹²、PEPATAC¹³或ChIP-AP¹⁴流程所示，数据处理的不同环节对计算资源和用户操作熟练程度的要求也存在较大差异。

总体而言，这导致了分析过程及分析结果报告的不一致，进而使得任何生物信息学知识有限的研究人员在研究的可重复性、可及性和便利性方面面临困难。我们通过CATCH-UP（完整的ATAC-seq和ChIP-seq上游分析流程）解决了上述所有问题，该流程是一个易于使用、灵活且模块化的ChIP-seq和ATAC/DNase-seq数据分析流程。CATCH-UP的实施仅需最少的生物信息学经验，可在多种计算基础设施上运行，并支持在不同研究团队内部及之间实现可重复的数据分析。

CATCH-UP 是一个基于 Python 的 Snakemake 分析流程，旨在标准化 ChIP-seq 和 ATAC-seq 数据的分析。该流程以原始测序数据（fastq.gz 文件）作为输入，输出各分析步骤结果，包括峰（.bed）文件。我们提供一个 YAML 格式的配置文件（config.yaml），用户可在其中编辑每个分析步骤的参数。Snakemake 内置的管理系统支持使用不同的计算基础设施（如服务器、集群、云系统或个人计算机），当用户提供大量数据时，还可实现并行处理。

以下，我们将详细描述工作流程的每个步骤（工作流程示意图见图1）。理解本说明对于后续遵循方案部分的逐步操作至关重要：

移动 fastq 文件： 流程的第一步是将原始 fastq 文件复制到指定的分析目录中。这样可以保留原始数据，避免对原始数据文件造成损坏或修改。

合并：如果原始测序数据包含多个测序通道（lane），则必须在分析前将这些通道的数据进行合并。默认情况下，分析流程将所有 fastq 文件视为单个样本，因此该合并步骤必须在配置文件中明确定义。

修剪：可选的数据清洗步骤。该步骤允许使用 trimmomatic¹⁵ 去除低质量的测序读段或接头序列。用户可提供自定义的接头序列 fasta 文件；接头目录中提供了示例文件。其他修剪参数可在配置文件中定义。默认情况下，工作流程跳过此步骤。

比对工具：默认使用 Bowtie2¹⁶ 进行序列比对；也可指定其他比对工具，例如 bwa-mem2¹⁷。选择 Bowtie2 作为默认比对工具，是因为它特别擅长将相对较短的测序读段比对到相对较大的基因组上，因此非常适合将 ChIP-seq 和 ATAC-seq 数据比对至哺乳动物基因组。为避免生成中间文件，比对过程通过管道直接传递至 samtools view，以将结果保存为输出的 bam 文件。在此规则中，用户必须指定用于比对读段的参考基因组版本，例如 hg19/hg38（人）、mm10/mm39（小鼠）。

过滤：保留正确比对的读段，过滤掉低质量的读段。默认使用 samtools view，参数为：-bShuF 4 -f 3 -q 30。

排序：将比对后的 reads 按照最左侧坐标顺序进行排序。默认：使用 samtools sort（snakemake 封装工具），参数为：-m 4G。

标记重复序列：识别并标记所有重复的测序读段。用户可通过修改配置文件参数决定是否将其去除。默认使用 Picard MarkDuplicates（snakemake 封装工具），参数设置为：--REMOVE_DUPLICATES False，以标记并保留重复读段。

合并 bam：如果测序数据包含重复样本或多个样本，用户可能希望将其合并为单个 bam 文件。在此情况下，用户可选择合并 bam 文件，或在分析过程中始终保持 bam 文件相互独立。若用户选择合并 bam 文件（使用 samtools merge 工具），则必须为合并后的 bam 文件指定一个共同的前缀。

索引：此步骤对已排序的坐标进行索引。默认：使用 samtools index（snakemake 封装器），采用 samtools 指定的默认参数。

BamCoverage：此规则通过比对后的测序读段生成 bigwig 覆盖度轨迹文件。该流程使用 deepTools 中的 bamCoverage 工具，覆盖度的计算方式为每个 bin 中的读段数量，其中 bin 表示指定大小的滑动窗口。在此分析流程中，bamCoverage 的默认参数设置如下：-bs 1 -normalizeUsing RPKM -extendReads。

峰检测（peak calling）：本分析流程默认选用 LanceOtron⁸ 作为峰检测工具。与传统主要基于统计检验的峰检测工具不同，LanceOtron 是一种基于深度学习的峰检测方法，它结合了基因组富集度测量与统计检验，已被证明性能优于行业标准工具 MACS2⁹。为了使 bigWig 文件与 LanceOtron 兼容，覆盖度必须以每个碱基对为单位进行计算，并进行 RPKM 标准化；这一点已在 BamCoverage 步骤的默认设置中体现。MACS2 可作为替代方案被选为峰检测工具。我们将持续关注新型峰检测工具的发布，并在适用情况下将其纳入本流程，以维持并优化该分析流程的性能。

TrackDb：该功能可创建 bigwig 文件的键值对关联，以便在 MLV⁶ 或 UCSC 基因组浏览器¹⁸ 等工具中加载和可视化这些文件。

除了输出数据外，分析流程的每个步骤都会输出一个日志文件，并提供适当的质量控制检查，以便用户跟踪分析进度。FastQC¹⁹ 应用于原始测序数据和经修剪的测序数据（如已选择）（步骤 1 - 快速移动 fastq 和 2- 修剪). Samtools stats 与 MultiQC 结合使用²⁰ 用于在第3步的输出中收集、生成和可视化bam文件的质量控制报告 矫治器, 6 - 标记重复序列，以及 7 - 合并 bam 文件. 有关上述步骤中所用每种工具的更多信息，请参见 表1.

方案

1. 运行 CATCH-UP 分析流程

1. 从 https://github.com/Genome-Function-Initiative-Oxford/UpStreamPipeline 克隆 UpStreamPipeline 仓库：
   导航至选定的工作目录，复制以下代码并在命令行中运行：
   git clone git@github.com:Genome-Function-Initiative-Oxford/UpStreamPipeline.git
2. 使用命令 cd UpStreamPipeline 进入下载的 UpStreamPipeline 文件夹。
3. 安装 anaconda 发行版（如尚未安装）：
   1. 使用命令 which conda 检查系统是否已安装 anaconda。如果该命令未显示任何 conda 发行版的路径，请从 https://github.com/conda-forge/miniforge#mambaforge 下载 mambaforge，并为系统选择正确的发行版和版本。例如，Linux 用户可使用 wget  https://github.com/conda-forge/miniforge/releases/latest/download/Mambaforge-Linux-x86_64.sh。访问此网页以获取不同操作系统的版本：https://github.com/conda-forge/miniforge/。
   2. 使用 sh Mambaforge-Linux-x86_64.sh 运行安装程序，并通过执行 conda init 在系统中初始化 conda。
4. 安装并激活 upstream conda 环境（upstream conda 环境的依赖项列于 表 2 中）：
   1. 使用命令 mamba env create -f envs/upstream.yml 安装环境。
   2. 使用命令 conda activate upstream 激活环境。
5. 成功安装 upstream conda 环境后，使用命令 conda activate upstream 激活环境，并使用 cd genetics/CATCH-UP 导航至 CATCH-UP 文件夹。
6. 编辑配置文件，该文件位于 config 文件夹中，使用命令 cd /config/analysis.yaml 打开，并使用文本编辑器根据分析需求修改参数。请按照文件内的逐行说明编辑各项参数。此文件将在分析完成后保留，用于记录运行参数，以促进结果的可重复性。
7. 使用文本编辑器（例如 Mac 上的 TextEdit 或 Windows 上的 Notepad）打开并编辑以下三个文件：
   1. 编辑 1_fastqfile_home_dir.txt 文件，使其包含所有待分析 fastq 文件的列表。
      注意：读段编号和扩展名（例如 _R1/_R2 和 .fastq.gz）必须排除。例如，若某项目包含以下 fastq 文件列表：
      Sample1_conditionA_L001_R1 . fastq . gz
      Sample1_conditionA_L001_R2 . fastq . gz
      Sample1_conditionA_L002_R1 . fastq . gz
      Sample1_conditionA_L002_R2 . fastq . gz
      Sample1_conditionB_L001_R1 . fastq . gz
      Sample1_conditionB_L001_R2 . fastq . gz
      Sample1_conditionB_L002_R1 . fastq . gz
      Sample1_conditionB_L002_R2 . fastq . gz
      Sample2_conditionA_L001_R1 . fastq . gz
      Sample2_conditionA_L001_R2 . fastq . gz
      Sample2_conditionA_L002_R1 . fastq . gz
      Sample2_conditionA_L002_R2 . fastq . gz
      Sample2_conditionB_L001_R1 . fastq . gz
      Sample2_conditionB_L001_R2 . fastq . gz
      Sample2_conditionB_L002_R1 . fastq . gz
      Sample2_conditionB_L002_R2 . fastq . gz
      在此情况下，1_fastqfile_home_dir.txt  文件内容应如下：
      Sample1_conditionA_L001
      Sample1_conditionA_L002
      Sample1_conditionB_L001
      Sample1_conditionB_L002
      Sample2_conditionA_L001
      Sample2_conditionA_L002
      Sample2_conditionB_L001
      Sample2_conditionB_L002
   2. 如果原始数据包含需要拼接的测序通道，请编辑 2_fastqfile_concat.txt 文件，定义需拼接文件名的前缀。若无需要拼接的测序通道，则无需编辑 2_fastqfile_concat.txt。确保 2_fastqfile_concat.txt 的每一行仅包含一个样本前缀，如下所示：
      Sample1_conditionA
      Sample1_conditionB
      Sample2_conditionA
      Sample2_conditionB
   3. 如果需要合并不同样本的数据，请编辑 3_merge_bams.txt 文件，填入待合并文件名的前缀。确保每一行仅包含一个样本前缀，如下所示：
      Sample1
      Sample2
      图 2 总结了如何整理这三个文件。本方案适用于单端或双端测序数据。除非另有指定，流程默认执行双端分析；此设置可在配置文件中修改（参见步骤 1.6）。
8. 完成所有必要文件的编辑后，使用 snakemake 按如下方式运行 CATCH-UP：snakemake --configfile=config/analysis_name.yaml all --cores 4。
   注意：如需更详细的说明和文档，请参阅 UpStreamPipeline GitHub 仓库中的 CATCH-UP 文件夹，链接为 此处。其中包含有关如何正确修改配置文件的详细说明，包括更改测序数据文件路径、结果存储路径以及各步骤参数的编辑方法。

结果

CATCH-UP 分析流程为每个步骤生成结果文件、日志文件以及质量控制（QC）输出。在配置文件中，用户可选择保留或删除输出文件，以减少所需的存储空间。所有输出内容说明如下：

00. fastq_home_dir：配置文件 fastqfile_home_dir.txt 和 merge_bams.txt 被复制到此文件夹中，以供参考和结果复现。
01. reads：为避免在工作流程中修改原始测序数据，fastq 文件被复制到此文件夹；若指定，可将不同测序通道（lanes）的数据合并。
02. trimming：若指定，对 fastq 文件中的测序读段（reads）和接头序列进行修剪。
03. aligner：将测序数据比对至选定的参考基因组。
04. filtering：进行质量控制过滤。
05. sorted：对 bam 文件进行排序。
06. duplicates：标记重复序列。
07. merge：若在 config.yaml 文件中指定，则合并 bam 文件。
08. bam_coverages：生成覆盖度的 bigwig 文件。
09. peak_calling：LanceOtron 峰值识别结果的 bed 文件。
10. track：生成格式化的文本文件，可根据需要直接用于基因组浏览器。

对于 01、02、03、06 和 07 的输出，提供了质量控制指标和 HTML 文件。此外，在图3中，我们展示了使用 CATCH-UP 处理数据的示例，并通过 MLV 平台对最终输出结果进行可视化。

图1：CATCH-UP工作流程。 给定一组fastq文件，CATCH-UP将并行处理所有样本的上游分析步骤。请点击此处查看该图的放大版本。

图2：图示说明如何正确修改 1_fastqfile_home_dir.txt、2_fastqfile_concat.txt 和 3_merge_bams.txt 文件，以运行 CATCH-UP。请点击此处查看该图的放大版本。

图 3：CATCH-UP 分析流程的示例输出结果。 原始测序数据（fastq 文件）从 ENCODE²¹ 下载。使用 CATCH-UP 分析流程对 DNase-seq 和 5 种类型的 ChIP-seq（H3K4me1、H3K4me3、H3K27ac、CTCF 和 POLR2A）的 fastq 文件进行处理。将生成的 bigwig 输出文件上传至 Multi Locus View，用于可视化和识别基因组调控元件。请点击此处查看该图的放大版本。

表1：文档资源。 本表列出了CATCH-UP工作流程中涉及的工具、其文档链接及相应参考文献。请点击此处下载该表格。

表2：上游conda环境的通道和依赖项要求列表。 请点击此处下载该表格。

表3：用于测试CATCH-UP的操作系统。 Ubuntu在高性能计算集群和本地计算机上进行了测试。请点击此处下载该表格。

讨论

随着高通量测序（NGS）技术在基因组数据生成中的应用日益广泛，用于分析这些数据的生物信息学工具也不断增多。在数据分析的每个步骤中，均可采用多种不同的工具，且每种工具内部还可设置大量不同的参数^{6,8,9,15,16,17,18,19,20,22,23,24}。这导致了数据分析策略组合的极大多样性，而每种组合都可能产生不同的分析结果。为了在不同实验之间进行准确比较，生物信息学分析的标准化至关重要。传统上，NGS数据由湿实验科学家生成，再由生物信息学专家进行数据分析。

下一代测序（NGS）数据分析可分为"上游"和"下游"分析流程，其中上游流程包括将测序仪输出的原始数据转换为研究人员可直观解读格式所必需的步骤。下游分析则包含针对具体研究问题和实验设计而定制的额外步骤。因此，上游流程具有通用性，适合标准化，以提高科学研究的可重复性；而下游流程则因依赖于具体的生物学问题，并需要研究者的专业判断，具有较强的定制性，因而较难实现标准化。我们开发了一个用户友好的上游分析流程，使湿实验研究人员无需具备生物信息学背景知识即可可重复地分析自身数据。本文介绍了CATCH-UP，该流程基于snakemake框架构建，旨在兼顾易用性，并解决ChIP-seq和ATAC-seq数据分析中的可重复性问题。该流程可处理ChIP-seq或ATAC-seq数据。用户在下载CATCH-UP后，需首先定义分析参数和样本命名，然后通过命令行使用一行代码运行流程。如何为ChIP-seq或ATAC-seq分析自定义参数的详细分步说明，已内置于配置文件中，并可在CATCH-UP的GitHub仓库提供的分步指南中找到。

目前已有用于 ChIP-seq 或 ATAC-seq 数据的分析流程，例如 PEPATAC 和 ChIP-AP。尽管这些流程具有一些优势，例如将上游和下游分析整合到单一工作流中，或支持图形用户界面（GUI），但这些工具主要面向具备一定计算技能的生物信息学研究人员和科学家^13,14。CATCH-UP 的设计旨在解决两个问题：一是使没有接受过生物信息学培训的湿实验科学家能够独立完成上游分析，二是通过便于报告和实现精确复现，促进实验室之间上游分析的标准化。CATCH-UP 故意仅限于上游分析，但其输出结果可与下游分析工具兼容，例如用于数据集的统计比较或推断转录因子结合的工具^25,26。

执行可重复的上游分析所需的所有关键步骤均已预先定义在 CATCH-UP 流程中，以确保分析的稳健性。该流程的详细输出特性使用户能够逐步跟踪其运行过程，这既有利于故障排查，也有助于实现分析流程的重复验证。鉴于高通量测序（NGS）技术的快速发展，该流程的模块化设计具有显著优势，可轻松适应各类工具版本更新以及新工具的引入与应用。CATCH-UP 已在以下操作系统中成功测试：Ubuntu、CentOS、macOS（Intel CPU）和 Windows（表 3）。该流程通过并行化处理工作流，能够应对包含数十个样本的大规模实验，因而可灵活适配不同的实验设计。总体而言，在 ChIP-seq 和 ATAC-seq 数据分析中应用 CATCH-UP，可实现用户友好、可重复且高度灵活的分析流程。

材料

本文使用的材料清单

姓名	公司	目录编号	评论
CATCH-UP	GitHub		https://github.com/Genome-Function-Initiative-Oxford/UpStreamPipeline/tree/main/genetics/CATCH-UP
CentOS	Linux	Version 7	此处列出的任何操作系统均可使用
macOS	Apple	Version 13 Ventura	此处列出的任何操作系统均可使用
Ubuntu	Ubuntu	Version 22.04 LTS	此处列出的任何操作系统均可使用
Windows	Microsoft	Version 11	此处列出的任何操作系统均可使用