摘要: SEACR（基于稀疏富集分析的染色质峰识别工具）是专为CUT&RUN等染色质分析数据设计的峰值检测软件，尤其适用于背景信号以大量零覆盖区域为主的稀疏数据。该工具依赖R语言和Bedtools，以双端测序生成的bedgraph文件作为输入。

在方法层面，SEACR通过分析信号块的总信号量与组成行数，区分真实富集区域与噪声。其核心改进包括：引入基于分位数的自适应密度截断阈值（取代原先硬编码的90%截止值），通过最大化正交距离确定最优阈值；将检测模式重新定义为"宽松模式"和"严格模式"，前者采用膝点与峰值中点作为替代阈值，以提升高测序深度下的性能；同时增加95%最优阈值容忍搜索，避免在某些数据集中出现过度筛选问题。

SEACR还新增了过滤机制，排除由极少bedgraph行组成的信号块（即便其总信号通过阈值），从而降低假阳性峰的识别率。该工具提供命令行和网页两种使用方式，已发表于《Epigenetics and Chromatin》期刊（Meers等，2019）。

SEACR：用于 CUT&RUN 染色质图谱的稀疏富集峰值识别工具¶

概述¶

SEACR（Sparse Enrichment Analysis for CUT&RUN）是一款专为稀疏 CUT&RUN 或染色质图谱数据设计的峰值识别（peak calling）工具。与传统 ChIP-seq 数据相比，CUT&RUN 数据的背景信号极低，大量基因组区域的读段覆盖度为零，这使得常规峰值识别算法（如 MACS2）在处理此类数据时往往表现欠佳。SEACR 的核心设计目标正是针对这种"背景以零值主导"的稀疏数据特征，通过分析总信号分布的拐点与峰值来自适应地确定富集区域阈值。

该工具由 Meers MP、Tenenbaum D 与 Henikoff S 于 2019 年发表于 Epigenetics and Chromatin 期刊，目前已成为 CUT&RUN 和 CUT&TAG 数据分析流程中的标准工具之一。SEACR 依赖 R 语言和 Bedtools 运行，接受双端测序（paired-end sequencing）生成的 bedgraph 文件作为输入，支持使用 IgG 对照数据生成经验阈值，也支持通过数值阈值直接筛选信号最强的富集区域。

该工具还提供了在线 Web 界面，方便无编程基础的研究人员直接上传数据完成分析，降低了使用门槛。如果您在研究中使用了 SEACR，请引用原始文献：

Meers MP, Tenenbaum D, Henikoff S. (2019). Peak calling by Sparse Enrichment Analysis for CUT&RUN chromatin profiling. Epigenetics and Chromatin 12(1):42.

核心原理与功能¶

适用数据类型¶

SEACR 专为以下场景设计：

• CUT&RUN（Cleavage Under Targets and Release Using Nuclease）：一种靶向染色质剪切技术，产生背景极低的稀疏信号数据。
• CUT&TAG（Cleavage Under Targets and Tagmentation）：类似技术，同样产生稀疏分布的富集信号。
• 其他背景信号被"零值"主导的染色质图谱（chromatin profiling）数据。

阈值确定策略：两种模式¶

SEACR 提供两种峰值识别模式，均基于对总信号曲线（total signal curve）形态的分析：

Stringent（严格）模式 - 使用总信号曲线的峰值点（peak of the curve）作为阈值。 - 对应文献中描述的"stringent"模式，识别的富集区域更为保守、可信度更高。

Relaxed（宽松）模式 - 使用总信号曲线拐点（knee）与峰值之间的中间点作为阈值。 - 对应文献中的"relaxed"模式，可识别更多潜在富集区域，适合信号较弱或读深（read depth）较高的数据集。 - 注意：v1.3 版本移除了宽松模式中的最大信号阈值上限，改用拐点与峰值中点，以改善高读深数据下的表现。

版本历史说明：在 v1.3 之前，这两种模式被分别称为"AUC"模式和"union"模式，v1.3 起统一更名为"stringent"和"relaxed"。

阈值自适应优化机制¶

SEACR v1.3 实现了多项关键算法改进：

