45. 《Bioinformatics Data Skills》

今天完成IRanges包最后一部分的学习。

---

slice函数

我们之前简单模拟了基因组测序并计算了各个位点的覆盖度：

> library(IRanges)
> set.seed(1234)
> rngs <- IRanges(start = sample(seq_len(60), 20), width = 7)
> rngs_cov <- coverage(rngs)

通过rngs_cov > 3只能看到每个run的深度是否大于3的逻辑值：

> rngs_cov > 3
logical-Rle of length 65 with 7 runs
  Lengths:     8     2    28     2     1     3    21
  Values : FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE

如果我们关心哪些区域的深度大于3，可以使用slice函数：

> min_cov3 <- slice(rngs_cov, lower = 3)
> min_cov3
Views on a 65-length Rle subject

views:
    start end width
[1]     6  11     6 [3 3 3 4 4 3]
[2]    15  20     6 [3 3 3 3 3 3]
[3]    28  28     1 [3]
[4]    38  45     8 [3 4 4 3 4 4 4 3]
[5]    47  48     2 [3 3]
[6]    59  63     5 [3 3 3 3 3]

slice的参数为Rle对象(rngs_cov)与深度的阈值(>3X)，结果为views对象（范围与Rle向量的结合），表明深度大于3的范围及各个位点深度。

可以提取views对象底层的范围：

> ranges(min_cov3)
IRanges object with 6 ranges and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]         6        11         6
  [2]        15        20         6
  [3]        28        28         1
  [4]        38        45         8
  [5]        47        48         2
  [6]        59        63         5

slice函数非常重要，可以帮助我们确定高深度的区域，进而关注这些区域是否拥有重要的生物学特征。

Views

我们已经知道slice函数会返回views对象，这里进一步地了解一下views。

views拥有一系列的总结函数，包括viewMeans, ViewMaxs, 甚至可以自定义的viewApply函数（通过help(viewMeans)查看所有相关函数），通过这些函数对每个区域信息进行总结，例如：

> viewMeans(min_cov3) # 查看各个区域深度平均值
[1] 3.333333 3.000000 3.000000 3.625000 3.000000 3.000000
> viewMaxs(min_cov3) # 查看各个区域深度最大值
[1] 4 3 3 4 3 3
> viewApply(min_cov3, median) # 查看各个区域深度中位值
[1] 3 3 3 4 3 3

使用Views还可以将基因组划分为小的窗口片段，这是因组分析常用的手段。例如我们将rngs_cov涉及的基因组划分成多个窗口，每个窗口包括5个位点：

> length(rngs_cov)
[1] 65
> bwidth <- 5L
> end <- floor(length(rngs_cov)/bwidth) * bwidth
> windows <- IRanges(start = seq(1, end, bwidth), width = bwidth)
> head(windows)
IRanges object with 6 ranges and 0 metadata columns:
          start       end     width
      <integer> <integer> <integer>
  [1]         1         5         5
  [2]         6        10         5
  [3]        11        15         5
  [4]        16        20         5
  [5]        21        25         5
  [6]        26        30         5

这里涉及到一些数学计算，rngs_cov包括65个位点，刚好划分为13个区域。通常会遇到最后一个区域元素不足5的情况，为了避免异常结果会忽略最后一个窗口的的元素（这是我们使用floor函数的原因）。

使用Views函数将序列深度信息映射到各个窗口上：

> cov_by_wnd <- Views(rngs_cov, windows)
> head(cov_by_wnd)
Views on a 65-length Rle subject

views:
    start end width
[1]     1   5     5 [0 0 0 1 2]
[2]     6  10     5 [3 3 3 4 4]
[3]    11  15     5 [3 2 1 2 3]
[4]    16  20     5 [3 3 3 3 3]
[5]    21  25     5 [2 2 1 1 1]
[6]    26  30     5 [2 2 3 2 2]

计算各个窗口平均的read深度（可以根据需要使用任意的总结函数）：

> viewMeans(cov_by_wnd)
 [1] 0.6 3.4 2.2 3.0 1.4 2.2 1.6 2.8 3.6 2.4 0.6 1.8 2.4

---

至此，我们学完了IRange包的基本思想与用法，后面会继续学习GenomicRanges包，它的函数与用法与IRange包基本一致，但是会进一步考虑多条染色体与链信息。