从bam文件中提取读取位置

Question

我有一个包含多个 SNP 的 vcf 文件，现在我想看看这些 SNP 是否均匀分布在我从中获取 SNP 的 bam 文件的读取中。具体来说，我想绘制读取位置上的 SNP 数量。 我想知道是否有一些工具可以执行此操作，或者我是否必须自己编写脚本。如果是这样，R 中是否有一个包可以让我做到这一点（我习惯了 R，但对 perl 没有太多经验）？

Answer 1

不确定“读取位置上的 SNP”是什么意思，但您可以使用 R / Bioconductor 包和函数 VariantAnnotation::readVcf 读取 VCF，并使用基因组坐标通过 Rsamtools::countBam 查询 bam 文件，使用

ScanBamParam

。未经测试，沿着

的思路

## first-time installation
source("http://bioconductor.org/biocLite.R")
biocLite(c("VariantAnnotation", "Rsamtools"))

安装相关软件包，然后

library(VariantAnnotation) # also loads Rsamtools
snps = readVcf("/some/file.vcf")
param = ScanBamParam(which=rowData(vcf))
reads = countBam("/some/file.bam", param=param)

实现这一点的最佳方法可能很大程度上取决于您感兴趣的 SNP 数量。我建议您使用预发布的 R-2.15 alpha，因为您将获得一组更新的 Bioconductor 软件包。这些包有大量的小插图（

vignette(package="VariantAnnotation")

和 Bioconductor 邮件列表上知识渊博的人，以及常用的帮助页面

?readVcf

。

Answer 2

@马丁·摩根， 谢谢你的解释。能否详细说明一下命令，

“参数= ScanBamParam（其中= rowData（vcf））”

此命令中的“vcf”是什么？我按照你的命令操作，但此时收到错误反馈。

谢谢你。