R专题：使用R读取，处理BED格式的文件

什么是BED?

BED的全称是Browser Extensible Data，是二代测序数据常见的数据格式，用来表述基因组区段位置及其附属信息，如果大家忘了，可以回顾 NGS数据格式之bed，需要注意两点：

1.bed文件一条染色体的起始位点是从0开始的，也就是txt文件转bed的时候起始位点要减1;

2.不管你把什么文件转成bed，都请严格按照bed格式转换，尤其是第六列代表链的正负，忘记的请点上面的链接；

使用R导入BED file

好了，理论看完了接下来就是实战了。既然BED文件很常见，那通往罗马的路肯定就不止一条，下面介绍几个可以读取bed文件的包：

1.BayesPeak

这是我们原始的测试文件test.bed：

下载BayesPeak并读取bed文件，读取的test.bed文件如下：

source('https://bioconductor.org/biocLite.R')
biocLite('BayesPeak')
library(BayesPeak)
bedR <>'/home/qians/test.bed')

可以看到BayesPeak读取bed文件只包含了bed文件的chr, start, end 和strand四列，并且自动把对象转换成了GRanges格式，这个格式也是一些R包的输入格式，省去了我们用GRanges函数转换其他格式为GRanges的步骤，但是有时候也会遇到问题，如果我们不想对象转换成GRanges格式，可以用data.frame转换为数据框。

2.rtracklayer

source('https://bioconductor.org/biocLite.R')
biocLite('rtracklayer')
library(rtracklayer)
rtracklayer <->'/home/qians/test.bed')
rtracklayer


下载rtracklayer包并读取test.bed，可以看到rtracklayer带的函数import读取也会自动转换成GRanges，但结果比BayesPeak更完整，除了chr, start, end 和strand，还包含了其他的信息，比如name, score。

如果大家足够细心，会发现rtracklayer和BayesPeak的GRanges对象strand属性是不一样的，BayesPeak中strand是character，rtracklayer中strand是Rle，也就是说：rtracklayer自动转换后的GRanges格式更标准。

3.TEQC

source('https://bioconductor.org/biocLite.R')
biocLite('TEQC')
library(TEQC)
teqc <->'/home/qians/test.bed',filetype = 'bed')
teqc


下载安装TEQC，读取bed文件，TEQC读取的结果里信息最少的，只包含了chr, start和end三列

小结一下：这里只列举了三个，个人觉得能处理bed文件的包都是带了读入bed文件函数的，R里面最简单的read.csv, read.table等也可以把bed文件导入R，只要设置好分隔符就OK。

按染色体号提取BED序列

我们就以BayesPeak为例讲解，仔细阅读一下我们发现BayesPeak的函数read.bed其实很简单，就两个参数，第一个是文件路径，第二个指定染色体号：

这里我们只选择2L号染色体，结果如下：

指定了染色体号，就会只读取相应染色体的条目，这一点就非常方便。

可视化BED文件中的reads分布

我们先看一下测试文件：

读取测试文件，选择一段区域plot：

dir <->'extdata', package='BayesPeak')
treatment <->path(dir, 'H3K4me3reduced.bed')
bed <->read.bed(treatment)
plot.bed(bed, 'chr16', 9.1E7, 9.4E7)
##选择区域chr16:91000000-94000000

可以看到bed文件中在chr16的91000000-94000000bp范围内reads的频率：

参考资料：

https://www.rdocumentation.org/packages/BayesPeak/versions/1.24.0