У меня есть большой файл формата Variant Call (VCF) ( > 4 ГБ), который имеет данные для нескольких образцов.
Я просмотрел Google, Stackoverflow, а также попробовал пакет VariantAnnotation в R, чтобы каким-то образом извлечь данные только для определенного образца, но не нашел никакой информации о том, как это сделать в R.
Кто-нибудь пробовал что-нибудь подобное или, может быть, знал о другом пакете, который бы включил это?
В VariantAnnotation используйте ScanVcfParam, чтобы указать данные, которые вы хотите извлечь. Использование образца файла VCF, включенного в пакет
library(VariantAnnotation)
vcfFile = system.file(package="VariantAnnotation", "extdata", "chr22.vcf.gz")
узнать информацию о файле
scanVcfHeader(vcfFile)
## class: VCFHeader
## samples(5): HG00096 HG00097 HG00099 HG00100 HG00101
## meta(1): fileformat
## fixed(0):
## info(22): LDAF AVGPOST ... VT SNPSOURCE
## geno(3): GT DS GL
Сформулировать запрос для информационных полей "LDAF", "AVGPOST", поля генотипа "GT" для образцов "HG00097", "HG00101" для вариантов на хромосоме 22 между координатами 50300000, 50400000
param = ScanVcfParam(
info=c("LDAF", "AVGPOST"),
geno="GT",
samples=c("HG00097", "HG00101"),
which=GRanges("22", IRanges(50300000, 50400000)))
Прочитайте запрошенные данные
vcf = readVcf(vcfFile, "hg19", param=param)
и извлечь из VCF соответствующие данные
head(geno(vcf)[["GT"]])
## HG00097 HG00101
## rs7410291 "0|0" "0|0"
## rs147922003 "0|0" "0|0"
## rs114143073 "0|0" "0|0"
## rs141778433 "0|0" "0|0"
## rs182170314 "0|0" "0|0"
## rs115145310 "0|0" "0|0"
head(info(vcf)[["LDAF"]])
## [1] 0.3431 0.0091 0.0098 0.0062 0.0041 0.0117
ranges(vcf)
## IRanges of length 1169
## start end width names
## [1] 50300078 50300078 1 rs7410291
## [2] 50300086 50300086 1 rs147922003
## [3] 50300101 50300101 1 rs114143073
## [4] 50300113 50300113 1 rs141778433
## [5] 50300166 50300166 1 rs182170314
## ... ... ... ... ...
## [1165] 50364310 50364312 3 22:50364310_GCA/G
## [1166] 50364311 50364313 3 22:50364311_CAT/C
## [1167] 50364464 50364464 1 rs150069372
## [1168] 50364465 50364465 1 rs146661152
## [1169] 50364609 50364609 1 rs184235324
Возможно, вас интересует только элемент генотипа "GS" как простая матрица R, а затем просто укажите образцы и/или диапазоны, которые вас интересуют, и используйте readGeno (или readGT или readInfo для аналогичные специализированные запросы).
В VariantAnnotation в документации и справочном руководстве имеется обширная документация; см. также ?ScanVcfParam; example(ScanVcfParam).
Я встретил ту же ситуацию. Следующая ссылка может предоставить возможное решение: http://biocluster.ucr.edu/~rsun/workshop/SNPINDEL/variantsAnalysis.pdf На странице, описывающей "Чтение файла VCF", они не используют набор данных в некотором пакете, а фактический путь для локального файла gz:
Библиотека (VariantAnnotation)
Библиотека (GenomicFeatures)
Пожалуйста, убедитесь, что точный data.zip под вашим текущим рабочим каталогом
vcf < - readVcf ( "data/var.raw.vcf.gz", "TAIR10" )
VCF
TAIR10 является эталонным геномом