1简介

VCFArray是一个Bioconductor包，将VCF文件表示为派生自DelayedArray包和DelayedArray类。它将数据项从VCF文件转换为一个DelayedArray-派生数据结构。后端VCF文件可以保存在本地磁盘上，也可以作为在线资源保存在远程磁盘上。可提取的数据项包括固定数据字段(REF, ALT, QUAL, FILTER)、信息字段(如AA, AF…)和单个格式字段(如GT, DP…)。由固定字段/信息字段生成的数组数据是一维的VCFArray，其中维度为变量的长度。由个体生成的数组数据格式字段总是返回第一个维度变体第二个维度是样品．这一特征与保存的化验数据一致SummarizedExperiment，并使VCFArray包可与其他已建立的包互操作Bioconductor数据基础设施。

2安装

  1. 下载软件包。
如果(!requireNamespace("BiocManager"， quiet = TRUE)) install.packages("BiocManager")::install("VCFArray")

开发版也可以从Github下载。

BiocManager:安装(“Bioconductor / VCFArray”)
  1. 将包加载到R会话。
库(VCFArray)

3.VCFArray

3．1VCFArray构造函数

要构造一个VCFArray对象，需要4个参数:文件，vindex而且名字,pfix．的文件参数可以取字符串(VCF文件名)，或者VcfFile对象，或RangedVcfStack对象。名字参数必须指定，以指示希望从输入文件中提取哪个数据项。它区分大小写，并且必须与VCF头文件中的名称一致。vindex参数如果不存在，则仅用于指示索引文件的文件路径。pfix用于指定类别名字字段属于。请注意这一pfix需要特别提供时，有相同的名字在多个类别中，否则将返回错误。

的vcfFields ()方法取VCF文件路径，VcfFile对象或RangedVcfStack对象作为输入，并返回一个包含特定类别中所有可用VCF字段的CharacterList。用户应查阅固定，信息而且基因族群可转换为的可用数据项的类别VCFArray实例。数据项名称可用作名字论点VCFArray构造函数。

args(VCFArray) #>函数(文件，vindex =字符()，name = NA, pfix = NULL) #> NULL fl <-系统。文件("extdata", "chr22.vcf.gz", package = "VariantAnnotation") library(VariantAnnotation) vcfFields(fl) #> CharacterList of length 4 #> [["fixed"]] REF ALT QUAL FILTER #> [["info"]] LDAF AVGPOST RSQ ERATE THETA ... ASN_AF AFR_AF EUR_AF VT SNPSOURCE #> [["geno"]] GT DS GL #> [["samples"]] HG00096 HG00097 HG00099 HG00100 HG00101

由于vcf文件的索引文件已经存在，所以vindex参数将不需要(这是磁盘上VCF文件最常见的情况)。所以我们可以构造VCFArray对象的GT提供的VCF文件中的数据项，参数为文件而且名字只有。

VCFArray(file = fl, name = "GT") #> <10376 × 5>矩阵类VCFMatrix和类型"character": #> HG00096 HG00097 HG00099 HG00100 HG00101 # b> rs7410291 "0|0" "0|0" "1|0" "0|0" "0|0" #> rs147922003 "0|0" "0|0" "0|0" "0|0" "0|0" " #> rs114143073 "0|0" "0|0" "0|0" "0|0" "0|0" #> ... ... ..# > rs5770892“| 0”“| 0”“| 0”“| 0”“0 | 1 # > rs144055359“| 0”“| 0”“| 0”“| 0”“| 0”# > rs114526001“| 0”“| 0”“| 0”“| 0”“| 0”

我们也可以构造aVCFArray对象的文件参数是VcfFile对象。

vcf <- VariantAnnotation::VcfFile(fl) VCFArray(file = vcf, name = "DS") #> <10376 x 5>矩阵类VCFMatrix和类型“double”:#> HG00096 HG00097 HG00099 HG00100 HG00101 #> rs7410291 000 000 0 #> rs147922003 000 000 #> rs114143073 000 000 0 #> ... ... ..#> rss5770892 00 001 #> rs144055359 00 00 00 #> rs114526001 00 00 00 0

