--- title: "VCFArrow" author: "Tomas Hrbek" date: "06/06/2026" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{VCFArrow} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ``` # Installation The released version can be installed from CRAN with `install.packages("VCFArrow")`, and the development version from GitHub with `pak::pak("legalLab/VCFArrow")` or `devtools::install_github("legalLab/VCFArrow")`. Two packages are optional and only needed for specific functions: `adegenet`, for `vcf2genlight()`, and `svglite`, for SVG plots (without it, SVG plots are written with R's own Cairo device). Install them with `install.packages(c("adegenet", "svglite"))`. The package is built around Apache Arrow (https://arrow.apache.org/), accessed through the `arrow` R package. # Introduction This package has six main functionalities: (1) reading in and writing out VCF data files, including precalculation of diverse metrics used in the filtering functions; (2) data filtering functions; (3) merging VCFs; (4) subsetting VCFs; (5) generation of summary statistics and graphs; and (6) data transformation and writing out of other data formats ( population genetic and phylogenetic formats) for downstream analyses. This package is a reimagination of the `vcf2others` package. Same as `vcf2others`, it uses the Apache Arrow library functions to manipulate SNP data via dplyr verb equivalents. However, unlike `vcf2others`, `VCFArrow` now reads in the VCF data file into an S4 R object constructed around an Apache Arrow data structure. It precalculate a number of metrics and stores them in the S4 object. The VCF genotype matrix is stored as a long tidy table, and the data in the S4 object are lazy loaded. Many common tasks have been delegated to shared, C++ backed APIs and functions, and these operations have been parallelized when possible. These modifications and optimizations now permit extremely fast data filtering, manipulation and transformation. It is possible to process VCFs comprising tens of millions of SNVs and hundreds of individuals in matter minutes with minimal RAM overhead. ## Meta * Please [report here any issues or bugs or suggestions](https://github.com/legalLab/VCFArrow/issues). * License: MIT. * Get citation information for `VCFArrow` in R by running `citation(package='VCFArrow')`. # How to use the functions of this package Following are examples of the usage of the functions of this package. ## Parallel processing (optional) VCFArrow stores the genotypes on disk in chunks of variants, and most functions process one chunk at a time. `vcf_set_workers()` starts a number of background R processes (workers) and from then on distributes the chunks over them, so that several chunks are processed at the same time. This applies to `read_vcf()`, the genotype-based filters, sample and group extraction, `vcf_stats()` and `vcf_theta()`, the assess functions, `vcf_bind()`, `write_vcf()`, the `vcf2*()` exporters and `vcf2gt_long()` with CSV output. Functions that only use the variant table (for example the quality, PASS, biallelic or indel filters) are fast and always run in the main session. ```{r, eval = FALSE} # process chunks on 4 background R processes vcf_set_workers(workers = 4) # ... read, filter, compute statistics and export as usual ... # stop the workers and return to processing in the main session vcf_set_workers(1) ``` The workers are started once and reused by every following call until `vcf_set_workers(1)` is called or the R session ends. Results do not depend on the number of workers. Each worker is a separate R process that holds about 1 GB of memory, so memory use grows with the number of workers. Workers pay off for large datasets: on 18.6 million variants by 56 samples, most functions ran about twice as fast with 4 workers, while more workers than that brought little further gain because reading and writing the disk becomes the limit. For small datasets such as the one in this tutorial, starting the workers costs more time than they save, so the code above is not run here. ## Load VCF and associated files of individuals and groups First we need to read in the VCF. If one has information on the grouping of individuals, this information can also be put into the VCFArrow object. This grouping information is stored in a 'strata' file (a tab-separated text file with two columns with an 'id' and 'pop' headers). The grouping of individuals can be used in filtering steps, and it is necessary for transformation of the VCFArrow object to other population genetic formats, many of which require it. The data used in this example are from ([Mota et al., 2026](https://dx.doi.org/10.1080/14772000.2026.2692965))---a population genomic analysis of the *Phyllomedusa vaillantii* species complex. The original dataset has 18579683 SNVs by 56 samples. For the purpose of this tutorial, the dataset was randomly subsampled to 10000 SNVs by 18 samples (five samples by three main lineages of *Phyllomedusa vaillantii* as ingroups, plus three samples of *Phyllomedusa bicolor* included as outgroups). ![white-lined leaf frog _Phyllomedusa vaillantii_ Photo: Albertina Lima](https://github.com/legalLab/VCFArrow/blob/main/man/figures/P_vaillantii.jpeg) ```{r setup, message = FALSE} library(VCFArrow) library(dplyr) # set path to example files and project name data_path <- system.file("extdata", package="VCFArrow") project <- "vaillantii_" postfix <- "discosnp_sub" # load vcf and assign individuals to groups based on 'strata' vcf <- read_vcf(file.path(data_path, paste0(project, postfix, ".vcf.gz"))) |> set_vcf_groups(data_path) # visualize the VCFArrow object vcf # read individuals to include # if indivs is blank, the default is to use all individuals indivs <- read.table(file.path(data_path, "indivs_b"), header = TRUE)$id # check if all samples are in VCF sample names if (any(!(indivs %in% vcf@samples))) stop(paste("Some individuals in list not in VCF")) ``` ## Assessment of missing data The function `assess_vcf_missing_data()` is used to generate a table and graph of missing data for each sample in a VCFArrow object, while the `assess_vcf_coverage()` generates a violin plot of read depths of each sample. Both are useful for visualizing data before and after filtering to evaluate the effect of filtering. The functions accept parameters which form part of the name of the output file. The idea is that the output file name contains information on the name of the project (project), how the VCF was extracted (postfix), and how it was filtered (fltr). The 'postfix' and 'fltr' can be left blank if the VCF name does not contain this information. Samples are ordered and colored by group assignment. The 'details' flag specifies whether or not filtering information is presented in the figure or if only the species name is reported. The default is to report species name and details. ```{r} # define results path res_path <- tempdir() # species for plot title species <- "Phyllomedusa vaillantii" # filter - filter parameters in file name fltr <- "" assess_vcf_missing_data(vcf, res_path, species, paste0(project, postfix, fltr)) assess_vcf_coverage(vcf, res_path, species, paste0(project, postfix, fltr)) ``` ## Basic stats for individuals in VCF The function `vcf_stats()` is used to generate a table with basic stats for each sample in the VCFArrow object. These stats include per sample average read depth, heterozygosity, number of heterozygotes, homozygotes, REF homozygotes and ALT homogygotes, percent missing SNVs, total missing SNVs, total non-missing SNVs and total SNVs, Watterson's Theta and Pi for the entire dataset, and per group Watterson's Theta and Pi, plus the number of individuals in each group. The calculation of Watterson's Theta and Pi can be included or excluded based on the 'theta' flag. ``` {r} # get a table of basic sample stats, including Watterson's Theta and Pi vcf_stats(vcf, res_path, paste0(project, postfix, fltr), theta = TRUE) ``` ## Filtering, subsetting, merging and otherwise wrangling VCF files Most but not all functions can be used by the user directly on the `VCFArrow` object. Internally, all filters use two common functions, `.vcf_filter_rows()` for variants and `.vcf_filter_columns()` for samples, which are not exported. Other functions, such as `vcf_filter_invariant()` will remove invariant SNPs from a VCFArrow object; however, VCF by definition should not have invariant SNPs. So this function is primarily called by other functions to remove loci that may have become invariant following filtering. ### Subsetting samples and variants There are two types of subsetting functions. The first type is used to subset the VCFArrow object by specific samples `vcf_extract_samples()` or groups (all samples of a group) of samples `vcf_extract_groups()`. The 'keep' flag in both functions determines whether the indicated samples should be extracted from the VCF ('keep = TRUE') or excluded from the VCF ('keep = FALSE'). Once subsetted, the invariant SNVs are removed by default, however, this behavior is controlled by the 'f_invar' flag. The second type subsets a random number of SNVs from the VCF. The subsetting is done either across the entire VCF `vcf_sub_SNVs()` or proportionately within each chromosome/linkage group `vcf_sub_SNVs_stratified()`. For repeatability, both functions accept a 'seed'; if 'seed' is not specified, it is randomly generated. ```{r} # subset VCF by individuals (default keep = TRUE, remove invariants f_invar = TRUE) indivs1 <- c("Pv14", "Pv27", "Pv28", "Pv78", "Pv126") vcf1 <- vcf_extract_samples(vcf, indivs1) # subset VCF by groups (default keep = TRUE, remove invariants f_invar = TRUE) groups1 <- c("GS", "BS", "WA") vcf1 <- vcf_extract_groups(vcf, groups1) # subset VCF by a random number of variants (default n_SNVs = 10000) vcf1 <- vcf_sub_SNVs(vcf, n_SNVs = 5000) # subset VCF by a random number of variants proportional to chromosome length (default n_SNVs = 10000) # should be used when SNVs are mapped to chromosomes or scaffolds vcf1 <- vcf_sub_SNVs_stratified(vcf, n_SNVs = 5000) ``` ### Filtering Filtering consists of removing samples and variants based on specific properties. The function `vcf_filter_missing()` is the only function that filters samples based on % missing SNV threshold. The other functions filter on properties of individuals SNVs across samples. The function `vcf_filter_missingness()` removes SNVs above a % missing threshold, `vcf_filter_quality()` removes SNVs below a quality threshold, `vcf_filter_pass()` removes SNVs that did not PASS filters, `vcf_filter_maf()` removes SNVs with minor allele frequency below a % threshold, `vcf_filter_hets()` removes SNVs with heterozigosity above a % threshold, `vcf_filter_rank()` removes SNPs with rank below a % threshold (rank is specific to VCFs generated by DiscoSnpRad and measures the likelihood of observed SNP variation being due to paralogs), `vcf_filter_rpt()` removes SNPs with a repeat score above a threshold (repeat is also specific to DiscoSnpRad VCFs; high values indicate paralogs), `vcf_filter_paralogs()` applies both the rank and the repeat filter, `vcf_filter_adr()` either removes (make genotype missing) or corrects (make genotype homozygous) genotypes with allele ratio imbalance, `vcf_filter_coverage()` removes SNVs below a specific coverage (make genotypes missing), and `vcf_filter_biallelic()` removes all variants that are not biallelic. Finally the `vcf_filter_invariant()` removes invariant SNVs (SNVs may become invariant after subsetting and filtering). The function `vcf_filter_indels()` removes indels from the VCF, retaining only SNPs. The functions `vcf_filter_oneSNV()` and `vcf_filter_multiSNV()` filter the VCF to remove linked SNVs (retain only unliked SNPs) or retain only linked SNVs, respectively. In general, all analyses assume that SNVs are unlinked, however, some analyses such as FineRadStructure work with linked SNVs. ```{r} # filter VCF for analyses (unlinked SNPs) vcf_oneSNP <- vcf_extract_samples(vcf, indivs) |> vcf_filter_indels() |> vcf_filter_biallelic() |> vcf_filter_maf(.03) |> vcf_filter_coverage(10) |> vcf_filter_missingness(.2) |> vcf_filter_missing(.3) |> vcf_filter_oneSNV() # see how many variants remained nrow(vcf_oneSNP@variants) # filter VCF for analyses (linked SNPs) vcf_multiSNP <- vcf_extract_samples(vcf, indivs) |> vcf_filter_indels() |> vcf_filter_biallelic() |> vcf_filter_maf(.03) |> vcf_filter_coverage(10) |> vcf_filter_missingness(.2) |> vcf_filter_missing(.3) |> vcf_filter_multiSNV() # see how many variants remained nrow(vcf_multiSNP@variants) ``` ### Extracting and binding The functions `vcf_bind()`, `vcf_extract_samples()` and `vcf_extract_groups()` are used to bind two or more VCFArrow objects, and to extract samples/groups of samples from a VCFArrow objects storing them in a separate VCFArrow object. Variants are matched by `CHROM`, `POS`, `REF` and `ALT`, so the objects can also come from different VCF files, for example from separate sequencing runs. Matching requires each locus to be described the same way in every dataset. This is assured when all datasets were called against the same reference genome assembly, which gives the same `CHROM`, `POS` and `REF` for the same locus. It is generally not possible for de novo (reference-free) SNP calling, where locus names and positions are generated anew in each calling run, so the same locus will usually not match between separately called datasets; such datasets should be called together in a single run instead. Multi-allelic sites can list their `ALT` alleles in a different order in different files, so it is best to normalize the VCFs first, for example by splitting multi-allelic sites with `bcftools norm -m -any`. ```{r} # extract individuals from VCFArrow object (keep all loci) indivs1 <- c("Pb2Jp", "Pb2Scx", "Pb1Rd") vcf_outgrp <- vcf_extract_samples(vcf, indivs1, f_invar = FALSE) vcf_ingrp <- vcf_extract_samples(vcf, indivs1, keep = FALSE, f_invar = FALSE) # extract groups of individuals from VCFArrow object (keep all loci) groups1 <- c("GS", "BS", "WA") vcf_outgrp <- vcf_extract_groups(vcf, groups1, keep = FALSE, f_invar = FALSE) vcf_ingrp <- vcf_extract_groups(vcf, groups1, f_invar = FALSE) # bind vcf_outgrp and vcf_ingrp (same variants) vcf1 <- vcf_bind(vcf_ingrp, vcf_outgrp) # extract groups of individuals from VCFArrow object (remove invariant loci) groups1 <- c("GS", "BS", "WA") vcf_outgrp <- vcf_extract_groups(vcf, groups1, keep = FALSE) vcf_ingrp <- vcf_extract_groups(vcf, groups1) # bind vcf_outgrp and vcf_ingrp (intersection - variable SNVs in all objects) vcf1 <- vcf_bind(vcf_ingrp, vcf_outgrp, mode = "intersect") # bind vcf_outgrp and vcf_ingrp (union - variable SNVs in any object) vcf1 <- vcf_bind(vcf_ingrp, vcf_outgrp, mode = "union", absent_as = "missing") ``` With `mode = "union"`, loci removed from an object stay removed for it: a locus that one object removed by filtering (for example with `vcf_filter_maf()`) but that another object still contains is included in the result, and the genotypes of the object that removed it are coded according to `absent_as` (missing, `./.`, by default). Binding in steps, `vcf_bind(vcf_bind(A, B, mode = "union"), C, mode = "union")`, and binding in one call, `vcf_bind(A, B, C, mode = "union")`, therefore give the same result. To join datasets that were filtered separately before the join, set `recover_loci = TRUE`: the real genotypes of removed loci are then restored wherever an object's genotype files still hold them. With `recover_loci = TRUE`, binding in steps can differ from binding in one call, because an intermediate object keeps only the loci it contains, so bind all objects in a single call. ### Filtering with outgroup Outgroups, which tend to have many fewer individuals than ingroups, also have more missing data than ingroups due to phylogenetic effects. Filtering the entire dataset often results in the removal of the outgroup taxa when `vcf_filter_missing()` is used during filtering to remove individuals with % missing data above some threshold. It also results in fewer SNVs being retained after filtering. The solution is to extract the outgroup taxa with `vcf_extract_samples()` or `vcf_extract_groups()`, filter the ingroup VCFArrow object, and add the outgroup taxa with `vcf_bind()` using `outgroup = TRUE`. With `outgroup = TRUE`, the result keeps exactly the loci of the filtered ingroup, which must be the first object: loci removed from the ingroup stay removed, and properties of the outgroup, such as its missing data, cannot remove ingroup loci. The outgroup genotypes are added for these loci wherever they exist, including loci the outgroup dropped as invariant, so it does not matter whether the outgroup was extracted with `f_invar = TRUE` or `FALSE`. By contrast, `mode = "union"` also keeps loci found only in the outgroup (see above). ```{r} # extract outgroup from VCFArrow groups1 <- "OG" vcf_outgrp <- vcf_extract_groups(vcf, groups1) # filter VCFArrow ingroup, then bind outgroups vcf1 <- vcf_extract_groups(vcf, groups1, keep = FALSE) |> vcf_filter_indels() |> vcf_filter_biallelic() |> vcf_filter_coverage(10) |> vcf_filter_missingness(.2) |> vcf_filter_missing(.3) |> vcf_filter_oneSNV() |> vcf_bind(vcf_outgrp, outgroup = TRUE, absent_as = "missing") # see how many variants remained nrow(vcf1@variants) # apply the same filter to the full dataset (ingroup + outgroup) vcf1 <- vcf |> vcf_filter_indels() |> vcf_filter_biallelic() |> vcf_filter_coverage(10) |> vcf_filter_missingness(.2) |> vcf_filter_missing(.3) |> vcf_filter_oneSNV() # see how many variants remained nrow(vcf1@variants) ``` ## Saving and copying VCFArrow objects The filtering steps never actually remove any sample or SNVs from the GT matrix which is a design decision. If you want a new VCFArrow object that only has the post filtering samples and SNVs, you can use the function `vcf_copy()` which creates a new object and removes any filtered samples and SNVs. After filtering the VCFArrow object you will also want to save it as a VCF. This is accomplished using the function `write_vcf()`, which, by default will save an uncompressed VCF, but setting the 'gzip' flag to TRUE will gzip compress the VCF. With large VCFArrow objects this has a large overhead, so it is more efficient to save an uncompressed VCF and then compress it later. ```{r} # write the VCFArrow object as a VCF (and compress it using gzip) write_vcf(vcf1, out_file = file.path(res_path, paste0(project, postfix, "_filtered.vcf.gz")), gzip = TRUE) ``` ## Converting a VCF file to other population genetic and phylogenetic formats All the following functions will take a VCFArrow object, and convert it other populations genetic and phylogenetic formats, exporting/writing a file of this format. Group/population information is extracted from the VCFArrow object for those formats that require it. The function `vcf2genlight()` automatically returns a genlight object and optionally can also save it as an .rds file given by `out_file`. Generally the `vcf2genlight()` function is called within a script using functions of the `adegenet` and `poppr` packages rather than importing the genlight object. `vcf2genlight()` requires the `adegenet` package. ```{r} res_path <- tempdir() project <- "vaillantii_" postfix <- "discosnp_sub" fltr <- "" ########## # datasets for analyses with unlinked SNPs vcf <- vcf_oneSNP ########## # export data formats # migrate-n https://peterbeerli.com/migrate-html5/ vcf2migrate(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_migrate.txt'))) # arlequin http://cmpg.unibe.ch/software/arlequin35/ vcf2arlequin(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.arp'))) # structure https://web.stanford.edu/group/pritchardlab/structure.html vcf2structure(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.str'))) # faststucture http://rajanil.github.io/fastStructure/ vcf2structure(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.fstr')), method = "F") # sNMF http://membres-timc.imag.fr/Olivier.Francois/snmf/index.htm vcf2snmf(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.geno'))) # Admixture https://dalexander.github.io/admixture/index.html # takes as input binary PLINK (.bed), ordinary PLINK (.ped), or EIGENSTRAT (.geno) formatted files # so a thin wrapper around vcf2plink_bed() plus optional generation of known-ancestry reference populations vcf2admixture(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_admixture'))) # PLINK .bed https://www.cog-genomics.org/plink/1.9/formats#bed vcf2plink_bed(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_plink'))) # PLINK .ped https://www.cog-genomics.org/plink/1.9/formats#ped vcf2plink_ped(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_plink'))) # eigenstrat https://github.com/DReichLab/EIG/tree/master vcf2eigenstrat(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_eigenstrat'))) # genepop https://gitlab.mbb.univ-montp2.fr/francois/genepop vcf2genepop(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.gen'))) # smartsnp https://github.com/ChristianHuber/smartsnp vcf2smartsnp(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.smartsnp'))) # bayescan https://github.com/mfoll/BayeScan vcf2bayescan(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.bayescan'))) # bayesass https://github.com/brannala/BA3 vcf2bayesass(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.bayesass'))) # treemix https://bitbucket.org/nygcresearch/treemix/wiki/Home vcf2treemix(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.treemix'))) # apparent https://github.com/halelab/apparent/tree/master vcf2apparent(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.apparent'))) # related https://github.com/timothyfrasier/related vcf2related(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.related'))) # long tidy dataframe of genotypes vcf2gt_long(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.csv')), format = 'csv') # snapp https://www.beast2.org/snapp/ vcf2snapp(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_snapp.nex'))) # nexus - only SNPs, meant for SVDq analyses https://www.asc.ohio-state.edu/kubatko.2/software/SVDquartets/ vcf2nexus(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_sdvq.nex'))) # fasta https://www.ncbi.nlm.nih.gov/genbank/fastaformat/ vcf2fasta(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.fna'))) # genlight object https://www.rdocumentation.org/packages/adegenet/versions/2.0.0/topics/genlight-class genlight <- vcf2genlight(vcf) # genlight object with an optional save genlight <- vcf2genlight(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_genlight.rds')), save = TRUE) ########## # datasets for analyses with linked SNPs vcf <- vcf_multiSNP ########## # export data formats # fineRadStructure - expects VCF of linked SNPs vcf2fineradstructure(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.finerad'))) ```