Найти все диапазоны за пределами определенного набора диапазонов

Мне интересно, что было бы лучшим способом определить все диапазоны, которые не охватываются данным набором диапазонов. Например, если у меня есть набор генов с известными координатами:

dtGenes <- fread(
  "id,start,end
 1,1000,1300
 2,1200,1500
 3,1600,2600
 4,3000,4000
")

Допустим, я знаю, что общая длина хромосомы (для простоты предположим, что они все находятся на одной хромосоме) составляет 10000. Итак, в конечном итоге я ожидаю получить следующий список межгенных областей:

"startR,endR
    0,1000
 1500,1600
 2600,3000
 4000,10000
"

может биокондуктор IRange быть полезным здесь? или есть другой хороший способ решить эту проблему?

2 ответа

Решение

Используя пакет Bioconductor GenomicRanges, преобразуйте свои исходные данные в GRanges

library(GenomicRanges)
gr <- with(dtGenes, GRanges("chr1", IRanges(start, end, names=id),
                            seqlengths=c(chr1=10000)))

Затем найдите промежутки между вашими генами

gaps <- gaps(gr)

GRanges знает о пряди. Вы не указали прядь в GRanges конструктор, так что нить была назначена *, Таким образом, в цепях +, - и * есть "пробелы", и вас интересуют только те, которые находятся в цепях *

> gaps[strand(gaps) == "*"]
GRanges with 4 ranges and 0 metadata columns:
      seqnames        ranges strand
         <Rle>     <IRanges>  <Rle>
  [1]     chr1 [   1,   999]      *
  [2]     chr1 [1501,  1599]      *
  [3]     chr1 [2601,  2999]      *
  [4]     chr1 [4001, 10000]      *
  ---
  seqlengths:
    chr1
   10000

Обратите внимание на соглашение Биокондуктора, что хромосомы начинаются с 1, и что диапазоны закрыты - start а также end координаты включены в диапазон. использование shift а также narrow на gr чтобы ваши диапазоны соответствовали соглашениям Bioconductor. Операции GRanges эффективны на десятках миллионов диапазонов.

Ты можешь использовать reduce от IRanges пакет

сначала уменьшите порядки диапазонов по x слева направо, затем объедините перекрывающиеся или смежные.

library(IRanges)
dat <- read.csv(text="id,start,end
 1,1000,1300
 2,1200,1500
 3,1600,2600
 4,3000,4000
")

ir <- IRanges(dat$start,dat$end)
rir <- reduce(ir)
IRanges of length 3
    start  end width
[1]  1000 1500   501
[2]  1600 2600  1001
[3]  3000 4000  1001
Другие вопросы по тегам