搜索
查看: 3521|回复: 0

[software] Rsamtools

[复制链接]

64

主题

138

帖子

681

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
681
发表于 2017-7-8 09:41:34 | 显示全部楼层 |阅读模式
这是一个查看比对文件很好用的包。坐台女看群主的视频,然后也有用到,感觉很有用。
可以直接查看文件内容
也可以通过设置参数只查看某几列内容。
赶快试试吧
## Overview

`r Biocpkg(“Rsamtools”)`包含读取和检查BAM格式的比对读段的功能。

## Rsamtools

Rsamtools包是广泛使用的`samtools` /`htslib`库的界面。 该软件包的主要功能是支持读取BAM文件。

### The BAM / SAM file format


SAM格式是一种基于文本的比对表示形式。 BAM格式是SAM的二进制版本,更小更快。

格式相当复杂,这意味着R表示也有点复杂。 由于文件格式的以下特征,会发生此并发症

- 它可能包含未对齐的序列。
- 一个序列可以对齐到多个位置。
- 支持拼接(分割)比对。
- 它可能包含来自多个样品的读段。


BAM文件可以以多种方式进行排序。 如果根据基因组位置进行排序,如果它也被“索引”,则可以检索映射到基因组位置的所有读段,这些可能非常有用。 在`r Biocpkg(“Rsamtools”)`这是由函数`sortBam()`和`indexBam()`完成的。
## scanBam

如何阅读BAM文件在概念上像这样
1.一个指向该文件的指针由`BamFile()`构造函数创建。
2.(可选)要报告的读取的参数由“ScanBamParams()”构建。
3.该文件正在根据这些参数被`scanBam()`读取。

首先我们设置一个`BamFile`对象:

[AppleScript] 纯文本查看 复制代码
```{r bamPath}	
bamPath <- system.file("extdata", "ex1.bam", package="Rsamtools")
bamFile <- BamFile(bamPath)
bamFile
```


一些高级信息可以在这里访问,像
```{r bamFileInfo}
seqinfo(bamFile)
```

(显然也支持`seqinfo()`和`seqlevels()`等)。

现在,我们使用`scanBam()`读取文件中的所有读取,忽略使用“ScanBamParams()”选择读取的可能性(我们将返回到下面)。

```{r scanBam}
aln <- scanBam(bamFile)
length(aln)
class(aln)
```

我们得到长度为1的列表; 这是因为`scanBam()`可以从多个基因组区域返回输出,这里我们只有一个 。 我们因此将输出分组 这再次给我们一个列表,我们显示第一个比对信息
[AppleScript] 纯文本查看 复制代码
```{r lookAtBam}
aln <- aln[[1]]
names(aln)
lapply(aln, function(xx) xx[1])
```



注意`scanBam()`函数如何返回一个基本的R对象,而不是一个S4类。将比对表示为S4对象(由`r Biocpkg(“GenomicAlignments”)`);这对于从RNA测序数据获得拼接的比对是特别有用的。

“aln”列表的名称基本上是BAM规范中使用的名称。以下是比较重要的列表:

- `qname`:读段的名称。
- `rname`:染色体/序列/ contig的名称。
- `strand`:比对的链。
- `pos`:比对的最左边部分的坐标。
- `qwidth`:读段的长度。
- “mapq”:对齐的映射质量。
- “seq”:对齐的实际顺序。
- `qual`:对齐的质量字符串。
- `cigar`:CIGAR字符串(下图)。
- “flag”:flag。

###在文件的部分阅读

BAM文件可能非常大,通常需要阅读文件的部分内容。你可以用不同的方法来做到这一点

1.读取设定数量的记录(比对)。
2.只读取符合特定标准的比对。
3.仅在某些基因组区域读取比对。

让我们从第一个开始吧。通过在使用`BamFile()'时指定`yieldSize`,每次调用`scanBam然后可以再次调用`scanBam()`来获取下一组比对方式;这需要你先打开()文件(否则你会保持读取相同的对齐方式)。

[AppleScript] 纯文本查看 复制代码
```{r yieldSize}
yieldSize(bamFile) <- 1
open(bamFile)
scanBam(bamFile)[[1]]$seq
scanBam(bamFile)[[1]]$seq
## Cleanup
close(bamFile)
yieldSize(bamFile) <- NA
```


在BAM文件的一部分读取的另外两种方法是使用“ScanBamParams()”,特别是`what`和`which`参数。
[AppleScript] 纯文本查看 复制代码
```{r ScanBamParams}
gr <- GRanges(seqnames = "seq2",
              ranges = IRanges(start = c(100, 1000), end = c(1500,2000)))
params <- ScanBamParam(which = gr, what = scanBamWhat())
aln <- scanBam(bamFile, param = params)
names(aln)
head(aln[[1]]$pos)
```


注意`pos`如何小于`which`参数中指定的值?这是因为比对与“which”参数重叠。 `what = scanBamWhat()`告诉函数读取所有内容。通常,您可能不会对阅读的实际序列或其质量得分感兴趣。这些占用了大量空间,因此您可能会考虑禁用阅读此信息。

#### CIGAR字符串


“CIGAR”是BAM格式如何表示拼接比对。例如,文件存储比对的最左边的坐标。要获得右边的坐标,您必须解析CIGAR字符串。在这个例子中,“36M”表示它没有插入或删除。如果您需要使用包含插入或删除的拼接比对或比对方式,则应使用“r Biocpkg(”GenomicAlignments“)包。
#### Flag

比对可能有一些“flag”设置或未设置。这些标志提供有关比对的信息。flag整数同时表示多个标志。flag的示例表示(对于配对的端比对)两个比对是否已经被正确比对。有关更多信息,请参阅BAM规范。

在`r Biocpkg(“Rsamtools”)`中有许多帮助函数仅处理读取某些flag;使用这些。

### BAM总结

有时需要快速总结BAM文件中的对齐方式:

[AppleScript] 纯文本查看 复制代码
```{r summary}
quickBamFlagSummary(bamFile)
```


## Other functionality from Rsamtools

### BamViews

不是阅读单个文件的时候,可以构造一个名为“BamViews”的东西,这是一个指向多个文件的链接。 在许多方面,它具有与其他视图相同的“view”功能。下例是一个很好的说明

[AppleScript] 纯文本查看 复制代码
```{r BamViews}
bamView <- BamViews(bamPath)
aln <- scanBam(bamView)
names(aln)
```


这给我们一个extra的列表级别的返回对象; 第一级是文件,二级是范围。

我们也可以在`BamViews`上设置`bamRanges()`来指定只读一些范围; 这与“ScanBamParams()”的`which`参数相似。
[AppleScript] 纯文本查看 复制代码
```{r BamViews2}
bamRanges(bamView) <- gr
aln <- scanBam(bamView)
names(aln)
names(aln[[1]])
```


### countBam

有时候,你想做的只是计数,使用`countBam()` 。

回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 08:28 , Processed in 0.029264 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.