Pipeline Detail

CUT&Tag表观调控实验

CUT&Tag T2T 基因组差异结合分析

面向 CUT&Tag 数据的 T2T 基因组比对、CPM 标准化、SEACR 搜峰、DiffBind 差异 Peak 与 GO 富集分析流程。

创建时间

2026/6/3

分析难度

入门

推荐场景

表观调控

预计耗时

0.5-1 天

Metadata

流程元数据

先看应用场景、输入输出和工具依赖，再进入正文命令细节。

Difficulty

入门

Scenario

表观调控

Estimated Time

0.5-1 天

Tools

DESeq2STAR

Inputs

FASTQBAM

Outputs

report

Workflow DAG

流程图

用步骤节点快速理解这个分析从原始数据到结果报告的流转关系。

STEP 1

Trim Galore/FastQC 质控

→

STEP 2

Bowtie2 构建 T2T 索引

→

STEP 3

T2T 基因组严格双端比对

→

STEP 4

Spike-in 诊断

→

STEP 5

CPM BedGraph 生成

→

STEP 6

SEACR 精确搜峰

→

STEP 7

DiffBind 差异 Peak

→

STEP 8

ChIPseeker/GO 注释

Protocol

流程文档

正文保留 Markdown 排版、代码语言标识和表格样式，适合边学边复现。

CUT&Tag T2T 基因组差异结合分析

适用场景

适用于 CUT&Tag 或 CUT&RUN 染色质结合数据，尤其是使用 T2T 级别参考基因组、需要比较处理组与对照组蛋白结合重分布的项目。

本流程来自 Zeocin vs Control 的 TOP1 CUT&Tag 分析场景：使用拟南芥 Col-PEK T2T 基因组进行 Bowtie2 比对，结合 Spike-in 诊断、CPM 标准化、SEACR 搜峰、DiffBind 差异结合分析、ChIPseeker 注释和 GO 富集，解释 DNA 损伤诱导的 TOP1 全基因组重分布。

样本设计

SampleID	Condition	Replicate	BAM	Peak
C1	Control	1	`C1_T2T.sorted.bam`	`C1_seacr.stringent.bed`
C2	Control	2	`C2_T2T.sorted.bam`	`C2_seacr.stringent.bed`
C3	Control	3	`C3_T2T.sorted.bam`	`C3_seacr.stringent.bed`
Z1	Zeocin	1	`Z1_T2T.sorted.bam`	`Z1_seacr.stringent.bed`
Z2	Zeocin	2	`Z2_T2T.sorted.bam`	`Z2_seacr.stringent.bed`
Z3	Zeocin	3	`Z3_T2T.sorted.bam`	`Z3_seacr.stringent.bed`

1. 原始 reads 质控与接头过滤

WORKDIR="/public/home/yhpeng/cut_tag"
RAW_DIR="${WORKDIR}/01_raw_data"
OUT_DIR="${WORKDIR}/03_trimmed_data"

for sample in C1 C2 C3 Z1 Z2 Z3
do
  trim_galore     --paired     --quality 20     --phred33     --fastqc     --cores 4     --gzip     -o ${OUT_DIR}     ${RAW_DIR}/${sample}_R1.fq.gz     ${RAW_DIR}/${sample}_R2.fq.gz
done

2. 构建 T2T 基因组 Bowtie2 索引

WORKDIR="/public/home/yhpeng/cut_tag"
REF_DIR="${WORKDIR}/ref_genenic"

bowtie2-build   --threads 4   ${REF_DIR}/Arabidopsis_thaliana_Col-PEK.genomic.fa   ${REF_DIR}/Col_PEK_T2T

3. 比对到 T2T 基因组并过滤多重比对

WORKDIR="/public/home/yhpeng/cut_tag"
TRIM_DIR="${WORKDIR}/03_trimmed_data"
ALN_DIR="${WORKDIR}/04_alignment"
REF_INDEX="${WORKDIR}/ref_genenic/Col_PEK_T2T"

for sample in C1 C2 C3 Z1 Z2 Z3
do
  bowtie2     --end-to-end     --very-sensitive     --no-mixed     --no-discordant     --phred33     -I 10     -X 700     -p 8     -x ${REF_INDEX}     -1 ${TRIM_DIR}/${sample}_R1_val_1.fq.gz     -2 ${TRIM_DIR}/${sample}_R2_val_2.fq.gz     -S ${ALN_DIR}/${sample}_T2T.sam

  samtools view -bS -q 20 -@ 8 ${ALN_DIR}/${sample}_T2T.sam |     samtools sort -@ 8 -o ${ALN_DIR}/${sample}_T2T.sorted.bam

  samtools index ${ALN_DIR}/${sample}_T2T.sorted.bam
  rm ${ALN_DIR}/${sample}_T2T.sam
done

4. Spike-in 诊断与归一化策略

