Извлечь позицию чтения из файла bam

Question

У меня есть файл vcf, который содержит несколько SNP, и теперь я хочу посмотреть, равномерно ли распределены эти SNP по чтению файла bam, из которого я получил SNP. В частности, я хочу нанести на график количество SNP по положению чтения. Мне интересно, есть ли какой-нибудь инструмент для этого или я должен написать сценарий самостоятельно. Если да, то есть ли в R пакет, с которым я могу это сделать (я привык к R, но не имею большого опыта работы с Perl)?

Answer 1

Не уверен, что означает «SNPs over read position», но вы можете прочитать VCF с помощью пакета R / Bioconductor и функции VariantAnnotation :: readVcf и использовать геномные координаты для запроса файла bam с помощьюRsamtools :: countBam, используя 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.