Read a VCF file and convert it to a
SeqArray GDS object
(SeqVarGDSClass; Zheng et al. 2017), with a
/radiator node containing radiator-specific metadata.
The function has an "advanced" mode (via ...) that allows a number of
filters to be applied on the fly (MAC, coverage, genotyping rate, LD, etc.),
as well as several VCF-specific clean-ups (Stacks, ipyrad, PLINK, DArT,
FreeBayes, GATK, samtools, etc.).
Used internally in radiator and might be of interest for users who want a fast and robust VCF → GDS import with basic QC.
Internal helper used by read_vcf() to synchronise radiator's
FILTERS field with the VCF FILTER annotation.
If the VCF contains more than one distinct FILTER value, markers with FILTER != "PASS" (e.g., low quality, failed filters, etc.) are tagged in radiator as:
"vcf filter column" in FILTERS
The function updates both markers.meta inside the GDS and the
filters.parameters tracking object.
read_vcf(
data,
strata = NULL,
filename = NULL,
vcf.stats = FALSE,
parallel.core = parallel::detectCores() - 1,
verbose = TRUE,
...
)
filter_vcf_filter_column(
gds,
markers.meta,
filters.parameters,
path.folder,
file.date,
verbose = TRUE
)(character)
Path to a VCF file (optionally bgzipped, *.vcf or *.vcf.gz).
Markers can be biallelic SNPs or haplotypes; downstream filters
(e.g. filter.haplotype.format) will normalise variants to the
desired representation.
(optional)
Strata definition, passed to read_strata.
Can be a path to a strata file or an object.
Default: strata = NULL.
(optional, character)
Base name of the GDS file to generate. Radiator will append
.gds.rad to the filename. If the chosen filename already exists
in path.folder, a timestamped default name is used instead.
Default: filename = NULL.
(logical, optional)
Generate basic statistics for individuals and markers (missingness, coverage,
etc.) and write them to disk.
Computational cost can be high for very large unfiltered VCF
Default: vcf.stats = FALSE.
(integer, optional)
Approximate number of cores to use where parallelism is supported
(e.g. SeqArray::seqApply, some filter helpers).
Default: parallel.core = parallel::detectCores() - 1.
(logical) Verbosity.
(optional) To pass further arguments for fine-tuning the function.
A SeqArray GDS object containing the imported VCF.
The current markers metadata tibble (from radiator).
The filters.parameters object to update.
Path to the filtering results folder.
Character string used for naming result files.
A SeqVarGDSClass object (radiator GDS) with a /radiator node
containing per-marker and per-individual metadata and a record of all
filters applied.
A list with:
markers.meta: updated tibble;
filters.parameters: updated parameters object.
Typical performance (rough order of magnitude):
a 35 GB VCF with ~4M SNPs: \(\sim\)7 minutes with 8 CPU;
a 21 GB VCF with ~2M SNPs: \(\sim\)5 minutes with 7 CPU.
The resulting GDS file can be reopened almost instantly in a later R session
with radiator::read_rad().
Heterozygosity and inbreeding (Fis):
The heterozygosity statistics generated here are global across strata and
primarily descriptive. For filtering on het/Fis/HWE, we recommend using
filter_het, filter_fis and
filter_hwe after obvious outlier individuals and
markers have been removed.
PLINK:
LOCUS is filled with an integer based on CHROM
(as.integer(factor(x = CHROM)));
COL is set to 1L (no within-read position available).
ipyrad:
the pattern "locus_" is stripped from CHROM;
COL is set to POS.
GATK / platypus / FreeBayes / samtools:
if the VCF ID column is ., it is replaced with
the position (POS);
short-read locus identity is assumed to be encoded in
CHROM + POS.
Stacks:
de novo: CHROM is typically "1";
LOCUS corresponds to "CHROM" in the Stacks VCF;
COL is POS - 1;
reference: ID is split into LOCUS, COL,
STRANDS.
DArT VCFs:
CHROM == "." is replaced by "denovo";
missing POS (NA) are set to 50;
COL is extracted from LOCUS (read position);
LOCUS is the first group of digits, then joined with POS
using "_", and POS is replaced by COL.
...)The ... lets you pass many additional arguments used by radiator’s
filtering framework, for example:
whitelist.markers, blacklist.id,
pop.select, pop.levels, pop.labels;
filter.strands – handle duplicate SNPs on opposite strands;
filter.common.markers – keep markers common to all strata;
filter.ma – global MAC/MAF/DP filtering;
filter.coverage, filter.genotyping;
filter.snp.position.read, filter.snp.number;
filter.short.ld, filter.long.ld, ld.method;
filter.individuals.missing, coverage filters by individual;
markers.info, vcf.metadata;
path.folder, random.seed, subsample.markers.stats;
filter.monomorphic, filter.haplotype.format, etc.
See the documentation of the individual filter functions
(e.g. filter_ma, filter_ld,
filter_monomorphic) for more detail.
Zheng X, Gogarten S, Lawrence M, Stilp A, Conomos M, Weir BS, Laurie C, Levine D (2017). SeqArray – A storage-efficient high-performance data format for WGS variant calls. *Bioinformatics*.
Danecek P, Auton A, Abecasis G et al. (2011) The variant call format and VCFtools. *Bioinformatics* 27:2156–2158.
if (FALSE) { # \dontrun{
# Simple import, no strata, defaults:
gds <- radiator::read_vcf(data = "populations.snps.vcf")
# With strata and a few filters:
gds <- radiator::read_vcf(
data = "populations.snps.vcf",
strata = "strata_salamander.tsv",
path.folder = "salamander",
filter.individuals.missing = "outliers",
filter.common.markers = TRUE,
filter.strands = "blacklist",
filter.ma = 4,
filter.genotyping = 0.3,
filter.snp.position.read = "outliers",
filter.short.ld = "mac",
filter.long.ld = NULL,
verbose = TRUE
)
# Later, in a new R session, reopen the GDS:
gds <- radiator::read_rad(data = "radiator_20200911@0748.gds")
} # }