太长不看系列
- IRanges提供最简单的注释信息,起点,终点以及宽度(终点位置-起点位置)
- Granges和IRanges类似,但是具有染色体编号和正负链信息
废话超多系列
Granges是genomic ranges的简称,它是bioconductor中的R包存储genomic intervals的数据结构,快速且有效,是想要深入分析数据所必须掌握的。它的重要性类似于bed注释文件,这么一比较这玩意的重要性就显而易见了
做基因组学分析时,最重要的就是intervals或者是intervals的集合。那么什么是intervals呢?它是基因组上的一些实体,我也想不出来一个合适的名词去翻译它,你可以粗鲁的称呼它为对象。
虽然翻译不出来,但是我举几个例子,一定能让你明白它究竟是个什么玩意
如果不能,我就罚自己今天不能吃肉~
如下图,这是UCSC browser上随便挑选出的一段序列的的截图,我们能从上面看到DNA clusters,SNPs ,promoters, genes,这些东西都可以称之为intervals,甚至一个测序read被map到基因组上的位置也可以称之为一个intervals
许多基因组工作涉及到各种intervals的各种相互关系,如哪些promoters含有SNP位点,某个TF结合位点又与哪些promoter重合,哪些基因被测序reads覆盖。说了这么多intervals,那么它和我们要介绍的GRanges有什么关系呢?
上图是一个GRanges的结构,其中的A1,A2和A3是它的三个intervals,所以我们可以理解,intervals是GRanges对象的横向最小单位。
从纵向上看,GRanges对象由三列组成,第一列是染色体名称,第二列是intervals所在的具体位置(这个位置信息是IRanges对象所提供的),第三列则是该位置所在链的方向
你肯定会疑惑什么是IRanges对象,不着急,下一部分就是IRange对象的介绍
除此之外,我们还能看到左下角是染色体的信息,比如seqlengths告诉了我们chr1的长度,若未定义则为NA
我们可以使用R来处理GRanges对象,但除此之外我们还可以使用bedtools处理intervals之间的相互关系。比如,我最常用的是intersectBed,它可以很快的帮助我知道我所需要的信号处于哪些位置。
IRanges基本用法
IRanges是包含intervals的向量, IRanges和GRanges类似,但是GRanges具有染色体编号和正负链信息,IRanges为GRanges提供了基础功能,可以进行intervals的简单代数运算
IRanges有三列,起始,终止以及宽度。当我们给定了任意两列信息就可以由此推断剩下一列的信息。
下面使用IRanges函数构建IRanges对象
library('IRanges')
ir1 <- IRanges(start = c(1,3,5), end = c(3,5,7))
ir1
## IRanges of length 3
## start end width
## [1] 1 3 3
## [2] 3 5 3
## [3] 5 7 3
这个个Iranges对象包含三个不同的intervals,每一个intervals都称之为一个range(范围区间),这里的ir1有三个ranges
start(), end(), width()这三个方法可以获得IRanges对象的起始列,终止列,以及宽度那一列的信息。IRanges对象的每个ranges也可以有自己的名字,如
names(ir1) <- paste("A", 1:3, sep = "")
ir1
## IRanges of length 3
## start end width names
## [1] 1 3 3 A1
## [2] 3 5 3 A2
## [3] 5 7 3 A3
因为是向量,因此IRanges只有一个维度,不能使用dim()方法,但是我们可以使用length()获取其长度。我们可以像操作向量一样,对IRanges对象进行操作,如
ir1[1]
## IRanges of length 1
## start end width names
## [1] 1 3 3 A1
ir1["A1"]
## IRanges of length 1
## start end width names
## [1] 1 3 3 A1
甚至可以利用c()将两个IRanges对象合并
c(ir1, ir2)
## IRanges of length 6
## start end width names
## [1] 1 3 3 A1
## [2] 3 5 3 A2
## [3] 5 7 3 A3
## [4] 1 1 1
## [5] 3 3 1
## [6] 5 5 1
有时候在阅读R包的资料时,你也许会看到NormalIRanges对象,它是由使用reduce函数对IRanges对象处理得到。一个NormalIRanges,他其中的ranges(就是每一行数据)有以下几个特点:
- 非空(即:它的宽度不能为0)
- 各个ranges没有重合
- 顺序是从左到右(即:第n+1个ranges的起始坐标一定比第n个ranges的终止坐标大)
- 两个ranges彼此不连接(即:第n+1个ranges的起始坐标和第n个ranges的终止坐标不能相等)
用文字可能比较抽象,用一张图应该能帮助大家更好的理解。上面的表示我们的其实Iranges对象,而下面的则是我们使用reduce函数处理后所得到的normal Ranges对象。
下面是reduce()函数的对象:
ir
## IRanges of length 4
## start end width
## [1] 1 4 4
## [2] 3 4 2
## [3] 7 8 2
## [4] 9 10 2
reduce(ir)
## IRanges of length 2
## start end width
## [1] 1 4 4
## [2] 7 10 4
应用场景:给定了一系列重叠的reads位置,有哪些碱基属于exxon,又有哪些碱基属于同一个外显子呢?
disjoin函数和reduce函数恰好相反,找出完全不重合的intervals,将重合的intervals进行分割。
还有一些方法可以用于批量转换IRanges对象的坐标位置,比如:
-
shift()
:帮助我们将IRanges对象位置整体迁移,而保持宽度不变 -
narrow()
:可以指定IRanges某一范围内的ranges的宽度 -
restric()
:可以帮我们提取某一个范围内的ranges
除此之外,还有一些方法对于上述三个方法做了进一步的拓展,由于篇幅原因就不在这里介绍,大家感兴趣可以自己搜搜看
reseize()函数可能是一个处理数据很有用的方法,它有一个fix参数,可以固定原始ranges中不变的列,具体用法如下:
ir
## IRanges of length 4
## start end width
## [1] 1 4 4
## [2] 3 4 2
## [3] 7 8 2
## [4] 9 10 2
resize(ir, width = 1, fix = "start")
## IRanges of length 4
## start end width
## [1] 1 1 1
## [2] 3 3 1
## [3] 7 7 1
## [4] 9 9 1
resize(ir, width = 1, fix = "center")
## IRanges of length 4
## start end width
## [1] 2 2 1
## [2] 3 3 1
## [3] 7 7 1
## [4] 9 9 1
既然IRanges是一个向量,那么肯定可以对IRanges进行集合操作,我们通常使用unnion取并集,intersect取交集:
ir1 <- IRanges(start = c(1, 3, 5), width = 1)
ir2 <- IRanges(start = c(4, 5, 6), width = 1)
union(ir1, ir2)
## IRanges of length 2
## start end width
## [1] 1 1 1
## [2] 3 6 4
intersect(ir1, ir2)
## IRanges of length 1
## start end width
## [1] 5 5 1
reduce(c(ir1, ir2))
## IRanges of length 2
## start end width
## [1] 1 1 1
## [2] 3 6 4
这些函数,还有一些变种,punion(), pintersect(), psetdiff(), pgap()。它们与R基础函数中的pmax很相似,很少被用到,因此不做讲解
其实是我自己不懂,没有用过
对测序数据进行下游分析时,经常需要寻找overlap,使用findOverlaps()函数寻找两个IRanges对象中的overlap。它的返回值是一个Hists对象,它可以描述两个Iranges之间重叠关系,是一个两列的矩阵,由IRanges的索引组成。可以由queryHist()和subjectHits()得到其索引信息。
ir1 <- IRanges(start = c(1,4,8), end = c(3,7,10))
ir2 <- IRanges(start = c(3,4), width = 3)
ov <- findOverlaps(ir1, ir2)
ov
## Hits object with 3 hits and 0 metadata columns:
## queryHits subjectHits
## <integer> <integer>
## [1] 1 1
## [2] 2 1
## [3] 2 2
## -------
## queryLength: 3
## subjectLength: 2
结果解释:query的第一个range和subject的第二个range有重合,query的第二个和subject的第一和第二两个range都有重合。
findOverlaps的参数很多,是否有最小重叠,是否完全重叠等。我自己也记不住,所以就不在这里介绍了
看到这里可以思考一下,findOverlaps()和union()有什么区别
countOverlaps()函数可以返回overlap的数目,该功能速度更快,占用内存更小,因为不需要记录ranges重叠的具体信息,仅仅考虑重叠的个数
nearest()函数用于寻找两个IRanges中,第一个IRanges对象中的ranges最接近第二个IRanges对象中的哪个ranges。
ir1
## IRanges of length 3
## start end width
## [1] 1 3 3
## [2] 4 7 4
## [3] 8 10 3
ir2
## IRanges of length 2
## start end width
## [1] 3 5 3
## [2] 4 6 3
nearest(ir1, ir2)
## [1] 1 1 2
如果你用过Y叔的R包ChIPseeker,你可能会对这个功能比较熟悉,它将peak注释到基因的一个方法就是通过找该peak与哪个基因最接近
类似的函数还有precede,follow等,实在是太多了,而且不常用,就不一一介绍了……
相信看到这里你们对于IRanges对象已经有了一个比较深的认识,对于GRanges也有了大致的了解。下一期不出意外,我会详细写写GRanges对象的操作,希望能带你进一步认识GRanges对象