《Bioinformatics Data Skills》之使用IRanges包存储基因组范围

R有一系列的包可以帮助我们去了解与操作基因组，例如GenomicRanges，GenomicFeatures，Biostrings，BSgenome，rtracklayer等。
今天我们先学习最基本的，使用IRanges包存储基因组范围。

---

IRanges的安装方式如下：

BiocManager::install("IRanges")

创建范围对象

使用IRanges函数创建范围对象只需要提供起始位点与终止位点的坐标：

> rng <- IRanges(start = 1, end = 5)
> rng
IRanges object with 1 range and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]         1         5         5

注意：IRange的坐标系统是从1开始的（这与R体系一致），两端都是闭区间。

除了使用起止位置，还可以指定宽度来生成范围，如下：

> rng <- IRanges(start = 1, width = 5)
> rng <- IRanges(end = 1, width = 5)

可以一次生成多组范围：

> x <- IRanges(start = c(1,3,33,8), end = c(5,10, 40, 20))
> x
IRanges object with 4 ranges and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]         1         5         5
  [2]         3        10         8
  [3]        33        40         8
  [4]         8        20        13

给每个范围命名：

> names(x) <- letters[1:4]
> x
IRanges object with 4 ranges and 0 metadata columns:
        start       end     width
    <integer> <integer> <integer>
  a         1         5         5
  b         3        10         8
  c        33        40         8
  d         8        20        13

我们生成的范围对象和R的dataframe很相似，不过其是IRanges对象：

> class(x)
[1] "IRanges"
attr(,"package")
[1] "IRanges"

> attributes(x)
$start
[1]  1  3 33  8
$width
[1]  5  8  8 13
$NAMES
[1] "a" "b" "c" "d"
$elementType
[1] "ANY"
$elementMetadata
`\001NULL\001`
$metadata
list()
$class
[1] "IRanges"
attr(,"package")
[1] "IRanges"

> str(x)
Formal class 'IRanges' [package "IRanges"] with 6 slots
  ..@ start          : int [1:4] 1 3 33 8
  ..@ width          : int [1:4] 5 8 8 13
  ..@ NAMES          : chr [1:4] "a" "b" "c" "d"
  ..@ elementType    : chr "ANY"
  ..@ elementMetadata: NULL
  ..@ metadata       : list()

基本操作

有一些方便的函数快速地获取起点，终点与宽度信息：

> start(x)
[1]  1  3 33  8
> end(x)
[1]  5 10 40 20
> width(x)
[1]  5  8  8 13

可以直接使用<- 更改范围的起止点或者宽度，例如：

> start(x) <- start(x) + 1
> x
IRanges object with 4 ranges and 0 metadata columns:
        start       end     width
    <integer> <integer> <integer>
  a         2         5         4
  b         4        10         7
  c        34        40         7
  d         9        20        12

通过range函数可以得到所有起止点所跨越的范围：

> range(x)
IRanges object with 1 range and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]         2        40        39

提取子集

判断一下各个范围的宽度是否大于5：

> width(x) > 5
[1] FALSE  TRUE  TRUE  TRUE

可以通过如下的方式提取范围大于5的子集：

> x[width(x)>5]
IRanges object with 3 ranges and 0 metadata columns:
        start       end     width
    <integer> <integer> <integer>
  b         4        10         7
  c        34        40         7
  d         9        20        12

同理，根据起点或者终点位置提取子集：

> x[start(x)>5]
IRanges object with 2 ranges and 0 metadata columns:
        start       end     width
    <integer> <integer> <integer>
  c        34        40         7
  d         9        20        12
> x[end(x)>20]
IRanges object with 1 range and 0 metadata columns:
        start       end     width
    <integer> <integer> <integer>
  c        34        40         7

也可以直接提取指定名字的范围：

> x["a"]
IRanges object with 1 range and 0 metadata columns:
        start       end     width
    <integer> <integer> <integer>
  a         2         5         4

合并范围

合并多个IRange对象只需要使用c函数：

> a <- IRanges(1,2)
> b <- IRanges(8,10)
> c(a,b)
IRanges object with 2 ranges and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]         1         2         2
  [2]         8        10         3