bioconductor入门第二弹--GRanges概述

太长不看系列

• 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

image.png

许多基因组工作涉及到各种intervals的各种相互关系，如哪些promoters含有SNP位点，某个TF结合位点又与哪些promoter重合，哪些基因被测序reads覆盖。说了这么多intervals，那么它和我们要介绍的GRanges有什么关系呢？

image.png

上图是一个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（就是每一行数据）有以下几个特点：

1. 非空（即：它的宽度不能为0）
2. 各个ranges没有重合
3. 顺序是从左到右（即：第n+1个ranges的起始坐标一定比第n个ranges的终止坐标大）
4. 两个ranges彼此不连接（即：第n+1个ranges的起始坐标和第n个ranges的终止坐标不能相等）

用文字可能比较抽象，用一张图应该能帮助大家更好的理解。上面的表示我们的其实Iranges对象，而下面的则是我们使用reduce函数处理后所得到的normal Ranges对象。

image.png
image.png

下面是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进行分割。

image.png
image.png

还有一些方法可以用于批量转换IRanges对象的坐标位置，比如：

1. shift()：帮助我们将IRanges对象位置整体迁移，而保持宽度不变
2. narrow()：可以指定IRanges某一范围内的ranges的宽度
3. 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对象

image.png