bedtools — 유전체 구간 연산: intersect·merge·coverage와 BED 포맷

CODE BENCH · NGS PIPELINES

bedtools — 유전체 구간 연산: intersect·merge·coverage와 BED 포맷

TL;DR — bedtools는 유전체 위의 좌표 구간(genomic interval)을 다루는 명령줄 도구 모음으로, '유전체 산술(genome arithmetic)의 스위스 아미 나이프'라 불립니다. BED·BAM·VCF·GFF에 담긴 구간을 교집합(intersect)·병합(merge)·커버리지(coverage) 같은 연산으로 한 줄에 처리하죠.

핵심은 BED 포맷과 그 좌표계입니다. BED는 chrom·start·end 세 칸이 기본인데, start가 0-based이고 end는 배제되는 반열림(half-open) 좌표라 chr1 0 100은 첫 100개 염기를 뜻합니다. 이 0-based가 1-based인 VCF·GFF와 자주 헷갈리는 지점이라 이 글에서 특히 강조합니다.

작은 예제 BED 파일로 intersect의 -wa·-wb·-v·-c, merge·sort, coverage·genomecov, subtract·closest·slop을 한 명령씩 따라간 뒤, ChIP-seq 피크와 프로모터의 교집합, 변이(VCF)와 엑손의 교집합을 samtools와 파이프(|)로 잇는 실전까지 만듭니다. 명령은 bedtools 2.31 기준입니다.

🔗 관련 글 · bedtools가 다루는 구간은 대개 다른 도구의 산출물입니다 — 정렬 결과는 SAM/BAM 파일이란?과 서열 정렬이란?, 변이 정보는 VCF 파일이란?, 그 BAM·VCF를 직접 다루는 samtools·bcftools 실전, 그리고 구간의 대표 사례인 ChIP-seq 피크에서 이어집니다.

이 글에서 만들 것

NGS 분석을 하다 보면 결국 좌표를 가진 '구간'을 비교하는 일이 반복됩니다. ChIP-seq 피크가 프로모터 안에 있나, 변이가 엑손에 걸리나, 이 유전자 구간의 커버리지는 얼마인가 같은 질문이죠. 스프레드시트로 좌표를 일일이 맞춰 보긴 어렵고, 직접 파이썬으로 겹침을 계산하자니 경계 조건에서 실수하기 쉽습니다.

bedtools는 이런 '구간 대 구간' 질문을 명령 한 줄로 풉니다. 두 파일의 구간이 겹치는지(intersect), 한 파일 안에서 인접한 구간을 하나로 합칠지(merge), 어떤 구간이 리드로 얼마나 덮이는지(coverage)를 각각 전용 서브명령으로 처리하죠. 입력은 BED뿐 아니라 BAM·VCF·GFF도 되고, 출력은 다시 BED라 파이프(|)로 samtools 같은 도구와 자연스럽게 이어집니다.

이 글은 그 bedtools를 처음부터 익힙니다. 구체적으로 ① BED 포맷과 0-based 좌표를 짚고, ② intersect로 교집합·겹침을 뽑고, ③ merge·sort로 구간을 병합하고, ④ coverage·genomecov로 커버리지를 계산하고, ⑤ subtract·closest·slop으로 빼기·최근접·확장을 다룬 뒤, ⑥ ChIP-seq 피크 ∩ 프로모터와 변이 ∩ 엑손을 실전으로 묶습니다.

# 핵심 흐름 한눈에
bedtools intersect -a peaks.bed -b promoters.bed -u   # ① 겹치는 구간(교집합)
bedtools merge -i sorted.bed                          # ② 인접·겹치는 구간 병합
bedtools coverage -a windows.bed -b reads.bed         # ③ 구간별 커버리지

1. 사전 준비 — 설치와 BED 포맷

bedtools는 리눅스·macOS용 명령줄 도구입니다. conda로 까는 게 가장 간편합니다.

# conda (bioconda) — 권장
conda install -c bioconda bedtools

# 설치 확인(버전 출력)
bedtools --version

