Rolling Join, чтобы найти геномные регионы - PullRequest
Новая
       51

Rolling Join, чтобы найти геномные регионы

0 голосов
/ 11 декабря 2018

Быстрая версия - это то, что я хочу найти все совпадения на определенном расстоянии от определенных значений.Пример: 

library(data.table)
reads <- structure(list(seqnames = c(1, 1, 1, 1, 1, 1, 1), start = c(100, 
100, 100, 130, 130, 132, 133), end = c(130, 130, 131, 160, 160, 
155, 160), strand = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "+", class = "factor")), row.names = c(NA, 
-7L), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x15af570>, sorted = c("start", 
"end", "strand"))

# take only the end position
p3 <- reads[, list(N3 = .N), by = c("seqnames", "end", "strand")]
setnames(p3, "end", "loc.3p")

# take only the start position
p5 <- reads[, list(N5 = .N), by = c("seqnames", "start", "strand")]
setnames(p5, "start", "loc.5p")

# find the start positions close to the end positions
p3[, loc := loc.3p]
p5[, loc := loc.5p]

p3[p5, roll=50, rollends=c(FALSE, TRUE), nomatch=0, on = c("seqnames", "strand", "loc")]

   seqnames loc.3p strand N3 loc loc.5p N5
1:        1    130      +  2 130    130  2
2:        1    131      +  1 132    132  1
3:        1    131      +  1 133    133  1

Это работает, но возвращается только первое совпадение.Моя цель сейчас найти все совпадения (в пределах указанного диапазона 50).Ожидаемый результат: 

   seqnames loc.3p strand N3 loc loc.5p N5
1:        1    130      +  2 130    130  2
2:        1    130      +  2 130    132  1
3:        1    130      +  2 133    133  1
4:        1    131      +  1 132    132  1
4:        1    131      +  1 133    133  1

Обратите внимание, что loc.3p 130 также должно совпадать с loc.5p 132 и 132.

Как это сделать?Мне нужно сохранить оценки для следующего набора операций.

Биоинформатическая версия , я пытаюсь найти все последующие 5 'чтения до определенного расстояния (в том жецепь).В этом примере я использую только «+», но мне также придется сделать это для «-».Поскольку это будет сделано для миллионов операций чтения, data.table представляется целесообразным.Я изучил GenomicRanges, но хранение метаданных для обоих наборов данных (5 'и 3' позиций) немного запутано.

1 Ответ

0 голосов
/ 11 декабря 2018

Когда вы говорите

в пределах указанного диапазона 50

Это заставляет меня думать, что соединение не равно.В будущем может оказаться полезным, если вы покажете желаемый результат:)

небольшое изменение p3, обратите внимание, что я пропустил создание столбца loc, который у вас был. 

p3[, "upper" := 50 + loc.3p]

Non-equi join: 

p3[p5, .(seqnames, x.loc.3p, strand,  N3, loc.5p, N5),
   on = .(seqnames, strand, loc.3p <= loc.5p, upper >= loc.5p), nomatch = 0
   ][order(x.loc.3p, loc.5p),
     .(seqnames, loc.3p = x.loc.3p, strand, N3, loc.5p, N5)
     ]

   seqnames loc.3p strand N3 loc.5p N5
1:        1    130      +  2    130  2
2:        1    130      +  2    132  1
3:        1    130      +  2    133  1
4:        1    131      +  1    132  1
5:        1    131      +  1    133  1

Последняя скобка для упорядочивания и переименования, я не совсем уверен, почему мне пришлось указать x.loc.3p в первых скобках, но я получил ошибку без нее...

...