GenomicRanges добавить покрытие

Question

Я работаю на данных Секвенирование РНК и пытается построить средние профили покрытия в зависимости от генотипа, подобно тому, как это делается здесь

РНК сл Расход на генотип (Источник: Pickrell и др, Nature, 2010)

Чтобы сделать этот сюжет, у меня есть шишка файлы из 100 человек, которые содержат информацию о зоне от РНК-сл данных (в конкретном регионе), и что я прочитал в R, так как GenomicRanges объектов.

Это дает мне Granges объекты, такие как те, полученное в следующем примере: игрушечной

GR1 = Granges (seqname = 1, диапазон = IRanges (Start = с (1,5,10,15,30 , 55), конец = с (4,9,14,29,39,60)))

gr1 $ COV = с (3,1,8,6,2,10)

gr2 = Границы (seqname = 1, range = IRanges (start = c (3,20,24), end = c (7,23,26)))

gr2 $ cov = c (3,5,3)

старт = (уникальные сортировки (C (диапазоны (gr1) @ начать, диапазоны (gr2) @start)))

gr1

GRanges object with 6 ranges and 1 metadata column: 
seqnames ranges strand |  cov 
    <Rle> <IRanges> <Rle> | <numeric> 
     1 [ 1, 4]  * |   3 
     1 [ 5, 9]  * |   1 
     1 [10, 14]  * |   8 
     1 [15, 29]  * |   6 
     1 [30, 39]  * |   2 
     1 [55, 60]  * |  10 
     ------- 
seqinfo: 1 sequence from an unspecified genome; no seqlengths 

gr2

GRanges object with 3 ranges and 1 metadata column: 
seqnames ranges strand |  cov 
    <Rle> <IRanges> <Rle> | <numeric> 
     1 [ 3, 7]  * |   3 
     1 [20, 23]  * |   5 
     1 [24, 26]  * |   3 
     ------- 
seqinfo: 1 sequence from an unspecified genome; no seqlengths 

Проблема в том, что у меня есть эти данные для каждого человека (gr1 и gr2 были бы двумя разными людьми), и я хотел бы объединить их для создания объект геномных диапазонов, что дает мне полное покрытие в каждом положении между индивидами, 1 и 2 , которые будут выглядеть следующим образом:

gr3

GRanges object with 6 ranges and 1 metadata column: 
seqnames ranges strand |  cov 
    <Rle> <IRanges> <Rle> | <numeric> 
     1 [ 1, 2]  * |   3 
     1 [ 3, 4]  * |   6 (=3+3) 
     1 [ 5, 7]  * |   4 (=1+3) 
     1 [ 8, 9]  * |   1 
     1 [10, 14]  * |   8 
     1 [15, 19]  * |   6 
     1 [20, 23]  * |   11 (=6+5) 
     1 [24, 26]  * |   9 (=6+3) 
     1 [27, 29]  * |   6 
     1 [30, 39]  * |   2 
     1 [55, 60]  * |  10 

Кто-нибудь знает простой способ сделать это ? или я обречен?

Спасибо за ваши ответы.

PS: мои данные не застряли, но если у вас есть данные о мельничных данных, это еще лучше.

PPS: В идеале я также хотел бы иметь возможность умножения или применять любую функцию с двумя аргументами x и y вместо простого добавления покрытия.

Answer 1

Прошло почти год, но вот мой ответ для дальнейшего использования.

Всякий раз, когда я не нахожу функцию для выполнения такой задачи, как эта, я просто расширяю объекты GRanges до разрешения одной-двух точек. Это позволяет мне выполнять любую требуемую операцию в столбцах метаданных, обрабатывая их как простые столбцы data.frame, так как IRanges теперь сопоставляются между двумя объектами Granges.

В конкретном случае этого вопроса выполняются следующие работы.

### Sort seqlevels 
# (not necessary here, but in real world examples, 
# with multiple sequences, you will want to do this) 
gr1 <- sort(GenomeInfoDb::sortSeqlevels(gr1)) 
gr2 <- sort(GenomeInfoDb::sortSeqlevels(gr2)) 

### Add seqlengths 
# (this corresponds to the actual sequence lengths; 
# here we use the highest position between the two objects: 60) 
seqlengths(gr1) <- 60 

### Make 1-bp tiles covering the genome 
# (using either one of gr1 and gr2 as a reference) 
bins <- GenomicRanges::tileGenome(GenomeInfoDb::seqlengths(gr1), 
            tilewidth=1, 
            cut.last.tile.in.chrom=TRUE) 

### Get coverage signal as Rle object 
gr1_cov <- coverage(gr1, weight="cov") 
gr2_cov <- coverage(gr2, weight="cov") 

### Get average coverage in each bin 
# (since the bins are 1-bp wide, this just keeps the original coverage value) 
gr1_bins <- GenomicRanges::binnedAverage(bins, gr1_cov, "binned_cov") 
gr2_bins <- GenomicRanges::binnedAverage(bins, gr2_cov, "binned_cov") 

### Make final object: 
# We can now sum the values in the metadata columns 
# Addressing the PPS, you could do any other operation or apply a function 
gr3 <- gr1_bins 
gr3$binned_cov <- gr1_bins$binned_cov + gr2_bins$binned_cov 

Это дает окончательный GRanges объект на одной п.н. разрешением.

> gr3 

GRanges object with 60 ranges and 1 metadata column: 
    seqnames ranges strand | binned_cov 
     <Rle> <IRanges> <Rle> | <numeric> 
[1]  1 [1, 1]  * |   3 
[2]  1 [2, 2]  * |   3 
[3]  1 [3, 3]  * |   6 
[4]  1 [4, 4]  * |   6 
[5]  1 [5, 5]  * |   4 
...  ...  ... ... .  ... 
[56]  1 [56, 56]  * |   10 
[57]  1 [57, 57]  * |   10 
[58]  1 [58, 58]  * |   10 
[59]  1 [59, 59]  * |   10 
[60]  1 [60, 60]  * |   10 
------- 
seqinfo: 1 sequence from an unspecified genome 

Чтобы сжать его и получить точную gr3 в вопросе, мы можем сделать следующее.

### Compress back to variable-width IRanges (by cov) 
gr3_Rle <- coverage(gr3, weight='binned_cov') 
gr3 <- as(gr3_Rle, "GRanges") 

### Drop 0-score rows 
gr3 <- gr3[gr3$score > 0] 

### Rename metadata column 
names(mcols(gr3)) <- 'cov' 

> gr3 

GRanges object with 11 ranges and 1 metadata column: 
     seqnames ranges strand |  cov 
      <Rle> <IRanges> <Rle> | <numeric> 
    [1]  1 [ 1, 2]  * |   3 
    [2]  1 [ 3, 4]  * |   6 
    [3]  1 [ 5, 7]  * |   4 
    [4]  1 [ 8, 9]  * |   1 
    [5]  1 [10, 14]  * |   8 
    [6]  1 [15, 19]  * |   6 
    [7]  1 [20, 23]  * |  11 
    [8]  1 [24, 26]  * |   9 
    [9]  1 [27, 29]  * |   6 
    [10]  1 [30, 39]  * |   2 
    [11]  1 [55, 60]  * |  10 
    ------- 
    seqinfo: 1 sequence from an unspecified genome