的文件论证也可以取RangedVcfStack对象。注意，一个普通的VcfStack对象没有范围这些信息不能用来构建一个VCFArray．

Extdata <- system. Extdata <- system. Extdata。文件(package = "GenomicFiles", "extdata") files <- dir(extdata, pattern="^CEUtrio.*bgz$", full=TRUE)[1:2] names(files) <- sub(".*_([0-9XY]+).*", "\\1", basename(files)) seqinfo <- as(readRDS(file.path(extdata, "seqinfo.rds")), "Seqinfo") stack <- GenomicFiles::VcfStack(files, seqinfo) gr <- as(GenomicFiles::seqinfo(stack)[rownames(stack)], "GRanges") ## RangedVcfStack rgstack <- GenomicFiles::RangedVcfStack(stack, rowRanges = gr) rgstack #> VcfStack object with 2 files and 3 samples #> GRanges object with 2 ranges and 0 metadata columns #> Seqinfo object with 25 sequences from hg19 genome #> use 'readVcfStack()' to extract VariantAnnotation VCF.

这里我们选择名字= SB，它返回一个三维图像VCFArray对象，其中前2个维度分别对应变体和样本。

vcfFields (rgstack) $基因族群# >[1]“GT”“广告”“DP”“GQ”“MIN_DP”“页面表”“PID”“PL”# >[9]“某人”VCFArray (name = "某人" rgstack) # > < 318 x 3 x 4 >类VCFArray和类型的数组“整数”:# >,,1 # > [1][2][3]# > [1]NA NA NA # > [2] NA NA NA  #> ... . . .#> [317，] na na na #> [318，] na na na #> #>…# > # >, 4 # > [1] [2] [3] # > [1] NA NA NA # > [2] NA NA NA  #> ... . . .#> [318，] na na na

正如小插图标题所示，后端VCF文件也可以是远程文件。在这里，我们包括了一个表示1000基因组计划(第3阶段)22号染色体VCF文件的例子。注意对于远程VCF文件，使用vindex参数必须指定。由于这个VCF文件比较大，而且花费的时间比较长，所以这里只显示代码，没有计算。

chr22url <- "https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz" chr22url。tbi <- paste0(chr22url， ".tbi") va <- VCFArray(chr22url, vindex =chr22url. tbi <- paste0(chr22url， ".tbi")tbi, name = "GT")

3.2VCFArray方法

VCFArray表示VCF文件为DelayedArray实例。它有这样的方法昏暗的，dimnames的类数组操作和方法DelayedArray的细分方法［．
请注意对于一维空间VCFArray由VCF文件的固定/信息数据字段生成的对象，drop = FALSE应该总是与［细分以确保VCFArray对象作为返回值。

3.2.1之上槽访问器

种子返回VCFArraySeed的VCFArray对象，其中包括有关后端VCF文件的信息，例如VCF文件路径、索引文件路径、数据项的名称(以类别为前缀)、维度等。

va <- VCFArray(fl, name = "GT") seed(va) #> VCFArraySeed #> VCF文件路径:/home/biocbuild/bb -3.15-bioc/R/library/VariantAnnotation/extdata/chr22.vcf.gz #> VCF索引路径:/home/biocbuild/bb -3.15-bioc/R/library/VariantAnnotation/extdata/chr22.vcf.gztbi #> array data: GT #> dim: 10376 × 5

vcffile返回VcfFile对象对应于后端VCF文件。

vcffile(va) #>类:vcffile #>路径:/home/biocbuild/bbs-3.15-bioc/R/library/VariantAnnotation/…/chr22.vcf.gz #>索引:/home/biocbuild/bbs-3.15-bioc/R/library/VariantAnnota…tbi #> isOpen: FALSE #> yieldSize: NA