bedtools v2.31.0

이제 포맷을 봅시다. BED (Browser Extensible Data)는 구간을 적는 표준 텍스트 포맷으로, 탭으로 나뉜 세 칸이 기본입니다.

• chrom — 염색체 이름(chr1)
• chromStart — 시작 좌표, 0-based
• chromEnd — 끝 좌표, 배제됨(half-open)

여기서 좌표계가 가장 중요합니다. BED는 0-based 반열림(half-open) 좌표라, 시작은 0부터 세고 끝은 포함하지 않습니다. 그래서 chr1 0 100은 1번째부터 100번째까지 100개 염기를 가리키고, 단일 염기 하나(1-based로 5번 위치)는 chr1 4 5로 적습니다. 반면 VCF·GFF·SAM은 1-based라 같은 위치를 5로 씁니다. bedtools를 쓸 때 좌표가 1씩 어긋나 보이면 대개 이 0-based를 잊은 경우죠.

칸을 더 붙이면 정보가 늘어납니다. BED3는 위 세 칸, BED6은 여기에 name·score·strand를 더한 여섯 칸, BED12는 exon 블록 구조(blockCount·blockSizes·blockStarts)까지 담아 유전자 하나를 통째로 표현합니다.

그림 1. BED 포맷과 0-based 반열림 좌표, 그리고 두 구간의 교집합(A∩B) — bedtools의 두 가지 기초.

예제로 쓸 작은 파일 두 개를 만듭니다. peaks.bed는 ChIP-seq 피크, genes.bed는 유전자 구간이라 치죠.

# peaks.bed — 피크 3개
printf "chr1\t100\t200\tpeak1\nchr1\t400\t500\tpeak2\nchr1\t700\t900\tpeak3\n" > peaks.bed

# genes.bed — 유전자 2개
printf "chr1\t150\t350\tgeneA\nchr1\t800\t1000\tgeneB\n" > genes.bed

2. intersect — 구간 교집합과 겹침

가장 많이 쓰는 서브명령이 intersect입니다. -a와 -b 두 파일을 받아, A의 각 구간이 B의 구간과 겹치는 부분을 찾습니다. 옵션 없이 돌리면 겹치는 구간(교집합) 자체를 출력합니다.

# A(peaks)와 B(genes)가 겹치는 부분만 출력
bedtools intersect -a peaks.bed -b genes.bed

chr1	150	200	peak1
chr1	800	900	peak3

peak1(100–200)은 geneA(150–350)와 150–200에서 겹치고, peak3(700–900)은 geneB(800–1000)와 800–900에서 겹칩니다. peak2(400–500)는 어느 유전자와도 안 겹쳐 빠졌죠. 출력의 좌표는 겹친 구간으로 잘려 나오고, peak 이름 같은 A의 나머지 칸은 그대로 따라옵니다.

원본 구간을 통째로 보고 싶으면 옵션을 붙입니다. 자주 쓰는 조합은 이렇습니다.

• -wa — 겹친 A 구간을 원본 그대로 출력(잘리지 않음)
• -wb — 짝이 된 B 구간도 함께 출력
• -v — B와 안 겹치는 A만(여집합)
• -c — A마다 겹친 B의 개수를 뒤에 붙임
• -u — 겹침이 하나라도 있으면 A를 한 번만 출력(중복 제거)

# -wa -wb: 원본 A와 짝이 된 B를 나란히
bedtools intersect -a peaks.bed -b genes.bed -wa -wb

chr1	100	200	peak1	chr1	150	350	geneA
chr1	700	900	peak3	chr1	800	1000	geneB

# -v: 어느 유전자와도 안 겹치는 피크(여집합)
bedtools intersect -a peaks.bed -b genes.bed -v

chr1	400	500	peak2

# -c: 피크마다 겹친 유전자 수를 마지막 칸에
bedtools intersect -a peaks.bed -b genes.bed -c

chr1	100	200	peak1	1
chr1	400	500	peak2	0
chr1	700	900	peak3	1

