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 アルファを使用することをお勧めします。そうすれば、より新しい Bioconductor パッケージのセットを取得できます。これらのパッケージには、通常のヘルプ ページと同様に、豊富なビネット ( vignette(package="VariantAnnotation")Bioconductorメーリング リスト?readVcfの知識豊富な人々) が含まれています。