ATAC-seq 데이터 분석: 정렬 후 처리 단계 (4)

  1. 데이터 처리 개요

이제 Greenleaf 연구팀의 ATAC-seq 쌍단-end 데이터를 처리했으니, 정렬된 결과를 분석해보자.

먼저 ATAC-seq 데이터의 예상 조각 길이 분포를 확인해야 한다. GenomicAlignments 패키지를 사용하여 새로 정렬된 데이터를 읽어들인다.

여기서는 정확하게 쌍을 이루는 읽기만 필요하므로, ScanBamParam()과 scanBamFlag() 함수를 활용하여 R로 읽어들일 내용을 제어한다.

scanBamFlag() 함수의 isProperPair 매개변수를 TRUE로 설정하여, 사전 정의된 최대 조각 길이(2000bp) 내에서 쌍을 이루는 읽기만 필터링한다.

library(GenomicAlignments)
pairFlags = scanBamFlag(isProperPair = TRUE)
  1. 정렬 데이터 읽기

이제 이 플래그를 ScanBamParam() 함수와 함께 사용하여 정확히 쌍을 이루는 읽기만 읽어들인다.

또한 what 매개변수를 통해 R로 읽어들일 정보를 지정한다. 특히 삽입 크기 정보인 isize를 지정하는 것이 중요하다. 메모리 사용량을 줄이기 위해 GRanges 객체 매개변수를 지정하여 20번染色体에서만 정보를 읽어들인다.

paramObj = ScanBamParam(flag = pairFlags, what = c("qname", "mapq", "isize"), which = GRanges("chr20",
    IRanges(1, 63025520)))
paramObj

ScanBamParam 객체를 설정했으니, readGAlignmentPairs() 함수를 사용하여 쌍단-end ATAC-seq 데이터를 읽어들인다. 이 함수는 readGAlignments() 함수로 단일-end ChIP-seq 데이터를 읽는 방식과 유사하게 동작한다. 결과는 GAlignmentPairs 객체이다.

atacData <- readGAlignmentPairs(sortedBAM, param = paramObj)
class(atacData)
  1. GAlignmentPairs 구조

GAlignmentPairs 객체에는 쌍을 이루는 읽기에 대한 정보가 포함되어 있다. 각 읽기 쌍의 정보가 병렬 GAlignments 객체에 쌍으로 저장된다.

atacData[1:2, ]

first()와 second() 접근자 함수를 사용하여 GAlignments 객체에 접근할 수 있다. 각각 첫 번째 읽기 또는 두 번째 읽기에 대한 정보를 얻을 수 있다.

firstRead <- first(atacData)
secondRead <- second(atacData)
secondRead[1, ]
  1. 매핑 품질 점수 분석

가장 먼저 수행할 수 있는 작업은 firstRead와 secondRead의 MapQ(매핑 품질) 점수 분포를 구하는 것이다. mcols() 함수를 사용하여 각 읽기의 GAalignments 객체 mapq 슬롯에 접근한다.

firstReadMapQ <- mcols(firstRead)$mapq
secondReadMapQ <- mcols(secondRead)$mapq
firstReadMapQ[1:2]
  1. 매핑 품질 빈도 계산

table() 함수를 사용하여 각 쌍 읽기의 점수 빈도를 요약할 수 있다.

firstReadMapQFreqs <- table(firstReadMapQ)
secondReadMapQFreqs <- table(secondReadMapQ)
firstReadMapQFreqs

secondReadMapQFreqs
  1. 시각화

마지막으로 ggplot2를 사용하여 각 쌍 읽기의 MapQ 분포를 시각화한다.

library(ggplot2)
plotData <- data.frame(MapQ = c(names(firstReadMapQFreqs), names(secondReadMapQFreqs)), Frequency = c(firstReadMapQFreqs,
    secondReadMapQFreqs), Read = c(rep("Read1", length(firstReadMapQFreqs)), rep("Read2",
    length(secondReadMapQFreqs))))
plotData$MapQ <- factor(plotData$MapQ, levels = unique(sort(as.numeric(plotData$MapQ))))
ggplot(plotData, aes(x = MapQ, y = Frequency, fill = MapQ)) + geom_bar(stat = "identity") +
    facet_grid(~Read)
  1. 삽입 크기 추출

이제 쌍을 이루는 정렬 데이터를 R로 읽어들였으니, 각 읽기 쌍에 연결된 GAlignments 객체의 elementMetadata()에서 삽입 크기를 검색할 수 있다. 올바르게 쌍을 이루는 읽기는 동일한 삽입 크기 길이를 가지므로, firstRead에서 삽입 크기를 추출한다.

atacData_first <- first(atacData)
insertLengths <- abs(elementMetadata(atacData_first)$isize)
head(insertLengths)
  1. 삽입 크기 분포 시각화

ATAC-seq는 무핵소체, 단핵소체 및 다핵소체 부분에 해당하는 조각 길이의 혼합물을 나타내야 한다. table() 함수를 사용하여 각 조각 길이 발생 빈도를 가진 벡터를 검색한다.

fragmentSizes <- table(insertLengths)
fragmentSizes[1:5]

이제 20번染色체의 삽입 길이를 사용하여 모든 조각 길이의 분포를 시각화할 수 있다.

library(ggplot2)
plotData <- data.frame(InsertSize = as.numeric(names(fragmentSizes)), Count = as.numeric(fragmentSizes))
fragmentPlot <- ggplot(plotData, aes(x = InsertSize, y = Count)) + geom_line()
fragmentPlot + theme_bw()

핵소체 패턴을 명확히 하기 위해 카운트에 log2 변환을 적용할 수 있다.

fragmentPlot + scale_y_continuous(trans = "log2") + theme_bw()

이제 Greenleaf 연구에서처럼 무핵소체(<100bp), 단핵소체(180bp-247bp) 및 이중핵소체(315-437)를 주석으로 표시할 수 있다.

fragmentPlot + scale_y_continuous(trans = "log2") + geom_vline(xintercept = c(180,
    247), colour = "red") + geom_vline(xintercept = c(315, 437), colour = "darkblue") +
    geom_vline(xintercept = c(100), colour = "darkgreen") + theme_bw()

7월 28일 01:57에 게시됨

인기 태그

©2026 괴물 클럽