rtracklayer - Data Import

学习最快乐

总结：
与多个基因组浏览器（目前是UCSC内置）进行交互，并以各种格式（目前为GFF，BED，BedGraph，BED15，WIG，BigWig和2bit内置）操作注释轨迹的可扩展框架。 用户可以向/从支持的浏览器导出/导入track，以及查询和修改浏览器状态，如当前视口。

[AppleScript] 纯文本查看 复制代码

 
## Dependencies

This document has the following dependencies:

```{r dependencies, warning=FALSE, message=FALSE}
library(rtracklayer)
library(AnnotationHub)
library(Rsamtools)
```

Use the following commands to install these packages in R.

```{r biocLite, eval=FALSE}
source("http://www.bioconductor.org/biocLite.R")
biocLite(c("rtracklayer", "AnnotationHub", "Rsamtools"))
```

## Corrections

Improvements and corrections to this document can be submitted on its [GitHub](https://github.com/kasperdanielh ... acklayer_Import.Rmd) in its [repository](https://github.com/kasperdanielhansen/genbioconductor).


## Overview

`r Biocpkg（“rtracklayer”）`包接口（UCSC）基因组浏览器。 它包含用于将数据导入和导出到此浏览器的功能。

这包括用于解析与UCSC基因组浏览器相关联的文件格式的功能，例如BED，Wig，BigBed和BigWig。

## Other Resources

- The vignette from the [rtracklayer webpage](http://bioconductor.org/packages/rtracklayer).

## The import function


解析数据格式的功能是`import（）`。 这个函数有一个`format`参数，取值为`BED`或`BigWig`。

请注意，有一个通用的“import（）”函数的帮助页面，但也有特定于文件格式的帮助页面。 找到这些帮助页面的最简单的方法是使用“XX”作为格式来查找“XXFile”。

[AppleScript] 纯文本查看 复制代码

```{r help, eval=FALSE}
?import
?BigWigFile
```

There are often format specific arguments.

## BED files


大多数BED文件很小，可以作为单个对象读取。 `import（format =“BED”）`的输出是`GRanges`格式。

你可以指定`genome`(例如`hg19`)，该函数将填充“GRanges”的“seqinfo”。

您还可以使用`which`参数来选择性地分析与GRanges重叠的文件的一个子集。 如果文件已经被tabix索引（见下文），这将变得更加有效。

## BigWig files



BigWig文件通常存储全基因组覆盖向量（或至少全基因组数据）。 因此，BigWig文件的R表示通常相当大，因此可能需要将其读入R中的小块。

对于BED文件，`import（format =“BigWig”）`支持`which` 它输出的数据类型是“GRanges”，默认情况下，使用`as` agurment可以使用`as =“Rle”`等几个选项。

`import（format =“BigWig”）`不支持 `genome` '参数。

##其他的文件形式

- GFF#序列注释的通用格式
- TwoBit就是Bed吧
- Wig #仅仅是为了追踪参考基因组的各个区域的覆盖度，测序深度的文件，UCSC提供
- bedGRaph #
https://www.plob.org/article/9505.html格式的网站
## Extensive example

Let us start an `AnnotationHub`:

[AppleScript] 纯文本查看 复制代码

```{r ahub}
library(AnnotationHub)
ahub <- AnnotationHub()
table(ahub$rdataclass)
```

以上可以看到多个文件格式，`GRanges`格式通常由BED文件直接建立，

```{r granges}

[AppleScript] 纯文本查看 复制代码

ahub.gr <- subset(ahub, rdataclass == "GRanges" & species == "Homo sapiens")
gr <- ahub.gr[[1]]
gr
seqinfo(gr)
```

BigWig格式的文件

[AppleScript] 纯文本查看 复制代码

```{r BigWig}
ahub.bw <- subset(ahub, rdataclass == "BigWigFile" & species == "Homo sapiens")
ahub.bw
bw <- ahub.bw[[1]]
bw
```

用 `import`来显示

[AppleScript] 纯文本查看 复制代码

```{r importBigWig}
gr1 <- gr[1:3]
out.gr <- import(bw, which = gr1)
out.gr
```

最后得到的格式 `GRanges`. 通常，一个`Rle` 可能是更恰当。

[AppleScript] 纯文本查看 复制代码

```{r importBigWig2}
out.rle <- import(bw, which = gr1, as = "Rle")
out.rle
```

得到chr22

[AppleScript] 纯文本查看 复制代码

```{r importBigWig3}
gr.chr22 <- GRanges(seqnames = "chr22",
                    ranges = IRanges(start = 1, end = seqlengths(gr)["chr22"]))
out.chr22 <- import(bw, which = gr.chr22, as = "Rle")
out.chr22[["chr22"]]
```

## LiftOver

We can re-use our `AnnotationHub`:

LiftOver是UCSC Genome浏览器中用于在不同基因组版本之间进行转换的流行工具。 `r Biocpkg（“rtracklayer”）`包也通过`liftOver'暴露了这个功能。 要使用“liftOver”，您需要一个所谓的“链”文件，描述如何从一个基因组转换到另一个基因组。 这可以从UCSC手中获取，或直接从`r Biocpkg（“AnnotationHub”）`获得。

我们可以重新使用我们的“AnnotationHub”：

[AppleScript] 纯文本查看 复制代码

```{r liftOver}
ahub.chain <- subset(ahub, rdataclass == "ChainFile" & species == "Homo sapiens")
query(ahub.chain, c("hg18", "hg19"))
chain <- ahub.chain[ahub.chain$title == "hg19ToHg18.over.chain.gz"]
chain <- chain[[1]]
gr.hg18 <- liftOver(gr, chain)
gr.hg18
```

这会将一个GRanges转换成一个GRangesList，为什么？ 这是因为单个范围（间隔）可以在另一个基因组中分裂成多个间隔。 因此，输出中的每个元素都对应于输入中的单个范围。 如果范围较小，则大多数范围应映射到单个范围。 我们来看看输出中的元素数量：

[AppleScript] 纯文本查看 复制代码

```{r liftOver2}
table(elementLengths(gr.hg18))
```

只有几个范围没有映射，只有少数分裂。


##直接从UCSC  import

使用`r Biocpkg（“rtracklayer”）`你可以直接从UCSC Genome浏览器导入表和tracks。然而，现在可以从`r Biocpkg（“AnnotationHub”）得到很多（全部）这个数据.

可能并非所有的轨迹/表格和/或来自UCSC的不同轨迹/表格的所有信息都暴露在`r Biocpkg（“AnnotationHub”）`中。


## Tabix索引

Tabix索引是使用染色体位置对文本文件进行索引的方法，用于随机访问。这将大大加快对这样一个文件的查询。 [tabix]（http://www.htslib.org/doc/tabix.html）功能在SAMtools库中引入;这个图书馆后来改名为htslib。

Tabix索引通常是您在命令行中执行的操作，但也可以使用`r Biocpkg（“Rsamtools”）包中的`indexTabix'从Bioconductor内部进行此操作。首先，文件需要被压缩bgzip2，你可以使用`bgzip2`函数。一个完整的管道，使用`r Biocpkg（“Rsamtools”）`的示例SAM文件

[AppleScript] 纯文本查看 复制代码

```{r tabixIndex}
library(Rsamtools)
from <- system.file("extdata", "ex1.sam", package="Rsamtools",
                    mustWork=TRUE)
from
to <- tempfile()
zipped <- bgzip(from, to)
idx <- indexTabix(zipped, "sam")
idx
```