1. 正交距离最大化截断法：之前版本使用硬编码的 90% 分位数截断数据框（dataframe），v1.3 改为构建分位数（list quantile，即行号/最大行号）与值分位数（value quantile，即值/最大值）的对应关系，通过寻找到 (0,0)→(1,1) 连线正交距离最大的点来确定最优截断阈值，使阈值选择更具数据自适应性。

2. 95% 最优阈值搜索：实现了备选阈值检验机制，搜索所有达到最优阈值 95% 水平的阈值点，避免在某些数据集中出现过度严格的虚假阈值。

3. 信号块行数过滤器（Bedgraph Line Counter）：新增计数器追踪每个信号块（signal block）由多少条 bedgraph 行组成，并计算每个信号块的行数最小阈值——当目标信号块的剩余比例低于对照信号块时触发过滤，用于剔除虽通过总信号阈值但由极少 bedgraph 行构成的信号块（这类信号块不太可能是真实峰值）。

4. 零信号行过滤：增加了对输入 bedgraph 中含零信号行的过滤检查，防止零值干扰阈值计算。

归一化（Normalization）支持¶

SEACR 支持对对照数据（通常为 IgG 对照）进行归一化处理，以校正目标数据与对照数据之间的测序深度差异。建议在目标数据与对照数据未经过精确归一化（如 spike-in 归一化）时使用 norm 参数。

输入与依赖¶

• R 语言：需在系统路径（PATH）中可调用
• Bedtools：需在系统路径中可调用
• 输入格式：来自双端测序的 bedgraph 文件，反映读段对（read pair）的覆盖密度，而非单条读段

安装与使用¶

环境依赖¶

运行 SEACR 前，请确保以下工具已安装并加入系统 PATH：

• R（https://www.r-project.org）
• Bedtools（https://bedtools.readthedocs.io/en/latest/）

命令行语法¶

bash SEACR_1.3.sh experimental_bedgraph [control_bedgraph | numeric_threshold] ["norm" | "non"] ["relaxed" | "stringent"] output_prefix

参数详解¶

准备输入 Bedgraph 文件¶

Bedgraph 文件应反映读段对（read pair）的覆盖密度，而非单条读段密度。推荐从双端比对的 BAM 文件按如下步骤转换：

# 步骤1：将 BAM 文件转换为双端 BED 格式（-bedpe 输出配对信息）
bedtools bamtobed -bedpe -i $sample.bam > $sample.bed

# 步骤2：过滤跨染色体的配对及片段长度 >= 1000 bp 的读段对
awk '$1==$4 && $6-$2 < 1000 {print $0}' $sample.bed > $sample.clean.bed

# 步骤3：提取片段的起始和终止坐标（第1列：染色体，第2列：5'端，第6列：3'端），并排序
cut -f 1,2,6 $sample.clean.bed | sort -k1,1 -k2,2n -k3,3n > $sample.fragments.bed

# 步骤4：使用 bedtools genomecov 生成 bedgraph 文件（-bg 输出 bedgraph 格式）
bedtools genomecov -bg -i $sample.fragments.bed -g my.genome > $sample.fragments.bedgraph

关键说明：awk 过滤步骤中，$1==$4 确保配对双端位于同一染色体，$6-$2 < 1000 过滤异常大片段（通常为噪声），my.genome 为基因组染色体大小文件（可由 samtools faidx 生成）。

输出文件格式¶

SEACR 生成以下输出文件（以输出前缀命名）：

<output prefix>.stringent.bed   # stringent 模式的富集区域 BED 文件
<output prefix>.relaxed.bed     # relaxed 模式的富集区域 BED 文件

每个输出 BED 文件包含 6 列：

实战示例¶

以下三个示例均来自 SEACR 官方文档，涵盖最常见的使用场景。

示例 1：使用归一化 IgG 对照 + 严格阈值¶

# 使用归一化的 IgG 对照轨迹，以 stringent 阈值识别 target 数据中的富集区域
bash SEACR_1.3.sh target.bedgraph IgG.bedgraph norm stringent output

适用场景：目标样本与 IgG 对照测序深度差异较大，需要归一化校正；希望获得高可信度的峰值结果。输出文件为 output.stringent.bed。

示例 2：使用非归一化 IgG 对照 + 宽松阈值¶

