2013-12-13 7 views
1

У меня есть список диапазонов геномных областей, полученных из серии пациентов.Поиск перекрывающих областей хромосомы с помощью GRanges

> head(dotoo) 
GRanges with 6 ranges and 3 metadata columns: 
    seqnames     ranges strand |  Id  CN Histology 
     <Rle>    <IRanges> <Rle> | <factor> <factor> <factor> 
[1]  3 [167946693, 168005541]  * |  9  3  MD 
[2]  3 [189907623, 189954633]  * |  9  3  MD 
[3]  6 [132274121, 132384438]  * |  9  3  MD 
[4]  11 [ 67685096, 70138399]  * |  9  4  MD 
[5]  12 [ 53859037, 53927595]  * |  9  3  MD 
[6]  15 [ 19830049, 20089383]  * |  9  1  MD 

Когда я сюжет геномные аберраций с помощью

autoplot(dotoo, aes(fill=as.factor(Id), color=as.factor(Id)))

Я вижу много перекрывающихся областей, см изображения

enter image description here

Как узнать, какие регионы перекрывается между, по меньшей мере, 3 пациента, и поделились CN?

В принципе, если вы посмотрите на изображение, как мне найти регионы, которые «стекают», и только часть, которая является общей? Есть ли вообще способ?

+0

Вы пробовали 'findOverlaps()'? – nacnudus

+0

У меня есть, и это не то, что мне нужно. –

ответ

3

Получить список регионов «непересекающиеся» (может быть, это не то, что вы хотите? Другие варианты reduce и только оригинальный dotoo объект без этого шага на всех)

d = disjoint(dotoo) 

находка перекрывается между оригиналом и каждый из непересекающихся областей

olap = findOverlaps(query=dotoo, subject=d) 

разделить индекс в наложения на основе на основе субъекта и CN

splt = split(seq_along(olap), list(subjectHits(olap), dotoo$CN[queryHits(olap)])) 

фильтра это для тех, удовлетворяющих условия вашего

filt = Filter(function(x) length(x) >= 3, splt) 

filt теперь список индексов в olap. Вы можете создать GRangesList с перекрывающимися элементами, как

idx = unlist(filt) 
grp = rep(seq_along(filt), sapply(filt, length)) 
splitAsList(dotoo[queryHits(olap)[idx]], grp) 

Задавайте вопросы о Bioconductor пакетов на Bioconductor mailing list (без абонентской платы не требуется).