r - 从 bam 文件中提取读取位置

Question

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

Answer 1

不确定“SNP超过读取位置”是什么意思，但您可以使用 R/Bioconductor 读取 VCF封装并调用 VariantAnnotation::readVcf 函数，并使用 ScanBamParam 使用基因组坐标通过 Rsamtools::countBam 查询 bam 文件。未经测试，沿着

## 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 mailing list 方面知识渊博的人员，以及常用的帮助页面 ?readVcf。