44.《Bioinformatics Data Skills》之

基因组由无数的碱基位点组成，每个位点都可以生成不同类型的注释信息包括测序深度，GC含量等等。这里了解一下IRanges包对基因组深度的计算与存储方式。

---

rung-length(rle）对象

长序列非常吃内存，基因组上每个位点都使用数字来表示深度的话，光1号染色体所有位点的深度信息就需要高达1.9G字节的内存。IRanges采用了一种更聪明的叫做run-length编码表示方式：将每组重复值压缩为1 run。原理很简单，例如下面这串数字：

4 4 4 3 3 2 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 4 4 4 4 4 4 4  

可以表示为：3个4，2个3，1个2，5个1，7个0，3个1，7个4。

在R里面采用Rle函数来生成run-length(Rle)对象：

> x <- as.integer(c(4, 4, 4, 3, 3, 2, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 4, 4, 4, 4, 4, 4, 4))
> xrln <- Rle(x)
> xrln
integer-Rle of length 28 with 7 runs
  Lengths: 3 2 1 5 7 3 7
  Values : 4 3 2 1 0 1 4

Rle对象也可以转回为正常的向量：

> as.vector(xrln)
 [1] 4 4 4 3 3 2 1 1 1 1 1 0 0 0 0 0 0 0 1 1 1 4 4 4 4 4 4 4

Rle对象支持大多基本的向量操作，例如四则运算，逻辑运算，取子集，总结与函数运算：

> xrln + 4L
integer-Rle of length 28 with 7 runs
  Lengths: 3 2 1 5 7 3 7
  Values : 8 7 6 5 4 5 8
> xrln / 2
numeric-Rle of length 28 with 7 runs
  Lengths:   3   2   1   5   7   3   7
  Values :   2 1.5   1 0.5   0 0.5   2
> xrln > 3
logical-Rle of length 28 with 3 runs
  Lengths:     3    18     7
  Values :  TRUE FALSE  TRUE
> xrln[xrln > 3]
integer-Rle of length 10 with 1 run
  Lengths: 10
  Values :  4
> sum(xrln)
[1] 56
> summary(xrln)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.
   0.00    0.75    1.00    2.00    4.00    4.00
> round(cos(xrln), 2)
numeric-Rle of length 28 with 7 runs
  Lengths:     3     2     1     5     7     3     7
  Values : -0.65 -0.99 -0.42  0.54     1  0.54 -0.65

使用runLength与runValue函数可以分别提取Rlen的length与value值：

> runLength(xrln)
[1] 3 2 1 5 7 3 7
> runValue(xrln)
[1] 4 3 2 1 0 1 4

深度

Rle对象只有在表示大规模的数据才能够带来可观的内存收益，常用于表示基因组测序深度。通过模拟测序数据来了解一下深度概念，假设在60多个位点上进行20条read测序（第4条命名为“A”）：

> set.seed(0)
> rngs <- IRanges(start = sample(seq_len(60), 20), width = 7)
> names(rngs)[4] <- "A"

使用coverage函数计算得所有位点的深度（示意图见图1）：

> rngs_cov <- coverage(rngs)
> rngs_cov
integer-Rle of length 66 with 31 runs
  Lengths: 3 3 1 1 1 1 4 1 1 1 4 1 2 3 2 3 1 5 1 1 1 1 5 2 3 2 1 1 3 4 3
  Values : 1 2 3 2 3 4 3 4 3 2 3 2 3 2 1 0 1 2 3 2 1 2 3 4 3 2 1 0 1 2 1

图1

返回值为Rle对象，观察一下深度大于3区域：

> rngs_cov > 3
logical-Rle of length 66 with 7 runs
  Lengths:     9     1     4     1    32     2    17
  Values : FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE

整合深度大于3的区域的数目：

> rngs_cov[as.vector(rngs_cov) > 3]
integer-Rle of length 4 with 1 run
  Lengths: 4
  Values : 4

同样可以观察某一个区域的深度信息：

> rngs_cov[rngs["A"]]
integer-Rle of length 7 with 5 runs
  Lengths: 1 1 1 1 3
  Values : 3 2 1 2 3

我们通常更关心一个区域的平均深度：

> mean(rngs_cov[rngs["A"]])
[1] 2.428571

---

非常多的分析会涉及到基因组序列的分析，抽象为序列与范围的概念后我们可以通过IRanges提供的方法方便地解决这些问题。