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

Question

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

Answer 1

不确定“SNPs over read position”是什么意思，但您可以使用 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。