Получите самые длинные диапазоны за последующие имена - PullRequest
0 голосов
/ 20 февраля 2020

У меня есть GRanges объект с разными геномными c диапазонами для каждого seqnames (например, хромосомы).
Как я могу получить GRanges, содержащий только самый длинный диапазон для каждого seqname / хромосомы?

Например, если gr является GRanges:

library(GenomicRanges)

# Make a GRanges object
set.seed(123)
gr <- GRanges(seqnames = rep(c("chr1", "chr2", "chr3"), times=2:4),
              ranges = IRanges(start=sample.int(10000, 9), 
                               width = c(3,5,50,20,10,500,100,500,200)))

# Add a column with the width for clarity:
mcols(gr)$width <- width(gr)

gr
#GRanges object with 9 ranges and 1 metadata column:
#      seqnames    ranges strand |     width
#         <Rle> <IRanges>  <Rle> | <integer>
#  [1]     chr1 2463-2465      * |         3
#  [2]     chr1 2511-2515      * |         5
#  [3]     chr2 8718-8767      * |        50
#  [4]     chr2 2986-3005      * |        20
#  [5]     chr2 1842-1851      * |        10
#  [6]     chr3 9334-9833      * |       500
#  [7]     chr3 3371-3470      * |       100
#  [8]     chr3 4761-5260      * |       500
#  [9]     chr3 6746-6945      * |       200
#  -------
#  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Тогда я хочу получить следующее GRanges:

#GRanges object with 3 ranges and 1 metadata column:
#    seqnames      ranges strand |     width
#       <Rle>   <IRanges>  <Rle> | <integer>
#  [1]     chr1 2511-2515      * |         5
#  [2]     chr2 8718-8767      * |        50
#  [3]     chr3 9334-9833      * |       500

Для моего Приложение Я в порядке, получив только первый самый длинный диапазон для chr3, но я был бы признателен за решение, которое также может выбрать все связи, если таковые имеются.

1 Ответ

0 голосов
/ 20 февраля 2020

Я нашел это решение:

library(GenomicRanges)

# Make the GRanges object
set.seed(123)
gr <- GRanges(seqnames = rep(c("chr1", "chr2", "chr3"), times=2:4),
              ranges = IRanges(start=sample.int(10000, 9), 
                               width = c(3,5,50,20,10,500,100,500,200)))

# Add a column with the width for clarity:
mcols(gr)$width <- width(gr)

# split gr by seqnames
grl <- split(gr, seqnames(gr))

# Get the width of each range organized as an IntegerList
wgr <- width(grl)
# if gr is ordered by seqnames this is equivalent to:
wgr <- splitAsList(width(gr), seqnames(gr))
wgr
#IntegerList of length 3
#[["chr1"]] 3 5
#[["chr2"]] 50 20 10
#[["chr3"]] 500 100 500 200

# Get the ranges with the largest width
# See ?'which.max,NumericList-method' for the global argument
longestRanges <- gr[which.max(wgr, global = TRUE)]

longestRanges
#GRanges object with 3 ranges and 1 metadata column:
#      seqnames      ranges strand |     width
#         <Rle>   <IRanges>  <Rle> | <integer>
#  [1]     chr1 2511-2515      * |         5
#  [2]     chr2 8718-8767      * |        50
#  [3]     chr3 9334-9833      * |       500
#  -------
#  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Обратите внимание, что это решение работает, только если диапазоны в gr изначально сгруппированы / упорядочены по seqnames.
Если это не так, это необходимо начать с сортировки gr с gr <- sort(gr)

Чтобы получить все связи, можно использовать:

longestRanges <- unlist(grl[which(wgr == max(wgr))])

longestRanges
#GRanges object with 4 ranges and 1 metadata column:
#       seqnames    ranges strand |     width
#          <Rle> <IRanges>  <Rle> | <integer>
#  chr1     chr1 2511-2515      * |         5
#  chr2     chr2 8718-8767      * |        50
#  chr3     chr3 9334-9833      * |       500
#  chr3     chr3 4761-5260      * |       500
#  -------
#  seqinfo: 3 sequences from an unspecified genome; no seqlengths
...