ChIP-seq分析：提取峰±1kb序列并生成Bed文件

一、生成±1kb窗口的BED文件 #

从你的HOMER注释输出中，可直接用awk处理制表符分隔文本，快速生成目标窗口的BED文件，分两种常见需求处理：

1. 基于峰中心扩展±1kb（推荐，聚焦峰核心区域） #

HOMER输出坐标通常为1-based，而BED格式要求0-based左闭右开，计算时需做转换：

# 跳过表头，计算峰中心并生成±1kb窗口的BED
awk 'NR>1 {
    # 转换为0-based坐标计算中心
    bed_start = $3 - 1
    bed_end = $4
    center = int((bed_start + bed_end) / 2)
    win_start = center - 1000
    win_end = center + 1000
    # 避免窗口起始小于0
    if (win_start < 0) win_start = 0
    # 输出BED格式列：染色体、起始、终止、峰ID、得分、链
    print $2 "\t" win_start "\t" win_end "\t" $1 "\t" $6 "\t" $5
}' homer_annotation.txt > peaks_center_plus_minus_1kb.bed

2. 基于峰全区间扩展±1kb（覆盖峰区域及上下游1kb） #

若需要保留峰的完整区间再扩展上下游，用以下命令：

awk 'NR>1 {
    bed_start = $3 - 1
    win_start = bed_start - 1000
    win_end = $4 + 1000
    if (win_start < 0) win_start = 0
    print $2 "\t" win_start "\t" win_end "\t" $1 "\t" $6 "\t" $5
}' homer_annotation.txt > peaks_region_plus_minus_1kb.bed

二、提取窗口内的序列 #

以下两种工具适配不同工作流，选你熟悉的即可：

1. 用bedtools（通用基因组工具） #

先确保参考基因组refGenome.fa已生成索引（执行samtools faidx refGenome.fa生成.fai文件），再运行：

bedtools getfasta -fi refGenome.fa \
                  -bed peaks_center_plus_minus_1kb.bed \
                  -fo peaks_plus_minus_1kb.fa \
                  -name  # 序列名用BED的第四列峰ID

2. 用HOMER自带工具（适配现有工作流） #

既然已经在用HOMER，直接用extractGenomeBed更便捷，无需额外索引：

extractGenomeBed -bed peaks_center_plus_minus_1kb.bed \
                 -genome refGenome.fa \
                 -fa peaks_plus_minus_1kb.fa

注意事项 #

• 坐标系统：确认HOMER输出坐标是1-based还是0-based（MACS3的peak文件为0-based，若直接传给HOMER，输出可能保留0-based），如果是0-based，命令中bed_start = $3 - 1这一步需去掉。
• 边界处理：命令中加入if (win_start < 0) win_start = 0避免窗口起始为负数，防止提取序列时出错。
• 链方向处理：若需提取负链峰的互补序列，bedtools可加-s参数，HOMER的extractGenomeBed会自动根据链方向处理。

内容的提问来源于stack exchange，提问作者LoganLee