## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) ## ----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) ## ----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")) ## ----------------------------------------------------------------------------- # 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)) ## ----------------------------------------------------------------------------- # get a table of basic sample stats, including Watterson's Theta and Pi vcf_stats(vcf, res_path, paste0(project, postfix, fltr), theta = TRUE) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- # 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") ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- # 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) ## ----------------------------------------------------------------------------- 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')))