bedtools区间操作：intersect/merge/coverage/closest

1. BED 格式——先搞懂你操作的对象#

BED（Browser Extensible Data）是 bedtools 的默认输入格式，最简单：

1
chr1  100  500  geneA  0  +
2
chr1  600  900  geneB  0  -

BED 最少 3 列（chrom, start, end），常见到 6 列：

最重要的一点——BED 是 0-based，跟人类直觉的 1-based 不同：

踩坑： 如果你从 GFF（1-based）手动转 BED，需要把 start 减 1：bed_start = gff_start - 1。

2. intersect——求交集，最核心的功能#

1
# 基因A（genes.bed）和peak（peaks.bed）在哪些区域重叠？
2
bedtools intersect -a genes.bed -b peaks.bed

1
# 要求基因区域至少50%被peak覆盖
2
bedtools intersect -a genes.bed -b peaks.bed -f 0.5
3

4
# 要求重叠>=100bp
5
bedtools intersect -a genes.bed -b peaks.bed -F 0.5 -e

1
# 先扩展promoter区域（TSS±2000bp作为promoter）
2
bedtools slop -i genes.bed -g chrom.sizes -b 2000 > promoters.bed
3

4
# 求peak和promoter的交集
5
bedtools intersect -a peaks.bed -b promoters.bed -wa -wb > peak_promoter.txt

1
chr1  150  450  peak1  100  chr1  100  500  geneA  0  +
2
chr1  580  850  peak2  80   chr1  600  900  geneB  0  -

1
bedtools intersect -a genes.bed -b peaks.bed -wa -wb -wo

3. merge——合并重叠区域#

ChIP-seq 的原始 peak 经常有多个相邻 peak 实际上来自同一个结合事件。merge 把它们合起来：

1
# 基础合并
2
bedtools merge -i peaks.bed
3

4
# 合并时要求重叠≥100bp才算连接
5
bedtools merge -i peaks.bed -d 100
6

7
# 合并并保留统计信息（如：合并了多少个原始peak）
8
bedtools merge -i peaks.bed -c 4,5 -o count,sum

-c 4,5 表示对第4列做 count（多少个原始peak）、第5列做 sum（总信号强度）。

4. coverage——计算覆盖度#

这在 RNA-seq 中用于计算每个基因的 reads 覆盖情况：

1
bedtools coverage -a genes.bed -b reads.bed

输出解读（每列含义）：

1
chr1  100  500  geneA  0  +  1  300  400  0.7500000
2
                               │   │    │      │
3
                reads数 ───────┘   │    │      │
4
           covered bases ──────────┘    │      │
5
             gene length ───────────────┘      │
6
        coverage ratio ────────────────────────┘

$coverage\_ratio = \frac{covered\_bases}{gene\_length}$

如果 coverage_ratio = 0.75，说明该基因 75% 的区域有 reads 覆盖。

生信场景： RNA-seq 的基因表达量化（老的 featureCounts 前身就是类似逻辑），或者 ChIP-seq 的 peak 富集度评估。

5. closest——找最近邻#

“每个基因 TSS 最近的 ChIP-seq peak 是哪个？”

1
bedtools closest -a genes.bed -b peaks.bed -d

-d 输出额外一列：两个区间之间的距离（bp）。如果重叠则为 0。

6. slop——扩展/收缩区间#

1
# 每个基因上下游各扩展200bp
2
bedtools slop -i genes.bed -g chrom.sizes -b 200
3

4
# 只向上游扩展（5'方向）
5
bedtools slop -i genes.bed -g chrom.sizes -l 500 -r 0
6

7
# 对于正链基因，-l 代表上游（5'端）；对于负链基因，自动翻转
8
# 用 -s 参数根据 strand 自动判断方向：
9
bedtools slop -i genes.bed -g chrom.sizes -l 2000 -r 500 -s

注意： -g chrom.sizes（染色体大小文件）是必须的——否则 slop 不知道染色体边界在哪里，可能导致区间超出染色体范围。

chrom.sizes 文件制作：

1
samtools faidx ref.fa
2
cut -f1,2 ref.fa.fai > chrom.sizes

7. genomecov——生成覆盖轨迹#

生成全基因组的覆盖度 bedGraph，之后可转为 bigWig 上传到 UCSC：

1
bedtools genomecov -i reads.bed -g chrom.sizes -bg > coverage.bedgraph
2

3
# -bg 输出 bedGraph 格式
4
# -bga 连零覆盖区域也输出（文件更大但更完整）
5
# -split 对包含N>CIGAR的BAM输入时，避免把intron区域也算进去

转 bigWig：

1
bedGraphToBigWig coverage.bedgraph chrom.sizes coverage.bw

8. 踩坑记录#

1
# 查看双方用了什么染色体名
2
cut -f1 genes.bed | sort -u
3
cut -f1 peaks.bed | sort -u

1
# 先排序再操作
2
sort -k1,1 -k2,2n genes.bed > genes_sorted.bed
3
bedtools intersect -a genes_sorted.bed -b peaks_sorted.bed -sorted

9. 小结#

本文于 2025-03-18 在 Debian 13 上实测完成。bedtools v2.31.1，测试数据均为模拟 BED 文件。