【Bioconductor系列】Bioconductor的地基--IRanges

Bioconductor的地基--IRanges

Bioconductor是一个开源项目,包括许多R生物信息学包。这里,首先介绍Bioconductor的核心包:

  • GenomicRanges: 用于表示和使用基因组范围,genomic ranges
  • GenomicFeatures: 用于表示和使用基因模型和基因组其他特性(genes, exons,UTRs,transcripts等)的范围
  • Biostrings和BSgenomes:用于在R中操作基因组序列(比如有些包可以从ranges中提取序列)
  • rtrcklayer: 用于读取常见的生物信息数据格式,BED,GTF/GFF,和WIG

安装和帮助说明

source("http://biocondutor.org/biocLite.R")
biocLite("GenomicRanges")
验证或升级版本:biocValid()、biocUpdatePackages()

使用IRanges储存通用数据

IRanges构建语法:IRanges(start=NULL, end=NULL, width=NULL, names=NULL)

library(IRanges)
rng <- IRanges(start=4, end=13)
x <- IRanges(start=c(4,7,2,20), end=c(13,7,5,23))
names(x) <- letters[1:4] #letter内置常量,小写字母,LETTER,大写字母,
class(x) # 查看类
#一些常用方法: start(), end(), width(),range()
# 算术、切片、逻辑操作
start(x) + 4
x[2:4]
x['a']
x[start(x) <4]
# merge
a <- IRanges(start=7, width=4)
b <- IRanges(start=2, end=5)
c(a,b)

基本Range操作,算术、变形、集合

IRanges对象可以通过算术运算增大或缩小范围

x <- IRanges(start=c(40,80), end=c(67,114)
x + 4L
x - 10 L
# 这些运算从两端同时进行
# 截取部分区域,restrict()
y <- IRanges(start=c(4,6,10,12), width=13)
restrict(y,5,10)
# flank()可以提取每个range的两端部分,
flank(x,width=7)
# reduce() 相当于read组装
set.seed(0)
alns <- IRanges(start=sample(seq_len(50),20),width=5) #sample随机取样,seq_len产生范围数据
head(alns,4)
reduce(alns)
# gaps找出不同序列间隔,默认不关注the beginning of the sequences to the start position of the first ranges, same as end
gaps(alns)
# 集合操作:交集intersect, 差级setdiff, 合集union, 逐对操作,pintersect,psetdiff,punion,pgap

寻找overlapping ranges

寻找overlap是许多基因组学分析任务的必要环节。RNA-seq,overlaps用于分析细胞活动
我们使用findOverlaps()函数寻找两个IRanges对象的overlaps.

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])
hts <- findOverlaps(qry, sbj)
hts
Show in New WindowClear OutputExpand/Collapse Output
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

Overlaps表示的是query和subject之间的映射关系。每个query根据我们寻找overlaps的方式不同,在不同的subjects可以有多个hits.单个subjects也可以有多个query的hits.举findOverlaps给出的结果[1] 1 1为例,表示为query的1和subject的1存在映射。我们可以使用accessor函数--queryHits(),和subjectHits()查看索引。

names(qry)[queryHits(hts)]  # 分解成index <- queryHts(hts); name <- names(qry); name[index]
names(suj)[subjectHits(hts)]

默认是当range任意部分与subject range有重叠,就被认为是overlap,这个默认type为"any"。也可以修改为type=within

hts_within <- findOverlaps(qry, sbj, type="within")

within和any是最常用的两种类型,其他的可以通过help(findOverlaps)了解。还有一个参数为select,用来处理多对多的映射关系,默认是select="all",first表示第一个,last表示最后一个,arbitrary表示是任意。
时刻留意Overlaps,因为有可能会造成结果的巨大差异
计算overlaps需要强大的计算能力,需要对query和subject进行一一对比。比较明智的方法是使用interval tree进行恰当的排序,使用help(IntervalTree)查看帮助,似乎这个功能好像取消了。
运行findOverlaps,我们需要提取Hits对象中的信息,例如之前用到的subjectHits。