3.2.2昏暗的()而且dimnames ()

的dimnames (VCFArray)返回一个未命名的列表，每个元素的长度与return from相同暗(VCFArray)．

va <- VCFArray(fl, name = "GT") dim(va) #> [1] 10376 5 class(dimnames(va)) #> [1] "list" length (dimnames(va)) #> [1] 10376 5

3.2.3［构造子集

VCFArray实例可以按照通常的方法进行细分R数字或逻辑向量的约定;逻辑向量被循环到适当的长度。

va[1:3, 1:3] #> <3 × 3>矩阵类DelayedMatrix和类型"character": #> HG00096 HG00097 HG00097 hgb15 0" #> rs147922003 "0|0" "0|0" "0|0" #> rs114143073 "0|0" "0|0" "0|0" va[c(TRUE, FALSE)，] #> <5188 × 5>矩阵类DelayedMatrix和类型"character": #> HG00096 HG00097 HG00099 HG00100 HG00101 #> rs7410291 "0|0" "0|0" "1|0" "0|0" "0|0" " #> rs114143073 "0|0" "0|0" "0|0" "0|0" "0|0" #> rs182170314 "0|0" "0|0" "0|0" "0|0" # b> ... ... ..# > rs9628212“| 0”“| 0”“| 0”“| 0”“| 0”# > rs9628178“| 0”“| 0”“| 0”“| 0”“| 0”# > rs144055359“| 0”“| 0”“| 0”“| 0”“| 0”

3.2.4一些数值计算

数值计算可以评估VCFArray对象。

ds <- VCFArray(fl, name = " ds ") log(ds+5) #> <10376 x 5>矩阵类DelayedMatrix和类型“double”:#> HG00096 HG00097 HG00099 HG00100 HG00101 #> rs7410291 1.609438 1.609438 1.791759 1.609438 1.609438 1.609438 1.609438 #> rs11414143073 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 #> ... ... ..#> rss5770892 1.609438 1.609438 1.609438 1.791759 #> rs144055359 1.609438 1.609438 1.609438 1.609438 1.609438 #> rs114526001 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438 1.609438

4内部:VCFArraySeed

的VCFArraySeed类的“种子”VCFArray对象。它不是从VCFArray包中。种子对象应包含VCF文件路径，并期望满足的“种子契约”DelayedArray，即支持昏暗的()而且dimnames ()．

- VCFArray:::VCFArraySeed(fl, name = "GT"， pfix = NULL) seed #> VCFArraySeed #> VCF文件路径:/home/biocbuild/bbs-3.15-bioc/R/library/VariantAnnotation/extdata/chr22.vcf.gz #> VCF索引路径:/home/biocbuild/bbs-3.15-bioc/R/library/VariantAnnotation/extdata/chr22.vcf.gztbi #> array data: GT #> dim: 10376 x 5 path(vcffile(seed)) #> [1] "/home/biocbuild/bb -3.15-bioc/R/library/VariantAnnotation/extdata/chr22.vcf.gz"

种子可以用来构建一个VCFArray实例。

(va <- VCFArray(seed)) #> <10376 × 5>矩阵类VCFMatrix和类型“character”:#> HG00096 HG00097 HG00099 HG00100 HG00101 # b> rs7410291“0|0”“0|0”“1|0”“0|0”“0|0”#> rs147922003“0|0”“0|0”“0|0”“0|0”“0|0”#> rs114143073“0|0”“0|0”“0|0”“0|0”“0|0”#> ... ... ..# > rs5770892“| 0”“| 0”“| 0”“| 0”“0 | 1 # > rs144055359“| 0”“| 0”“| 0”“| 0”“| 0”# > rs114526001“| 0”“| 0”“| 0”“| 0”“| 0”

