그레인이있는 여러 간격의 평균 신호가

Question

입니다. bedGraph 파일에서 GRanges 객체로 가져온 게놈 전체 ChIP-seq 신호가 있습니다. 모든 피크를 덮고있는 고정 폭 간격에 대한 평균 신호를 플로팅하고 싶습니다. 신호를 숫자 벡터로 추출하여 평균을 낼 수 있습니까?

예로서 고려 :

는

library(GenomicRanges) 
set.seed(1) 

signal <- GRanges(
    seqnames = Rle(c("chr1"), c(10)), 
    ranges = IRanges(1:10*10, end = 1:10*10+5), 
    score = runif(10)) 

intervals <- GRanges(
    seqnames = Rle(c("chr1"), c(5)), 
    ranges = IRanges(1:5*20 + floor(runif(5)*4), width = 10)) 

때문에 신호는 다음과 같습니다

GRanges with 10 ranges and 1 metadata column: 
     seqnames  ranges strand |    score 
      <Rle> <IRanges> <Rle> |   <numeric> 
    [1]  chr1 [ 10, 15]  * | 0.2655086631421 
    [2]  chr1 [ 20, 25]  * | 0.37212389963679 
    [3]  chr1 [ 30, 35]  * | 0.572853363351896 
    [4]  chr1 [ 40, 45]  * | 0.908207789994776 
    [5]  chr1 [ 50, 55]  * | 0.201681931037456 
    [6]  chr1 [ 60, 65]  * | 0.898389684967697 
    [7]  chr1 [ 70, 75]  * | 0.944675268605351 
    [8]  chr1 [ 80, 85]  * | 0.660797792486846 
    [9]  chr1 [ 90, 95]  * | 0.62911404389888 
    [10]  chr1 [100, 105]  * | 0.0617862704675645 
    --- 
    seqlengths: 
    chr1 
    NA 

와 간격는 다음과 같다 : 그래서 평균 싶습니다

GRanges with 5 ranges and 0 metadata columns: 
     seqnames  ranges strand 
     <Rle> <IRanges> <Rle> 
    [1]  chr1 [ 20, 29]  * 
    [2]  chr1 [ 40, 49]  * 
    [3]  chr1 [ 62, 71]  * 
    [4]  chr1 [ 81, 90]  * 
    [5]  chr1 [103, 112]  * 
    --- 
    seqlengths: 
    chr1 
    NA 

벡터 :

Rle(c(0.372, 0), c(6, 4))   # [ 20, 29] 
Rle(c(0.908, 0), c(6, 4))   # [ 40, 49] 
Rle(c(0.898, 0, 0.945), c(4, 4, 2)) # [ 62, 71] 
Rle(c(0.661, 0, 0.629), c(5, 4, 1)) # [ 81, 90] 
Rle(c(0.061, 0), c(3, 7))   # [103,112] 

루프 및 많은 지루한 오류가 발생하기 쉬운 간격 계산없이이 작업을 수행하려면 어떻게해야합니까? GenomicRanges 패키지에 이런 종류의 기능이 포함되기를 기대했지만 수동으로 볼 수는 없었습니다. 나는 subsetByOverlaps를 사용하려고 노력해 왔지만 이것은 신호 스코어를 결과로 가져 오는 것 같지 않고 위의 Rle 벡터를 추출하는 데 도움이되지 않습니다.

Answer 1

나는 그것을 알아 냈을 수도 있습니다. 나는 간격으로 각 범위에 아래의 getScores() 기능을 적용 할 수 있습니다. 기능이 대답 https://stackoverflow.com/a/9913411/959926에서 적응으로 findOverlaps 사용 : 지금까지 작동하는 것 같다 있지만 개선은 환영받을

getScores <- function(interval) { 
    scores <- Rle(0, width(interval)) 
    bases <- GRanges(
     seqnames = seqnames(interval), 
     ranges = IRanges(start(interval):end(interval), width = 1)) 
    overlaps <- findOverlaps(signal, bases) 
    scores[start(bases)[subjectHits(overlaps)] - start(interval) + 1] <- score(signal)[queryHits(overlaps)] 
    scores 
} 
Reduce('+', sapply(split(intervals, 1:length(intervals)), getScores))/length(intervals) 

. 예를 들어 신호 및/또는 간격이 길면 상당히 느립니다.

Answer 2

이 솔루션은 어떻게됩니까?

overlaps <- findOverlaps(signal, intervals) 
signal <- signal[overlaps@queryHits] 
averagedSignal <- aggregate(score(sites), list(overlaps@subjectHits), mean) 

Answer 3

overlaps <- findOverlaps(signal, intervals) 
sites <- signal[queryHits(overlaps)] 
intervals$averagedSignal <- aggregate(score(sites), list(subjectHits(overlaps)), mean)