# 使用非归一化的 IgG 对照轨迹，以 relaxed 阈值识别 target 数据中的富集区域
bash SEACR_1.3.sh target.bedgraph IgG.bedgraph non relaxed output

适用场景：目标数据与 IgG 对照已经过 spike-in 等严格归一化处理，无需 SEACR 内部归一化；或希望识别更多潜在富集区域（包括信号较弱的峰）。输出文件为 output.relaxed.bed。

示例 3：使用数值阈值（无对照）¶

# 不使用 IgG 对照，直接选取总信号（AUC）排名前 1% 的区域作为峰值
bash SEACR_1.3.sh target.bedgraph 0.01 non stringent output

适用场景：没有可用的 IgG 对照数据；希望基于信号强度直接筛选最富集的 1% 区域（0.01 即前 1% 比例）。该参数值可根据实验需求调整，例如 0.05 表示前 5%。

常见问题¶

Q1：SEACR 与 MACS2 有什么区别？为什么 CUT&RUN 数据推荐使用 SEACR？

MACS2 等传统峰值识别工具是为 ChIP-seq 数据设计的，假设背景信号服从泊松分布，在低背景稀疏数据中容易产生大量假阳性或假阴性结果。SEACR 专为 CUT&RUN/CUT&TAG 的稀疏信号特征设计，通过分析总信号曲线的形态特征自适应确定阈值，无需依赖背景信号分布假设，更适合背景以零值主导的数据。

Q2：什么时候应该选择 stringent 模式，什么时候选择 relaxed 模式？

• 选择 stringent：当优先追求峰值的可信度和精确性时，或数据信噪比较高时。
• 选择 relaxed：当信号较弱、需要识别更多潜在富集区域时，或处理高读深数据时（v1.3 改进了宽松模式在高读深情况下的表现）。
• 实践中建议两种模式均运行，对比结果后根据下游分析目标选择。

Q3：Field 2 使用数值阈值（如 0.01）和使用 IgG 对照有何区别？

使用 IgG 对照时，SEACR 根据对照数据的信号分布生成经验阈值，能够更准确地区分真实富集信号与背景噪声，是推荐的首选方法。使用数值阈值（0~1 之间）时，SEACR 直接按总信号（AUC）排名返回前 n 比例的区域，不依赖对照数据，适用于无对照样本的情况，但阈值的生物学意义较弱。

Q4：为什么 bedgraph 文件需要反映读段对密度而非单条读段密度？如何避免错误？

CUT&RUN 技术产生的是片段（fragment）信号，使用读段对（read pair）的 5' 和 3' 坐标定义完整片段，能更准确地反映蛋白质结合位点的实际范围。如果直接使用单条读段，可能产生偏差。官方推荐的 BAM→BED→bedgraph 转换流程（见"安装与使用"章节）已内置片段长度过滤（< 1000 bp）和同染色体配对过滤，应严格遵循该流程以确保输入数据格式正确。

Q5：运行 SEACR 报错"找不到 R 或 Bedtools"，如何解决？

SEACR 要求 R 和 Bedtools 均在系统 PATH 中可直接调用。可在终端分别运行 R --version 和 bedtools --version 验证是否已安装。若未安装，请参考 R 官网（https://www.r-project.org）和 Bedtools 文档（https://bedtools.readthedocs.io/en/latest/）进行安装。如使用 conda 环境，可通过 conda install -c bioconda bedtools r-base 安装两者，并激活对应环境后再运行 SEACR。

总结¶

SEACR 是专为 CUT&RUN 和 CUT&TAG 稀疏染色质图谱数据设计的峰值识别工具，其核心优势在于针对背景零值主导的数据特征进行了专门的算法优化。工具通过分析总信号曲线的拐点与峰值，提供 stringent（严格）和 relaxed（宽松）两种识别模式，支持使用 IgG 对照数据或数值阈值进行峰值筛选，并支持对照数据的内部归一化。v1.3 版本引入了正交距离最大化截断法、95% 最优阈值搜索、信号块行数过滤等多项算法改进，显著提升了高读深数据下的识别准确性。使用时需注意正确准备 bedgraph 输入文件（反映读段对密度），并根据实验设计选择合适的归一化方式和阈值模式。该工具还提供 Web 界面，适合不同技术背景的研究人员使用。