피크 탐지 개요
전사 인자 Myc의 게놈 내 결합 부위를 규명하기 위해선 차등적으로 리드가 축적된 영역을 식별하는 피크 콜러(peak caller) 도구가 필요하다. 다양한 알고리즘이 존재하지만, 현재까지도 가장 신뢰성 있고 널리 사용되는 도구는 MACS2이다. MACS2는 다음 단계를 통해 유의미한 피크를 도출한다:
- 시퀀싱 리드로부터 프래그먼트 길이 추정
- 각 리드를 추정된 프래그먼트 중심으로 이동
- 대조군(입력 DNA 등) 대비 실험군에서 통계적으로 유의하게 축적된 영역 탐색
MACS2 설치 방법
MACS2는 R 패키지로 제공되지 않으며, 파이썬 기반 커맨드라인 도구이므로 일반적으로 Conda를 통해 설치한다. macOS 및 Linux 환경에서 가능하며 Windows는 공식 지원하지 않는다. Bioconductor에 포함된 Herper 패키지를 활용하면 R 내부에서 Conda 환경과 상호작용할 수 있다.
BiocManager::install("Herper")
library(Herper)
다음 명령어로 전용 환경에 MACS2를 설치할 수 있다:
install_CondaTools(tools = "macs2", env = "ChIPseq_analysis")
R 내부에서 MACS2 실행
Conda 환경이 준비되면, with_CondaEnv() 함수를 사용해 R 스크립트 내에서 직접 MACS2를 호출할 수 있다. 주요 입력은 다음과 같다:
-t: 처리 샘플 BAM 파일 (ChIP 샘플)-c: 대조 샘플 BAM 파일 (Input DNA)-n: 결과 출력 이름--outdir: 결과 저장 디렉터리
chip_bam <- "Sorted_Myc_MEL_1.bam"
input_bam <- "Sorted_Input_MEL.bam"
with_CondaEnv("ChIPseq_analysis",
system2("macs2",
args = c("callpeak",
"-t", chip_bam,
"-c", input_bam,
"-n", "Mel_Rep1",
"--outdir", "PeakDirectory")),
stdout = TRUE)
결과 파일 구조 이해
MACS2는 여러 출력 형식을 생성하며, 그 중 대표적인 것은 _peaks.xls와 _peaks.narrowPeak이다. XLS 파일은 실제 엑셀 파일이 아니라 탭 구분 텍스트이며, 앞부분에 명령행과 파라미터 정보가 주석 형태로 포함된다:
# Command line: callpeak -t Sorted_Myc_MEL_1.bam -n Mel1 -c Sorted_Input_MEL.bam
# effective genome size = 2.70e+09
# band width = 300
chr start end length fold_enrichment log10.pvalue name
chr1 4785371 4785642 272 5.33590 10.66553 Mel1_peak_1
피크 데이터 불러오기
R에서 XLS 형식 피크 파일을 로드할 때는 주석 행을 무시해야 하며, 이를 위해 comment.char = "#" 옵션을 설정한다:
peak_file <- "PeakDirectory/Mel_Rep1_peaks.xls"
peak_df <- read.delim(peak_file, comment.char = "#")
head(peak_df[, c("chr", "start", "end", "fold_enrichment", "log10.pvalue", "name")])
GRanges 객체로 변환
게놈 범위 기반 분석을 위해선 GenomicRanges 패키지의 GRanges 객체로 변환하는 것이 효율적이다:
library(GenomicRanges)
peak_gr <- GRanges(
seqnames = peak_df$chr,
ranges = IRanges(start = peak_df$start, end = peak_df$end)
)
메타데이터 추가
피크의 생물학적 중요도를 평가하기 위해 fold enrichment, p-value 등의 값을 메타데이터 열로 추가할 수 있다:
mcols(peak_gr) <- DataFrame(
summit = peak_df$abs_summit,
fold_change = peak_df$fold_enrichment,
pval = peak_df$X.log10.pvalue.,
qval = peak_df$X.log10.qvalue.
)
narrowPeak 형식 활용
.narrowPeak 파일은 UCSC에서 정의한 표준 형식으로, rtracklayer 패키지로 직접 임포트 가능하다:
library(rtracklayer)
narrow_peak_gr <- import("PeakDirectory/Mel_Rep1_peaks.narrowPeak", format = "narrowPeak")
이 형식은 signalValue, pValue, qValue, peak 위치 등 10개의 필드를 포함하며, 시각화 및 필터링에 용이하다.
블랙리스트 영역 제거
기술적 아티팩트가 발생하기 쉬운 게놈 영역(예: 반복 서열, 복제 기원)은 블랙리스트로 관리되며, 피크 분석 전 반드시 제외해야 한다. %over% 연산자를 사용해 중첩 여부를 판단하고 제거한다:
blacklist <- import.bed("hg38-blacklist.v2.bed")
filtered_peaks <- narrow_peak_gr[!overlapsAny(narrow_peak_gr, blacklist)]