ChIP-seq 데이터 분석 과정에서 여러 샘플 또는 반복실험 간의 피크를 비교하고, 공통되거나 특정 조건에만 나타나는 피크를 식별하는 것은 매우 중요합니다. 이 글에서는 MACS2로 호출된 피크 데이터를 R 환경으로 가져와 GenomicRanges 객체로 처리하고, 최종적으로 합의(consensus) 피크를 정의하는 방법을 다룹니다.
1. 피크 데이터 로드 및 초기 처리
먼단, MACS2에서 생성된 피크 파일을 R로 불러와야 합니다. 여기서는 'data/peaks/' 디렉토리에 있는 '.peaks' 확장자를 가진 모든 파일을 찾아 목록화합니다.
peak_파일_경로 <- list.files("data/peaks/", pattern = "*.peaks", full.names = TRUE)
print(peak_파일_경로)
이제 이 파일들을 순회하면서 각각의 탭으로 구분된 피크 파일을 데이터 프레임 형태로 R로 가져옵니다. MACS2 출력 파일 중 '.xls' 확장자가 실제로는 탭으로 구분된 텍스트 파일인 경우가 많습니다.
macs_피크_DF_목록 <- lapply(peak_파일_경로, function(path) read.delim(path, comment.char = "#"))
cat("로드된 데이터 프레임 수:", length(macs_피크_DF_목록), "\n")
로드된 데이터 프레임 목록을 바탕으로, 각 피크 세트를 GenomicRanges (GRanges) 객체로 변환합니다. GRanges 객체는 유전체 데이터를 다루는 데 최적화된 Bioconductor의 핵심 데이터 구조입니다.
library(GenomicRanges)
macs_피크_GR_목록 <- lapply(macs_피크_DF_목록, function(df) {
GRanges(seqnames = df[, "chr"], IRanges(df[, "start"], df[, "end"]))
})
print("첫 번째 GRanges 객체:")
print(macs_피크_GR_목록[[1]])
각 피크 세트에 의미 있는 이름을 부여하기 위해 파일 이름에서 확장자를 제거하여 샘플 이름을 생성합니다. basename() 함수는 전체 경로에서 파일 이름만 추출하고, gsub() 함수는 특정 패턴을 다른 문자열로 대체합니다.
원본_파일_이름 <- basename(peak_파일_경로)
print("원본 파일 이름:")
print(원본_파일_이름)
샘플_이름 <- gsub("_peaks.xls", "", 원본_파일_이름)
print("추출된 샘플 이름:")
print(샘플_이름)
이제 GRanges 객체의 목록을 GRangesList() 함수를 사용하여 하나의 GRangesList 객체로 변환하고, 위에서 생성한 샘플 이름을 할당합니다.
macs_피크_GRL <- GRangesList(macs_피크_GR_목록)
names(macs_피크_GRL) <- 샘플_이름
cat("GRangesList 객체 클래스:", class(macs_피크_GRL), "\n")
cat("GRangesList 객체 이름:", names(macs_피크_GRL), "\n")
2. GRangesList 객체 활용
GRangesList는 표준 R 리스트와 유사하게 동작하면서도 GRanges 객체에 특화된 여러 함수를 직접 적용할 수 있는 장점이 있습니다. 예를 들어, 각 샘플의 피크 수를 확인하려면 lengths() 함수를 사용할 수 있습니다.
cat("각 샘플별 피크 수:\n")
print(lengths(macs_피크_GRL))
GRangesList 객체는 GRanges 접근자 및 연산자 함수를 직접 적용할 수 있습니다. 예를 들어, 피크의 너비를 조절하는 resize() 함수를 GRangesList에 직접 적용할 수 있습니다.
library(rtracklayer) # resize 함수는 rtracklayer 패키지에서 제공될 수 있습니다.
# (GenomicRanges 패키지에 포함되어 있으므로 rtracklayer는 필수는 아닙니다.)
중심_조정된_피크_GRL <- resize(macs_피크_GRL, width = 10, fix = "center")
cat("중심 조정된 피크의 너비 (첫 번째 샘플의 첫 번째 피크):\n")
print(width(중심_조정된_피크_GRL[[1]])[1])
이제 특정 샘플의 피크를 추출하여 개수를 확인해 봅니다. 예를 들어, 'Mel_1' 및 'Mel_2' 반복실험의 피크를 분리할 수 있습니다.
멜_반복1_피크 <- macs_피크_GRL$Mel_1
멜_반복2_피크 <- macs_피크_GRL$Mel_2
cat("Mel_1 피크 수:", length(멜_반복1_피크), "\n")
cat("Mel_2 피크 수:", length(멜_반복2_피크), "\n")
3. 고유(Unique) 및 공통(Common) 피크 식별
두 반복실험에서 서로 겹치지 않는 고유한 피크를 식별할 수 있습니다. %over% 연산자는 두 GRanges 객체 간의 중첩 여부를 확인합니다.
멜_반복1_고유_피크 <- 멜_반복1_피크[!멜_반복1_피크 %over% 멜_반복2_피크]
멜_반복2_고유_피크 <- 멜_반복2_피크[!멜_반복2_피크 %over% 멜_반복1_피크]
cat("Mel_1에만 존재하는 고유 피크 수:", length(멜_반복1_고유_피크), "\n")
cat("Mel_2에만 존재하는 고유 피크 수:", length(멜_반복2_고유_피크), "\n")
# 식별된 고유 피크를 BED 파일로 저장
export.bed(멜_반복1_고유_피크, "Mel_1_Unique_peaks.bed")
export.bed(멜_반복2_고유_피크, "Mel_2_Unique_peaks.bed")
마찬가지로, 두 반복실험에 모두 나타나는 공통 피크도 식별할 수 있습니다.
멜_반복1_공통_피크 <- 멜_반복1_피크[멜_반복1_피크 %over% 멜_반복2_피크]
멜_반복2_공통_피크 <- 멜_반복2_피크[멜_반복2_피크 %over% 멜_반복1_피크]
cat("Mel_1과 Mel_2에 중복되는 Mel_1 피크 수:", length(멜_반복1_공통_피크), "\n")
cat("Mel_1과 Mel_2에 중복되는 Mel_2 피크 수:", length(멜_반복2_공통_피크), "\n")
# 공통 피크를 BED 파일로 저장 (두 목록의 길이가 다를 수 있음)
export.bed(멜_반복1_공통_피크, "Mel_1_Common_peaks.bed")
export.bed(멜_반복2_공통_피크, "Mel_2_Common_peaks.bed")
위 결과에서 공통 피크의 개수가 반복1과 반복2에서 다르게 나타날 수 있습니다. 이는 하나의 샘플에서 두 개의 피크가 다른 샘플의 하나의 피크와 겹칠 수 있기 때문입니다. 진정한 '합의' 피크를 정의하기 위해서는 다른 접근 방식이 필요합니다.
4. 비중복 통합 피크 및 합의 피크 정의
모든 샘플에서 정의되는 비중복 피크 세트를 생성하는 것이 일반적입니다. 먼저, 모든 샘플의 피크를 하나의 GRanges 객체로 통합하여 중복되는 피크를 포함하는 집합을 만듭니다.
모든_피크_중복_집합 <- unlist(macs_피크_GRL)
print("모든 피크 통합 (중복 포함):")
print(모든_피크_중복_집합)
다음으로, reduce() 함수를 사용하여 중복되는 피크 영역을 병합하여 비중복적이고 고유한 피크 집합을 생성합니다. 이 집합은 어떤 샘플에서든 피크가 존재하는 모든 영역을 대표합니다.
비중복_통합_피크 <- reduce(모든_피크_중복_집합)
print("비중복 통합 피크 집합:")
print(비중복_통합_피크)
export.bed(비중복_통합_피크, "all_non_redundant_peaks.bed")
이 비중복 통합 피크 세트를 기반으로, 특정 조건(예: 두 반복실험 모두)에서 나타나는 합의 피크를 정의할 수 있습니다.
합의_피크 <- 비중복_통합_피크[비중복_통합_피크 %over% 멜_반복1_피크 & 비중복_통합_피크 %over% 멜_반복2_피크]
print("Mel_1과 Mel_2에 모두 존재하는 합의 피크:")
print(합의_피크)
export.bed(합의_피크, "consensus_peaks.bed")
마찬가지로, 특정 반복실험에만 고유하게 나타나는 피크를 비중복 통합 피크 세트에서 식별할 수 있습니다.
멜1_전용_피크 <- 비중복_통합_피크[비중복_통합_피크 %over% 멜_반복1_피크 & !비중복_통합_피크 %over% 멜_반복2_피크]
멜2_전용_피크 <- 비중복_통합_피크[!비중복_통합_피크 %over% 멜_반복1_피크 & 비중복_통합_피크 %over% 멜_반복2_피크]
cat("Mel_1에만 존재하는 비중복 피크 수:", length(멜1_전용_피크), "\n")
cat("Mel_2에만 존재하는 비중복 피크 수:", length(멜2_전용_피크), "\n")
export.bed(멜1_전용_피크, "mel1_only_peaks.bed")
export.bed(멜2_전용_피크, "mel2_only_peaks.bed")
5. 복합적인 피크 중복 분석
많은 수의 샘플을 다룰 때는 비중복 통합 피크가 각 샘플에 존재하는지 여부를 나타내는 논리 행렬을 생성하는 것이 유용합니다. 먼저, 각 샘플에 대해 비중복 통합 피크가 중첩되는지 여부를 나타내는 논리 벡터 목록을 생성합니다.
중복_여부_목록 <- lapply(macs_피크_GRL, function(gr) 비중복_통합_피크 %over% gr)
cat("첫 번째 샘플의 중복 여부 (상위 2개):\n")
print(중복_여부_목록[[1]][1:2])
이 논리 벡터 목록을 do.call(cbind, ...)를 사용하여 열로 묶어 중복 행렬을 만듭니다. 이 행렬의 각 열은 샘플을, 각 행은 비중복 통합 피크를 나타내며, 해당 피크가 샘플에 존재하는지(TRUE) 여부를 표시합니다.
중복_행렬 <- do.call(cbind, 중복_여부_목록)
colnames(중복_행렬) <- names(macs_피크_GRL)
cat("중복 행렬 (상위 2행):\n")
print(중복_행렬[1:2, ])
이 중복 행렬을 비중복 통합 피크 GRanges 객체의 메타데이터 열(mcols())에 추가합니다. 이제 각 비중복 피크에 대해 어떤 샘플에 존재하는지에 대한 정보를 직접 확인할 수 있습니다.
mcols(비중복_통합_피크) <- 중복_행렬
cat("메타데이터가 추가된 비중복 통합 피크 (상위 2개):\n")
print(비중복_통합_피크[1:2, ])
limma 패키지는 RNA-seq 및 마이크로어레이 데이터 분석에 자주 사용되지만, 논리 행렬의 중복을 시각화하는 데 유용한 vennDiagram() 함수도 제공합니다. 이 함수를 사용하여 샘플 간의 피크 중복 관계를 벤 다이어그램으로 나타낼 수 있습니다.
library(limma)
# (이 코드를 실행하려면 X11 또는 기타 그래픽 장치가 필요할 수 있습니다.)
# vennDiagram(mcols(비중복_통합_피크))
limma 패키지의 vennCounts() 함수는 벤 다이어그램에 표시되는 각 영역의 피크 수를 데이터 프레임 형태로 반환합니다.
cat("벤 다이어그램 카운트:\n")
print(vennCounts(mcols(비중복_통합_피크)))
6. 고신뢰도 피크 식별
비중복 통합 피크 세트와 피크 존재 행렬을 활용하여, 특정 조건 내에서 반복적으로 나타나는 '고신뢰도' 피크를 정의할 수 있습니다. 예를 들어, 'ch12' 샘플의 두 반복실험 모두에 나타나는 피크를 추출합니다. 논리 행렬은 1(TRUE) 또는 0(FALSE)으로 간주될 수 있으므로, rowSums() 함수를 사용하여 특정 반복실험에서 피크가 나타나는 횟수를 계산할 수 있습니다.
ch12_고신뢰도_피크 <- 비중복_통합_피크[rowSums(as.data.frame(mcols(비중복_통합_피크)[, c("ch12_1", "ch12_2")])) >= 2]
export.bed(ch12_고신뢰도_피크, "ch12_high_confidence_peaks.bed")
cat("Ch12의 고신뢰도 피크 (상위 2개):\n")
print(ch12_고신뢰도_피크[1:2, ])
또한, 'Ch12'에서는 고신뢰도로 나타나지만 'Mel' 샘플에서는 전혀 나타나지 않는 피크와 같이 더 복합적인 조건도 정의할 수 있습니다.
ch12_고신뢰도_전용_피크 <- 비중복_통합_피크[rowSums(as.data.frame(mcols(비중복_통합_피크)[, c("ch12_1", "ch12_2")])) >= 2 &
rowSums(as.data.frame(mcols(비중복_통합_피크)[, c("Mel_1", "Mel_2")])) == 0]
export.bed(ch12_고신뢰도_전용_피크, "ch12_high_confidence_unique_peaks.bed")
cat("Ch12 고신뢰도 전용 피크 (첫 번째):\n")
print(ch12_고신뢰도_전용_피크[1, ])