的DelayedArray ()构造函数与VCFArraySeed对象作为输入返回的内容与VCFArray ()同样的构造函数VCFArraySeed．

da <- DelayedArray(seed) class(da) #> [1] "VCFMatrix" #> attr(，"package") #> [1] "VCFArray" all。= (da, va) #> [1] TRUE

5sessionInfo ()

sessionInfo() #> R version 4.2.0 RC (2022-04-19 r82224) #>平台:x86_64-pc-linux-gnu (64-bit) #>运行在:Ubuntu 20.04.4 LTS #> #>矩阵产品:默认#> BLAS: /home/biocbuild/bbs-3.15-bioc/R/lib/libRblas。所以#> LAPACK: /home/biocbuild/bbs-3.15-bioc/R/lib/libRlapack。so #> #> locale: #> [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# > [3] LC_TIME=en_GB LC_COLLATE= c# > [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 #> [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# > [9] LC_ADDRESS=C LC_TELEPHONE= c# > [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION= c# > #>附加基础包:#> [1]stats4 stats graphics grDevices utils datasets methods #>[8]基础#> #>其他附加包:#> [1] VariantAnnotation_1.42.0 Rsamtools_2.12.0 #> [3] Biostrings_2.64.0 XVector_0.36.0 #> [5] SummarizedExperiment_1.26.0 Biobase_2.56.0 #> [9] VCFArray_1.12.0 DelayedArray_0.22.0 #> [11] IRanges_2.30.0 S4Vectors_0.34.0 #> [13] MatrixGenerics_1.8.0 matrixStats_0.62.0 #> [15] Matrix_1.4-1 BiocGenerics_0.42.0 #> [17] BiocStyle_2.24.0 #> #>通过命名空间加载(并且没有附加):#> [13] yaml_2.3.5 progress_1.2.2 pillar_1.7.0 #> [16] RSQLite_2.2.12 lattice_0.20-45 glue_1.6.2 #> [19] digest_0.6.29 htmltools_0.5.2 XML_3.99-0.9 #> [22] pkgconfig_2.0.3 biomaRt_2.52.0 bookdown_0.26 #> [25] zlibbioc_1.42.0 purrr_0.3.4 BiocParallel_1.30.0 #> [28] tibble_3.1.6 # [7] assertthat_1.2.1 BiocManager_1.30.17 BiocFileCache_2.4.0 #> [10] blob_1.2.3 BSgenome_1.64.0KEGGREST_1.36.0 generics_0.1.2 #> [31] ellipsis_0.3.2 cachem_48.0 #> [34] cli_3.3.0 magrittr_2.0.3 crayon_1.5.1 #> [37] memoise_2.0.1 evaluate_0.15 fansi_1.0.3 #> [40] xml2_1.3.3 tools_4.2.0 prettyunits_1.1.1 #> [43] hms_1.1.1 BiocIO_1.6.0 lifecycle_1.0.1 #> [46] string_1 .4.0 AnnotationDbi_1.58.0 compiler_4.2.0 #> [52] RCurl_1.98-1.6 rjson_0.2.21 rmarkdown_2.14 restfulr_0.0.13 #> > > > > [49] jquerylib_0.1.4 rlang_1.0.2 grid_4.2.0 #> [52] RCurl_1.98-1.6 rjson_1.0 -7 rmarkdown_2.14 restfulr_0.0.13 #> > [58]curl_4.3.2 DBI_1.1.2 R6_2.5.1 #> [61] GenomicAlignments_1.32.0 rtracklayer_1.56. 6.0 knitr_1.38 #> [64] dplyr_1.0.8 utf8_1.2.2 fastmap_1.1.0 #> [67] bit_4.0.4 filelock_1.0.2 stringi_1.7.6 #> [70] parallel_4.2.0 Rcpp_1.0.8.3 vctrs_0.4.1 #> [73] png_0.1-7 tidyselect_1.1.2 dbplyr_2.1.1 #> [76] xfun_0.30