-a에는 파일 여럿을, 겹침 기준에는 -f(A의 최소 겹침 비율)·-r(양쪽 상호 비율)을 더할 수 있습니다. 예컨대 -f 0.5는 A가 절반 이상 겹칠 때만 셈에 넣죠.

3. merge·sort — 인접·겹치는 구간 병합

merge는 한 파일 안에서 겹치거나 맞닿은 구간을 하나로 합칩니다. 단, merge는 입력이 정렬돼 있어야 하므로 먼저 sort로 좌표순(염색체→시작)으로 맞춥니다.

# regions.bed — 겹치거나 인접한 구간들
printf "chr1\t100\t200\nchr1\t180\t300\nchr1\t400\t500\nchr1\t480\t600\nchr1\t900\t1000\n" > regions.bed

# 좌표순 정렬(염색체별, 시작 위치 오름차순)
bedtools sort -i regions.bed > regions.sorted.bed

# 겹치거나 맞닿은 구간을 하나로 병합
bedtools merge -i regions.sorted.bed

chr1	100	300
chr1	400	600
chr1	900	1000

100–200과 180–300은 겹쳐 100–300으로, 400–500과 480–600은 겹쳐 400–600으로 합쳐지고, 멀리 떨어진 900–1000은 그대로 남습니다. 조금 떨어진 구간까지 합치려면 -d로 허용 간격을 줍니다. -d 350이면 600과 900 사이(300 간격)도 이어 붙죠.

병합하면서 값을 요약할 수도 있습니다. -c로 대상 칸을, -o로 연산(count·sum·mean·collapse 등)을 지정합니다.

# 병합된 각 구간이 원래 몇 개 구간에서 왔는지 세기
bedtools merge -i regions.sorted.bed -c 1 -o count

chr1	100	300	2
chr1	400	600	2
chr1	900	1000	1

4. coverage·genomecov — 커버리지 계산

coverage는 A의 각 구간이 B의 구간(리드 등)으로 얼마나 덮이는지를 셉니다. 기본 출력은 각 A 구간 뒤에 겹친 B의 개수, 덮인 염기 수, 구간 길이, 덮인 비율을 덧붙입니다.

# windows.bed의 각 구간을 reads.bed가 얼마나 덮는지
bedtools coverage -a windows.bed -b reads.bed

chr1	0	500	12	480	500	0.9600000
#                ↑개수 ↑덮인bp ↑길이 ↑비율

마지막 열의 비율(0.96)은 이 창(window)의 96%가 리드로 덮였다는 뜻입니다. 창 단위 커버리지를 훑을 때 유용하죠. 구간별이 아니라 염기별(per-base) 커버리지가 필요하면 genomecov를 씁니다. BAM을 바로 받아 BedGraph로 뽑을 수 있습니다.

# BAM에서 염기별 커버리지를 BedGraph로(-bg: 0인 구간은 생략)
bedtools genomecov -ibam aln.sorted.bam -bg > coverage.bedgraph

chr1	0	120	3
chr1	120	340	7
chr1	340	500	2

각 줄은 '이 구간의 모든 염기가 N개 리드로 덮였다'는 뜻입니다. genomecov는 BAM이 좌표순으로 정렬돼 있어야 하고, 리드가 아닌 BED 입력이면 -i와 함께 염색체 크기 파일(-g)을 줍니다.

5. subtract·closest·slop — 빼기·최근접·확장

나머지 자주 쓰는 서브명령 셋을 묶어 봅니다.

subtract는 A에서 B와 겹치는 부분을 잘라 냅니다. 블랙리스트 구간이나 반복서열을 제거할 때 씁니다.

# peaks에서 genes와 겹치는 부분을 제거
bedtools subtract -a peaks.bed -b genes.bed

chr1	100	150	peak1
chr1	400	500	peak2
chr1	700	800	peak3

peak1은 150부터 geneA와 겹쳐 100–150만 남고, peak3은 800부터 겹쳐 700–800만 남습니다. peak2는 겹침이 없어 통째로 유지되죠.

