GenomicRangesパッケージにおけるreduceとdisjoinの扱い #R

データセットの準備

library(GenomicRanges)
ir <- IRanges(c(1, 8, 14, 15, 19, 34, 40), width=c(12, 6, 6, 15, 6, 2, 7))
strand <- rep(c("+", "-"), c(4,3))
grngs <- GRanges(seqnames = "chr1", ranges = ir, strand = strand, seqlengths = c("chr1" = 50))

> grngs
GRanges with 7 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1  [ 1, 12]      +
  [2]     chr1  [ 8, 13]      +
  [3]     chr1  [14, 19]      +
  [4]     chr1  [15, 29]      +
  [5]     chr1  [19, 24]      -
  [6]     chr1  [34, 35]      -
  [7]     chr1  [40, 46]      -
  ---
  seqlengths:
   chr1
     50

青、赤の色分けはそれぞれゲノム上の+鎖、－鎖に対応させている。

reduce

＋鎖、－鎖に存在するGRangeオブジェクトを個別にそれぞれつなげてくれる。

disjoin

各GRangeオブジェクトについて、他のGRangeオブジェクトと重ならない部分を抽出する。

なお、上の図は http://www.bioconductor.org/help/course-materials/2009/GenentechNov2009/Module7/IRanges.R を参考に、以下のコードで描画しました。
条件分岐による色分けをするには、col= 引数にifelse()を与えてあげると便利です。

plotRanges <- function(x, xlim = x, main = deparse(substitute(x)), sep = 0.5, ...) 
{
  height <- 1
  if (is(xlim, "Ranges"))
    xlim <- c(min(start(xlim)), max(end(xlim)))
  bins <- disjointBins(IRanges(start(x), end(x) + 1))
  plot.new()
  plot.window(xlim, c(0, max(bins)*(height + sep)))
  ybottom <- bins * (sep + height) - height
  rect(start(x)-0.5, ybottom, end(x)+0.5, ybottom + height,
       col = ifelse(as.data.frame(x)$strand=="+", "blue", "red"), ...)
  title(main)
  axis(1)
}
png("qiita-20130404-default.png",width=640,height=320)
plotRanges(grngs,xlim=c(1,50))
dev.off()

png("qiita-20130404-reduce.png",width=640,height=320)
plotRanges(reduce(grngs),xlim=c(1,50), main="reduce")
dev.off()

png("qiita-20130404-disjoin.png",width=640,height=320)
plotRanges(disjoin(grngs),xlim=c(1,50), main="disjoin")
dev.off()