IRanges包可以用来确定两个range（范围）的重叠，此操作可以用于比对read，寻找变异位点，以及定量RNA-seq中基因的表达等等。
今天我们就了解一下如何寻找两个range的重叠。
基本操作
首先，我们建立两个range：qry(query range)与sbj(subject range)：
> qry <- IRanges(start=c(1, 26, 19, 11, 21, 7), end=c(16, 30, 19, 15, 24, 8), names=letters[1:6])
> sbj <IRanges(start=c(1, 19, 10), end=c(5, 29, 16), names=letters[24:26])
通过findOverlaps函数寻找两个range的重叠(hts)：
> hts <- findOverlaps(qry, sbj)
> hts
Hits object with 6 hits and 0 metadata columns:
queryHits subjectHits
<integer> <integer>
[1] 1 1
[2] 1 3
[3] 2 2
[4] 3 2
[5] 4 3
[6] 5 2
-------
queryLength: 6 / subjectLength: 3
hts的两列数字对应qry与sbj对象的坐标(坐标分别通过queryHits与subjectHits函数提取)：
> names(qry)[queryHits(hts)]
[1] "a" "a" "b" "c" "d" "e"
> names(sbj)[subjectHits(hts)]
[1] "x" "z" "y" "y" "z" "y"
这个结果实际上是query range与subject range之间的映射关系（图1，图2）：
图1
图2
query range选项
咋一看range重叠操作非常简单直接，细思的话这个问题其实并不是那么简单。我们现在得到的结果是基于认为qry和sbj存在任意的重合即认为发生了重叠(默认type = "any")，重叠的判定方式还可以是认为qry只有完全被sbj容纳才认为是重叠，可以修改参数type = "within"来匹配这种判定方式的结果（如图3）：
> hts_within <- findOverlaps(qry, sbj, type = "within")
> hts_within
Hits object with 3 hits and 0 metadata columns:
queryHits subjectHits
<integer> <integer>
[1] 3 2
[2] 4 3
[3] 5 2
-------
queryLength: 6 / subjectLength: 3
图3
除了这两个最常用的，type还可以取其它值（使用help(findOverlaps)了解）。
subject range选项
通过select参数修改subject所选取的结果，默认是取与当前query range重叠的所有subject range（即select = “all”）。我们可以取与qry重叠的第一个，最后一个或者任意一个range结果(参考图2)：
> findOverlaps(qry, sbj, select = "first")
[1] 1 2 2 3 2 NA
> findOverlaps(qry, sbj, select = "last")
[1] 3 2 2 3 2 NA
> findOverlaps(qry, sbj, select = "arbitrary")
[1] 1 2 2 3 2 NA
寻找重叠的操作非常耗时，假设query range有Q个对象而subject range有S个对象，每个query range与所有subject range进行比对的话需要计算QxS次。但是，明显有更好的做法：我们可以将subject range进行排序，构建一个NCList对象，只要比对query range的起点与终点之间的subject range就可以了。IRanges包提供了这样的做法，实现起来也很简单：
> sbj_ncl <- NCList(sbj)
> findOverlaps(qry, sbj_ncl)
Hits object with 6 hits and 0 metadata columns:
queryHits subjectHits
<integer> <integer>
[1] 1 1
[2] 1 3
[3] 2 2
[4] 3 2
[5] 4 3
[6] 5 2
-------
queryLength: 6 / subjectLength: 3
虽然这里没有体现性能上的优势，但是当我们有一个固定的subject range要与非常多的query range进行重叠匹配的话将节省大量时间。
hits对象
我们可以将hits转化为matrix:
> as.matrix(hts)
queryHits subjectHits
[1,] 1 1
[2,] 1 3
[3,] 2 2
[4,] 3 2
[5,] 4 3
[6,] 5 2
通过countQueryHits提取query range发生重叠的次数：
> countQueryHits(hts)
[1] 2 1 1 1 1 0
> setNames(countQueryHits(hts), names(qry))
a b c d e f
2 1 1 1 1 0
类似，使用countSubjectHits提取subject range发生重叠的次数：
> countSubjectHits(hts)
[1] 1 3 2
> setNames(countSubjectHits(hts), names(sbj))
x y z
1 3 2
使用overlapsRanges函数查看qry与sbj间发生重叠的具体区域：
> overlapsRanges(qry, sbj, hts)
IRanges object with 6 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 1 5 5
[2] 10 16 7
[3] 26 29 4
[4] 19 19 1
[5] 11 15 5
[6] 21 24 4
注意：这个函数与intersect函数（见上节）存在区别，因为intersect函数会整合范围（即qry的a与d合并为一个range），而overlapsRanges返回了所有qry range的结果：
> intersect(qry, sbj)
IRanges object with 5 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 1 5 5
[2] 10 16 7
[3] 19 19 1
[4] 21 24 4
[5] 26 29 4
通过R我们可以总览一下重叠区域长度的分布（甚至可以使用ggplot进行可视化）：
> summary(width(overlapsRanges(qry, sbj, hts)))
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 4.000 4.500 4.333 5.000 7.000
最后补充两个方便的函数：
- 使用
countOverlaps可以直接查看qry与sbj重叠的数目（等同于setNames(countQueryHits(hts)), names(qry)）：
> countOverlaps(qry, sbj)
a b c d e f
2 1 1 1 1 0
- 使用
subsetByOverlaps可以只保留在sbj中出现的qry range（等同于qry[unique(queryHits(hts))]）:
> subsetByOverlaps(qry, sbj)
IRanges object with 5 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
a 1 16 16
b 26 30 5
c 19 19 1
d 11 15 5
e 21 24 4
网友评论