closest는 A의 각 구간에 가장 가까운 B 구간을 찾습니다. 겹치는 게 없어도 제일 가까운 것을 돌려주죠. sort된 입력이 필요하고, -d를 붙이면 거리(bp)를 마지막 칸에 적습니다.

# 각 피크에서 가장 가까운 유전자와 그 거리
bedtools closest -a peaks.sorted.bed -b genes.sorted.bed -d

chr1	100	200	peak1	chr1	150	350	geneA	0
chr1	400	500	peak2	chr1	150	350	geneA	51
chr1	700	900	peak3	chr1	800	1000	geneB	0

거리 0은 겹친다는 뜻이고, peak2는 geneA에서 51bp 떨어져 있습니다(0-based라 400과 349의 차이).

slop은 각 구간을 양옆으로 넓힙니다. TSS 주변 프로모터를 만들 때처럼 자주 쓰죠. 염색체 밖으로 삐져나가지 않게 크기 파일(-g)이 필요합니다.

# genome.txt: 염색체 크기(chrom <TAB> size)
printf "chr1\t248956422\n" > genome.txt

# 각 구간을 양쪽으로 50bp씩 확장(-b), 한쪽만 넓히려면 -l/-r
bedtools slop -i genes.bed -g genome.txt -b 50

chr1	100	400	geneA
chr1	750	1050	geneB

한쪽으로만 크기 0인 flanking 구간을 만들려면 flank를 씁니다. slop이 원본을 포함해 넓힌다면, flank는 원본 양옆의 새 구간만 돌려주죠.

6. 실전 — ChIP-seq 피크 ∩ 프로모터, 변이 ∩ 엑손

그림 2. 피크 ∩ 프로모터 워크플로 — TSS → slop → 프로모터, peaks와 intersect로 프로모터 속 피크만. 곁가지는 변이 ∩ 엑손과 파이프 연결.

이제 조각을 모아 실전 질문에 답합니다. 첫째, 어떤 ChIP-seq 피크가 프로모터에 있나? 전사 시작점(TSS) 좌표에서 slop으로 프로모터(TSS ±2kb)를 만든 뒤, intersect의 -u로 프로모터에 걸친 피크만 추립니다.

# 1) TSS ±2kb로 프로모터 구간 만들기
bedtools slop -i tss.bed -g genome.txt -b 2000 > promoters.bed

# 2) 프로모터에 겹치는 피크만(한 번씩)
bedtools intersect -a peaks.bed -b promoters.bed -u > peaks_in_promoters.bed

둘째, 어떤 변이가 엑손에 있나? intersect는 VCF도 그대로 받습니다. -header를 주면 VCF 헤더를 살려 결과가 다시 유효한 VCF가 되죠.

# 엑손에 걸친 변이만(VCF 헤더 유지)
bedtools intersect -a variants.vcf -b exons.bed -header > exonic_variants.vcf

마지막으로, bedtools의 출력은 다시 BED라 파이프로 이어 붙일 수 있습니다. 예컨대 특정 구간에 정렬된 리드만 BAM으로 뽑아 바로 커버리지를 볼 수 있죠.

# 표적 구간(-L)의 리드만 뽑아 그대로 커버리지 계산
samtools view -b -L targets.bed aln.bam \
  | bedtools coverage -a targets.bed -b - > target_coverage.txt

-b -의 -는 '앞 명령의 출력을 표준입력으로 받으라'는 뜻입니다. 중간 파일 없이 samtools와 bedtools가 곧장 연결되죠.

그림 3. 세 핵심 명령의 입출력 — intersect는 교집합, merge는 병합, coverage는 구간별 덮인 개수·비율을 낸다.

7. 자주 발생하는 에러 & 해결

