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

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


r bioinformatics perl bioconductor vcf-variant-call-format
person UUU    schedule 13.03.2012    source источник
comment
вы должны спросить биостар: biostar.stackexchange.com   -  person Pierre    schedule 16.03.2012

Ответы (1)


arrow_upward
2
arrow_downward

Не знаю, что означает «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, так как вы получите более свежий набор пакетов Bioconductor. Эти пакеты имеют обширные виньетки (vignette(package="VariantAnnotation") и знающие люди из списка рассылки компании Bioconductor, а также обычные страницы справки ?readVcf.

person Martin Morgan    schedule 14.03.2012
comment
Спасибо за вашу помощь! Я попробую это. Под «количеством SNP над прочитанной позицией» я имею в виду, что я хочу иметь по оси X все основания чтения (в случае чтения Illumina 100 п.н.), а по оси Y — кумулятивное количество SNP, которые были находится на этой базовой позиции. Пример показан в этой статье, рисунок 2: biomedcentral.com/1471-2164/12/150. Можно ли делать такие вещи с указанными вами пакетами? - person UUU; 14.03.2012