as.matrix(hts) # 直接强迫转换成matrix
countQueryHits(hts) # 计算每个query的命中数
setNames(countQueryHits(hts),names(qry)) # 给对象命名
countSubjectHits(hts)
setNames(countSubjectHits(hts),names(sbj)
ranges(hts,qry,sbj)  # 找出qry,sbj的共同之处

ranges所创建的对象可以使用R的所有向量数据分析方法。例如summary.
subsetByOverlaps()仅保留的query的子集,countOverlaps()则计算重叠数。

寻找邻近Ranges以及计算距离

有三种方法可以完成邻近ranges的寻找:neatest(), precede(), fellow()

qry <- IRanges(start = 6, end = 13, names = 'query')
sbj <- IRanges(start = c(2,4,18,19), end=c(4,5,21,24),names=1:4)
nearest(qry, sbj) #最近的距离
precede(qry,sbj) #寻找前面的ranges
follow(qry,sbj) #寻找后面的ranges

这些操作允许vertorization.·qry2 <- IRanges(start=c(6,7), width=3); nearest(qry2,sbj)`.
使用distanceToNearest()和distance()确定两个ranges之间的距离

qry <- IRanges(sample(seq_len(1000),5),width=10)
sbj <- IRanges(sample(seq_len(1000),5),width=10)
distanceToNearest(qry,sbj)
distance(qry,sbj) #仅返回距离

Run Lenth Enconding and Views

IRanges可以操作任意类型的序列的ranges。如果是基因组数据,那么ranges的坐标依赖于核酸和特定的染色体。当然爱有许多其他类型的基因组数据,例如:

  • Coverage, 一定长度序列中ranges的重叠深度
  • Conservation tracks, 两个物种间,碱基与碱基之间的进化保守计分
  • Per-base pair estimates of population genomics summary statistics like nucleotide diversity.

Run-length encoding and coverage()

序列越长所占的容量越大,为了能在R中存放如此大的数据,IRanges使用一个小技巧来压缩这些序列。例如444333222就可以表示为3个4,3个3,3个2。这个压缩方法称为run-length encoding, RLE 。并且压缩后的数据支持原来的所有常规的R向量操作,如截取、算数、summary,和数学函数

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))
xrle <- Rle(x) # 压缩为run coding length

integer-Rle of length 28 with 7 runs
  Lengths: 3 2 1 5 7 3 7
  Values : 4 3 2 1 0 1 4

as.vector(xrle) #解压缩
runLength(xrle) # 仅获取lengths
runValue(xrle) # 仅获取values

我们可以在converage()中遇到run-length encoded values. converage()输入一组ranges, 返回以Rle对象表示的覆盖度.

set.seed(0)
rngs <- IRanges(start =sample(seq_len(60), 10), width=7)
names(rngs)[9] <- 'A'
rngs_cov <- coverage(rngs)

rngs_cov > 3 # 找到覆盖度大于3的区域
rngs_cov[as.vector(rngs_cov) >3] # 提取覆盖度大于3的值
rngs_cov[rngs['A']] #标记为A的区域的覆盖度
mean(rngs_cov[rngs['A'])

在基因组学中有大量的分析工作,例如对一组ranges(重复区域,蛋白编码区域,低重组区等)计算一些序列的描述性统计量(coverage, GC含量,核酸多样性等)。这些操作在我们使用IRanges中的methods会选择特别琐碎,不过可以用GenomicRanges中相似方法来处理。

使用slice()从RLE序列创建ranges

define new ranges: taking coverage and defining ranges correspnding to extremely high-coverage peaksm or a map of per-base pair recombination estimates and defining a recombinationa hotspot region.
例如找出覆盖率低于2的ranges.

min_cov2 <- slice(rngs_cov, lower=2)
min_cov2
Views on a 63-length Rle subject

views:
    start end width
[1]    16  18     3 [2 2 2]
[2]    22  22     1 [2]
[3]    35  39     5 [2 2 2 2 2]
[4]    51  62    12 [2 2 2 3 3 3 4 3 3 3 2 2]

rangs(min_cov2) # 提取views中的range

IRanges进阶:Views

Views通过整合sequence vector和ranges,使得对根据特定ranges的sequence vector的聚合操作更加简单。这与之前所提到分组描述性统计分析类似,只不过这里的组别为ranges。
例如,我们可以对之前使用slice()创建的views进行描述性统计分析,使用的函数为viewMeans(), viewMaxs(),和viewApply(可以使用任意函数)

viewMeans(min_cov2)
viewMaxs(min_cov2)
viewApply(min_cov2, median)

还有其他类似的函数,见help(viewMeans)。
同样还可以根据window/bin对序列进行描述性分析。首先是在序列中创建一组windows ranges,然后在分析.

length(rngs_cov)
bwidth <- 5L
end <- bwidth*floor(length(rngs_cov / bwidth))
windows <- IRanges(start= seq(1,end, bwidth), width = bwidth)
head(windows)
cov_by_wnd <- View(rngs_cov, windows)
head(cov_by_wnd)
viewMeans(cov_by_wnd)

命令行有一个工具叫做bedtools,能够实现如上代码所做的事情。
顺便放一下自己的知识星球,如果你觉得我对你有帮助的话。


知识星球