• ***** ERROR: Requested column ... , but the number of columns ... — 입력이 탭이 아니라 공백으로 나뉘었을 때 흔합니다. BED는 탭 구분이 원칙이라, 스프레드시트에서 붙여넣은 공백을 sed 's/ /\t/g'나 printf "...\t..."로 탭으로 바꿔야 합니다.
• Error: Sorted input specified, but the file ... has ... out of order record — merge·closest처럼 정렬을 요구하는 명령에 정렬 안 된 파일을 넣은 경우입니다. 먼저 bedtools sort -i in.bed(또는 sort -k1,1 -k2,2n)로 정렬하세요.
• 좌표가 1씩 어긋나 보임 — 0-based를 잊은 경우입니다. BED start는 0-based, VCF·GFF는 1-based라, 같은 위치라도 BED가 1 작습니다. 다른 포맷과 맞대볼 때는 좌표계를 먼저 확인하세요.
• Error: Unable to open file ... genome file — slop·genomecov에 크기 파일(-g)을 안 줬거나 경로가 틀린 경우입니다. samtools faidx ref.fa로 만든 ref.fa.fai의 앞 두 칸을 쓰거나, chrom<TAB>size 형식으로 직접 만듭니다.
• 염색체 이름이 안 맞음 — 한쪽은 chr1, 다른 쪽은 1이면 겹침이 0으로 나옵니다. 이름 규칙(UCSC식 chr1 vs Ensembl식 1)을 한쪽으로 통일해야 합니다.

8. 확장 — 더 해 볼 것

• complement로 빈 구간 찾기 — bedtools complement는 크기 파일 기준으로 구간이 없는 영역을 돌려줍니다. 캡처 패널 밖 영역이나 갭을 찾을 때 유용하죠.
• makewindows로 유전체 창 나누기 — bedtools makewindows -g genome.txt -w 10000은 유전체를 10kb 창으로 잘라, coverage와 짝지어 창 단위 커버리지·밀도 프로파일을 만듭니다.
• getfasta로 서열 뽑기 — bedtools getfasta -fi ref.fa -bed peaks.bed는 구간에 해당하는 실제 염기서열을 FASTA로 추출합니다. 모티프 분석의 입력으로 바로 이어지죠.
• pybedtools로 파이썬에서 — 같은 연산을 파이썬 스크립트 안에서 하려면 pybedtools가 bedtools를 감싸 줍니다. pandas와 엮어 대량 구간을 프로그램적으로 다루기 좋습니다.
• 대안 도구와 비교 — 아주 큰 데이터에서 속도가 급하면 bedops가 정렬 기반으로 더 빠를 때가 있고, R에서는 GenomicRanges가 같은 구간 연산을 객체로 제공합니다. 셸에서 빠르게 훑기엔 bedtools가 여전히 표준이죠.

다음 학습

🎓 이어서 읽기
• (정렬 결과 포맷) SAM/BAM 파일이란? — bedtools가 입력으로 받는 BAM의 구조
• (변이 포맷) VCF 파일이란? — intersect로 걸러 낸 변이가 담기는 표준 포맷
• (정렬 원리) 서열 정렬이란? — 구간의 출발점인 리드 매핑과 좌표
• (BAM·VCF 실전) samtools·bcftools 실전 — 파이프로 bedtools와 잇는 짝 도구
• (구간의 사례) ChIP-seq란? — 피크라는 구간이 어떻게 만들어지나

References

1. Quinlan, A. R., & Hall, I. M. (2010). BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics, 26(6), 841–842. doi:10.1093/bioinformatics/btq033 (bedtools 원논문)
2. Quinlan, A. R. (2014). BEDTools: The Swiss-Army Tool for Genome Feature Analysis. Current Protocols in Bioinformatics, 47, 11.12.1–11.12.34. doi:10.1002/0471250953.bi1112s47
3. Quinlan lab. (2023). bedtools documentation (v2.31.0) — a powerful toolset for genome arithmetic. bedtools.readthedocs.io
4. UCSC Genome Browser. Frequently Asked Questions: Data File Formats (BED · 0-based half-open 좌표). genome.ucsc.edu/FAQ/FAQformat.html

Pipette & Pipeline · A bio portfolio journal