Group	Samples	Spike-in rate
Control	C1, C2, C3	~6.63%
Zeocin	Z1, Z2, Z3	~1.45%

Zeocin 处理组 Spike-in 比对率相较对照组出现约 4.5 倍下降，同时拟南芥基因组比对 reads 明显上升。这说明处理组中真实靶序列发生大规模扩增，外源 Spike-in 被物理稀释。

关键判断： 对于这种全局剧烈变化的 CUT&Tag 实验，不应直接使用外源 Spike-in scale factor 强行放大处理组信号，否则容易造成全基因组背景假阳性。这里采用 CPM 或 DESeq2 内部相对定量，重点寻找特异性重分布靶点。

5. 生成 CPM 标准化 BedGraph

WORKDIR="/public/home/yhpeng/cut_tag"
ALN_DIR="${WORKDIR}/04_alignment"
PEAK_DIR="${WORKDIR}/05_peaks"
mkdir -p ${PEAK_DIR}

for sample in C1 C2 C3 Z1 Z2 Z3
do
  bamCoverage     -b ${ALN_DIR}/${sample}_T2T.sorted.bam     -o ${PEAK_DIR}/${sample}.bedgraph     --outFileFormat bedgraph     --normalizeUsing CPM     --binSize 10     -p 8
done

6. SEACR 无 IgG 模式搜峰

SEACR 适合 CUT&Tag 的窄峰/稀疏信号场景。对于无 IgG 对照的实验，可以使用 stringent 模式获得更干净的候选峰。

wget https://github.com/FredHutch/SEACR/raw/master/SEACR_1.3.sh
wget https://github.com/FredHutch/SEACR/raw/master/SEACR_1.3.R

bash SEACR_1.3.sh C1.bedgraph 0.01 non stringent C1_seacr

7. DiffBind 差异 Peak 定量

library(DiffBind)

cut_tag <- dba(sampleSheet = "sample_sheet.csv")
cut_tag <- dba.count(cut_tag, bUseSummarizeOverlaps = TRUE)
cut_tag <- dba.contrast(cut_tag, categories = DBA_CONDITION, minMembers = 2)
cut_tag <- dba.analyze(cut_tag)

res_deseq <- dba.report(cut_tag, method = DBA_DESEQ2, contrast = 1, th = 0.05)
write.csv(as.data.frame(res_deseq), file = "Zeocin_vs_Control_DiffPeaks.csv")

dba.plotPCA(cut_tag, DBA_CONDITION, label = DBA_ID)
dba.plotVolcano(cut_tag)

8. T2T 注释、TAIR 名称转换与 GO 富集

library(ChIPseeker)
library(GenomicFeatures)
library(GenomicRanges)
library(txdbmaker)

txdb <- makeTxDbFromGFF("Arabidopsis_thaliana_Col-PEK.genomic.gff", format = "gff3")
diff_peaks <- read.csv("Zeocin_vs_Control_DiffPeaks.csv")

peak_gr <- GRanges(
  seqnames = diff_peaks$seqnames,
  ranges = IRanges(start = diff_peaks$start, end = diff_peaks$end),
  strand = diff_peaks$strand
)
mcols(peak_gr) <- diff_peaks[, c("Fold", "FDR")]

peak_anno <- annotatePeak(peak_gr, tssRegion = c(-3000, 3000), TxDb = txdb)
write.csv(as.data.frame(peak_anno), "Zeocin_vs_Control_Annotated_DiffPeaks.csv", row.names = FALSE)

library(clusterProfiler)
library(org.At.tair.db)

diff_genes <- read.csv("Zeocin_vs_Control_TAIR_Ready.csv")
diff_genes <- diff_genes[grep("^AT[1-5]G\d{5}$", diff_genes$tairId), ]

genes_up <- unique(diff_genes$tairId[diff_genes$Fold > 0])
genes_down <- unique(diff_genes$tairId[diff_genes$Fold < 0])

go_up <- enrichGO(gene = genes_up, OrgDb = org.At.tair.db, keyType = "TAIR", ont = "BP")
go_down <- enrichGO(gene = genes_down, OrgDb = org.At.tair.db, keyType = "TAIR", ont = "ALL")

生物学解释

正常状态下，TOP1 更集中地占据启动子或活跃转录区，形成清晰的局部高峰；Zeocin 诱导 DNA 损伤后，TOP1 发生全基因组重分布，更多 reads 映射到拟南芥基因组，但局部孤立峰可能减少，呈现“高原化”的广泛结合模式。

主要输出

Trim Galore/FastQC 质控报告
T2T 基因组比对 BAM 与索引
Spike-in 诊断统计
CPM 标准化 BedGraph
SEACR peak BED 文件
DiffBind 差异 Peak 表
ChIPseeker 注释结果
TAIR ID 转换表与 GO 富集图