在

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

制定样品“HG00097”,“HG00101”的“LDAF”,“AVGPOST”信息字段,“GT”基因型字段的请求,用于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(或类似的专门查询的readGen或readInfo).

VariantAnnotation小插曲和参考手册中有大量文档;另请参见？ScanVcfParam;示例(ScanVcfParam).