VCFArrow
Fast and Memory-Efficient Manipulation and Conversion of VCF Files
Reads Variant Call Format (VCF) files into an S4 object whose genotypes are stored on disk as an 'Apache Arrow' dataset and loaded lazily, allowing large datasets to be processed with little memory. Provides functions to filter variants and samples, subset and merge datasets, compute summary statistics and diagnostic plots, write VCF files, and convert data to input formats of population genetic and phylogenetic software such as 'STRUCTURE', 'ADMIXTURE', 'PLINK', 'EIGENSTRAT', 'Arlequin', 'GENEPOP', 'BayeScan', 'TreeMix', 'SNAPP' and 'adegenet'.
VCFArrow 
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.
- License: MIT.
- Get citation information for
VCFArrowin R by runningcitation(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.
# 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)—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).
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
#>
#> An object of class "VCFArrow"
#>
#> Dimensions:
#> Variants: 10000
#> Samples: 18
#>
#> Quick stats:
#> Non-missing variants: 10000
#>
#> Phased genotypes: FALSE
#>
#> Storage:
#> Path: /tmp/RtmpOHmUZa/arrow_vcf_7e345bc28489
#>
#> Genotype storage (Arrow):
#> FileSystemDataset with 1 Feather file
#> 9 columns
#> .row_id: int32
#> sample: string
#> a1: int32
#> a2: int32
#> phased: bool
#> fmt: string
#> DP: double
#> GQ: double
#> ADR: double
#>
#> See $metadata for additional Schema metadata
#>
#> Variants (first 5 rows):
#> CHROM POS ID REF ALT QUAL FILTER Rk RPT n_alt
#> 1 SNP_higher_path_9994239 41 9994239_1 C G . . 1 NA 1
#> 2 SNP_higher_path_9984432 105 9984432 C T . . 1 NA 1
#> 3 SNP_higher_path_9967574 50 9967574 A C . . 1 NA 1
#> 4 SNP_higher_path_993510 88 993510_4 A G . . 1 NA 1
#> 5 SNP_higher_path_9803974 33 9803974 A G . . 1 NA 1
#> is_biallelic is_indel .row_id
#> 1 TRUE FALSE 1
#> 2 TRUE FALSE 2
#> 3 TRUE FALSE 3
#> 4 TRUE FALSE 4
#> 5 TRUE FALSE 5
#> ... 9995 more
#>
#> INFO (first 5):
#> [1] "Ty=SNP;Rk=1.0;UL=11;UR=8;CL=11;CR=27;Genome=.;Sd=.;Cluster=20022615;ClSize=3"
#> [2] "Ty=SNP;Rk=1.0;UL=3;UR=1;CL=75;CR=7;Genome=.;Sd=.;Cluster=17296599;ClSize=1"
#> [3] "Ty=SNP;Rk=1.0;UL=20;UR=4;CL=20;CR=4;Genome=.;Sd=.;Cluster=18266172;ClSize=3"
#> [4] "Ty=SNP;Rk=1.0;UL=1;UR=0;CL=1;CR=0;Genome=.;Sd=.;Cluster=16400006;ClSize=2"
#> [5] "Ty=SNP;Rk=1.0;UL=3;UR=7;CL=3;CR=7;Genome=.;Sd=.;Cluster=10880986;ClSize=5"
#>
#> FORMAT (first 5):
#> FORMAT .row_id
#> 1 GT:DP:PL:AD:HQ 1
#> 2 GT:DP:PL:AD:HQ 2
#> 3 GT:DP:PL:AD:HQ 3
#> 4 GT:DP:PL:AD:HQ 4
#> 5 GT:DP:PL:AD:HQ 5
#>
#> Samples (first 5):
#> [1] "Pv120" "Pv126" "Pv13" "Pv14" "Pv27"
#> ... 13 more
# 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.
# 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.
# get a table of basic sample stats, including Watterson's Theta and Pi
vcf_stats(vcf, res_path, paste0(project, postfix, fltr), theta = TRUE)
#> ℹ Computing per-sample stats: 10000 variants x 18 samples, reading 1 chunk directly
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.
# 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)
#> ℹ Compacting GT: 18 -> 5 samples across 1 chunk
#> ℹ Applying invariant filter
#> ℹ Removed 7220 invariant variants; 2780 retained.
#> ℹ Removed samples: Pv120, Pv13, Pv2, Pv31, Pv56, Pv62, Pv68, Pv73, Pv79, Pv93, Pb2Jp, Pb2Scx, and Pb1Rd
#> ℹ Variants retained: 2780 | Samples retained: 5
# subset VCF by groups (default keep = TRUE, remove invariants f_invar = TRUE)
groups1 <- c("GS", "BS", "WA")
vcf1 <- vcf_extract_groups(vcf, groups1)
#> ℹ Compacting GT: 18 -> 15 samples across 1 chunk
#> ℹ Applying invariant filter
#> ℹ Removed 1705 invariant variants; 8295 retained.
#> ℹ Removed samples: Pb2Jp, Pb2Scx, and Pb1Rd
#> ℹ Variants retained: 8295 | Samples retained: 15
# 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.
# 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()
#> ℹ No samples to keep or remove - keeping all samples
#> ℹ Applying indel filter
#> ℹ Retained 9313 / 10000 variants (non-Indels)
#> ℹ Applying biallelic filter
#> ℹ Retained 9313 / 9313 variants (biallelic SNVs)
#> ℹ Applying MAF filter
#> ℹ Retained 9069 / 9313 variants (MAF >= 0.03)
#> ℹ Applying read coverage filter
#> ℹ Retained 5171 / 9069 variants (polymorphic with DP >= 10)
#> ℹ Applying locus missingness filter
#> ℹ Retained 1734 / 5171 variants (per-variant missingness <= 0.2)
#> ℹ Applying sample missingness filter
#> ℹ Variants retained: 1734 | Samples retained: 18
#> ℹ Applying unlinked SNV filter
#> ℹ Retained 1734 / 1734 variants (unlinked SNVs)
# see how many variants remained
nrow(vcf_oneSNP@variants)
#> [1] 1734
# 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()
#> ℹ No samples to keep or remove - keeping all samples
#> ℹ Applying indel filter
#> ℹ Retained 9313 / 10000 variants (non-Indels)
#> ℹ Applying biallelic filter
#> ℹ Retained 9313 / 9313 variants (biallelic SNVs)
#> ℹ Applying MAF filter
#> ℹ Retained 9069 / 9313 variants (MAF >= 0.03)
#> ℹ Applying read coverage filter
#> ℹ Retained 5171 / 9069 variants (polymorphic with DP >= 10)
#> ℹ Applying locus missingness filter
#> ℹ Retained 1734 / 5171 variants (per-variant missingness <= 0.2)
#> ℹ Applying sample missingness filter
#> ℹ Variants retained: 1734 | Samples retained: 18
#> ℹ Applying linked SNV filter
#> ℹ Retained 0 / 1734 variants (linked SNVs)
# see how many variants remained
nrow(vcf_multiSNP@variants)
#> [1] 0
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.
# extract individuals from VCFArrow object (keep all loci)
indivs1 <- c("Pb2Jp", "Pb2Scx", "Pb1Rd")
vcf_outgrp <- vcf_extract_samples(vcf, indivs1, f_invar = FALSE)
#> ℹ Compacting GT: 18 -> 3 samples across 1 chunk
#> ℹ Removed samples: Pv120, Pv126, Pv13, Pv14, Pv27, Pv28, Pv2, Pv31, Pv56, Pv62, Pv68, Pv73, Pv78, Pv79, and Pv93
#> ℹ Variants retained: 10000 | Samples retained: 3
vcf_ingrp <- vcf_extract_samples(vcf, indivs1, keep = FALSE, f_invar = FALSE)
#> ℹ Compacting GT: 18 -> 15 samples across 1 chunk
#> ℹ Removed samples: Pb2Jp, Pb2Scx, and Pb1Rd
#> ℹ Variants retained: 10000 | Samples retained: 15
# 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)
#> ℹ Compacting GT: 18 -> 3 samples across 1 chunk
#> ℹ Removed samples: Pv120, Pv126, Pv13, Pv14, Pv27, Pv28, Pv2, Pv31, Pv56, Pv62, Pv68, Pv73, Pv78, Pv79, and Pv93
#> ℹ Variants retained: 10000 | Samples retained: 3
vcf_ingrp <- vcf_extract_groups(vcf, groups1, f_invar = FALSE)
#> ℹ Compacting GT: 18 -> 15 samples across 1 chunk
#> ℹ Removed samples: Pb2Jp, Pb2Scx, and Pb1Rd
#> ℹ Variants retained: 10000 | Samples retained: 15
# bind vcf_outgrp and vcf_ingrp (same variants)
vcf1 <- vcf_bind(vcf_ingrp, vcf_outgrp)
#> ℹ Binding 2 VCFArrow objects (intersect): 10000 variants, 18 total samples.
# extract groups of individuals from VCFArrow object (remove invariant loci)
groups1 <- c("GS", "BS", "WA")
vcf_outgrp <- vcf_extract_groups(vcf, groups1, keep = FALSE)
#> ℹ Compacting GT: 18 -> 3 samples across 1 chunk
#> ℹ Applying invariant filter
#> ℹ Removed 8904 invariant variants; 1096 retained.
#> ℹ Removed samples: Pv120, Pv126, Pv13, Pv14, Pv27, Pv28, Pv2, Pv31, Pv56, Pv62, Pv68, Pv73, Pv78, Pv79, and Pv93
#> ℹ Variants retained: 1096 | Samples retained: 3
vcf_ingrp <- vcf_extract_groups(vcf, groups1)
#> ℹ Compacting GT: 18 -> 15 samples across 1 chunk
#> ℹ Applying invariant filter
#> ℹ Removed 1705 invariant variants; 8295 retained.
#> ℹ Removed samples: Pb2Jp, Pb2Scx, and Pb1Rd
#> ℹ Variants retained: 8295 | Samples retained: 15
# bind vcf_outgrp and vcf_ingrp (intersection - variable SNVs in all objects)
vcf1 <- vcf_bind(vcf_ingrp, vcf_outgrp, mode = "intersect")
#> ℹ Binding 2 VCFArrow objects (intersect): 421 variants, 18 total samples.
# bind vcf_outgrp and vcf_ingrp (union - variable SNVs in any object)
vcf1 <- vcf_bind(vcf_ingrp, vcf_outgrp, mode = "union", absent_as = "missing")
#> ℹ Binding 2 VCFArrow objects (union): 10000 variants, 18 total samples.
#> ℹ Absent genotypes will be filled 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).
# extract outgroup from VCFArrow
groups1 <- "OG"
vcf_outgrp <- vcf_extract_groups(vcf, groups1)
#> ℹ Compacting GT: 18 -> 3 samples across 1 chunk
#> ℹ Applying invariant filter
#> ℹ Removed 8904 invariant variants; 1096 retained.
#> ℹ Removed samples: Pv120, Pv126, Pv13, Pv14, Pv27, Pv28, Pv2, Pv31, Pv56, Pv62, Pv68, Pv73, Pv78, Pv79, and Pv93
#> ℹ Variants retained: 1096 | Samples retained: 3
# 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")
#> ℹ Compacting GT: 18 -> 15 samples across 1 chunk
#> ℹ Applying invariant filter
#> ℹ Removed 1705 invariant variants; 8295 retained.
#> ℹ Removed samples: Pb2Jp, Pb2Scx, and Pb1Rd
#> ℹ Variants retained: 8295 | Samples retained: 15
#> ℹ Applying indel filter
#> ℹ Retained 7706 / 8295 variants (non-Indels)
#> ℹ Applying biallelic filter
#> ℹ Retained 7706 / 7706 variants (biallelic SNVs)
#> ℹ Applying read coverage filter
#> ℹ Retained 4483 / 7706 variants (polymorphic with DP >= 10)
#> ℹ Applying locus missingness filter
#> ℹ Retained 2409 / 4483 variants (per-variant missingness <= 0.2)
#> ℹ Applying sample missingness filter
#> ℹ Variants retained: 2409 | Samples retained: 15
#> ℹ Applying unlinked SNV filter
#> ℹ Retained 2409 / 2409 variants (unlinked SNVs)
#> ℹ Binding 2 VCFArrow objects (outgroup): 2409 variants, 18 total samples.
#> ℹ Absent genotypes will be filled as: missing
# see how many variants remained
nrow(vcf1@variants)
#> [1] 2409
# 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()
#> ℹ Applying indel filter
#> ℹ Retained 9313 / 10000 variants (non-Indels)
#> ℹ Applying biallelic filter
#> ℹ Retained 9313 / 9313 variants (biallelic SNVs)
#> ℹ Applying read coverage filter
#> ℹ Retained 5388 / 9313 variants (polymorphic with DP >= 10)
#> ℹ Applying locus missingness filter
#> ℹ Retained 1951 / 5388 variants (per-variant missingness <= 0.2)
#> ℹ Applying sample missingness filter
#> ℹ Variants retained: 1951 | Samples retained: 18
#> ℹ Applying unlinked SNV filter
#> ℹ Retained 1951 / 1951 variants (unlinked SNVs)
# see how many variants remained
nrow(vcf1@variants)
#> [1] 1951
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.
# 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)
#> ℹ VCF is being written in 1 chunk
#> ✔ VCFArrow object successfully written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_filtered.vcf.gz'
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.
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')))
#> ℹ Accumulating Migrate-N (S): 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing Migrate-N (S) file...
#> ✔ Migrate-N file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_migrate.txt'
# arlequin http://cmpg.unibe.ch/software/arlequin35/
vcf2arlequin(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.arp')))
#> ℹ Accumulating Arlequin: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing Arlequin file...
#> ✔ Arlequin file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.arp'
# structure https://web.stanford.edu/group/pritchardlab/structure.html
vcf2structure(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.str')))
#> ℹ Accumulating Structure: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing Structure file...
#> ✔ Structure file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.str'
# faststucture http://rajanil.github.io/fastStructure/
vcf2structure(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.fstr')), method = "F")
#> ℹ Accumulating Structure: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing Structure file...
#> ✔ Structure file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.fstr'
# sNMF http://membres-timc.imag.fr/Olivier.Francois/snmf/index.htm
vcf2snmf(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.geno')))
#> ℹ Formatting sNMF: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing sNMF file...
#> ✔ sNMF file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.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')))
#> ℹ Building PLINK: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing PLINK files...
#> ℹ Chromosome names were recoded as integers for PLINK compatibility; mapping written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_admixture.chrommap'
#> ✔ PLINK binary fileset written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_admixture.bed', '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_admixture.bim', '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_admixture.fam'
# 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')))
#> ℹ Building PLINK: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing PLINK files...
#> ℹ Chromosome names were recoded as integers for PLINK compatibility; mapping written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_plink.chrommap'
#> ✔ PLINK binary fileset written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_plink.bed', '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_plink.bim', '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_plink.fam'
# 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')))
#> ℹ Accumulating PLINK .ped: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing PLINK file...
#> ℹ Chromosome names were recoded as integers for PLINK compatibility; mapping written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_plink.chrommap'
#> ✔ PLINK text fileset written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_plink.ped', '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_plink.map'
# eigenstrat https://github.com/DReichLab/EIG/tree/master
vcf2eigenstrat(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_eigenstrat')))
#> ℹ Building EIGENSTRAT: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing EIGENSTRAT files...
#> ✔ EIGENSTRAT fileset written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_eigenstrat.geno', '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_eigenstrat.ind', '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_eigenstrat.snp'
# genepop https://gitlab.mbb.univ-montp2.fr/francois/genepop
vcf2genepop(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.gen')))
#> ℹ Accumulating Genepop: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing Genepop file...
#> ✔ Genepop file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.gen'
# smartsnp https://github.com/ChristianHuber/smartsnp
vcf2smartsnp(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.smartsnp')))
#> ℹ Building SmartSNP: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing SmartSNP file...
#> ✔ SmartSNP file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.smartsnp'
# bayescan https://github.com/mfoll/BayeScan
vcf2bayescan(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.bayescan')))
#> ℹ Accumulating BayesScan: 1734 variants x 4 pops (0 MiB raw storage, vs 0 MiB with integer matrices)
#> ℹ Writing BayesScan file...
#> ✔ BayesScan file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.bayescan'
# bayesass https://github.com/brannala/BA3
vcf2bayesass(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.bayesass')))
#> ℹ Accumulating BayesAss: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing BayesAss file...
#> ✔ BayesAss file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.bayesass'
# treemix https://bitbucket.org/nygcresearch/treemix/wiki/Home
vcf2treemix(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.treemix')))
#> ℹ Building Treemix: 1734 variants x 4 pops (0 MiB raw storage)
#> ℹ Writing Treemix file...
#> ✔ Treemix file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.treemix'
# apparent https://github.com/halelab/apparent/tree/master
vcf2apparent(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.apparent')))
#> ℹ Accumulating Apparent: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing Apparent file...
#> ✔ Apparent file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.apparent'
# related https://github.com/timothyfrasier/related
vcf2related(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.related')))
#> ℹ Accumulating Related: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing Related file...
#> ✔ Related file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.related'
# long tidy dataframe of genotypes
vcf2gt_long(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.csv')), format = 'csv')
#> ℹ Writing gt_long: 1734 variants x 18 samples
#> ✔ gt_long table written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.csv'
# snapp https://www.beast2.org/snapp/
vcf2snapp(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_snapp.nex')))
#> ℹ Accumulating SNAPP: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing SNAPP file...
#> ✔ SNAPP file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_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')))
#> ℹ Accumulating Nexus: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing Nexus file...
#> ✔ Nexus file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_sdvq.nex'
# fasta https://www.ncbi.nlm.nih.gov/genbank/fastaformat/
vcf2fasta(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '.fna')))
#> ℹ Accumulating FASTA: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing FASTA file...
#> ✔ FASTA file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.fna'
# genlight object https://www.rdocumentation.org/packages/adegenet/versions/2.0.0/topics/genlight-class
genlight <- vcf2genlight(vcf)
#> ℹ Accumulating Genlight: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Building Genlight object...
# genlight object with an optional save
genlight <- vcf2genlight(vcf, out_file = file.path(res_path, paste0(project, postfix, fltr, '_genlight.rds')), save = TRUE)
#> ℹ Accumulating Genlight: 1734 variants x 18 samples (0 MiB raw storage)
#> ℹ Building Genlight object...
#> ℹ Writing Genlight object...
#> ✔ Genlight object written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub_genlight.rds'
##########
# 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')))
#> ℹ Building fineRADstructure: 0 variants x 18 samples (0 MiB raw storage)
#> ℹ Writing fineRADstructure file...
#> ✔ fineRADstructure file written to '/tmp/RtmpOHmUZa/vaillantii_discosnp_sub.finerad'